算法逻辑:

将方程组Ax = b 中的A分解为A=LU,其中L为单位下三角矩阵,U为上三角矩阵,则方程组Ax=b化为解两个方程组Ly=b,Ux=y

matlab 代码如下

function test

A = [1 2 -12 8;
    5 4 7 -2;
    -3 7 9 5;
    6 -12 -8 3];

b = [27;4;11;49];
[n m] = size(A);
Ab = [A b];
mb = m + 1;
Ab(2:n, 1) = Ab(2:n, 1) / Ab(1, 1);
for k = 2 : n
    for j = k : mb
        Ab(k, j) = Ab(k, j) - Ab(k, 1:k-1)*Ab(1:k-1, j);
    end
    
    for i = k+1 : n
        Ab(i, k) = ( Ab(i, k) - Ab(i, 1:k-1)*Ab(1:k-1, k) ) / Ab(k, k);
    end
end
Ab
x = zeros(n, 1);
x(n,1) = Ab(n, mb) / Ab(n, n);
for i = n-1 : -1 : 1
    x(i, 1) = (Ab(i, mb) - Ab(i, i+1:n)*x(i+1:n, 1)) / Ab(i, i);
end
x