我数学就没及过格,爱莫能助
代码是:
A=[1 2 3 4;2 3 1 2;1 1 1 -1;1 0 -2 -6];
n=length(A);
L=eye(n);
U=A;
for k=1:n-1
for j=k+1:n
L(j,k)=U(j,k)/U(k,k);
end
for i=k+1:n
for j=k+1:n
U(i,j)=U(i,j)-L(i,k)*U(k,j);
end
end
end
U=triu(U);
L
U
norm(A-L*U)
运行结果是:
L =
1.0000 0 0 0
2.0000 1.0000 0 0
1.0000 1.0000 1.0000 0
1.0000 2.0000 1.6667 1.0000
U =
1.0000 2.0000 3.0000 4.0000
0 -1.0000 -5.0000 -6.0000
0 0 3.0000 1.0000
0 0 0 0.3333
ans =
0