doolittle.m

来自「线性方程组的数值解法」· M 代码 · 共 52 行

function Doolittle(A,B)
%LU-Doolittle分解法，参数说明(系数矩阵A,常数项B)
error(nargchk(2,2,nargin));%判断输入矩阵大小
d=max(size(A));%获取系数矩阵的阶
p=max(size(B));
if d==p
    dimention=d;
else
    fprintf('输入有误，不能构成线性方程组');
end
  
for k=1:d
    for j=k:d
        for r=1:k-1
            A(k,j)=A(k,j)-A(k,r)*A(r,j);
        end
    end
    for j=k:d
        if(j==k)
            continue
        else  
            for i=1:k-1        
                A(j,k)=A(j,k)-A(j,i)*A(i,k);
            end
            A(j,k)=A(j,k)/A(k,k);          
        end
    end
end

U=triu(A);
L=A-U;
for n=1:d
    L(n,n)=1;
end

L=L
U=U

for i=1:d
    for k=1:i-1
        B(i)=B(i)-A(i,k)*B(k);
    end
end
Y=B

for i=d:-1:1
    for k=i+1:d
        B(i)=B(i)-A(i,k)*B(k);
    end
    B(i)=B(i)*1/A(i,i);
end
X=B