算法逻辑:
将方程组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
Comments