matlab中矩阵LDLT分解与Cholesky分解.docx

矩阵LDLT分解与Cholesky分解：矩阵的LDLT消去函数的程序代码：%矩阵的LDLT分解function [s,l,d]=ldlt(a)s=1;l=0;d=0;%判断矩阵是否对称if a~=a' %矩阵不对称，输出错误信息 s=0;else b=diag(a); %列向量b存放矩阵a的对角元素，矩阵D的元素也放在该向量 n=size(a,1); %矩阵a维数nfor k=1:n b(k)=b(k)-(a(k,1:k-1).^2)*b(1:k-1);if ~b(k) %如果矩阵D的对角元素出现0，出现错误，停止计算 s=0;break else %进行递推 a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*(b(1:k-1).*a(k,1:k-1)'))/b(k);endendif s l=tril(a)-diag(diag(a))+diag(b); d=diag(b);endend矩阵的Cholesky消去函数程序代码：%矩阵的cholesky分解function [s,g]=cholesky(a)s=1;g=0;%判断矩阵是否对称正定if a~=a'|min(eig(a))<=0 %矩阵不是对称正定，输出错误信息 s=0;else %矩阵对称正定 n=size(a,1); %矩阵a维数na(1,1)=a(1,1)^(1/2);for i=2:nfor j=1:i-1 a(i,j)=(a(i,j)-a(i,1:j-1)*a(j,1:j-1)')/a(j,j); %对角线下方元素计算end a(i,i)=(a(i,i)-a(i,1:i-1)*a(i,1:i-1)')^(1/2); %对角线元素计算end g=tril(a); %输出结果end主程序为：%计算方法上机第二题clear;clc;%输入矩阵am=1:20;[mm,nn]=meshgrid(m);a=min(mm,nn);%首先进行ldlt分解[s,l,d]=ldlt(a); %调用ldlt分解函数if ~s %分解不能进行，输出错误信息disp('Error!Theldlt decomposition cannot go!');elsedisp('The Matrix L is: ');disp(l); %输出矩阵Ldisp('The Matrix D is: ');disp(d); %输出矩阵Ddisp('The Matrix LT is: ');disp(l'); %输出矩阵LTend%进行cholesky分解[s,g]=cholesky(a); %调用cholesky分解函数if ~s %分解不能进行，输出错误信息disp('Error!Thecholesky decomposition cannot go!');elsedisp('The Matrix G is: ');disp(g); %输出矩阵Gdisp('The Matrix GT is: ');disp(g'); %输出矩阵Gend矩阵LDLT分解与Cholesky分解：%用追赶法解方程组function x=chase(a,b,c,d)n=length(b);for k=2:n %进行"追"a(k-1)=a(k-1)/b(k-1);b(k)=b(k)-a(k-1)*c(k-1);d(k)=d(k)-a(k-1)*d(k-1);endd(n)=d(n)/b(n);for k=n-1:-1:1 %进行"赶"d(k)=(d(k)-c(k)*d(k+1))/b(k);endx=d';主程序为：%计算方法上机第三题clear;clc;n=20;a=-1+zeros(1,n-1); %对角线下方的元素b=4+zeros(1,n); %对角线上的元素c=-1+zeros(1,n-1); %对角线上方的元素d=[3,2+zeros(1,n-2),3]; %右端向量d=chase(a,b,c,d);disp('The solution is:');disp(d); %输出解向量