矩阵论上机作业.doc

矩阵论上机作业 班级：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，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; 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;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;案例1、 crout分解%本函数将一个满秩方阵按Crout方式分解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); %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,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)=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)=sqrt(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 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.0769 0.0090 0.0840 0.0615 0.6805 1.0769 4.4924请验证Cholesky分解正确（是/否）？4、 吉文斯方法的QR分解function Q,R=qrgv(A)% 基于Givens变换，将方阵A分解为A=QR，其中Q为正交矩阵，R为上三角阵% 参数说明% A：需要进行QR分解的方阵% Q：分解得到的正交矩阵% R：分解得到的上三角阵% Q,R=qr(A) % 调用MATLAB自带的QR分解函数进行验证% q,r=qrgv(A) % 调用本函数进行QR分解% q*r-A % 验证 A=QR% q*q % 验证q的正交性% norm(q) % 验证q的标准化，即二范数等于1% 线性代数基础知识% 1.B=P*A*inv(P)，称A与B相似，相似矩阵具有相同的特征值% 2.Q*Q=I，称Q为正交矩阵，正交矩阵的乘积仍为正交矩阵% by dynamic of Matlab技术论坛% see also % contact me % 2010-01-17 22:51:18%n=size(A,1);R=A;Q=eye(n);for i=1:n-1 for j=2:n-i+1 x=R(i:n,i); rt=givens(x,1,j); r=blkdiag(eye(i-1),rt); Q=Q*r; R=r*R; endendfunction R,y=givens(x,i,j)% 求解标准正交的Given变换矩阵R，使用Rx=y，其中y(j)=0，y(i)=sqrt(x(i)2+x(j)2)% 参数说明% x：需要进行Givens变换的列向量% i：变为sqrt(x(i)2+x(j)2)的元素下标% j：变为0的元素的下标% R：Givens变换矩阵% y：Givens变换结果% 实例说明% x=1 3 5 9 6; % 将3等效到9上% R,y=givens(x,4,2) % 注意3的下标为2，9的下标为4% R*x-y % 验证Rx=y% R*R % 验证正交性% norm(R) % 验证标准性，就是范数为1% 关于Givens变换说明% 1.Givens矩阵是标准正交矩阵，也叫平面旋转矩阵，它是通过坐标旋转的原理将元素j的数值等效到元素i上% 2.Givens变换每次只能将一个元素变为0，而Householder变换则一次可以将任意个元素变为0% 3.Givens变换常用于将矩阵A变为对角阵%xi=x(i);xj=x(j);r=sqrt(xi2+xj2);cost=xi/r;sint=xj/r;R=eye(length(x);R(i,i)=cost;R(i,j)=sint;R(j,i)=-sint;R(j,j)=cost;y=x(:);y(i,j)=r,0;案例5、 豪斯霍尔德方法function Q,R=qrhs(A)% 基于Householder变换，将方阵A分解为A=QR，其中Q为正交矩阵，R为上三角阵% 参数说明% A：需要进行QR分解的方阵% Q：分解得到的正交矩阵% R：分解得到的上三角阵% 实例说明% A=-12 3 3;3 1 -2;3 -2 7;% Q,R=qr(A) % 调用MATLAB自带的QR分解函数进行验证% q,r=qrhs(A) % 调用本函数进行QR分解% q*r-A % 验证 A=QR% q*q % 验证q的正交性% norm(q) % 验证q的标准化，即二范数等于1% 线性代数基础知识% 1.B=P*A*inv(P)，称A与B相似，相似矩阵具有相同的特征值% 2.Q*Q=I，称Q为正交矩阵，正交矩阵的乘积仍为正交矩阵% 2010-01-17 22:49:51%n=size(A,1);R=A;Q=eye(n);for i=1:n-1 x=R(i:n,i); y=1;zeros(n-i,1); Ht=householder(x,y); H=blkdiag(eye(i-1),Ht); Q=Q*H; R=H*R;endfunction H,rho=householder(x,y)% 求解正交对称的Householder矩阵H，使Hx=rho*y，其中rho=-sign(x(1)*|x|/|y|% 参数说明% x：列向量% y：列向量，x和y必须具有相同的维数% 实例说明% x=3 5 6 8;% y=1 0 0 0;% H,rho=householder(x,y);% H*x-rho*y % 验证Hx=rho*y% H*H % 验证正交% 关于HouseHolder变换% 1.H=I-2vv，其中|v|=1，我们称H为反射(HouseHolder)矩阵，易证H对称且正交% 2.如果|x|=|y|，那么存在HouseHolder矩阵H=I-2vv，其中v=(x-y)/|x-y|，使Hx=y % 3.如果|x|y|，那么存在HouseH