matlab - 高斯消元法求解 A*x = b 线性系统 (MATLAB)

Question

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

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

有更好/更简单的方法吗？

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

带旋转的高斯消除如下。

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