《矩阵论上机作业.doc》由会员分享，可在线阅读，更多相关《矩阵论上机作业.doc（13页珍藏版）》请在金锄头文库上搜索。
1、 矩阵论上机作业 班级：070921 学号：07092002 姓名：黄涛 完成时间： 时间：2011年10月21号 一、 doolittle 方法 定义：设A Cn n . 如果A 可分解成A = LR,其中L 是对角元素为1 的下三角矩阵(称为单位下三角矩阵) , R 是上三角矩阵, 则称之为A 的Doolit tle分解5 .根据定义获得的计算公式r1 j = a1 j ( j = 1 ,2 , , n) ,li1 = ai1 / r11 ( i = 2 ,3 , , n) ,rkj = akj - k- 1t = 1lkt3 rtj ( j = k , k + 1 , , n; k =2
2、，3 , , n) ,lik = 1/ rkk3 aik - k- 1t = 1lit3 rtk ( i = k + 1 , ,n; k = 2 ,3 , , n) .由定义知,Doolit tle 分解是在n 行n 列的矩阵上进行的,由公式知道当i , j , k 较大时, lik , rkj 的计算量是比较大的,当矩阵的行数与列数不相同时是无法进行Doolittle 分解的, 为克服这些缺陷, 从分块的角度出发研究Doolit tle 分解. 使Doolit tle 分解能够适合可分任何矩阵。 程序代码#include#include#define N 20int main() int n
3、; coutn; static double aNN; int i,j; for(i=0;in;i+) for(j=0;jaij; /input A coutA=endl; for(i=0;in;i+) for(j=0;jn;j+) printf(%-8.2f,aij); coutendl;/output A static double lNN,uNN; int k; /step number for(i=0;in;i+) lii=1; for(k=0;kn;k+) /用k步实现，第k步计算：U的第k行，L的第k列 for(j=k;jn;j+) /U的第k行 ukj=akj; for(i=0;
4、i=k-1;i+) ukj-=lki*uij; for(i=k+1;in;i+) /L的第k列 lik=aik; for(j=0;j=k-1;j+) lik-=lij*ujk; lik/=ukk; coutL=endl; for(i=0;in;i+) for(j=0;jn;j+) printf(%-8.2f,lij); coutendl; /output L coutU=endl; for(i=0;in;i+) for(j=0;jn;j+) printf(%-8.2f,uij); coutendl; /output U return 0;案例一、 crout分解%本函数将一个满秩方阵按Crou
5、t方式分解function L,U=Crout(A) b=size(A); %b(1)行 %b(2)列 n=b(1);%这里只处理n*n的非奇异矩阵 %错误检查 if b(1)=b(2)%非方阵错误 error(MATLAB:Crout:Input Matrix should be a Square matrix. See Crout.); end if n=rank(A)%非满秩矩阵错误 error(MATLAB:Crout:Input Matrix should be FULL RANK. See Crout.); end %初始化 L=zeros(n,n); U=zeros(n,n);
6、%U中的主对角线元素均为1 for i=1:n U(i,i)=1; end %开始计算 for k=1:n for i=k:n %L中计算的方式是行为外循环，列为内循环 temp_sum=0; %临时和 for t=1:k-1 temp_sum=temp_sum+L(i,t)*U(t,k); end L(i,k)=A(i,k)-temp_sum; %计算L的第k列元素 end for j=k+1:n %U中计算的方式是列为外循环，行为内循环 temp_sum=0; %临时和 for t=1:k-1 temp_sum=temp_sum+L(k,t)*U(t,j); end U(k,j)=(A(k
7、,j)-temp_sum)/L(k,k);%计算U的第k行元素 end endend %end of function三、cholesky分解A=input(A=)m,n=size(A);if m=n %判断输入的矩阵是不是方阵 disp(输入的矩阵不是方阵，请重新输入); return;endd=eig(A); %根据方阵的特征值判定是不是正定矩阵for i=1:n if d(i)=0 disp(输入的矩阵不是正定矩阵，请重新输入); return; else break; endenddisp(输入的矩阵可以进行Cholesky分解); %如果是正定矩阵，可以进行下面的分解操作S(1,1)
8、=sqrt(A(1,1); %确定第1列元素for i=2:n S(i,1)=A(i,1)/S(1,1);endsum1=0;for j=2:n-1 %确定第j列元素 for k=1:j-1 sum1=sum1+S(j,k)2; end S(j,j)=sqrt(A(j,j)-sum1); sum1=0; for i=j+1:n for k=1:j-1 sum1=sum1+S(i,k)*S(j,k); end S(i,j)=(A(i,j)-sum1)/S(j,j); end sum1=0;endfor k=1:n-1 %确定第n列元素 sum1=sum1+S(n,k)2;endS(n,n)=sq
9、rt(A(n,n)-sum1);SP=S*S%验证分解是否正确disp(请验证Cholesky分解正确（是/否）？);end任取A=，则其仿真结果为：A = 1 6 9 6 5 3 3 5 6 4 2 4 3 7 9 0 4 3 5 4 6 5 4 9 3 8 5 7 6 8 9 3 6 4 9 6输入的矩阵可以进行Cholesky分解ans =Columns 1 through 4 1.0000 3.0000 3.0000 5.0000 0 0 - 2.0000i 0 - 1.0000i 0 -10.0000i 0 0 1.0000 1.0000 0 0 0 8.8882 0 0 0 0 0
10、 0 0 0 Column 5 through 63.0000 9.0000 0-12.5000 i 0 -28.5000i 3.5000 6.0000 12.7697 38.7593 0-4.6795i 0 -37.8279i0 25.0981 P = 1.0e+003 *0.0010 0.0030 0.0030 0.0050 0.0030 0.0090 0.0030 0.0130 0.0110 0.0350 0.0340 0.0840 0.0030 0.0110 0.0110 0.0260 0.0250 0.0615 0.0050 0.0350 0.0260 0.2050 0.2570 0.6805 0.0030 0.0340 0.0250 0.2570 0.3626 1
