用MATLAB语言实现Cholesky分解求解SPD（对称正定）线性方程组

函数部分

function [L,X]= Cholesky(A,b)
N=size(A);%记录矩阵的大小
m=N(1);%矩阵的行数
n=N(2);%矩阵的列数
%判断输入矩阵是否为正定矩阵
if m~=n
    disp('输入的矩阵不是方阵，请重新输入')
    return;
end
%判断输入的矩阵是否为正定矩阵
d=eig(A);
for i=1:n
    if d(i)<0
        disp('输入的矩阵不是正定矩阵，请重新输入')
        return;
    else
        break;
    end
end
L=zeros(m,n);%初始化L矩阵
L(1,1)=sqrt(A(1,1));
for k=2:n
    L(k,1)=A(k,1)/L(1,1);
end
    for k=2:n
    L(k,k)=sqrt(A(k,k)-L(k,1:k-1)*L(k,1:k-1)');
    L(k+1:m,k)=(A(k+1:m,k)-L(k+1:m,1:k-1)*L(k,1:k-1)')/L(k,k);
    if L(k,k)==0
        break
    end
    end



% 前代法
y=L\b;
X=L'\y;

## 测试代码部分
M=diag(rand(100,1));
U=orth(rand(100,100));
A=U'*M*U;
b=rand(100,1);
[L,X]=Cholesky(A,b)
eps=norm(A*X-b);
display(eps)