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 = 1lkt 3 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;iaij; /input Acout“A=“endl;for(i=0;in;i+)for(j=0;jn;j+)printf(“%-8.2f“,aij);coutendl;/output Astatic double lNN,uNN;int k; /step numberfor(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-=
4、lki*uij;for(i=k+1;in;i+) /L 的第 k 列lik=aik;for(j=0;j=k-1;j+)lik-=lij*ujk;lik/=ukk;cout“L=“endl; for(i=0;in;i+)for(j=0;jn;j+)printf(“%-8.2f“,lij);coutendl; /output Lcout“U=“endl; for(i=0;in;i+)for(j=0;jn;j+)printf(“%-8.2f“,uij);coutendl; /output Ureturn 0;案例1、 crout分解%本函数将一个满秩方阵按Crout方式分解function L,U=
5、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.);endif n=rank(A)%非满秩矩阵错误error(MATLAB:Crout:Input Matrix should be FULL RANK. See Crout.);end %初始化L=zeros(n,n);U=zeros(n,n); %U中的主对角线元素均为1for i=1:nU(i,i)=1;
6、end%开始计算for k=1:nfor i=k:n %L中计算的方式是行为外循环，列为内循环temp_sum=0; %临时和for t=1:k-1temp_sum=temp_sum+L(i,t)*U(t,k);endL(i,k)=A(i,k)-temp_sum; %计算L的第k列元素endfor j=k+1:n %U中计算的方式是列为外循环，行为内循环temp_sum=0; %临时和for t=1:k-1temp_sum=temp_sum+L(k,t)*U(t,j); endU(k,j)=(A(k,j)-temp_sum)/L(k,k);%计算U的第k行元素endendend %end of
7、 function三、cholesky 分解A=input(A=)m,n=size(A);if m=n %判断输入的矩阵是不是方阵disp(输入的矩阵不是方阵，请重新输入);return;endd=eig(A); %根据方阵的特征值判定是不是正定矩阵for i=1:nif d(i)=0disp(输入的矩阵不是正定矩阵，请重新输入);return;elsebreak;endenddisp(输入的矩阵可以进行 Cholesky 分解); %如果是正定矩阵，可以进行下面的分解操作S(1,1)=sqrt(A(1,1); %确定第 1 列元素for i=2:nS(i,1)=A(i,1)/S(1,1);e
8、ndsum1=0;for j=2:n-1 %确定第 j 列元素for k=1:j-1sum1=sum1+S(j,k)2;endS(j,j)=sqrt(A(j,j)-sum1);sum1=0;for i=j+1:nfor k=1:j-1sum1=sum1+S(i,k)*S(j,k);endS(i,j)=(A(i,j)-sum1)/S(j,j);endsum1=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 分解正确（是/否）？);en
9、d任取 A= ，则其仿真结果为：A =1 6 9 6 5 33 5 6 4 2 43 7 9 0 4 35 4 6 5 4 93 8 5 7 6 89 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.5000i3
10、.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.00900.0030 0.0130 0.0110 0.0350 0.0340 0.08400.0030 0.0110 0.0110 0.0260 0.0250 0.06150.0050 0.0350 0.0260 0.2050 0.2570 0.68050.0030 0.0340 0.0250 0.2570 0.3626 1.07690.0090 0.0840 0.0615 0
11、.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 的标准化，即二范
12、数等于 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-1for j=2:n-i+1x=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
13、,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% 关于 Given
14、s 变换说明% 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、豪斯霍尔德方法functi
15、on 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*i
16、nv(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-1x=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|
17、% 参数说明% 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|，那么存在 HouseHolder 矩阵 H，使 Hx=rho*y，其中 rho=-sign(x(1)*|x|/|y|% 4.Householder 矩阵，常用于将一个矩阵 A 通过正交变换对角阵%x=x(:);y=y(:);if length(x)=length(y)error(The Column Vectors X and Y Must Have The Same Length!);endrho=-sign(x(1)*norm(x)/norm(y);y=rho*y;v=(x-y)/norm(x-y);I=eye(length(x);H=I-2*v*v案例
