Устранение Гаусса для решения линейной системы A * x = b (MATLAB)

Question

Я пытаюсь создать код, который решает линейные системы A * x = b.

Я сделал приведенный ниже код, используя процесс исключения Гаусса, и он работает каждый раз, когда в A нет нулей. Если в А есть нули, то иногда это работает, иногда нет. В основном я пробую альтернативу "A \ b" в MATLAB.

Есть ли лучший / более простой способ сделать это?

A = randn(5,5);
b = randn(5,1);

nn = size(A);
n = nn(1,1);
U = A;
u = b;

for c = 1:1:n
    k = U(:,c);
    for r = n:-1:c
        if k(r,1) == 0
            continue;
        else
            U(r,:) = U(r,:)/k(r,1);
            u(r,1) = u(r,1)/k(r,1);
        end
    end
    for r = n:-1:(c+1)
        if k(r,1) == 0
            continue;
        else
            U(r,:) = U(r,:) - U(r-1,:);
            u(r,1) = u(r,1) - u(r-1,1);
        end
    end
end
x = zeros(size(b));
for r = n:-1:1
    if r == n
        x(r,1) = u(r,1);
    else
        x(r,1) = u(r,1);
        x(r,1) = x(r,1) - U(r,r+1:n)*x(r+1:n,1);
    end
end
error = A*x - b;
for i = 1:1:n
    if abs(error(i)) > 0.001
        disp('ERROR!');
        break;
    else
        continue;
    end
end
disp('x:');
disp(x);

Рабочий пример с 0:

A = [1, 3, 1, 3;
    3, 4, 4, 1;
    3, 0, 3, 9;
    0, 4, 0, 1];

b = [3;
    4;
    5;
    6];

Пример, который терпит неудачу (A * x-b не [0])

A = [1, 3, 1, 3;
    3, 4, 4, 1;
    0, 0, 3, 9;
    0, 4, 0, 1];
b = [3;
    4;
    5;
    6];

Объяснение моего алгоритма: Допустим, у меня есть следующая матрица A:

|4, 1, 9|
|3, 4, 5|
|1, 3, 5|

Для первого столбца я делю каждую строку на первое число в строке, поэтому каждая строка начинается с 1

|1, 1/4, 9/4|
|1, 4/3, 5/3|
|1,   3,   5|

Затем я вычитаю последнюю строку с той, что над ней, и затем сделаю то же самое для строки выше и так далее.

|1,     1/4,     9/4|
|0, 4/3-1/4, 5/3-9/4|
|0,   3-4/3,   5-5/3|

|1,    0.25,   2.250|
|0,   1.083, -0.5833|
|0,   1.667,   3.333|

Затем я повторяю то же самое для остальных столбцов.

|1,    0.25,   2.250|
|0,       1, -0.5385|
|0,       1,   1.999|

|1,    0.25,    2.250|
|0,       1,  -0.5385|
|0,       0,  -8.7700|

|1,    0.25,    2.250|
|0,       1,  -0.5385|
|0,       0,        1|

Те же самые операции, которые я выполняю в A, я делаю в b, поэтому система остается эквивалентной.

повторно UPDATE:

Я добавил это сразу после "для c = 1: 1: n"

Поэтому, прежде чем что-либо делать, он сортирует строки A (и b), чтобы в столбце «c» были элементы с уменьшением (0 будут оставлены в нижних строках A). Прямо сейчас это, кажется, работает для любой обратимой квадратной матрицы, хотя я не уверен, что это будет.

r = c;
a = r + 1;
while r <= n
    if r == n
        r = r + 1;
    elseif a <= n
        while a <= n
            if abs(U(r,c)) < abs(U(a,c))
                UU = U(r,:);
                U(r,:) = U(a,:);
                U(a,:) = UU;
                uu = u(r,1);
                u(r,1) = u(a,1);
                u(a,1) = uu;
            else
                a = a+1;
            end
        end
    else
        r = r+1;
        a = r+1;
    end
end

Answer 1

Устранение по Гауссу с поворотом выглядит следующим образом.

function [L,U,P] = my_lu_piv(A)
n = size(A,1);
I = eye(n);
O = zeros(n);
L = I;
U = O;
P = I;
    function change_rows(k,p)
        x = P(k,:); P(k,:) = P(p,:); P(p,:) = x;
        x = A(k,:); A(k,:) = A(p,:); A(p,:) = x;
        x = v(k); v(k) = v(p); v(p) = x;
    end

    function change_L(k,p)
        x = L(k,1:k-1); L(k,1:k-1) = L(p,1:k-1);
        L(p,1:k-1) = x;
    end
for k = 1:n
    if k == 1, v(k:n) = A(k:n,k);
    else
        z = L(1:k-1,1:k -1)\ A(1:k-1,k);
        U(1:k-1,k) = z;
        v(k:n) = A(k:n,k)-L(k:n,1:k-1)*z;
    end
    if k<n
        x = v(k:n); p = (k-1)+find(abs(x) == max(abs(x))); % find index p
        change_rows(k,p);
        L(k+1:n,k) = v(k+1:n)/v(k);
        if k > 1, change_L(k,p); end
    end
    U(k,k) = v(k);
end
end

Для того, чтобы решить систему ..

% A x = b (1) исходная система% L U = P A
(2) факторизация P A или A (p, :) в продукт L U% P A x = P b (3) умножить обе части (1) на P% L U x = P b
(4) подставить (2) в (3)%, пусть y = U x (5) определить y как U x% let c = P b (6) определяют c как P b% L y = c
(7) подузлы (5) и (6) в (4)% U * x = y (8) a переписать (5)

Для этого ..

% [L U p] = lu (A); % факторизовать% y = L \ (P * b); % вперед решить из (7) нижнюю треугольную систему% x = U \ y; % обратный ход (8), верхняя треугольная система

Answer 2

Алгоритм Гаусса предполагает, что матрица преобразуется в верхнюю треугольную матрицу. Это не происходит в вашем примере. Результат вашего алгоритма

A =

   1   3   1   3
   3   4   4   1
   0   0   3   9
   0   4   0   1

U =

   1.00000   3.00000   1.00000   3.00000
  -0.00000   1.00000  -0.20000   1.60000
   0.00000   0.00000   1.00000   3.00000
   0.00000   4.00000  -0.00000   1.00000

Как видите, это не верхняя треугольная фигура. Вы пропускаете строки, если элемент разворота равен нулю. Это не работает. Чтобы это исправить, вам нужно поменять местами столбцы в матрице и строки в векторе, если элемент разворота равен нулю. В конце вы должны поменять местами строки в вашем результате b соотв. u.

Гауссовский алгоритм:

1 Set n = 1
2 Take pivot element (n, n)
3 If (n, n) == 0, swap column n with column m, so that m > n and (n, m) != 0 (swap row m and n in vector b)
4 Divide n-th row by pivot element (divide n-th row in vector b)
5 For each m > n
6   If (m, n) != 0
7      Divide row m by m and subtract element-wise row n (same for vector b)
8 n = n + 1
9 If n <= number of rows, go to line 2

С точки зрения числовой устойчивости лучше всего использовать максимум каждой строки в качестве основного элемента. Также вы можете использовать максимум матрицы как элемент поворота, меняя столбцы и строки. Но не забудьте поменять местами b и вернуться обратно в свое решение.