lltpingfan.m

来自「这是线性议程组的直接解法」· M 代码 · 共 30 行

function [L,LT,x]=LLTpingfan(A,b)
% A为对称正定矩阵
n=length(A);
for k=1:n
    if k==1
        L(k,k)=sqrt(A(k,k));
        L(k+1:n,k)=A(k+1:n,k)/L(k,k);
    else
        L(k,k)=sqrt(A(k,k)-L(k,1:k-1)*L(k,1:k-1)');
        if k~=n
            L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*L(k,1:k-1)')/L(k,k);
        end
    end
end
LT=L';
for k=1:n
    if k==1
        y(k,1)=b(k,1)/L(k,k);
    else
        y(k,1)=(b(k,1)-L(k,1:k-1)*y(1:k-1,1))/L(k,k);
    end
end
for k=n:-1:1
    if k==n
        x(k,1)=y(k,1)/LT(k,k);
    else
        x(k,1)=(y(k,1)-LT(k,k+1:n)*x(k+1:n,1))/LT(k,k);
    end
end
end