《数值分析计算实习作业二(共20页).docx》由会员分享，可在线阅读，更多相关《数值分析计算实习作业二(共20页).docx（20页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上数 值 分 析计算实习题二学号：姓名： 院系：2015年11月28日专心-专注-专业一、算法采用拟上三角化的双步位移的QR分解法计算特征值，并利用列主元素Guass消去法求实数特征值对应的特征向量，定义全局变量n=10。1 主程序在定义一维和二维数组时将数组的大小定为n+1（大小为11），这样在使用单个元素时的角标可以与矩阵或向量元素的角标相对应。先输入矩阵A的各个元素；保留数组A，定义数组B的元素与A相同，即B=A，后面的矩阵变换用数组B来进行；对B进行拟上三角化；输出拟上三角化后的矩阵B；对拟上三角化后的B进行QR分解，输出矩阵Q、R、RQ；对拟上三角化的数组矩阵
2、B进行QR分解；输出A的全部特征值；判断所有特征值，如果特征值为实数（即储存特征值虚部的数组中元素为0），利用列主元素高斯消去法求其特征向量；输出实特征值对应的特征向量。2 拟上三角化程序对于r=1,2，n-2循环执行（1）求第r列的第r+2个元素到第n个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r的循环。如果和不为0，则继续执行。（2）计算第r列从第r+1到第n个元素的平方和的平方根当Br+1r小于等于0时，取cr=dr，否则取cr=-drhr=cr*cr-cr*Br+1r（3）令（4）计算（5）执行下一次循环。3 拟上三角化后的矩阵的QR分解初始令Q=I。对于r=1,2，n-1
3、循环执行（1）求第r列的第r+1个元素到第n个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r的循环。如果和不为0，则继续执行。（2）计算第r列从第r到第n个元素的平方和的平方根当Br+1r小于等于0时，取cr=dr，否则取cr=-drhr=cr*cr-cr*Br+1r（3）令（4）计算（5）执行下一次循环。循环结束后的Q和B就是分解后的Q和R，之后再求RQ。4 QR分解过程给定精度水平sigma=e-12，给定迭代的次数限制L，令k=1，m=n。利用goto和if判断语句在下面几个步骤中跳转。loop3：如果fabs(Bmm-1)sigma，则得到一个特征值Bmm，将其放入特征值数组
4、中，并令m=m-1，跳转到loop4上；否则跳转到loop5。loop4：如果m=1，则得到特征值B11放入特征值数组中，跳转到loop9，如果m=2，则求二阶子阵特征值放入特征值数组中，跳转到loop9，否则跳转到loop3。loop5：如果fabs(Bm-1m-2)=sigma，则对右下角相邻的一个二阶子阵执行求特征值函数得到两个特征值放入数组中，并且令m=m-2，跳转到loop4；否则跳转到loop6。loop6：如果k=L，则跳出整个循环，输出”未得到全部特征值；否则跳转到loop7。loop7：计算之后对M执行QR分解，同时得到新的经过QR分解后的矩阵B，跳转到loop8。loop8
5、：令k=k+1，跳转到loop3。loop9：结束整个循环。最后输出特征值。5 对M的QR分解对于r=1,2，m循环执行（1）将M矩阵的第r列的第r+1个元素到第m个元素的绝对值相加，如果和为0，则跳出该次循环直接进行下一次r的循环。如果和不为0，则继续执行。（2）计算矩阵M第r列从第r到第m个元素的平方和的平方根当Mrr小于等于0时，取cr=dr，否则取cr=-drhr=cr*cr-cr*Mrr（3）令（4）计算（5）执行下一次循环。最终得到的矩阵B便是对M进行QR分解之后随之变动的矩阵B。6 二阶子阵的特征值此函数的输入量为m，计算之后判断一元二次函数的是否小于零，从而得到两个特征值的实部
6、和虚部分别放入特征值数组的相应位置。7 列主元素高斯消去法参照课本的步骤，需要注意的是因为为奇异，所以在求特征向量的元素时令最后一个元素为1，再利用课本上的步骤倒推其他元素的值。二、程序#include#include#include#define L 2500/定义迭代次数#define n 10double An+1n+1=0,Bn+1n+1=0,Mn+1n+1=0,Xn+1n+1=0; double lamrn+1=0,lamcn+1=0;/lamr储存特征值的实部，lamc储存特征值的虚部double sigma=1e-12;/定义精度void hess();/拟上三角化函数void
7、AQR();/A的QR分解函数void QRmeth();/带双步位移的QR分解函数void MQRfenjie(int m);/M的QR分解函数void gauss(int y);/列主元素gauss消去法求特征向量函数/*以下为主程序*/void main() int i,j,k; for(i=1;i=n;i+)/输入矩阵A for(j=1;j=n;j+) if(i!=j) Aij=sin(0.5*i+0.2*j); else if(i=j) Aij=1.52*cos(i+1.2*j); for(i=1;i=n;i+)/令矩阵B=A，对B做拟上三角化 for(j=1;j=n;j+) Bij
8、=Aij; hess();/拟上三角化矩阵B AQR();/将拟上三角化后的矩阵B进行QR分解 QRmeth();/带双步位移的QR分解 for(k=1;k=n;k+)/判断特征值是否为实数，并求实数的特征向量 if(lamck!=0) continue;else for(i=1;i=n;i+)/求特征向量方程的矩阵系数 for(j=1;j=n;j+) if(i=j) Aij=Aij-lamrk; else Aij=Aij; gauss(k);/gauss消去法求特征向量 /*以下为拟上三角化程序*/void hess()/拟上三角化函数int r=0,i=0,j=0; double temp
9、1,sum=0;double dr=0,cr=0,hr=0,tr=0;double urn+1=0,prn+1=0,qrn+1=0,wrn+1=0;for(r=1;r=n-2;r+) for(temp1=0.0,i=r+2;i0) cr=-dr; else cr=dr; hr=cr*cr-cr*Br+1r; for(i=1;i=r;i+) uri=0.0; urr+1=Br+1r-cr; for(i=r+2;i=n;i+) uri=Air; for(j=1;j=n;j+) for(prj=0.0,i=1;i=n;i+) prj+=Bij*uri/hr; for(tr=0.0,i=1;i=n;i
10、+) tr+=pri*uri/hr; for(i=1;i=n;i+) for(qri=0.0,j=1;j=n;j+) qri+=Bij*urj/hr; for(i=1;i=n;i+) wri=qri-tr*uri; for(i=1;i=n;i+) for(j=1;j=n;j+) Bij-=(wri*urj+uri*prj); for(i=1;i=n;i+) for(j=1;j=n;j+) if(fabs(Bij)sigma) Bij=0.0;printf(拟上三角化后A(n-1):n);for(i=1;i=n;i+)/输出拟上三角化函数 for(j=1;j=n;j+) printf(%1.12
11、e,Bij); if(j%3=0) printf(n); printf(n);/*以下为QR分解程序*/void AQR()/QR分解函数double urn+1,wrn+1,prn+1,cr,dr,hr,Rn+1n+1,Qn+1n+1=0,rqn+1n+1=0;int i,j,r,k;for(i=1;i=n;i+)/令矩阵R等于矩阵B，对R进行QR分解 for(j=1;j=n;j+) Rij=Bij;for(i=1;i=n;i+) Qii=1;/Q初始为单位阵for(r=1;r=n-1;r+) for(dr=0.0,i=r;i0) cr=-dr; else cr=dr; hr=cr*cr-c
12、r*Rrr; for(i=1;ir;i+) uri=0; urr=Rrr-cr; for(i=r+1;i=n;i+) uri=Rir; for(i=1;i=n;i+) for (pri=0.0,j=1;j=n;j+) pri+=Rji*urj/hr; for(i=1;i=n;i+) for(j=1;j=n;j+) Rij=Rij-uri*prj; for(i=1;in;i+) for(wri=0.0,j=1;j=n;j+) wri+=Qij*urj; for(i=1;i=n;i+) for(j=1;j=n;j+) Qij-=wri*urj/hr; for(i=1;i=n;i+)/求RQ的矩阵值
13、 for(j=1;j=n;j+) for(rqij=0.0,k=1;k=n;k+) rqij+=Rik*Qkj; printf(生成的Q阵：n);for(i=1;i=n;i+) for(j=1;j=n;j+) if(fabs(Qij)sigma) Qij=0.0;for(i=1;i=n;i+)/输出矩阵Q for(j=1;j=n;j+) printf(%1.12e,Qij);if(j%3=0) printf(n); printf(n);printf(生成的R阵：n);for(i=1;i=n;i+) for(j=1;j=n;j+) if(fabs(Rij)sigma) Rij=0.0;for(i
14、=1;i=n;i+)/输出矩阵R for(j=1;j=n;j+) printf(%1.12e,Rij);if(j%3=0) printf(n); printf(n);printf(生成的RQ阵：n);for(i=1;i=n;i+) for(j=1;j=n;j+) if(fabs(rqij)sigma) rqij=0.0;for(i=1;i=n;i+)/输出矩阵RQ for(j=1;j=n;j+) printf(%1.12e,rqij);if(j%3=0) printf(n); printf(n);/*以下为带双步位移的QR分解的程序*/void QRmeth()/带双步位移的QR分解int i
15、,j,k=1,m=n,s;int nR=0,nC=0;double s1,t,x,sum,In+1n+1=0;for(i=1;i=n;i+) Iii=1;loop3:if(fabs(Bmm-1)=0) lamrm-1=(s1+sqrt(x)/2; lamrm=(s1-sqrt(x)/2; else lamrm-1=s1/2; lamcm-1=sqrt(-x)/2; lamrm=s1/2; lamcm=-sqrt(-x)/2; goto loop9; else goto loop3;loop5: if(fabs(Bm-1m-2)=0) lamrm-1=(s1+sqrt(x)/2; lamrm=(
16、s1-sqrt(x)/2; else lamrm-1=s1/2; lamcm-1=sqrt(-x)/2; lamrm=s1/2; lamcm=-sqrt(-x)/2; m-; m-; goto loop4; else goto loop6; loop6: if(k=L) printf(计算终止，未能得到全部特征值n); else goto loop7; loop7: s1=Bm-1m-1+Bmm; t=Bm-1m-1*Bmm-Bmm-1*Bm-1m; for(i=1;i=m;i+)/生成M for(j=1;j=m;j+) for(sum=0.0,s=1;s=m;s+) sum+=Bis*Bsj
17、; Mij=sum-s1*Bij+t*Iij; MQRfenjie(m);/M进行QR分解，求出新的矩阵B goto loop8; loop8:k+;goto loop3;loop9: ; printf(矩阵的全部特征值为：n);for(i=1;i=n;i+)/对于实数特征值和复数特征值两种不同的输出方法 if(lamci=0) printf(lam%d=%13.11len,i,lamri); else printf(lam%d=(%13.11le,%13.11le)n,i,lamri,lamci); /*以下为对M的QR分解的程序*/void MQRfenjie(int m)/QR分解函数i
18、nt r=0,i=0,j=0; double temp1,sum=0;double dr=0,cr=0,hr=0,tr=0;double urn+1=0,vrn+1=0,prn+1=0,qrn+1=0,wrn+1=0; for(r=1;r=m-1;r+) for(temp1=0.0,i=r+1;i0) cr=-dr; else cr=dr; hr=cr*cr-cr*Mrr; for(i=1;ir;i+) uri=0.0; urr=Mrr-cr; for(i=r+1;i=m;i+) uri=Mir; for(j=1;j=m;j+) for(vrj=0.0,i=1;i=m;i+) vrj+=Mij
19、*uri/hr; for(i=1;i=m;i+) for(j=1;j=m;j+) Mij-=uri*vrj; for(j=1;j=m;j+) for(prj=0.0,i=1;i=m;i+) prj+=Bij*uri/hr; for(i=1;i=m;i+) for(qri=0.0,j=1;j=m;j+) qri+=Bij*urj/hr; for(tr=0.0,i=1;i=m;i+) tr+=pri*uri/hr; for(i=1;i=m;i+) wri=qri-tr*uri; for(i=1;i=m;i+) for(j=1;j=m;j+)Bij-=(wri*urj+uri*prj); /*以下为
20、列主元素高消去法求特征向量程序*/void gauss(int y)/列主元的高斯消元法求解特征向量 int k,i,j,p; double temp1,temp2,bn+1=0,max,m; for(k=1;k=n-1;k+)/换行过程 max=Akk;p=k;for(i=k;imax) p=i;for(j=k;j=n;j+) temp1=Akj; Akj=Apj; Apj=temp1;for(i=k+1;i=n;i+)/消元过程 m=Aik/Akk; for(j=k+1;j=1;k-) temp2=0;for(j=k+1;j=n;j+) temp2+=Akj*Xyj;Xyk=(bk-tem
21、p2)/Akk; printf(对应特征值%1.12e的特征向量为:n,lamry); for(i=1;i=n;i+) printf(%1.12en,Xyi);三、结果拟上三角化的结果QR分解后的Q矩阵QR分解后的R矩阵QR分解后的RQ矩阵矩阵的全部特征值（带括号的表示复数，逗号前表示实部，逗号后表示虚部）实数特征值对应的特征向量四、讨论1、对于复数特征值的储存方法，采用两个数组，一个储存实部，一个储存虚部，因为复数特征值只在二阶子阵的特征值求解中产生，很容易判断特征值是否为复数，也很容易分别求出复数特征值的实部和虚部，因此可以将复数特征值储存在两个数组中。2、关于特征向量的求解，需要求解这个方程组，对于之前讲的迭代法，高斯消去法等线性方程组的解法前提都是系数矩阵非奇异，而这个方程组的系数矩阵显然是奇异的（否则X只有零解了），所以在采用高斯消去法时在求解过程中先规定最后一个元素的值为1，然后再倒推出其余各元素。