cholesky分解

cholesky分解

-

2023年3月16日发(作者：干洗行业)

functionx=cholesky(A,b)

%求正定对称矩阵A的楚列斯基分解

%b为方程组的右端项

%x为方程组的Ax=b的解

[m,n]=size(A);

ifm~=n

disp('ArgumentmatrixAmustbesquare!');

return;

elseifm~=length(b)

disp('ThedimentionsofAandbdonotagree!');

return;

end

L=zeros(n);x=zeros(n,1);y=zeros(n,1);%分别对L，x，y赋初

始值

fork=1:n;

if(A(k,k)-L(k,1:k-1)*L(k,1:k-1)')<=0

disp('error');

return;

elseL(k,k)=sqrt(A(k,k)-L(k,1:k-1)*L(k,1:k-1)');%

求下三角矩阵L对角线的元素

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);%求L对角线下每列的元素

end

y(k)=(b(k)-L(k,1:k-1)*y(1:k-1))/L(k,k);%

求下三角方程组Ly=b的解y；

end

x(n)=y(n)/L(n,n);

fork=n:-1:1

x(k)=(y(k)-L(k+1:n,k)'*x(k+1:n))/L(k,k);%

求上三角方程组L'x=y的解x

end

a=10*ones(1,100);

c=ones(1,99);

d=ones(1,99);

A=diag(a)+diag(c,1)+diag(d,-1);

b=12*ones(100,1);

b(1)=11;

b(100)=11;

x=cholesky(A,b)