《北航数值分析大作业第二题精解.pdf》由会员分享，可在线阅读，更多相关《北航数值分析大作业第二题精解.pdf（11页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、目标目标：使用带双步位移的QR 分解法求矩阵A aij10*10的全部特征值，并对其中的每一个实特征值求相应的特征向sin(0.5i0.2 j)(i j)a 量。已知：ij1.5cos(i1.2 j)(i j)(i,j=1,2,10)算法算法：开始输入矩阵 A用函数“nishangsanjiadiv”将矩阵 A 拟上三角化即为 A用函数“characteristic”求解矩阵 A（n-1）即 A的所有特征值用函数“qrdiv”将 A（n-1）QR 分解用函数“characteristicvector”求解矩阵 A（n-1）即 A 的所有实输出 A 的所有特征值、A 的所有实特征值对应的特征向量
2、、拟上三角矩阵 A（n-1）、及其 Q、R 和 R*Q以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵 A 的实特征值对应的一个特征向量的求法：In1n1BAiI选主元的gause消元Bu 0 AiIu 0检查知无重特征值01n1Q1n10u 0由于AiI=0，因此在经过选主元的高斯消元以后，AiI即B的最后一行必然为零，左上方变为 n-1 阶单位矩阵In1n1，右上方变为 n-1 阶向量Q1n1，然后令un 1，则ujQjj 1,2,n1。这样即求出所有 A 所有实特征值对应的一个特征向量。#include#include#include#def
3、ine N 10#define E 1.0e-12#define MAX 10000/以下是符号函数double sgn(double a)double z;if(aE)z=1;else z=-1;return z;/以下是矩阵的拟三角分解void nishangsanjiaodiv(double ANN)int i,j,k;int m=0;double d,c,h,t;double uN,pN,qN,wN;for(i=0;iN-2;i+)for(j=i+2;jN;j+)if(Aji=E)m=m+1;if(m=(N-2-i)continue;for(j=i+1,d=0;jN;j+)d=d+Aj
4、i*Aji;d=sqrt(d);c=-1*sgn(Ai+1i)*d;h=c*c-c*Ai+1i;for(j=i+2;jN;j+)uj=Aji;for(j=0;ji+2;j+)uj=0;ui+1=Ai+1i-c;for(j=0;jN;j+)for(k=i+1,pj=0;kN;k+)pj=Akj*uk+pj;pj=pj/h;for(j=0;jN;j+)for(k=i+1,qj=0;kN;k+)qj=Ajk*uk+qj;qj=qj/h;for(j=0,t=0;jN;j+)t=t+pj*uj;t=t/h;for(j=0;jN;j+)wj=qj-t*uj;for(j=0;jN;j+)for(k=0;kN
5、;k+)Ajk=Ajk-wj*uk-uj*pk;/以下是矩阵的QR分解void qrdiv(double ANN,double QNN,double RNN)int i,j,k;/int m=0;double d,c,h;double uN,wN,pN;for(i=0;iN;i+)for(j=0;jN;j+)if(i=j)Qij=1;else Qij=0;for(i=0;iN;i+)for(j=0;jN;j+)Rij=Aij;for(i=0;iN-1;i+)/for(j=i+1;jN;j+)if(Rji=E)m=m+1;/if(m=(N-1-i)continue;for(j=i,d=0;jN;
6、j+)d=d+Rji*Rji;d=sqrt(d);c=-1*sgn(Rii)*d;h=c*c-c*Rii;for(j=i+1;jN;j+)uj=Rji;for(j=0;ji;j+)uj=0;ui=Rii-c;for(j=0;jN;j+)for(k=0,wj=0;kN;k+)wj=Qjk*uk+wj;for(j=0;jN;j+)for(k=0;kN;k+)Qjk=Qjk-wj*uk/h;for(j=0;jN;j+)for(k=i,pj=0;kN;k+)pj=Rkj*uk+pj;pj=pj/h;for(j=0;jN;j+)for(k=0;kN;k+)Rjk=Rjk-uj*pk;/矩阵的QR分解/以
7、下是二次多项式求根double root(double b,double c)double m;m=b*b-4*c;return m;/二次多项式求根/以下是求解矩阵的所有特征值void characteristic(double ANN,double chaRN,double chaIN)int k=0,m=N-1;int i,j;int L;double s,t,x;double MNN,BNN;int f=0;double d,c,h;double uN,wN,pN;double QNN,RNN;for(L=0;LMAX;L+)next:if(m=0)chaR0=A00;chaI0=0;
8、break;if(fabs(Amm-1)=E)x=sqrt(x);chaRm=s/2+x/2;chaRm-1=s/2-x/2;chaIm=0;chaIm-1=0;elsex=sqrt(fabs(x);chaRm=s/2;chaRm-1=s/2;chaIm=x/2;chaIm-1=-x/2;break;if(fabs(Am-1m-2)=E)x=sqrt(x);chaRm=s/2+x/2;chaRm-1=s/2-x/2;chaIm=0;chaIm-1=0;elsex=sqrt(fabs(x);chaRm=s/2;chaRm-1=s/2;chaIm=x/2;chaIm-1=-x/2;m=m-2;go
9、to next;for(i=0;i=m;i+)for(j=0;j=m;j+)if(i=j)for(k=0,Mij=0;k=m;k+)Mij=Aik*Akj+Mij;Mij=Mij-s*Aij+t;elsefor(k=0,Mij=0;k=m;k+)Mij=Aik*Akj+Mij;Mij=Mij-s*Aij;/以下是M的QR分解for(i=0;i=m;i+)for(j=0;j=m;j+)if(i=j)Qij=1;else Qij=0;for(i=0;i=m;i+)for(j=0;j=m;j+)Rij=Mij;for(i=0;im;i+)for(j=i+1;j=m;j+)if(Rji=E)f=f+1
10、;if(f=(m-i)continue;for(j=i,d=0;j=m;j+)d=d+Rji*Rji;d=sqrt(d);c=-1*sgn(Rii)*d;h=c*c-c*Rii;for(j=i+1;j=m;j+)uj=Rji;for(j=0;ji;j+)uj=0;ui=Rii-c;for(j=0;j=m;j+)for(k=0,wj=0;k=m;k+)wj=Qjk*uk+wj;for(j=0;j=m;j+)for(k=0;k=m;k+)Qjk=Qjk-wj*uk/h;for(j=0;j=m;j+)for(k=i,pj=0;k=m;k+)pj=Rkj*uk+pj;pj=pj/h;for(j=0;j
11、=m;j+)for(k=0;k=m;k+)Rjk=Rjk-uj*pk;for(j=0;j=m;j+)for(k=0;k=m;k+)Mjk=Qjk;/以上是M的QR分解for(i=0;i=m;i+)for(j=0;j=m;j+)for(k=0,Bij=0;k=m;k+)Bij=Mki*Akj+Bij;for(i=0;i=m;i+)for(j=0;j=m;j+)for(k=0,Aij=0;k=m;k+)Aij=Bik*Mkj+Aij;/以下是求矩阵的所有特征值的特征向量void eigenvector(double VNN,double TN)double ANN,baozNN,guodN;dou
12、ble c;int i,j,k,m,t;int W=0;for(i=0;iN;i+)for(j=0;jN;j+)baozij=Vij;for(t=0;t6;t+)for(i=0;iN;i+)for(j=0;jN;j+)Aij=baozij;for(i=0;iN;i+)Aii=Aii-Tt;for(i=0;iN-1;i+)for(j=i;jE)k=j;break;for(j=i;jN;j+)guodj=Aij;Aij=Akj;Akj=guodj;for(j=i;jE)for(m=i;mN;m+)Ajm=Ajm/c;for(j=0;jN;j+)c=Aji;if(j!=i)for(m=i;m=0;i
13、-)Vti=AiN-1;/以下是主函数void main()double aNN,bNN,chaRN,chaIN;double qNN,rNN,qrNN;double shiyanN;double f,g;int i,j,k;for(i=0;iN;i+)for(j=0;jN;j+)if(i!=j)aij=sin(0.5*(i+1)+0.2*(j+1);elseaij=1.5*cos(i+1)+1.2*(j+1);nishangsanjiaodiv(a);printf(矩阵A的拟上三角分解：n);for(i=0;iN;i+)for(j=0;jN-5;j+)if(fabs(aij)E)aij=0;
14、printf(%22.11e,aij);printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+)if(fabs(aij)E)aij=0;printf(%22.11e,aij);printf(n);printf(n);qrdiv(a,q,r);printf(n);printf(n);printf(n);printf(拟上三角矩阵A的QR分解：n);printf(上三角矩阵R：n);for(i=0;iN;i+)for(j=0;jN-5;j+)if(fabs(rij)E)rij=0;printf(%22.12e,rij);printf(n);printf(
15、n);for(i=0;iN;i+)for(j=N-5;jN;j+)if(fabs(rij)E)rij=0;printf(%22.12e,rij);printf(n);printf(n);printf(正交矩阵Q：n);for(i=0;iN;i+)for(j=0;jN-5;j+)if(fabs(qij)E)qij=0;printf(%22.12e,qij);printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+)if(fabs(qij)E)qij=0;printf(%22.12e,qij);printf(n);printf(n);for(i=0;iN;
16、i+)for(j=0;jN;j+)for(k=0,qrij=0;kN;k+)qrij=qrij+rik*qkj;printf(n);printf(n);printf(n);printf(R*Q：n);for(i=0;iN;i+)for(j=0;jN-5;j+)if(fabs(qrij)E)qrij=0;printf(%22.12e,qrij);printf(n);printf(n);printf(n);printf(n);for(i=0;iN;i+)for(j=N-5;jN;j+)if(fabs(qrij)E)qrij=0;printf(%22.12e,qrij);printf(n);prin
17、tf(n);printf(n);printf(n);characteristic(a,chaR,chaI);for(i=1;iN;i+)if(i3)f=chaRi;g=chaIi;chaRi=chaR7+i;chaIi=chaI7+i;chaR7+i=f;chaI7+i=g;if(i=5)f=chaRi;g=chaIi;chaRi=chaR7;chaIi=chaI7;chaR7=f;chaI7=g;printf(矩阵A所有特征值：n);for(j=0;jN;j+)if(fabs(chaIj)E)printf(%2d=%18.11e+%18.11ein,j+1,chaRj,chaIj);else
18、 printf(%2d=%19.11e%19.11ein,j+1,chaRj,chaIj);printf(n);printf(n);printf(n);for(i=0;iN;i+)for(j=0;jN;j+)if(i!=j)aij=sin(0.5*(i+1)+0.2*(j+1);elseaij=1.5*cos(i+1)+1.2*(j+1);/重新输入矩阵Aeigenvector(a,chaR);printf(相应实特征值对应的特征向量：n);for(i=0;i6;i+)printf(%d的一个特征向量为：n,(i+1);for(j=0;jN-5;j+)printf(%22.11e,aij);printf(n);for(j=N-5;jN;j+)printf(%22.11e,aij);printf(n);getch();