北航数值分析大作业 第二题 QR分解
北航 数值分析第二次大作业(带双步位移的QR方法)
一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。
总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。
为了防止矩阵在后面的计算中被破坏保存A[][]。
2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。
由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。
这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。
3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。
用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。
4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。
为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。
若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。
这里用函数intTwoStepDisplacement_QR()来实现。
这是用来存储计算得到的特征值的二维数组。
考虑到特征值可能为复数，因此将所有特征值均当成复数处理。
此函数中，QR分解部分用子函数void QR_decompositionMk()实现。
这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。
5）计算特征向量得到特征值后，计算实特征值相应的特征向量。
这里判断所得特征值的虚数部分是否为零。
求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。
因此先初始化矩阵Array，计算(A-λI)，再求解方程组。
北航数值分析大作业二(纯原创,高分版)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令 以及最大迭代步数L；
2.若m≤0，则结束计算，已求出A的全部特征值，判断 或 或m≤2是否成立，成立则转3，否则转4；
3.若 ，则得一个特征值 ，m=m-1，降阶；若 ，则计算矩阵：
的特征值得矩阵A的两个特征值，m=m-2，降阶，转2.；
4.若k≤L，成立则令
k=k+1，转2，否则结束计算，为计算出矩阵A的全部特征值；
北航数理统计大作业2-聚类与判别分析讲解
应用数理统计作业二学号：姓名：电话：二〇一四年十二月对NBA球队的聚类分析和判别分析摘要：NBA联盟作为篮球的最高殿堂深受广大球迷的喜爱，联盟的30支球队大家也耳熟能详，本文选取NBA联盟30支球队2013-2014常规赛赛季场均数据。
利用spss软件通过聚类分析对27个地区进行实力类型分类，并利用判断分析对其余3支球队对分类结果进行验证。
可以看出各球队实力类型与赛季实际结果相吻合。
关键词：聚类分析，判别分析，NBA目录1. 引言 (4)2、相关统计基础理论 (5)2.1、聚类分析 (5)2.2，判别分析 (6)3.聚类分析 (7)3.1数据文件 (7)3.2聚类分析过程 (9)3.3 聚类结果分析 (11)4、判别分析 (12)4.1 判别分析过程 (12)4.2判别检验 (17)5、结论 (20)参考文献 (21)致谢 (22)1. 引言1896年，美国第一个篮球组织"全国篮球联盟（简称NBL）"成立，但当时篮球规则还不完善，组织机构也不健全，经过几个赛季后，该组织就名存实亡了。
1946年4月6日，由美国波士顿花园老板沃尔特.阿.布朗发起成立了“美国篮球协会”（简称BAA）。
1949年在布朗的努力下，美国两大篮球组织BAA和NBL合并为“全国篮球协会”（简称NBA）。
NBA季前赛是 NBA各支队伍的热身赛，因为在每个赛季结束后，每支球队在阵容上都有相当大的变化，为了让各队磨合阵容，熟悉各自球队的打法，确定各队新赛季的比赛阵容、同时也能增进队员、教练员之间的沟通，所以在每个赛季开始之前，NBA就举办若干场季前赛，使他们能以比较好的状态投入到漫长的常规赛的比赛当中。
为了扩大NBA在全球的影响，季前赛有约三分之一的球队在美国以外的国家举办。
从总体上看，NBA的赛程安排分为常规赛、季后赛和总决赛。
常规赛采用主客场制，季后赛和总决赛采用七场四胜制的淘汰制。
[31]NBA常规赛从每年的11月的第一个星期二开罗，到次年的4月20日左右结束。
北航数值分析计算实习题目二 矩阵QR分解
数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。
2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。
Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
北航数值分析第二次大作业--QR分解
《数值分析A》计算实习题目二姓名学号联系方式班级指导教师2012年10月一、算法设计方案整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。
因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。
1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。
2.对经过拟上三角化处理的矩阵进行QR分解。
3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。
计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。
二、源程序代码#include<stdio.h>#include<math.h>#include<string.h>int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s;double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10]; double X[10][10],Y[10][10],Qt[10][10],M[10][10];double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0}; double epsilon=1e-12;void main(){void Quasiuppertriangular(double A[][10]);void QRdecomposition(double A[][10]);void DoublestepsQR(double A[][10]);int i,j;for(i=0;i<10;i++){for(j=0;j<10;j++){A[i][j]=sin(0.5*(i+1)+0.2*(j+1));Q[i][j]=0;AA[i][j]=A[i][j];}A[i][i]=1.5*cos(2.2*(i+1));AA[i][i]=A[i][i];}Quasiuppertriangular(A); //调用拟上三角化函数printf( "\n A经过拟上三角化矩阵为：\n\n");for(i=0;i<10;i++) //输出拟上三角化矩阵{for(j=0;j<10;j++){printf("%.12e ",A[i][j]); //输出拟上三角化矩阵}printf( "\n\n");}QRdecomposition(A); //调用QR分解函数printf( " 进行QR分解后，R矩阵为：\n\n"); //输出R矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",R[i][j]);}printf( "\n\n");}printf( " Q矩阵为:\n\n"); //输出Q矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",Q[i][j]);}printf( "\n\n");}printf( " RQ矩阵为：\n\n"); //输出RQ矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",RQ[i][j]);}printf( "\n\n");}DoublestepsQR(A); //调用双步位移函数printf( "\n\n 特征值实部依次为：\n\n"); //输出特征值实部for(j=0;j<10;j++){printf("%.12e ",Re[j]);}printf("\n\n 特征值虚部依次为：\n\n "); //输出特征值虚部for(j=0;j<10;j++){printf("%.12e ",Im[j]);}//按行输出特征向量printf( "\n\n 按行输出实特征根相应特征向量为：\n\n");for(i=0;i<10;i++){if(i==1||i==2||i==5||i==6){continue;}for(j=0;j<10;j++){printf("%.12e ",X[i][j]);}printf( "\n\n");}getchar();}//拟上三角化函数void Quasiuppertriangular(double A[][10]) {for(j=0;j<8;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;T[i]=0;W[i]=0;}m=0;for(i=j+2;i<10;i++){if(A[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j+1;i<10;i++){d=d+pow(A[i][j],2);}d=sqrt(d);c=-d;if(A[j+1][j]<=0){c=d;}h=c*(c-A[j+1][j]);U[j+1]=A[j+1][j]-c;for(i=j+2;i<10;i++){U[i]=A[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*A[k][i];}P[i]=P[i]/h;}t=0;for(i=0;i<10;i++){for(k=0;k<10;k++){T[i]=T[i]+U[k]*A[i][k];}T[i]=T[i]/h;t=t+P[i]*U[i];}t=t/h;for(i=0;i<10;i++){W[i]=T[i]-t*U[i];for(k=0;k<10;k++){A[i][k]=A[i][k]-W[i]*U[k]-U[i]*P[k];if(abs(A[i][k])<1e-12){A[i][k]=0;}}}}}//QR分解函数void QRdecomposition(double A[][10]) {for(i=0;i<10;i++){for(j=0;j<10;j++){RQ[i][j]=0;Q[i][j]=0;R[i][j]=A[i][j];}Q[i][i]=1;}for(j=0;j<9;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;W[i]=0;}m=0;for(i=j+1;i<10;i++){if(R[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j;i<10;i++){d=d+pow(R[i][j],2);}d=sqrt(d);c=-d;if(R[j][j]<=0){c=d;}h=c*(c-R[j][j]);U[j]=R[j][j]-c;for(i=j+1;i<10;i++){U[i]=R[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){W[i]=W[i]+U[k]*Q[i][k];}}for(i=0;i<10;i++){for(k=0;k<10;k++){Q[i][k]=Q[i][k]-((W[i]*U[k])/h);}}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*R[k][i];}P[i]=P[i]/h;}for(i=0;i<10;i++){for(k=0;k<10;k++){R[i][k]=R[i][k]-U[i]*P[k];if(abs(R[i][k])<epsilon){R[i][k]=0;}}}}for(i=0;i<10;i++) //计算A(n+1)=RQ {for(j=0;j<10;j++){for(k=0;k<10;k++){RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];}}}}//双步位移法计算特征值特征向量函数void DoublestepsQR(double A[][10]){int L=1000,m=9; //定义最大循环次数for(i=0;i<L;i++){for(;m>-1;){if(abs(A[m][m-1])<=epsilon){Re[m]=A[m][m];m=m-1; //降阶if(m==0) //4{Re[0]=A[0][0];break;}if(m==-1){break;}if(m>1){continue;}}b=-A[m][m]-A[m-1][m-1]; //5c=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];if(m==1) //6{if((b*b-4*c)>=0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-1; //循环出口条件break;}if((m>1)&&(abs(A[m-1][m-2])>epsilon)) //8{if(i==L-1){printf("No results! \n");m=0; //循环出口条件break;}break;}if((m>1)&&(abs(A[m-1][m-2])<=epsilon)) //7 {if((b*b-4*c)>0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-2; //降阶if(m>0){continue;}if(m==0){Re[0]=A[0][0];break;}}}if(m<=0){break;}s=A[m-1][m-1]+A[m][m]; //9t=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=0;Q[j][k]=0;M[j][k]=0;X[j][k]=0;Y[j][k]=0;}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){M[j][k]=M[j][k]+A[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){M[j][k]=M[j][k]-s*A[j][k];}M[j][j]=M[j][j]+t;}//调用QR分解函数对M矩阵进行分解并传递参数矩阵QQRdecomposition(M);for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=Q[k][j];}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){X[j][k]=X[j][k]+Qt[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];}}}for(j=0;j<10;j++){{A[j][k]=Y[j][k];}}}//应用列主元高斯消元法计算实部特征向量for(l=0;l<10;l++){if(l==1||l==2||l==5||l==6){continue;}for(k=0;k<10;k++){for(m=0;m<10;m++){A[k][m]=AA[k][m];}A[k][k]=A[k][k]-Re[l];}for(j=0;j<9;j++){m=j;for(i=j+1;i<10;i++){if(abs(A[i][j])>abs(A[m][j])){m=i;}}{Y[j][k]=A[j][k];A[j][k]=A[m][k];A[m][k]=Y[j][k];}for(k=j+1;k<10;k++){b=A[k][j]/A[j][j];for(i=j;i<10;i++){A[k][i]=A[k][i]-A[j][i]*b;}}}X[l][9]=1;for(i=8;i>=0;i--){c=0;for(j=i+1;j<10;j++){c=c+A[i][j]*X[l][j];}X[l][i]=-c/A[i][i];}}}三、程序输出结果1819。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
数值分析大作业之2
数值分析大作业二一、算法设计方案(一)算法流程1.将要求解的矩阵进行Householder变换，进行拟上三角化，并输出拟上三角化的结果；2.将拟上三角化后的矩阵进行带双步位移的QR分解(最大迭代次数1000次)，逐个求出特征值，输出QR法结束后得到的Q、R、RQ阵，输出求出的特征值（用实部和虚部表示）3.对于求出来的实特征值，使用带原点平移的反幂法求出对应的特征向量，输出这些特征向量。
（二）程序设计流程1. 定义精度和最大迭代次数；2. 使用容器vector进行编程（方便增减元素），使用传引用调用（提高执行效率）；3. 将各个步骤用到的数学算法封装成函数，方便调用。
具体需要的函数如下：double VectMod(vector< double > &p)：求向量的模vector< double > NumbMultiVect(vector< double > &vect, double a)：数乘向量double VectMultiVect(vector< double > &y, vector< double > &u)：求两个向量的积vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector<double > &u)：矩阵的转置乘向量vector< double > ArraMultiVect(vector< vector< double > > &B, vector< double > &u)：矩阵乘向量vector< vector< double > > VectMultiCovVect(vector< double > &a, vector< double >&b)：向量乘向量转置得到矩阵void ArraSubs(vector< vector< double > > &C, vector< vector < double > > D) ：两个矩阵相减Vector< double > VectSubs(vector< double > &a, vector< double > &b)：两个向量相减vector<vector<double>> GetM(vector<vector<double>> &A)：求解矩阵M （用于双部位移QR迭代用）double GetMode(vector < vector < double > > &B, const int r)：求解矩阵B的r列向量的模double GetModeH(vector < vector < double > > &B, const int r)：求解矩阵B的第r列向量的模，用于拟对角化vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)：一个实数乘矩阵bool IsBirZeroH(vector< vector< double > > &B, const int r)：判断B[i][r]对角线下是否为零void GausElim(vector< vector< double > > a)：列主元高斯消元法求齐次方程解向量void Stop(vector< vector< double > > &Ar)：停止，结束程序void SolutS1S2(complex< double > &s1, complex< double > &s2, vector< vector<double > > &A)：求解二阶子阵的特征值s1，s2；void Save2(complex< double > &s1, complex< double > &s2)：保存两个特征值void Save1(complex< double > &s)：保存一个特征值void JudgemBelow2(vector< vector< double > > &A, vector< vector< double > > Abk)：对于m == 1 及m == 0 的处理void Hessenberg(vector< vector< double > > &A)：矩阵拟上三角化void QRMethod(vector< vector< double > > A)：矩阵QR分解void CalculatAk(vector< vector < double > > &Ak)：带双步位移QR迭代法二、源程序#include "stdafx.h"#include <vector>#include <iostream>#include <math.h>#include <complex>#include <fstream>using namespace std;const double epsion = 1e-12;const int L = 1000;int m,n;int k = 1;vector< vector< double > > I;vector< complex< double > > Lambda;///////////////////////////以下为自定义的算法流程中用到的函数double VectMod(vector< double > &p) //求向量的模{double value = 0.0;vector< double >::size_type i,j;j = p.size();for (i=1; i<j; i++){value += p[i] * p[i];}value = sqrt(value);return value;}vector< double > NumbMultiVect(vector< double > &vect, double a) //数乘向量{int j = vect.size();vector< double > b(j, 0);for (int i=1; i<j; i++){b[i] = a * vect[i];}return b;}double VectMultiVect(vector< double > &y, vector< double > &u)//两个向量相乘{vector< double >::size_type a = y.size();double value = 0;for (vector< double >::size_type i=1; i<a; i++){value += y[i] * u[i];}return value;}vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector< double > &u)//矩阵的转置乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[j][i] * u[j];}}return vec;}}vector< double > ArraMultiVect(vector< vector< double > > &B,vector< double > &u) //矩阵乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[i][j] * u[j];}}return vec;}}vector< vector<double> > GetM(vector< vector <double> > &A){int a = A.size();double s = A[a-2][a-2] + A[a-1][a-1];double t = A[a-2][a-2] * A[a-1][a-1] - A[a-1][a-2] * A[a-2][a-1];vector<vector<double>> D(a, vector< double >(a, 0));for (int i=1; i<a; i++){{double sum = 0;for (int k=1; k<a; k++){sum += A[i][k] * A[k][j];}D[i][j] = sum - s * A[i][j] + t *(i==j ? 1.0 : 0);}}return D;}double GetMode(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}double GetModeH(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r+1; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)//数乘//向量{int b = D.size();vector< vector< double > > U(b, vector< double >(b, 0));for (int i=1; i<b; i++){{U[i][j] = a *D[i][j];}}return U;}bool IsBirZero(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+1; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}bool IsBirZeroH(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+2; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}vector< vector< double > > VectMultiCovVect(vector< double > &a,vector< double > &b)//向量乘向量的转置得到矩阵{int s1 = a.size();int s2 = b.size();if (s1 != s2){cerr << "Vectors not match in size ! ";}else{vector< vector< double > > U(s1, vector< double >(s1, 0));for (int i=1; i<s1; i++){for (int j=1; j<s1; j++){U[i][j] = a[i] * b[j];}}return U;}}void ArraSubs(vector< vector< double > > &C,vector< vector < double > > D){int a = C.size();int b = D.size();int c = C[0].size();int d = D[0].size();if (a!=b || c!=d){cerr << "Vectors not match in size !";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){C[i][j] -= D[i][j];}}}}vector< double > VectSubs(vector< double > &a,vector< double > &b){int s1 = a.size();int s2 = b.size();if(s1 != s2){cerr << "Vectors not match in size !";}else{vector< double > value(s1,0);for (int i=1;i<s1;i++){value[i] = a[i] - b[i];}return value;}}void GausElim(vector< vector< double > > a) //高斯消元{int sz = a.size();vector< int > fx, ufx,ufxp, record;vector< double > x(sz, 0);vector< vector< double > > ret;//*****消元过程********for (int k = 1; k < sz-1; k++){double max = a[k][k];int p=0;for (int i=k+1; i<sz; i++){if (abs(a[i][k]) > abs(max)){p =i;max = a[i][k];}}//选出主元行if (abs(max) >= epsion){if (p != 0){for (int j=k; j<sz; j++){double temp = a[k][j];a[k][j] = a[p][j];a[p][j] = temp;}//交换主元行}for (int i=k+1; i<sz; i++){double tt = a[i][k] / a[k][k];for (int j=k+1; j<sz; j++){a[i][j] = a[i][j] - tt * a[k][j];}}}// end of if (abs(max) >= epsion)}// end of for (int k = 1; k < sz; k++)for (int i=1; i<sz; i++){int p=0;for (int j=1; j<=i;j++){if(abs(a[j][i]) >= 10*epsion)//不为零{p=j;}}record.push_back(p);}for (int i=1; i<sz; i++){int p=0;for (int j=sz-2; j>=0; j--){if (record[j] == i){p=j + 1;}}if (p != 0){ufxp.push_back(p);ufx.push_back(i);}}for (int i=1; i<sz; i++){int p = 0;for (int j=0; j < ufxp.size(); j++){if (ufxp[j] == i){p = 1;}}if (p == 0){fx.push_back(i);}}//end of for (int i=1; i<sz; i++)int c = fx.size();//***************************************************** for (int i=0; i<c; i++){for (int j=0; j<c; j++){if (i == j){x[fx.at(j)] = 1.0 + epsion;}else{x[fx.at(j)] = 0;}}int b = ufxp.size();for (int s=b-1; s>=0; s--){double temp=0;for(int t=ufxp.at(s)+1; t<sz; t++){temp += x[t] * a[ufx.at(s)][t];}x[ufxp.at(s)] = (0 - temp) / a[ufx.at(s)][ufxp.at(s)];}ret.push_back(x);}//end of for (int i=0; i<c; i++)//******************************************************* int sz1 = ret.size();int sz2 = ret[0].size();ofstream result2("CharactVector .txt", ios::app);result2.precision(12);for (int i=0; i<sz1; i++){for (int j=1;j<sz2; j++){result2 << scientific << ret[i][j] << endl;}result2 << endl << endl;}result2 << "下一个特征值的特征向量！" << endl;}void Stop(vector< vector< double > > &Ar){int a = Lambda.size();for (int i=0; i<a; i++){if (abs(Lambda[i].imag())<=epsion){vector< vector< double > > An=Ar;ArraSubs(An,NumbMultiArra(I,Lambda[i].real()));GausElim(An);}}ofstream result("result.txt");result.precision(12);vector<complex< double > >::iterator i,j;j = Lambda.end();for (i = Lambda.begin(); i<j; i++){result << scientific <<(*i) << endl;}exit(0);}void SolutS1S2(complex< double > &s1, complex< double > &s2,vector< vector< double > > &A){double s = A[m-1][m-1] + A[m][m];double t = A[m-1][m-1] * A[m][m] - A[m][m-1] * A[m-1][m];double det = s *s - 4.0 * t;if (det > 0){s1.imag(0);s1.real ((s + sqrt(det)) / 2);s2.imag (0);s2.real((s - sqrt(det)) / 2);}else{s1.imag(sqrt(0 - det) / 2);s1.real(s / 2);s2.imag(0 - sqrt(0 - det) / 2);s2.real(s / 2);}}void Save2(complex< double > &s1, complex< double > &s2)//存储特征值s1，s2 {Lambda.push_back(s1);Lambda.push_back(s2);}void Save1(complex< double > &s){Lambda.push_back(s);}void JudgemBelow2(vector< vector< double > > &A,vector< vector< double > > Abk){if (m == 1){complex< double > lbd;lbd.imag(0);lbd.real(A[1][1]);Save1(lbd);Stop(Abk);}else if (m == 0){Stop(Abk);}}//拟上三角化void Hessenberg(vector< vector< double > > &A){int a = A.size();vector< double > u(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;for (int r=1; r<a-2; r++){if (!IsBirZeroH(A, r)){d = GetModeH(A,r);c = (A[r+1][r] > 0) ?(0-d) : d;h = c * c - c * A[r+1][r];for (int j=1; j<a; j++){if (j < r+1){u[j] = 0;}else if (j == r+1){u[j] = A[j][r] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)p = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);q = NumbMultiVect(ArraMultiVect(A, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(A, VectMultiCovVect(w, u));ArraSubs(A, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("NISHANGDANJIAO.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}}void QRMethod(vector< vector< double > > A){int a = A.size();vector< vector< double > > Q(a,vector < double > (a,0));for (int i=1; i<a; i++){Q[i][i] = 1.0;}vector< vector< double > > RQ(A);//vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);vector< double > l(a,0);double d, c, h, t;for (int r=1; r<a-1; r++){if (!IsBirZero(A, r)){d = GetMode(A,r);c = (A[r][r] > 0) ?(0-d) : d;h = c * c - c * A[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = A[j][j] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)l = ArraMultiVect(Q, u);ArraSubs(Q, VectMultiCovVect(l,NumbMultiVect(u, 1/h)));v = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);ArraSubs(A, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(RQ, u), 1 / h);q = NumbMultiVect(ArraMultiVect(RQ, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(RQ, VectMultiCovVect(w, u));ArraSubs(RQ, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("QR.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << Q[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << RQ[i][j] << " ";}result << endl;}}void CalculatAk(vector< vector < double > > &Ak){int a = Ak.size();vector< vector < double > > B(a,vector < double > (a,0));vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;B = GetM(Ak);for (int r=1; r<a-1; r++){if (!IsBirZero(B, r)){d = GetMode(B,r);c = (B[r][r] > 0) ?(0-d) : d;h = c * c - c * B[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = B[j][j] - c;}else{u[j] = B[j][r];}}// end of for (int j=1; j<=m; j++)v = NumbMultiVect(ConveArraMultiVect(B, u), 1 / h);ArraSubs(B, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(Ak, u), 1 / h);q = NumbMultiVect(ArraMultiVect(Ak, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(Ak, VectMultiCovVect(w, u));ArraSubs(Ak, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)}int _tmain(int argc, _TCHAR* argv[]){vector< vector< double > > A(11, vector< double >(11,0));vector< vector< double > > Abk;I=A;for (int i=1; i<11; i++){vector< double > temp;for (int j=1; j<11; j++){if(i != j){A[i][j] = sin(0.5 * i + 0.2 * j);I[i][j] = 0;}else{A[i][j] = 1.5 * cos(i + 1.2 *j);I[i][j] = 1.0;}}}Abk = A;n = A.size() - 1;m = n;//初始化问题Hessenberg(A);QRMethod(A);while (1){if (abs(A[m][m-1]) <= epsion){complex< double > lbdm;lbdm.imag(0);lbdm.real(A[m][m]);Save1(lbdm);m -= 1;//对A进行降维处理！！！！A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();}JudgemBelow2(A, Abk);}else{complex< double > va1, va2;SolutS1S2(va1, va2, A);if (m == 2){Save2(va1,va2);Stop(Abk);}//end of if (m == 2)if ( abs(A[m-1][m-2]) <= epsion){Save2(va1,va2);m = m - 2;//矩阵降维A.pop_back();A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();A[i].pop_back();}JudgemBelow2(A, Abk);}else{if (k == L){cerr << "Stop without solution";exit(-1);}else{CalculatAk(A);k += 1;}}}// end of if (abs(A[m][m-1]) >= epsion)}}三、实验结果（1）A经过拟上三角化后得到的矩阵-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000-7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000-8.782436382928e-002 -9.255889387184e-001 6.032599440534e-0011.518860956469e-001-2.347878362416e+000 2.372370104937e+000 1.819290822208e+0003.237804101546e-001 2.205798440320e-001 2.102692662546e+0001.816138086098e-001 1.278839089990e+000 -6.380578124405e-001-4.154075603804e-0011.0556********e-016 1.728274599968e+000 -1.171467642785e+000-1.243839262699e+000 -6.399758341743e-001 -2.002833079037e+0002.924947206124e-001 -6.412830068395e-001 9.783997621285e-0022.557763574160e-001-5.393383812774e-017 -3.165873865181e-017 -1.291669534130e+000-1.111603513396e+000 1.171346824096e+000 -1.307356030021e+0001.803699177750e-001 -4.246385358369e-001 7.988955239304e-0021.608819928069e-0011.533464996622e-017 5.963321406181e-017 0.000000000000e+0001.560126298527e+000 8.125049397515e-001 4.421756832923e-001-3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-0021.300562737229e-016 -3.097060010889e-017 0.000000000000e+0000.000000000000e+000 -7.707773755194e-001 -1.583051425742e+000-3.042843176799e-001 2.528712446035e-001 -6.709925401449e-0012.544619929082e-0011.610216724767e-016 -2.211571837369e-016 0.000000000000e+0000.000000000000e+000 6.483933100712e-017 -7.463453456938e-001-2.708365157019e-002 -9.486521893682e-001 1.195871081495e-0011.929265617952e-0021.368550186199e-016 7.151513190805e-017 0.000000000000e+0000.000000000000e+000 -2.522454384246e-017 1.072074718358e-016-7.701801374364e-001 -4.697623990618e-001 4.988259468008e-0011.137691603776e-001-2.780851300718e-017 -6.708630788363e-017 0.000000000000e+0000.000000000000e+000 -3.282796821065e-017 4.879774287475e-0171.854829286220e-016 7.013167092107e-001 1.582180688475e-0013.862594614233e-001-2.124604440055e-017 -1.707979758930e-016 0.000000000000e+0000.000000000000e+000 1.013496750890e-016 -4.153758325241e-0171.222621691887e-016 -5.551115123126e-017 4.843807602783e-0013.992777995177e-001（2）Q-3.519262579534e-001 4.427591982245e-001 -6.955982513607e-0016.486200753651e-002 3.709718861896e-001 1.855847143605e-001-1.628942319628e-002 -1.181053169648e-001 -5.255375383724e-002 -5.486582943568e-002-9.360277287361e-001 -1.664679186545e-001 2.615299548560e-001 -2.438671728934e-002 -1.394774360893e-001 -6.977585391241e-0026.124472142963e-003 4.440505443139e-002 1.975907909728e-0022.062836970533e-002-4.208697095111e-017 -8.810520554692e-001 -3.989762796959e-0013.720308728479e-002 2.127794064090e-001 1.064463557221e-001-9.343171079758e-003 -6.774200464527e-002 -3.014340698675e-002 -3.146955080444e-002-2.150178169911e-017 4.009681353156e-017 -5.371806806439e-001 -1.234945854205e-001 -7.063151608719e-001 -3.533456368496e-0013.101438948264e-002 2.248676491598e-001 1.000601783527e-0011.044622748702e-0016.113458775639e-018 -3.721194648197e-0177.068175586628e-0189.892235468621e-001 -1.239411731211e-001 -6.200358589825e-0025.442272839461e-003 3.945881637235e-002 1.755813350011e-0021.833059462907e-0025.184948268593e-017 -4.198303559531e-017 3.255071298412e-017-1.527665297025e-017 5.323610690264e-001 -6.733900344896e-0015.910581205868e-002 4.285425323867e-001 1.906901343193e-0011.990794495295e-0016.419444583601e-017 4.121668945107e-017 3.096964582901e-017-5.050360323651e-017 -7.078737686239e-017 -6.0597********e-001 -9.165783032818e-002 -6.645586508974e-001 -2.957110877580e-001 -3.0872********e-0015.455993559780e-017 -9.724896332186e-017 3.956603694870e-0171.913857001710e-018 -4.316583456405e-0172.797376665557e-0179.933396625117e-001 -9.690440311939e-002 -4.311990584470e-002-4.501694411183e-002-1.108640877071e-017 4.655237056115e-017 -1.0722********e-017 -9.470236121745e-018 4.277993227986e-017 8.866601870176e-017 -2.516725033065e-016 5.410088006061e-001 -5.817540838226e-001 -6.0734********e-001-8.470152033092e-018 9.650816729410e-017 -1.429362211392e-017 -2.789222052040e-017 -3.690793770141e-017 -8.090765462521e-017 -1.964050295697e-016 -6.772274683437e-017 -7.221591336735e-0016.917269588876e-001（3）R2.508342744917e+000 -2.185646885493e+000 -1.314609070786e+000-3.558787493836e-002 -2.609857850388e-001 -1.283121847090e+000 -1.390878610606e-001 -8.712897972161e-001 3.849367902971e-0013.353802899665e-0012.100627753398e-016 -1.961603277854e+000 2.407523727633e-0017.054714572823e-001 5.957204318279e-001 5.526978774676e-001-3.268209924413e-001 -5.769498668364e-002 2.871129330189e-001 -8.895128754189e-002-3.300197935770e-016 -1.872873410421e-016 2.404534601993e+0001.706758096328e+000 -4.239566704091e-001 3.405332305815e+000-1.050017655852e-001 1.462257102734e+000 -6.684487469283e-001 -4.027*********e-0013.0773********e-017 1.746386351950e-017 0.000000000000e+0001.577122080722e+000 6.399535133956e-001 3.468127872427e-001-5.701786649768e-002 4.014788054433e-001 -2.222476176311e-001 -6.317059236442e-0021.760039865880e-016 9.988285339980e-017 0.000000000000e+0000.000000000000e+000 -1.447846997770e+000 -1.415724007744e+000-2.806139044665e-001 -2.817910521892e-001 -4.611434881851e-0021.996629079956e-0018.804885435596e-017 4.996802050991e-017 0.000000000000e+0000.000000000000e+000 8.831633340975e-017 1.231641451542e+0001.619701003419e-001 1.962638275504e-001 5.350035621760e-001-1.509273424767e-001-7.728357669400e-018 -4.385868928757e-018 0.000000000000e+0000.000000000000e+000 -7.751835246841e-018 -1.279231078565e-017-7.753441914209e-001 -3.464514508821e-001 4.312226803504e-0011.234643696237e-001-5.603391361152e-017 -3.179943413314e-017 0.000000000000e+0000.000000000000e+000 -5.620413613517e-017 -9.274974944455e-0170.000000000000e+000 1.296312940612e+000 -4.288053318338e-0012.737334158165e-001-2.493361499851e-017 -1.414991023727e-017 0.000000000000e+0000.000000000000e+000 -2.500935953597e-017 -4.127119443934e-0170.000000000000e+000 0.000000000000e+000 -6.707396440648e-001-4.842320121884e-001-2.603055667460e-017 -1.477242832192e-017 0.000000000000e+0000.000000000000e+000 -2.610963355436e-017 -4.308689959101e-0170.000000000000e+000 0.000000000000e+000 -1.110223024625e-0167.168323926323e-002（4）RQ1.163074414164e+0002.632670934508e+000 -1.772796003272e+000-8.668899138521e-002 3.300503471047e-001 1.455162371214e+000 -9.730650448593e-001 -4.873031174655e-001 -7.756411630489e-001 -3.249201979113e-0011.836115060851e+000 1.144286420080e-001 -9.880381403133e-0015.589725694767e-001 4.694190067101e-002 -2.978478237007e-0011.617130577649e-002 6.936977702522e-001 1.367670571405e-0011.419099231519e-0025.598200524418e-016 -2.118520153533e+000 -1.876189745783e+000-5.407071940597e-001 1.171538359721e+000 -2.550323020223e+0001.691577936540e+000 1.229951613262e+000 1.387947777212e+0008.667502917242e-001-3.179059303049e-017 -5.261714527977e-017 -8.471995127808e-0014.382910468318e-001 -1.008632199185e+000 -7.959374261495e-0014.769258865577e-001 4.072683083890e-001 4.096390493527e-0013.363378940862e-001-3.514195999269e-016 3.391949386044e-017 -2.160938214545e-016 -1.432244342447e+000 -5.742284908055e-001 1.213151477723e+000 -3.457508625575e-001 -4.749853573124e-001 -3.176158274191e-001 -4.294507015032e-002-3.704735750656e-017 -1.0998********e-016 -1.953241363386e-0178.269089741494e-017 6.556779598004e-001 -9.275250974463e-0012.529079844053e-001 6.905949216976e-001 -2.374430675823e-002-2.429781119781e-001-6.420051822287e-017 3.865817713597e-017 -3.806718328740e-0172.129734126613e-017 7.853*********e-017 4.698400884876e-001-2.730776009527e-001 7.821296259798e-001 -9.580964936399e-0027.846239841323e-0021.478496020372e-016 -8.383952267131e-017 9.577540382744e-017-7.120338911078e-017 -1.220876056602e-016 -1.854471832043e-0161.287679058937e+000 -3.576058900348e-001 -4.116725408807e-0033.914268216423e-0014.477158378610e-017 -6.204017427568e-017 3.360322998507e-017-1.133882337819e-017 -2.759056827456e-017 -9.217582575819e-0172.214639994444e-016 -3.628760503545e-001 7.398980975354e-0017.241608309576e-0023.408891482791e-017 2.353495494255e-017 1.932283926688e-017-3.456910153081e-017 -2.015177337156e-017 -8.084422691634e-017 -5.839476980893e-017 5.551115123126e-017 -5.176670596524e-0024.958522909877e-002（5）特征值(6.360627875745e-001,0.000000000000e+000)(5.650488993501e-002,0.000000000000e+000)(9.355889078188e-001,0.000000000000e+000)(-9.805309562902e-001,1.139489127430e-001)(-9.805309562902e-001,-1.139489127430e-001)(1.577548557113e+000,0.000000000000e+000)(-1.484039822259e+000,0.000000000000e+000)(-2.323496210212e+000,8.930405177200e-001)(-2.323496210212e+000,-8.930405177200e-001)(3.383039617436e+000,0.000000000000e+000)（6）实特征值的特征向量4.745031936539e+0003.157868541750e+0001.729946912420e+001-1.980049331455e+000-3.187521973524e+0017.794009603201e+000-1.004255685845e+0011.670757770480e+0011.310524273054e+0011.000000000001e+000下一个实特征值对应的特征向量：-5.105003830625e+000-4.886319842360e+0009.505161576151e+000-6.788331788214e-001-9.604334110499e+000-3.0457********e+0001.574873885602e+001-7.395037126277e+000-7.109654943661e+0001.000000000001e+000下一个实特征值对应的特征向量：2.792418944529e+0001.598236841512e+000-5.207507040911e-001-1.667886451702e+000-1.225708535859e+0017.241214790777e+000-5.398214501433e+0002.841008912974e+001-1.216518754416e+0011.000000000001e+000下一个实特征值对应的特征向量：6.217350824581e-001-1.115111815226e-001-2.483773580804e+000-1.306860840421e+000-3.815605442533e+0008.117305509395e+000-1.239170883675e+000-6.800309586197e-0012.691900144840e+0001.000000000001e+000下一个实特征值对应的特征向量：-1.730784582112e+0012.408192534967e+0014.133852365119e-001-8.572044074538e+0009.287334657928e-002-7.832726042776e-002-6.374274025716e-001-3.403204760832e-001。
北航数值分析大作业 第二题 QR分解
《数值分析B》课计算实习第一题设计文档与源程序姓名：杨彦杰学号：SY10171341 算法的设计方案（1）运行平台操作系统：Windows XP；开发平台：VC6.0++；工程类型：文档视图类；工程名：Numanalysis；（2）开发描述首先新建类CMetrix，该类完成矩阵之间的相关运算，包括相乘、加减等，以主程序方便调用；题目的解算过程在视图类CNumanalysisView中实现，解算结果在视图界面中显示；（3）运行流程（4）运行界面2、全部源代码（1）类CMetrixMetrix.h文件：class CMetrix{public:double** MetrixMultiplyConst(double**A,int nRow,int nCol,double nConst);//矩阵乘常数double** MetrixMultiplyMetrix(double**A,double**mA,int nRow,int nCol);//矩阵相乘double** MetrixSubtractMetrix(double **A, double **subA, int nRow,int nCol);//矩阵减矩阵double VectorMultiplyVector(double*V,double*mulV,int nV);//向量点积double** VectorMultiplyVectortoMetrix(double*V,double*VT,int nV);//向量相乘为矩阵double* VectorSubtractVector(double*V,double*subV,int nV);//向量相减double* VectorMultiplyConst(double *V, int nV, double nConst);//向量乘常数double LengthofVector(double *V,int nV);//求向量的长度double* MetrixMultiplyVector(double**A,int nRow,int nCol,double*V,int nV);//矩阵与向量相乘double** AtoAT(double **A,int Row,int Col);//矩阵转置运算void FreeMem();CMetrix(int nRow,int nCol);uCMetrix();virtual ~CMetrix();double* vector; //过渡向量double** B; //过渡矩阵};Metrix.cpp文件：CMetrix::CMetrix(int nRow, int nCol){B = new double*[nRow];for (int i = 0;i < nCol;i++){B[i] = new double[nCol];}vector = new double[nRow];}CMetrix::~CMetrix(){delete vector;B = NULL;delete B;}double** CMetrix::AtoAT(double **A, int nRow, int nCol){for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[col][row] = A[row][col];}}return B;}double* CMetrix::MetrixMultiplyVector(double **A, int nRow, int nCol, double *V, int nV) {if (nCol != nV){AfxMessageBox("矩阵列数和向量维数不等,不能相乘!");return 0;}double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){sum += A[row][col]*V[col];}vector[row] = sum;sum = 0.0;}return vector;}double CMetrix::LengthofVector(double *V, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*V[col];}return length;}double* CMetrix::VectorMultiplyConst(double *V, int nV, double nConst){for (int col = 0;col < nV;col++){vector[col] = V[col]*nConst;}return vector;}double* CMetrix::VectorSubtractVector(double *V, double *subV, int nV){for (int col = 0;col < nV;col++){vector[col] = V[col]-subV[col];}return vector;}double** CMetrix::VectorMultiplyVectortoMetrix(double*V, double *VT, int nV){for (int row = 0;row < nV;row++){for (int col = 0;col < nV;col++){B[row][col] = V[row]*VT[col];}}return B;}double CMetrix::VectorMultiplyVector(double *V, double *mulV, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*mulV[col];}return length;}double** CMetrix::MetrixSubtractMetrix(double **A, double **subA, int nRow, int nCol) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]-subA[row][col];}}return B;}double** CMetrix::MetrixMultiplyMetrix(double **A, double **mA, int nRow, int nCol) {double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){for(int n = 0;n < nCol;n++){sum += A[row][n]*mA[n][col];}B[row][col] = sum;sum = 0.0;}}return B;}double** CMetrix::MetrixMultiplyConst(double **A, int nRow, int nCol, double nConst) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]*nConst;}}return B;}（2）类CNumanalysisViewNumanalysisview.hclass CNumanalysisView : public CEditView{…………public:double Sign(double x);void DisplayVector(double*V,int nV); // 显示向量数据void DisplayMetrix(double **A,int Row,int Col); //显示矩阵void DisplayText(CString str); //显示文本protected://{{AFX_MSG(CNumanalysisView)afx_msg void OnQRanalyze(); //运行主函数…………};Numanalysisview.cppvoid CNumanalysisView::OnQRanalyze(){//开辟空间int nRow = 10;int nCol = 10;CString str;CMetrix Metrix(nRow,nCol);double tempa = 0.0;double *V = new double[nCol]; //分配10*10矩阵空间double *ur = new double[nCol];double *pr = new double[nCol];double *qr = new double[nCol];double *wr = new double[nCol];double *tempV = new double[nCol];double **Ar = new double*[nRow];double **C = new double*[nRow];double **Cr = new double*[nRow];double **tempA = new double*[nRow];double **A = new double*[nRow];double **R = new double*[nRow];for (int col = 0;col < nRow;col++){A[col] = new double[nCol];Ar[col] = new double[nCol];C[col] = new double[nCol];Cr[col] = new double[nCol];tempA[col] = new double[nCol];R[col] = new double[nCol];}//矩阵A求解for (int i = 0;i < nRow;i++){for (int j = 0;j < nCol;j++){if(i == j)A[i][j] = 1.5*cos((i+1.0)+1.2*(j+1.0));elseA[i][j] = sin(0.5*(i+1.0)+0.2*(j+1.0));}}//--------------------拟上三角化-------------------------// double dr = 0.0,cr = 0.0,hr = 0.0,tr = 0.0;for (int r = 0;r < nCol - 2;r++){dr = 0.0;for (i = r+1;i < nCol;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+2;i < nCol;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r+1][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r+1][r])*dr;hr = cr*cr - cr*A[r+1][r]; //hrstr.Format("dr = %.6e, cr = %.6e, hr = %.6e",dr,cr,hr);for (int row = 0;row < nRow;row++) //ur{if (row < r+1)ur[row] = 0.0;else if (row == r+1)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,nRow,nCol,ur,nCol); //pr memcpy(pr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(pr,nCol,1.0/hr);memcpy(pr,tempV,nCol*8);tempV = Metrix.MetrixMultiplyVector(A,nRow,nCol,ur,nCol); //qr memcpy(qr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(qr,nCol,1.0/hr);memcpy(qr,tempV,nCol*8);tr = Metrix.VectorMultiplyVector(pr,ur,nCol)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,nCol,tr); //wr memcpy(wr,tempV,nCol*8);tempV = Metrix.VectorSubtractVector(qr,wr,nCol);memcpy(wr,tempV,nCol*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,nCol); //Arfor (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)A[row][col] = tempA[row][col];}tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}}DisplayText("矩阵A拟上三角化后所得的矩阵为：");DisplayMetrix(A,nRow,nCol);for (int row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)R[row][col] = A[row][col];}// -------------------------------------------------////--------------------带双步位移的QR分解-------------------------// int m = nCol;struct EigenVal //定义特征值结构，实数和虚数{double Realnum;double Imagnum;};EigenVal *eigenvalue = new EigenVal[m];EigenVal tmpEigen1,tmpEigen2;double b = 0.0,c = 0.0,delta = 0.0,s = 0.0,t = 0.0;double *vr = new double[m];for (int k = 1;k < 100; k++){//m代表矩阵阶数，判断式中直接用，运算中需要-1while (m > 1 && fabs(A[m-1][m-2]) <= IPSLEN)//第三步和第四步{eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;m = m - 1;}if (m == 1){eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;DisplayText("已求出A的全部特征值：");break;}b = -(A[m-2][m-2]+A[m-1][m-1]); //第五步求一元二次方程式的根s1,s2c = A[m-2][m-2]*A[m-1][m-1]-A[m-2][m-1]*A[m-1][m-2];delta =b*b - 4*c;if (delta >= 0.0){tmpEigen1.Realnum = (-b-sqrt(delta))/2;tmpEigen1.Imagnum = 0.0;tmpEigen2.Realnum = (-b+sqrt(delta))/2;tmpEigen2.Imagnum = 0.0;}else{tmpEigen1.Realnum = -b/2;tmpEigen1.Imagnum = -sqrt(fabs(delta))/2 ;tmpEigen2.Realnum = -b/2;tmpEigen2.Imagnum = sqrt(fabs(delta))/2;}if (m == 2) //第六步 m=2时结束运算{eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;DisplayText("已求出A的全部特征值：");break;}else //第七步 m > 1{if (fabs(A[m-2][m-3]) <= IPSLEN){eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;m = m - 2;continue;}}for (int row = 0;row < m;row++) //Mk求之前需要把A付给C{for (int col = 0;col < m;col++)C[row][col] = A[row][col];}double **I = new double*[m]; //第九步求Mk和Mk的QR分解for (int i = 0;i < m;i++) //求单位矩阵I,分配m*m矩阵空间{I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}s = A[m-2][m-2]+A[m-1][m-1];t = A[m-2][m-2]*A[m-1][m-1] - A[m-2][m-1]*A[m-2][m-1];tempA = Metrix.MetrixMultiplyMetrix(A,A,m,m);//A*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixMultiplyConst(A,m,m,s);//s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(Ar,A,m,m);//A*A-s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col]; }tempA = Metrix.MetrixMultiplyConst(I,m,m,-1*t);//-t*Ifor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m);//A*A - s*A + r*I for (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}delete I;//Mk的QR分解for (int r = 0;r < m - 1;r++){dr = 0.0;for (i = r;i < m;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+1;i < m;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r][r])*dr;hr = cr*cr - cr*A[r][r]; //hrfor (int row = 0;row < m;row++) //ur{if (row < r)ur[row] = 0.0;else if (row == r)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,m,m); //Btfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,m,m,ur,m); //Bt*ur memcpy(vr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(vr,m,1.0/hr); //vr = Bt*ur/hr memcpy(vr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(ur,vr,m);//Ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m); //Br-ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}tempA = Metrix.AtoAT(C,m,m); //Ctfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempV = Metrix.MetrixMultiplyVector(Cr,m,m,ur,m); //pr memcpy(pr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(pr,m,1.0/hr);memcpy(pr,tempV,m*8);tempV = Metrix.MetrixMultiplyVector(C,m,m,ur,m); //qr memcpy(qr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(qr,m,1.0/hr);memcpy(qr,tempV,m*8);tr = Metrix.VectorMultiplyVector(pr,ur,m)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,m,tr); //wr memcpy(wr,tempV,m*8);tempV = Metrix.VectorSubtractVector(qr,wr,m);memcpy(wr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,m);//Cr+1for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)C[row][col] = tempA[row][col]; }tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++){C[row][col] = tempA[row][col];if (fabs(C[row][col]) < IPSLEN){C[row][col] = 0.0;}}}}str.Format("矩阵A%d QR分解结束后所得到的矩阵为：",m);//计算结果输出DisplayText(str);DisplayMetrix(A,m,m);for (row = 0;row < m;row++) //Mk的QR分解后需要把C付给A{for (col = 0;col < m;col++)A[row][col] = C[row][col];}str.Format("迭代完成后的矩阵A%d = ",k);DisplayText(str);DisplayMetrix(A,m,m);}DisplayText("矩阵A的全体特征值如下: ");for (i = 0;i<nCol;i++){str.Format("%.6e + j%.6e",eigenvalue[i].Realnum,eigenvalue[i].Imagnum);DisplayText(str);}// -------------------------------------------------//求实特征值的特征向量，在拟上三角矩阵基础上直接求解即可////(A-egiI)X = 0.0;m = nRow;for (row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)A[row][col] = R[row][col];}double **I = new double*[m]; //求单位矩阵I,分配m*m矩阵空间double sum = 0.0;for (i = 0;i < m;i++){I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}for (i = 0;i < nRow;i++){if (eigenvalue[i].Imagnum != 0.0){str.Format("特征值%.6e+j%.6e为虚数，不需要求特征向量。
北航数值分析大作业 第二题 QR分解
数值分析第二题 梁进明SY0906529算法设计方案。
一．矩阵的QR 分解。
把矩阵A 分解为一个正交矩阵Q 与一个上三角矩阵R 的乘积，称为矩阵A 的正三角分解，简称QR 分解。
QR 分解的算法如下：记1A A =，并记[]rij n n Ar a ⨯=，令1Q I =（n 阶单位矩阵） 对于r=1,2,…,n-1执行(1) 若(1,2,...,)rir a i r r n =++全为零，则令1r r Q Q +=，1r r A A +=转(5);否则转(2)(2) 计算2r d =()sgn()r r rr r c a d =-(若()0r rr a =,则取r r c d =) 2()r r r r rrh c c a =- (3) 令()()()1,(0,...,0,,,...,)r r r T nr rr r r r nr u a c a a R +=-∈ (4) 计算11//r r rTr r r r rT r rr rT r r r rQ u Q Q u h p A u h A A u p ωω++==-==-(5) 继续当此算法执行完后就得到正交矩阵n Q Q =和上三角矩阵n R A =且有A QR =。
二．矩阵的 拟上三角化。
对实矩阵A 的拟上三角化具体算法如下：记(1)AA =，并记()r A 的第r 列到第n 列的元素为(1,2,...,;,1,...,)rij a i n j r r n ==+。
对于1,2,...,2r n =-执行（1） 若()(2,3,...,)r ir a i r r n =++全为零，则令(1)()r r AA +=,转(5);否则转(2)。
（2） 计算2()()1,1,2()1,sgn()(0,)r r r r r r r r r r r r r r r r rd c a d a c d h c c a +++==-===-若则取 （3） 令()()()1,2,(0,...,0,,,...,)r r r T nr r r r r r nr u a c a a R ++=-∈ （4） 计算()()(1)()///r T r r r r r r rTr r r rr r r rr r T T r r r rp A u h q A u h t p u h q t u A A u u p ωω+====-=--（5） 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵(1)n A -。
数值分析矩阵的正交分解(QR分解)
§10 矩阵的正交分解(QR 分解)设nm RA ⨯∈，则存在初等反射阵s H H 1使得)1(2+=s s A A H H (上梯形)[]nmn m m n n a a a a a a a a a a a a A ,,,21212222111211=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=(按列分块) (1)(1)第1步：当01=a 时，取I H =1这一步不需约化，不妨设01≠a ，于是有初等反射阵1H 使1111e a H σ-=，其中Tu u I H 11111--=β。
于是],,,[21)1(1n Ha Ha Ha A H =⎥⎦⎤⎢⎣⎡-=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-=222)2(121)2()2(2)2(2)2(22)2(2)2(121000D c B a a a a a a a mn m n n σσ)2(A =其中)2()1(21)2(2)2(222,),,(-⨯--∈∈=n m m T m R D Ra a c (2)第k 步：设已完成对A 上述第1步~第k-1步约化，即存在初等反射阵11,,-k H H 使)(121k k A A H H H =-其中 ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=-)()()()(12)2(1)2(1)2(121)(k mn k mkk knk kkk n k k a a a a a a a Aσσσ ⎥⎦⎤⎢⎣⎡=k k k k k D c B r R 0其中)()1(1)()(,],,[k n k m k k n T k m k k kk k R D Ra c c -⨯+-=+-∈∈= ，为 EMBED Equation.3阶上三角阵。
如果0=k c ，这一步不需约化，取I H k =。
不妨设0≠k c ，于是存在初等反射阵kH '使 1e c H k k kσ-='计算kH '的公式： T k k k ku u I H ''-='-1β ⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧+='=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡='=+=∑)(21)()()(22)(,)(,1)(2)()(k kk k k k k k k m k k k k kk k m k i k ik k kk k a u a a a u a a sign σσβσ ………………(2) 令mm k k m k k k R H I H ⨯-+--∈⎥⎦⎤⎢⎣⎡'=111第k 步约化：)1(1)(+==k k k k A A H H A H⎥⎦⎤⎢⎣⎡''⎥⎦⎤⎢⎣⎡'=-k k k k k k kk k D H c H B r R H I 01 )1(121+-=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡'----=k k k k k k k AD H B r σσσσ方框内为第k 步约化需要计算的部分，其中)1(+k A 左上角子阵，1+k R 为k阶上三角阵，这样就使A 三角化过程前进了一步。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北航2009级研究生《数值分析B》计算实习题目（第二题：QR分解）设计文档与源程序2009年11月16日打印内容1 算法的设计方案（1）运行平台（2）算法设计2 全部源代码3 输出结果：3.1 矩阵A经过拟上三角化后所得的矩阵A(n-1)；3.2 对矩阵A(n-1)进行QR分解方法后所得的矩阵；3.3 矩阵A的全部特征值λｉ；3.4 A的相对应于实特征值的特征向量。
1 算法的设计方案1.1 运行与开发平台操作系统：Windows 7；开发平台：VC++ 6.0；工程类型：Win32 Console Application；工程名：QR_EigenValue；1.2 算法设计设计思想：本题目要求采用带双步位移的QR分解法求矩阵A={a ij}10×10的全部特征值。
首先由已知公式计算出A各元素值；然后对矩阵A进行拟上三角化得到矩阵A(n-1)；然后对得出的A(n-1)进行QR分解，得到一个分块上三角矩阵以及矩阵的各个特征值λi；最后，采用高斯消去法，可以求得对应于各实特征值的特征向量。
具体算法如下：（精度eps=le-12，最大迭代次数L=50，n=10）(一)、计算矩阵A对i,j=0到10，执行：if i=j,a[i][j]=1.5*cos*(i+1.2j);else a[i][j]=sin*(0.5*i+0.2j);(二)、矩阵拟上三角化函数矩阵拟上三角化函数为：void Hessenberg_A(double a[][n],double a_n1[][n])。
定义矩阵ar = a，对于r=1，2，……，n-2执行以下算法：第一步：判断a[i-1][r-1]是否全为零（i=r+2,……,n）,如果是，则令a[i-1][r-1] = a[i-1][r-1]，并执行循环；否则，跳出第一步的循环。
第二步：对i=r+1，……，n计算d_r = sqrt(ar[i-1][r-1] * ar[i-1][r-1])cr=--sgn(d_r)h_r = c_r*c_r - c_r*ar[r][r-1]第三步：对 i=r+1，……，n，执行for(i=r+1;i<=n;i++)ur[i-1] = ar[i-1][r-1];ur[r] = ur[r] - c_r;第四步：计算pr=A’*ur/hrqr=A*ur/hrtr=pr’*ur/hrwr=qr-tr*urA=A-wr*ur-ur*pr’(三)、对A(n-1)做带双步位移的QR分解定义的带双步位移的QR分解函数为int Db_QR_Method(double Ak[n][n], double lamda[n][2], FILE *pFile) 记Ak[i][j] = A(n-1)=a_n1[i][j]，令k=0,m=n;(2) if | A[m][m-1]<=eps|,则得到A的一个特征值lamda[m][0]=A[m][m]，置m=m-1,转(3);否则转(4);(3) if m=1,则得到A的一个特征值lamda[m][0]=A[1][1]，转(10);if m=0，则直接转(10)；if m>1,则转(2);(4) 调用了解二次方程函数，用来求二阶子阵的两个特征值，即计算二次方程的两个根，x*x- s*x+ t=0 （*）其中s = A[m-1][m-1]+A[m][m]; t]= A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]；(5) if m=2,则得到A的两个特征值，转(10)；否则转(7)；(6) if |A[m-1][m-2]) |<=eps,则得到A的两个特征根，置m=m-2，转(3);否则转(7)；(7) if k=L,跳出程序，未能得到结果；否则转(8)；(8) 对i,j=1,…,m，计算s= A[m-1][m-1]+A[m][m];t= A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1];Mk[i][j] = Mk[i][j]-s*A[i][j]+t*I[i][j]Mk=Qk*Rk（对Mk作QR分解，可调用QR分解函数）Ak+1=Qk’*Ak*Qk(9) 置k=k+1，转(2)；(10) 得到所有的特征值，停止计算。
Mk的QR分解算法：对r=1到m-1执行(1) if |B[i][r]|>=eps，转(2);(2) 计算dr=sqrt(A[i][r]* A[i][r])cr=--sgn(A[r+1][r])hr = cr*(cr-A[r+1][r]);(3) 对i=r+1到m执行ur[i] = A[i][r]; ur[r+1] = ur[r+1]-cr；(4)计算vr=Br’*ur/hrB r+1=Br-ur*vr’pr=Crur/hrqr=Crr/hrtr=pr’*ur/hrwr=qr-tr*urC r+1=Cr-wr*ur-ur*pr’求二次方程的根的算法：对方程（*）(1)计算det = s*s-4*t;sqrt_det = sqrt(fabs(det));(2) if det>=0，则s1[0] = (s+ sqrt_det)/2.0，s2[0] = (s- sqrt_det)/2.0;否则转(3);(3) s1[0]=s/2; s1[1]= sqrt_det /2; s2[0]=s/2; s2[1]=- sqrt_det /2;(四)、用列主元Gauss消去法求解实特征根对应的特征向量(1) if lamda[m][1]=0，转(2);(2) 调用GaussOptimal子程序；(3) printf “When the eigenvalue is %4.12e, the corresponding eigenvector is\n",Lamda[i][0]”GaussOptimal算法记a1=aij，令b[N]={0}。
对i=1到10执行A1[i][i] = A1[i][i]-Lamda[p][0]消元过程：对于k=1到n-1执行(1) 选行号ik,使得|aik,k|=max|aik|;(2) 交换akj与aik,j(j=1到n)；(3) 对于i=k+1到n计算mik = A1[i][k]/A1[k][k]A1[i][j] = A1[i][j]-mik*A1[k][j]( j=k+1到n)回代过程(1) if fabs(A1[n][n])<=1e-9则if fabs(A1[n-1][n-1])<=1e-9，令xm=(0,0,…,0,1), xm-1=(*,*,…,*,1,0);否则转(3)；(2) 对k=n-2到1执行1）置temp=0.0;2）对j=k+1到n执行： temp = temp+A1[k][j]*x[p][j];3）x[p-1][k] = -temp/A1[k][k];(3) l令x[p][n]=1.0;对k=n-1到1执行1）temp=0.0;2）对j=k+1到n执行temp = temp+A1[k][j]*x[p][j];3）x[p][k] = -temp/A1[k][k];(4) 特征向量归一化.2、全部源代码#include <stdio.h>#include <math.h>#define n 10#define L 50#define eps 1e-12void Matrix_A(double a[][n]);void Hessenberg_A(double a[][n],double a_n1[][n]);void Solve_Quadratic(double s,double t,double s1[2],double s2[2]);void QR_Mk(int m, double Mk[][n], double Ak[][n]);int Db_QR_Method(double A_n1[][n], double lamda[n][2], FILE *pFile);void ReVagin_Vector(double A[n][n], double x[][n], double lamda[n][2], FILE *pFile);void GaussOptimal(double a1[][n],double x[][n],double Lamda[][2],int p);void main(){int i,j,flag = 0;double a[n][n],a_n1[n][n],lamda[n][2];double x[n][n];FILE *pFile;pFile=fopen("QR_EigenValue.txt","wt"); /*以只写方式打开 Power_EigenValue 文本文档*/ printf("**********************************************\n");printf("用带双步位移的QR分解法求解矩阵特征值及特征向量\n"); /*输出标题*/printf("**********************************************\n");fprintf(pFile,"**********************************************\n");fprintf(pFile,"用带双步位移的QR分解法求解矩阵特征值及特征向量\n"); /*输出标题*/ fprintf(pFile,"**********************************************\n");//输出矩阵Aprintf("\n已知矩阵A：\n");fprintf(pFile,"\n已知矩阵A：\n"); /*输出标题*/Matrix_A(a);for(i=0;i<n;i++){for(j=0;j<n;j++){printf(" %4.12e ",a[i][j]);fprintf(pFile," %4.12e ",a[i][j]);}printf("\n");fprintf(pFile,"\n");}//矩阵A的拟上三角化printf("\n拟上三角化后的矩阵A(n-1)：\n");fprintf(pFile,"\n拟上三角化后的矩阵A(n-1)：\n"); /*输出标题*/Hessenberg_A(a,a_n1);for(i=0;i<n;i++){for(j=0;j<n;j++){a[i][j] = a_n1[i][j];printf(" %4.12e ",a_n1[i][j]);fprintf(pFile," %4.12e ",a_n1[i][j]);}printf("\n");fprintf(pFile,"\n");}//带双步位移的QR分解函数调用flag = Db_QR_Method(a_n1, lamda, pFile);//对矩阵A(n-1)进行QR分解后的矩阵输出printf("\n输出对矩阵A(n-1)进行QR分解后所得的矩阵：\n");fprintf(pFile,"\n输出对矩阵A(n-1)进行QR分解后所得的矩阵：\n");for(i=0;i<n;i++){for(j=0;j<n;j++){printf(" %4.12e ",a_n1[i][j]);fprintf(pFile," %4.12e ",a_n1[i][j]);}printf("\n");fprintf(pFile,"\n");}//输出矩阵A的全部特征值printf("\n输出特征值\n");fprintf(pFile,"\n输出特征值\n");for(i=0;i<n;i++){printf("第%d 个特征值λ(%d) = %4.12e + %4.12e*i\n",i+1,i+1,lamda[i][0],lamda[i][1]);fprintf(pFile,"第%d 个特征值λ(%d) = %4.12e + %4.12e*i\n",i+1,i+1,lamda[i][0],lamda[i][1]);}//求A的相对于实特征值的特征向量printf("\n\n");fprintf(pFile,"\n\n");ReVagin_Vector(a, x, lamda, pFile);fclose(pFile);return;}//原始矩阵A的计算void Matrix_A(double a[][n]){int i,j;for(i=1;i<=n;i++){for(j=1;j<=n;j++){if(i == j)a[i-1][j-1] = 1.5*cos(i + 1.2*j);elsea[i-1][j-1] = sin(0.5*i + 0.2*j);}}return;}//矩阵的拟上三角化void Hessenberg_A(double a[][n],double a_n1[][n]){int i,j,r;double d_r,c_r,h_r,temp,t_r;double ur[n],pr[n],qr[n],wr[n];double ar[n][n];//准备：令Ar=Afor(i=0;i<n;i++)for(j=0;j<n;j++)ar[i][j] = a[i][j];//进入r循环for(r=1;r<=n-2;r++){d_r=0.0; c_r=0.0; h_r=0.0; temp=0.0; t_r=0.0;for(i=1;i<=n;i++)ur[i-1]=0.0; pr[i-1]=0.0; qr[i-1]=0.0; wr[i-1]=0.0;//step 1:for(i=r+2;i<=n;i++){if(fabs(a[i-1][r-1]) <= eps)a[i-1][r-1] = a[i-1][r-1];elsebreak;}//step 2:for(i=r+1;i<=n;i++)temp += ar[i-1][r-1] * ar[i-1][r-1];d_r = sqrt(temp);if(ar[r][r-1]>eps)c_r = -d_r;elsec_r = d_r;h_r = c_r*c_r - c_r*ar[r][r-1];//step 3:for(i=r+1;i<=n;i++)ur[i-1] = ar[i-1][r-1];ur[r] = ur[r] - c_r;//step 4:for(i=1;i<=n;i++){for(j=1;j<=n;j++){pr[i-1] += ar[j-1][i-1] * ur[j-1];qr[i-1] += ar[i-1][j-1] * ur[j-1];}pr[i-1] = pr[i-1]/h_r;qr[i-1] = qr[i-1] /h_r;}temp = 0;temp += pr[i-1]*ur[i-1];t_r = temp/h_r;for(i=1;i<=n;i++)wr[i-1] = qr[i-1] - t_r*ur[i-1];for(i=1;i<=n;i++){for(j=1;j<=n;j++){if((j==r)&&(i>j+1))ar[i-1][j-1] = 0.0;elsear[i-1][j-1] = ar[i-1][j-1] - wr[i-1]*ur[j-1] - ur[i-1]*pr[j-1];}}//得出矩阵A(n-1)各元素值for(i=0;i<n;i++)for(j=0;j<n;j++)a_n1[i][j] = ar[i][j];}}//Mk的QR分解与A(k+1)的计算void QR_Mk(int m, double Mk[][n], double Ak[][n]){double Br[n][n],Cr[n][n];double ur[n],vr[n],pr[n],qr[n],wr[n];double dr,cr,hr,tr,temp;int i,j,r;//初始化for(i=0;i<n;i++){for(j=0;j<n;j++){Br[i][j] = Mk[i][j];Cr[i][j] = Ak[i][j];}}//进入r循环for(r=1;r<=m-1;r++){dr=0;cr=0;hr=0;tr=0;temp=0;for(i=0;i<n;i++)ur[i]=0; vr[i]=0; pr[i]=0; qr[i]=0; wr[i]=0;//step 1:{if(fabs(Br[i-1][r-1]) <= eps)Br[i-1][r-1] = Br[i-1][r-1];elsebreak;}//step 2:for(i=r;i<=m;i++)temp += Br[i-1][r-1] * Br[i-1][r-1]; dr = sqrt(temp);if(Br[r-1][r-1]>eps)cr = -dr;elsecr = dr;hr = cr*cr - cr*Br[r-1][r-1];//step 3:for(i=r;i<=m;i++)ur[i-1] = Br[i-1][r-1];ur[r-1] = ur[r-1] - cr;//step 4:for(i=1;i<=m;i++){for(j=1;j<=m;j++){vr[i-1] += Br[j-1][i-1] * ur[j-1];pr[i-1] += Cr[j-1][i-1] * ur[j-1];qr[i-1] += Cr[i-1][j-1] * ur[j-1];}vr[i-1] = vr[i-1]/hr;pr[i-1] = pr[i-1]/hr;qr[i-1] = qr[i-1]/hr;}temp = 0;for(i=1;i<=m;i++)temp += pr[i-1]*ur[i-1];tr = temp/hr;for(i=1;i<=m;i++){wr[i-1] = qr[i-1] - tr*ur[i-1];}for(i=1;i<=m;i++){for(j=1;j<=m;j++){Br[i-1][j-1] = Br[i-1][j-1] - ur[i-1]*vr[j-1];Cr[i-1][j-1] = Cr[i-1][j-1] - wr[i-1]*ur[j-1] - ur[i-1]*pr[j-1];if(fabs(Cr[i-1][j-1])<=eps)Cr[i-1][j-1] = 0;}}//得出矩阵 A(k+1)for(i=0;i<m;i++)for(j=0;j<m;j++)Ak[i][j] = Cr[i][j];}}//对二阶阵Dk，求解二次方程,定义实部和虚部，在函数调用的时候得出的特征值要注意变换void Solve_Quadratic(double s,double t,double s1[2],double s2[2]){double det,sqrt_det;det = s*s - 4*t;sqrt_det = sqrt(fabs(det));if(det>=0){s1[0] = (s+sqrt_det)/2;//实根s2[0] = (s-sqrt_det)/2;//实根s1[1] = 0;s2[1] = 0;}else{s1[0] = s/2;//实根s1[1] = sqrt_det/2;//虚部s2[0] = s/2;//实根s2[1] = -sqrt_det/2;//虚部}return;}//带双步位移的QR方法求实矩阵A的全部特征值int Db_QR_Method(double Ak[n][n], double lamda[n][2], FILE *pFile){double I[n][n]={0}, Ak_2[n][n], s1[2],s2[2],Mk[n][n];double s,t;int flag = 1;int i,j,k,r;int m;//单位矩阵Ifor(i=0;i<n;i++)I[i][i] = 1;for(i=0;i<n;i++)for(j=0;j<2;j++)lamda[i][j] = 0;//step 1:调用拟上三角矩阵A(n-1)//step 2:for(i=0;i<n;i++)for(j=0;j<n;j++)Ak_2[i][j] = 0;k=1; m=n;while(1){step3: if( fabs(Ak[m-1][m-2]) <=eps){lamda[m-1][0] = Ak[m-1][m-1];m = m-1;step4: if(m==1){lamda[m-1][0] = Ak[m-1][m-1];printf("\n最大迭代次数k：%d\n",k);fprintf(pFile,"\n最大迭代次数k：%d\n",k);return(flag);}else if(m==0){printf("\n最大迭代次数k：%d\n",k);fprintf(pFile,"\n最大迭代次数k：%d\n",k);return(flag);}elsegoto step3;}//step5:s = Ak[m-2][m-2] + Ak[m-1][m-1];t = Ak[m-2][m-2]*Ak[m-1][m-1] - Ak[m-2][m-1]*Ak[m-1][m-2];//调用解二次方程函数Solve_Quadratic(s,t,s1,s2);//step6if(m ==2){lamda[m-2][0] = s1[0];lamda[m-2][1] = s1[1];lamda[m-1][0] = s2[0];lamda[m-1][1] = s2[1];printf("\n最大迭代次数k：%d\n",k);fprintf(pFile,"\n最大迭代次数k：%d\n",k);return(flag);}//step7if(fabs(Ak[m-2][m-3])<=eps){lamda[m-2][0] = s1[0];lamda[m-2][1] = s1[1];lamda[m-1][0] = s2[0];lamda[m-1][1] = s2[1];m = m-2;goto step4;}//step8if(k==L){printf("\n计算终止，未得到A的全部特征值！\n");fprintf(pFile,"\n计算终止，未得到A的全部特征值！\n");return(flag);}//step9s = Ak[m-2][m-2]+Ak[m-1][m-1];t = Ak[m-2][m-2]*Ak[m-1][m-1] - Ak[m-1][m-2]*Ak[m-2][m-1];//计算Ak*Akfor(i=0;i<m;i++)for(j=0;j<m;j++)Mk[i][j] = 0;for(i=0;i<m;i++)for(j=0;j<m;j++)for(r=0;r<m;r++)Mk[i][j] += Ak[i][r]*Ak[r][j];for(i=0;i<m;i++)for(j=0;j<m;j++)Mk[i][j] = Mk[i][j] - s*Ak[i][j] + t*I[i][j];//调用函数，对Mk进行QR分解，计算A(k+1)QR_Mk(m, Mk, Ak);//step10k = k+1;}printf("迭代次数：%d",k);}//高斯消去法子程序void GaussOptimal(double a1[][n],double x[][n],double Lamda[][2],int p){int i,j,k,maxindex;double max,mik,temp=0;double A1[n][n]={0},b[n]={0}; //解齐次方程//构造线性方程组for(i=1;i<=n;i++)for(j=1;j<=n;j++)A1[i-1][j-1] = a1[i-1][j-1];for(i=1;i<=n;i++)A1[i-1][i-1] = A1[i-1][i-1]-Lamda[p][0];//解方程for(k=1;k<=n-1;k++){//消元过程max=A1[k-1][k-1];maxindex=k-1; //最大元素下标for(i=k+1;i<=n;i++)if(fabs(A1[i-1][k-1])>fabs(max)) //找出第k列主元素{maxindex=i-1; //交换行号max=A1[i-1][k-1]; //交换元素}if(maxindex>k-1){for(j=k;j<=n;j++) //作第k行与列主元素所在行元素互换初等变换{temp = A1[k-1][j-1];A1[k-1][j-1] = A1[maxindex][j-1];A1[maxindex][j-1] = temp;}}for(i=k+1;i<=n;i++){mik = A1[i-1][k-1]/A1[k-1][k-1];for(j=k+1;j<=n;j++)A1[i-1][j-1] = A1[i-1][j-1]-mik*A1[k-1][j-1];}}//回代过程if(fabs(A1[n-1][n-1])<=1e-9)//RA(n-1)=n-1.注意要选取合适精度，否则会误判根{if(fabs(A1[n-1][n-2])<=1e-9){//二重根x[p][n-1]=1.0; x[p][n-2]=0.0; /*令xm=(0,0,…,0,1)*/for(k=n-1;k>=1;k--){for(j=k+1;j<=n;j++)temp = temp+A1[k-1][j-1]*x[p][j-1];x[p][k-1] = -temp/A1[k-1][k-1]; /*b[N]=0*/}x[p-1][n-1]=0.0; x[p-1][n-2]=1.0; /*令xm-1=(*,*,…,*,1,0)*/for(k=n-2;k>=1;k--){temp=0.0;for(j=k+1;j<=n;j++)temp = temp+A1[k-1][j-1]*x[p][j-1];x[p-1][k-1] = -temp/A1[k-1][k-1];}}else{//单根x[p][n-1]=1.0;for(k=n-1;k>=1;k--){temp=0.0;for(j=k+1;j<=n;j++)temp = temp+A1[k-1][j-1]*x[p][j-1];x[p][k-1] = -temp/A1[k-1][k-1];}}}//特征向量归一化temp=0.0;for(i=1;i<=n;i++)temp = temp+pow(x[p][i-1],2);temp = sqrt(temp);for(i=1;i<=n;i++)x[p][i-1]=x[p][i-1]/temp;}//求实根特征值对应的特征向量void ReVagin_Vector(double A[n][n], double x[][n], double lamda[n][2], FILE *pFile) {int i,j;for(i=1;i<=n;i++){if(lamda[i-1][1] == 0){GaussOptimal(A,x,lamda,i-1);printf("\n当实根特征值为λ = %4.12e 时，对应的特征向量为：\n",lamda[i-1][0]);fprintf(pFile,"\n当实根特征值为λ= %4.12e 时，对应的特征向量为：\n",lamda[i-1][0]);for(j=0;j<n;j++){printf(" %4.12e ",x[i-1][j]);fprintf(pFile," %4.12e ",x[i-1][j]);}printf("%\n");fprintf(pFile,"\n");}}}3 输出结果：**********************************************用带双步位移的QR分解法求解矩阵特征值及特征向量**********************************************已知矩阵A：-8.827516758830e-001 7.833269096275e-001 8.912073600614e-001 9.635581854172e-001 9.974949866041e-001 9.916648104525e-001 9.463000876874e-001 8.632093666489e-001 7.457052121767e-001 5.984721441040e-0019.320390859672e-001 -4.609993049676e-001 9.995736030415e-001 9.738476308782e-001 9.092974268257e-001 8.084964038196e-001 6.754631805512e-001 5.155013718215e-001 3.349881501559e-001 1.411200080599e-0019.916648104525e-001 9.463000876874e-001 1.425348887938e+000 7.457052121767e-001 5.984721441040e-001 4.273798802338e-001 2.392493292140e-001 4.158066243329e-002 -1.577456941432e-001 -3.507832276896e-0018.084964038196e-001 6.754631805512e-001 5.155013718215e-001 -1.216639521092e+000 1.411200080599e-001 -5.837414342758e-002 -2.555411020268e-001 -4.425204432949e-001 -6.118578909427e-001 -7.568024953079e-0014.273798802338e-001 2.392493292140e-001 4.158066243329e-002 -1.577456941432e-001 6.638546982076e-003 -5.298361409085e-001 -6.877661591840e-001 -8.182771110644e-001 -9.161659367495e-001 -9.775301176651e-001-5.837414342758e-002 -2.555411020268e-001 -4.425204432949e-001 -6.118578909427e-001 -7.568024953079e-001 1.208825936461e+000 -9.516020738895e-001 -9.936910036335e-001 -9.961646088358e-001 -9.589242746631e-001-5.298361409085e-001 -6.877661591840e-001 -8.182771110644e-001 -9.161659367495e-001 -9.775301176651e-001 -9.999232575641e-001 -1.429429375331e+000 -9.258146823277e-001-8.322674422239e-001 -7.0554********e-001-8.715757724136e-001 -9.516020738895e-001 -9.936910036335e-001 -9.961646088358e-001 -9.589242746631e-001 -8.834546557202e-001 -7.727644875560e-001 4.736156323789e-001 -4.646021794138e-001 -2.794154981989e-001-9.999232575641e-001 -9.824526126243e-001 -9.258146823277e-001 -8.322674422239e-001 -7.0554********e-001 -5.506855425976e-001 -3.738766648302e-001 -1.821625042721e-001 8.719827177217e-001 2.151199880878e-001-8.834546557202e-001 -7.727644875560e-001 -6.312666378723e-001 -4.646021794138e-001 -2.794154981989e-001 -8.308940281750e-002 1.165492048505e-001 3.115413635134e-001 4.941133511386e-001 -1.499941239592e+0003.1 矩阵A经过拟上三角化后所得的矩阵A(n-1)拟上三角化后的矩阵A(n-1)：-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000 -7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000 -8.782436382927e-002 -9.255889387184e-001 6.032599440534e-001 1.518860956469e-001-2.347878362416e+000 2.372370104937e+000 1.819290822208e+000 3.237804101546e-001 2.205798440320e-001 2.102692662546e+000 1.816138086098e-001 1.278839089990e+000 -6.380578124405e-001 -4.154075603804e-0010.000000000000e+000 1.728274599968e+000 -1.171467642785e+000 -1.243839262699e+000 -6.399758341743e-001 -2.002833079037e+000 2.924947206124e-001 -6.412830068395e-001 9.783997621285e-002 2.557763574160e-0010.000000000000e+000 0.000000000000e+000 -1.291669534130e+000 -1.111603513396e+0001.171346824096e+000 -1.307356030021e+000 1.803699177750e-001 -4.246385358369e-001 7.988955239304e-002 1.608819928069e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 1.560126298527e+000 8.125049397515e-001 4.421756832923e-001 -3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-0020.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -7.707773755194e-001 -1.583051425742e+000 -3.042843176799e-001 2.528712446035e-001 -6.709925401449e-001 2.544619929082e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+0000.000000000000e+000 -7.463453456938e-001 -2.708365157019e-002 -9.486521893682e-0011.195871081495e-001 1.929265617952e-0020.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -7.701801374364e-001 -4.697623990618e-001 4.988259468008e-001 1.137691603776e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+0000.000000000000e+000 0.000000000000e+000 0.000000000000e+000 7.013167092107e-0011.582180688475e-001 3.862594614233e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 4.843807602783e-001 3.992777995177e-001最大迭代次数k：153.2 对矩阵A(n-1)进行QR分解方法后所得的矩阵输出对矩阵A(n-1)进行QR分解后所得的矩阵：3.383039617436e+000 8.948776654951e-001 -8.956760706755e-001 8.465153677933e-002 2.612677876847e-001 1.610337179819e+000 -1.022*********e+000 9.371886210517e-002 -1.002578700816e+000 -4.086260812400e-0010.000000000000e+000 -2.118477462187e+000 -2.361527908150e+000 -3.455612566063e-002 -4.736641031445e-002 1.816383623683e+000 -2.320066714159e-001 -1.435516012783e-001 -6.537076947738e-001 3.227152127219e-0020.000000000000e+000 3.555130771199e-001 -2.528514958236e+000 -6.375262961974e-001 2.023*********e-002 1.838645450412e+000 1.867660529118e-001 -2.932577853277e-001 1.987065844229e+000 1.004631189796e+0000.000000000000e+000 0.000000000000e+000 0.000000000000e+000 1.577548557113e+000 -1.396212144884e-002 -6.971850695189e-001 1.556103289616e-001 8.405054265525e-003 -8.154391705226e-002 -1.086119035380e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -1.484039822259e+000 -1.005266147971e-001 4.250492249735e-002 2.623603124719e-002 1.040075636780e-001 -1.180720843081e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 -6.948931346337e-001 -5.289549237019e-001 2.679156722088e-001 -5.962368505102e-001 -4.911592278261e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 1.787927773015e-001 -1.266168777947e+000 4.715000043044e-002 2.904780036865e-001 -3.570390800282e-0020.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+0000.000000000000e+000 0.000000000000e+000 0.000000000000e+000 9.355889078188e-0011.877411050997e-001 1.360256519312e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 6.360511809656e-001 2.737034413516e-0010.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 0.000000000000e+000 2.457611471065e-005 5.651649654397e-0023.3 矩阵A的全部特征值λｉ输出特征值第 1 个特征值λ(1) = 3.383039617436e+000 + 0.000000000000e+000*i第 2 个特征值λ(2) = -2.323496210212e+000 + 8.930405177200e-001*i第 3 个特征值λ(3) = -2.323496210212e+000 + -8.930405177200e-001*i第 4 个特征值λ(4) = 1.577548557113e+000 + 0.000000000000e+000*i第 5 个特征值λ(5) = -1.484039822259e+000 + 0.000000000000e+000*i第 6 个特征值λ(6) = -9.805309562902e-001 + 1.139489127430e-001*i第 7 个特征值λ(7) = -9.805309562902e-001 + -1.139489127430e-001*i第 8 个特征值λ(8) = 9.355889078188e-001 + 0.000000000000e+000*i第 9 个特征值λ(9) = 6.360627875745e-001 + 0.000000000000e+000*i第 10 个特征值λ(10) = 5.650488993501e-002 + 0.000000000000e+000*i3.4 A的相对应于实特征值的特征向量当实根特征值为λ = 3.383039617436e+000 时，对应的特征向量为：-1.052215734361e-001 9.032057866994e-001 3.851012511178e-001 -1.353018197378e-001 -7.981862313658e-002 1.258084765899e-002 -2.915794543278e-003 6.007606769995e-004 1.332409759711e-004 2.163020012305e-005当实根特征值为λ = 1.577548557113e+000 时，对应的特征向量为：-6.249625599173e-002 -3.448156083143e-001 1.479393070867e-004 -4.538473787534e-001 -7.818522244103e-001 2.014926926704e-001 -1.235718198662e-001 5.457885803205e-002 3.036557783871e-002 1.248312544887e-002当实根特征值为λ = -1.484039822259e+000 时，对应的特征向量为：5.601181168002e-001 2.218928269229e-001 -3.451922928463e-001 -3.485819655265e-001 9.829530157885e-002 2.571357665263e-001 3.880798023772e-001 3.732217910296e-001 -1.696444052682e-001 4.363177249054e-002当实根特征值为λ = 9.355889078188e-001 时，对应的特征向量为：8.0565********e-002 2.046856173180e-001 1.265257503741e-001 5.941904802900e-002 1.807775549629e-001 -1.010*********e-001 -2.0754********e-001 3.798709310692e-001 6.217085576446e-001 5.615092791522e-001当实根特征值为λ = 6.360627875745e-001 时，对应的特征向量为：1.070374683972e-001 1.807102811285e-001 4.050115674185e-001 1.531478777025e-0012.488923958328e-001 -1.845539115249e-001 4.633916211656e-001 -1.298460949311e-001 2.915801247322e-001 5.964727901835e-001当实根特征值为λ = 5.650488993501e-002 时，对应的特征向量为：-2.0901********e-001 -8.235701107749e-002 -3.216570398280e-001 1.974260440778e-002 -1.472246723965e-001 3.149431253984e-001 -2.871918085228e-003 -2.924742077289e-001 -4.618510426992e-001 6.526529750275e-001。