求解A*x = b线性系统的高斯消元法（MATLAB）

Question

我正在尝试编写求解 A*x = b 线性系统的代码。

我使用高斯消去法编写了下面的代码，如果 A 中没有任何 0，它每次都有效。如果 A 中有零，则有时有效，有时无效。基本上我正在尝试替代 MATLAB 中的 "A\b"。

有没有 better/simpler 的方法？

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 中执行的操作相同，因此系统保持等效。

重新更新：

我在 "for 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

高斯算法假设矩阵转换为上三角矩阵。在您的示例中不会发生这种情况。你的算法的结果是

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 并换回您的解决方案。

Answer 2

旋转高斯消元法如下

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

为了解决系统..

% Ax = b (1) original system % LU = PA
(2) factorization of PA or A(p,:) into the product LU % PAx = Pb (3) multiply both sides of (1) by P % LUx = Pb
(4) substitute (2) into (3) % let y = Ux (5) define y as Ux % let c = Pb (6) define c as Pb % Ly = c
(7) subsitute (5) and (6) into (4) % U*x = y (8) a rewrite of (5)

要做到这一点..

% [L U p] = lu (A) ; % factorize % y = L \ (P*b) ; % forward solve of (7), a lower triangular system % x = U \ y ; % backsolve of (8), an upper triangular system

Answer 3

试试这个：

Ab = [A,b] % Extended matrix of the system of equations
rref(Ab)   % Result of applying the Gauss-Jordan elimination to the extended matrix

有关更多详细信息和示例，请参阅 rref documentation。

求解A*x = b线性系统的高斯消元法（MATLAB）

Gauss elimination to solve A*x = b linear system (MATLAB)

matlab

linear-algebra