北航数值分析大作业第二题
北航数值分析大作业二(纯原创,高分版)
(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的全部特征值；
北航数值分析实验报告
北航数值分析实验报告篇一：北航数值分析报告第一大题《数值分析》计算实习报告第一大题学号：DY1305姓名：指导老师：一、题目要求已知501*501阶的带状矩阵A，其特征值满足?1?2...?501。
试求：1、?1，?501和?s的值；2、A的与数?k??1?k?501??140最接近的特征值?ik（k=1，2，...，39）；3、A的（谱范数）条件数c nd(A)2和行列式de tA。
二、算法设计方案题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵A1如下所示。
?0?0?A1??a1??b??c0b a2bcc bb c............c bb ccb a500b0a 3...a499c?b??a501??0?0??由矩阵A的特征值满足的条件可知?1与?501之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为?501；否则为?1。
在求得?1与?501其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。
用反幂法求得的特征值必为?s。
由条件数的性质可得，c nd(A)2为模最大的特征值与模最小的特征值之比的模，因此，求出?1，?501和?s的值后，则可以求得c nd(A)2。
北航数值分析计算实习题目二 矩阵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 λλ---++=的两个根。
北航数理统计大作业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日左右结束。
北航研究生数值分析试题
∗⎞ ⎟的 A1 ⎠
矩阵。
三、（12 分）试用高斯列主元素法求解线性方程组
⎡ 1 3 −2 −4 ⎤ ⎡ x1 ⎤ ⎡3 ⎤ ⎢ 2 6 −7 −10 ⎥ ⎢ x ⎥ ⎢ −2 ⎥ ⎢ ⎥⎢ 2⎥ = ⎢ ⎥ ⎢ −1 −1 5 9 ⎥ ⎢ x3 ⎥ ⎢14 ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ x4 ⎦ ⎥ ⎣ −6 ⎦ ⎣ −3 −5 0 15 ⎦ ⎣ 四、（12 分）利用矩阵 A 的三角分解 A = LU 求解下列方程组 ⎛ 1 2 1 ⎞ ⎛ x1 ⎞ ⎛ 0 ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ 2 2 3 ⎟ ⎜ x2 ⎟ = ⎜ 3 ⎟ ⎜ −1 −3 0 ⎟ ⎜ x ⎟ ⎜ 2 ⎟ ⎝ ⎠⎝ 3 ⎠ ⎝ ⎠
第一章
1、近似数 x = 0.231 关于真值 x = 0.229 有（ （1）1；（2）2；（3）3；（4）4。
∗
绪论
一、选择题（四个选项中仅有一项符合题目要求，每小题 3 分，共计 15 分) ）位有效数字。
2、取 3 ≈ 1.732 计算 x = ( 3 − 1) ，下列方法中哪种最好？（
4
）
Ax
∞和
A ∞ 的值分别为（
）
3
(1) 8 , 8 ；
(2) 8 , 7 ；
(3) 8 , 6 ；
(4) 7 , 7 。
5 、若解线性代数方程组的 Gauss 部分选主元方法第二步得到的系数矩阵的第三列向量为
(2
6 3 2 −5 4 2 ) ，则第三步主行是（
T
） (4) 第 6 行。
(1) 第 2 行；
1 − cos x , sin x
x ≠ 0且 x << 1 ；
（2）
1 1− x , − 1+ 2x 1+ x
北航数值分析第二次大作业--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。
北航数理统计大作业 聚类分析
应用数理统计聚类分析与判别分析（第二次作业）学院：姓名：学号：2015年12月目录我国部分城市经济发展水平的聚类分析和判别分析................................. - 1 - 摘要：................................................................... - 1 -1. 引言 ................................................................ - 1 -2. 相关统计基础理论 .................................................... - 1 -2.1 聚类分析......................................................... - 1 -2.2 判别分析......................................................... - 2 -3. 模型建立 ............................................................ - 3 -3.1 设置变量......................................................... - 3 -3.2 数据收集和整理................................................... - 3 -4. 数据结果及分析 ...................................................... - 5 -4.1 聚类分析......................................................... - 5 -4.2 判别分析......................................................... - 7 -5. 结论 ............................................................... - 11 -参考文献................................................................ - 12 -我国部分城市经济发展水平的聚类分析和判别分析摘要：本文基于《中国统计年鉴》（2014年版）统计数据，统计全国各省市居民消费情况，包括各地区农村居民人均纯收入、农村居民人均现金消费、城镇居民人均可支配收入、城镇居民人均现金消费情况共4个指标，利用统计软件SPSS综合考虑各指标，对所选地区进行K-Means 聚类分析，利用Fisher 线性判别待判地区类型，进一步验证所建模型的有效性。
北航数值分析大作业 第二题 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为虚数，不需要求特征向量。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析第二次大作业史立峰SY1505327一、 方案（1）利用循环结构将sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的矩阵A ；（2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n -1)。
对A 拟上三角化，得到拟上三角矩阵A (n -1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)( +==。
对于2,,2,1-=n r 执行 1. 若()n r r i a r ir,,3,2)( ++=全为零，则令A(r+1) =A(r),转5；否则转2。
2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,12r rr r r r a c c h +-=3. 令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0 。
4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。
（3）使用带双步位移的QR 方法计算矩阵A (n -1)的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。
2. 记n n ij n a A A ⨯-==][)1()1()1(，令n m k ==,1。
3. 如果ε≤-)(1,k m m a ，则得到A 的一个特征值)(k mm a ，置1:-=m m （降阶），转4；否则转5。
4. 如果1=m ，则得到的一个特征值)(11k a ，转11；如果1>m ，则转3。
5. 求2阶子阵⎥⎥⎦⎤⎢⎢⎣⎡=----)()(1,(k),1)(1,1 k mm k m m m m k m m k a a a a D的两个特征值1s 和2s ，即计算二次方程det )()()(1.12=++---k k mm k m m D a a λλ的两个根1s 和2s 。
6. 如果2=m ，则得到A 的两个特征值1s 和2s ，转11；否则转7。
7. 如果ε≤--)(2,1k m m a ，则得到A 的两个特征值1s 和2s ，置2:-=m m （降阶），转4；否则转88. 如果L k =，则计算终止，未得到A 的全部特征值；否则转9。
9. 记),1(][)()(m j i a A m m k ij k ≤≤=⨯，计算kk T k k k k k k k k k k mm k m m k mm k m m k mmk m m Q A Q A QR M R Q M I tI sA A M a a a a t a a s ==+-=-=+=+------12)(,1)(1,)()(1,1)()(1,1)分解作对阶单位矩阵是( )m ( 10. 置1:+=k k ，转3。
11. A 的全部特征值已计算完毕，停止计算。
其中，k M 的QR 分解与1+k A 的计算用下列算法实现：记k m m r ij r m m ij k A C b B b M B ====⨯⨯1)()1(1,][,][。
对于1,,2,1-=m r 执行 1. 若()m r r i b r ir,,2,1)( ++=全为零，则令r r r r C C B B ==++11,,转5；否则转2。
2. 计算()∑==mri r irr b d 2)(()()r r r r r r r rr r d c a d b c ==-=+则取,0sgn )(,1)(若)(2r rrr r r b c c h -=3. 令()m Tr mr r r r r r rr r R b b c b u ∈-=+)()(,1)(,,,,0,,0 。
4. 计算r r T r r h u B v /=Trr r r v u B B -=+1r r T r r h u C p /=r r r r h u C q /=r r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u C C --=+ω15. 继续。
此算法执行完后，就得到m k C A =+1。
（4）计算Q,R ，一边求R*Q 矩阵。
记阶单位阵）n （令,][，并记A 1)(r 1I Q a A A n n r ij ===⨯ 对于1n ,,2,1-= r 执行 1.若()n r r i a r ir,,3,2)( ++=全为零，则令rr r r A A Q Q ==++11转5；否则转 2.2.计算()∑==mri r irr b d 2)(()()r r r r r r r rrr d c a d b c ==-=+则取,0sgn )(,1)(若)(2r rrr r r b c c h -=3.令()m Tr mr r r r r r rr r R b b c b u ∈-=+)()(,1)(,,,,0,,0 。
4.计算Trr r r rrTr rT r r r r r r p u A A h u A p h u Q Q u Q -==-==++1r 1r ωω5.继续当此算法执行完毕后，就得到正交矩阵n Q Q =和上三角矩阵R 6.然后计算出Q R 矩阵（5）用列主元素Gauss 消去法计算矩阵A 对应于实特征值的特征向量，具体算法如下：记)( )1(的实特征值为A I A Aλλ-=1. 消元过程对于1,,2,1-=n k 执行 (1) 选行号k i ，使)()(max k ikni k k k i aa k ≤≤=。
(2) 交换)(k kj a 与),,1,()(n k k j a k j i k +=所含的数值。
(3) 对于n k k i ,,2,1 ++=计算),,2,1( /)()()1()()(n k k j a m a a a a m k ki ik k ij k ij k kkk ik ik ++=-==+2. 回代过程()1,,2,1 1)(1)( --=⎪⎪⎭⎫ ⎝⎛-==∑+=n n k a x ax x k kk nk j j k kjk n最终得到的向量Tn x x x )1,,,(21= 的即为A 对应于实特征值λ的特征向量。
二、源程序#include<iostream.h>#include<math.h>#include<iomanip.h>#define n 10#define E 1.0e-12void Householder(double a[n][n]);//拟上三角化函数double sgn(double a);//符号函数void QRfenjie(double a[n][n],double Q[n][n],double R[n][n]);void QR(double a[n][n],double L[n][2]);//带双步位移的QR分解void MxM(double M[n][n],double A[n][n],double B[n][n],int m);//矩阵相乘void solve(double a[n][n],double s1[2],double s2[2],int m);//解方程函数void Gauss(double a[n][n], double x[n]);//定义列主元高斯消去法函数void main(){int i,j,k;double a[n][n],qr[n][n],q[n][n],r[n][n],L[n][2],x[n];for(i=0;i<n;i++)//录入要进行变换的矩a，并对q初始赋值。
{for(j=0;j<n;j++){if(i==j)a[i][j]=1.5*cos(i+1+1.2*(j+1));elsea[i][j]=sin(0.5*(i+1)+0.2*(j+1));}}// cout<<endl;Householder(a);//调用拟上三角化函数cout<<endl<<endl<<"对矩阵A拟上三角化前七列的结果:"<<endl;for(i=0;i<n;i++){//k=0;//为了显示方便，每行显示三个元素，使用k来实现//cout<<"矩阵A的第"<<i+1<<"行元素为："<<endl;for(j=0;j<7;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";//k++;//cout<<endl;//if(k%3==0)//判断每行是否到达三个元素}cout<<endl;}cout<<"对矩阵A拟上三角化后三列的结果:"<<endl;for(i=0;i<n;i++){//k=0;//为了显示方便，每行显示三个元素，使用k来实现//cout<<"矩阵A的第"<<i+1<<"行元素为："<<endl;for(j=7;j<n;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";//k++;//if(k%3==0)//判断每行是否到达三个元素//cout<<endl;}cout<<endl;}cout<<endl;QRfenjie(a,q,r);QR(a,L);cout<<endl<<"对矩阵A进行QR分解后所得矩阵前七列的结果:"<<endl; for(i=0;i<n;i++){for(j=0;j<7;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";}cout<<endl;}cout<<"对矩阵A进行QR分解后所得矩阵三列的结果:"<<endl;for(i=0;i<n;i++){for(j=7;j<n;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";}cout<<endl;}/*cout<<"对矩阵A进行QR分解后所得矩阵:"<<endl;for(i=0;i<n;i++){k=0;cout<<"矩阵A的第"<<i+1<<"行元素为："<<endl;for(j=0;j<n;j++){if (fabs(a[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<a[i][j]<<"";k++;if(k%3==0)cout<<endl;}cout<<endl;}*/for(i=0;i<n;i++){for(j=0;j<n;j++)for(k=0,qr[i][j]=0;k<n;k++)qr[i][j]=qr[i][j]+r[i][k]*q[k][j];}cout<<endl<<"R*Q相乘的前七列:"<<endl;for(i=0;i<n;i++){//k=0;//cout<<"R*Q的第"<<i+1<<"行元素为："<<endl;for(j=0;j<7;j++){if (fabs(qr[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<qr[i][j]<<"";//k++;//if(k%3==0)//cout<<endl;}cout<<endl;}cout<<endl<<"R*Q相乘的后三列:"<<endl;for(i=0;i<n;i++){for(j=7;j<n;j++){if (fabs(qr[i][j])<E)a[i][j]=0;cout<<setprecision(12)<<setiosflags(ios::scientific)<<qr[i][j]<<"";}cout<<endl;}cout<<endl<<"矩阵A的全部特征值为"<<endl;for(i=0;i<n;i++){if(L[i][1]==0)cout<<"λ"<<i+1<<"="<<L[i][0]<<endl;elsecout<<"λ"<<i+1<<"="<<L[i][0]<<"+"<<L[i][1]<<"i"<<endl;}cout<<endl<<"实特征值对应的特征向量分别是："<<endl;for(k=0;k<n;k++){if(L[k][1]==0)//判断特征值是否为实特征值{for(i=0;i<n;i++)//重新录入矩阵A{for(j=0;j<n;j++){if(i==j)a[i][j]=1.5*cos(i+1+1.2*(j+1))-L[k][0];elsea[i][j]=sin(0.5*(i+1)+0.2*(j+1));}}Gauss(a,x);cout<<"λ"<<k+1<<"="<<L[k][0]<<""<<"对应的特征向量是:"<<endl;for(j=0;j<n;j++){cout<<x[j]<<endl;}}else continue;}}void Householder(double a[n][n])//拟上三角化函数{int i,j,k;int m=0;double d,c,h,t;double u[n],p[n],q[n],w[n];for(i=0;i<n-2;i++){for(j=i+2;j<n;j++){if(a[j][i]<=E)m=m+1;}if(m==(n-2-i))continue;for(j=i+1,d=0;j<n;j++){d=d+a[j][i]*a[j][i];}d=sqrt(d);c=-1*sgn(a[i+1][i])*d;h=c*c-c*a[i+1][i];for(j=i+2;j<n;j++){u[j]=a[j][i];}for(j=0;j<i+2;j++){u[j]=0;}u[i+1]=a[i+1][i]-c;for(j=0;j<n;j++){for(k=i+1,p[j]=0;k<n;k++){p[j]=a[k][j]*u[k]+p[j];}p[j]=p[j]/h;}for(j=0;j<n;j++){for(k=i+1,q[j]=0;k<n;k++){q[j]=a[j][k]*u[k]+q[j];}q[j]=q[j]/h;}for(j=0,t=0;j<n;j++){t=t+p[j]*u[j];}t=t/h;for(j=0;j<n;j++){w[j]=q[j]-t*u[j];}for(j=0;j<n;j++){for(k=0;k<n;k++){a[j][k]=a[j][k]-w[j]*u[k]-u[j]*p[k];}}}} //以上是拟上三角化的结果double sgn(double a)//符号函数{if(a>0)return 1;else if(a<0)return -1;return 0;}//以上是符号函数void QRfenjie(double a[n][n],double Q[n][n],double R[n][n])//矩阵的QR分解{int i,j,k;double d,c,h;double u[n],w[n],p[n];for(i=0;i<n;i++){for(j=0;j<n;j++){if (i==j) Q[i][j]=1; else Q[i][j]=0;}}for(i=0;i<n;i++){for(j=0;j<n;j++) R[i][j]=a[i][j];}for(i=0;i<n-1;i++){for(j=i,d=0;j<n;j++)d=d+R[j][i]*R[j][i];d=sqrt(d);c=-1*sgn(R[i][i])*d;h=c*c-c*R[i][i];for(j=i+1;j<n;j++) u[j]=R[j][i];for(j=0;j<i;j++) u[j]=0;u[i]=R[i][i]-c;for(j=0;j<n;j++)for(k=0,w[j]=0;k<n;k++) w[j]=Q[j][k]*u[k]+w[j];}for(j=0;j<n;j++){for(k=0;k<n;k++) Q[j][k]=Q[j][k]-w[j]*u[k]/h;}for(j=0;j<n;j++){for(k=i,p[j]=0;k<n;k++)p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}for(j=0;j<n;j++){for(k=0;k<n;k++)R[j][k]=R[j][k]-u[j]*p[k];}}}//以上是矩阵的QR分解void MxM(double M[n][n],double A[n][n],double B[n][n],int m)//矩阵相乘{double C[n][n],D[n][n];int i,j,g;for(i=0;i<m;i++)for(j=0;j<m;j++){C[i][j]=A[i][j];D[i][j]=B[i][j];M[i][j]=0;}for(i=0;i<m;i++)for(j=0;j<m;j++)for(g=0;g<m;g++)M[i][j]=M[i][j]+C[i][g]*D[g][j];}void QR(double a[n][n],double L[n][2])//带双步位移的QR方法{double M[n][n],s,t,s1[2],s2[2],u[n],v[n],p[n],q[n],w[n],sum,h,c,d;int i,j,m,r,k;for(i=0;i<n;i++){for(j=0;j<2;j++){L[i][j]=0;}}m=n;for(k=1;k<=100;k++){if(fabs(a[m-1][m-2])<=E){L[m-1][0]=a[m-1][m-1];m=m-1;if(m==1){L[0][0]=a[0][0];break;}else if(m==0)break;elsecontinue;}else{solve(a,s1,s2,m);//调用解方程函数，将求解结果存到s1和s2中if(m==2){L[0][0]=s1[0];L[0][1]=s1[1];L[1][0]=s2[0];L[1][1]=s2[1];break;}else{if(fabs(a[m-2][m-3])<=E){L[m-2][0]=s1[0];L[m-2][1]=s1[1];L[m-1][0]=s2[0];L[m-1][1]=s2[1];m=m-2;if(m==1){L[0][0]=a[0][0];break;}else if(m==0)break;elsecontinue;}else{if(k==100){cout<<"未能得到全部特征值"<<endl;break;}else{t=a[m-1][m-1]*a[m-2][m-2]-a[m-1][m-2]*a[m-2][m-1];s=a[m-1][m-1]+a[m-2][m-2];MxM(M,a,a,m); //调用矩阵相乘函数构造矩阵Mfor(i=0;i<m;i++)for(j=0;j<m;j++)M[i][j]=M[i][j]-s*a[i][j];for(i=0;i<m;i++)M[i][i]=M[i][i]+t;for(r=0;r<(m-1);r++){for(i=(r+1);i<m;i++){if(M[i][r]==0)continue;else{sum=0;for(i=r;i<m;i++){sum+=pow(M[i][r],2);}d=sqrt(sum);if(M[r][r]!=0)c=-sgn(M[r][r])*d;else c=d;h=c*(c-M[r][r]);for(i=0;i<m;i++){if(i<r)u[i]=0;if(i==r)u[i]=M[r][r]-c;if(i>r)u[i]=M[i][r];}for(i=0;i<m;i++){sum=0;for(j=r;j<m;j++){sum+=M[j][i]*u[j];}v[i]=sum/h;}for(i=0;i<m;i++){for(j=0;j<m;j++){M[i][j]=M[i][j]-u[i]*v[j];} }for(i=0;i<m;i++){sum=0;for(j=r;j<m;j++){sum+=a[j][i]*u[j];}p[i]=sum/h;sum=0;for(j=r;j<m;j++){sum+=a[i][j]*u[j];}q[i]=sum/h;}sum=0;for(j=r;j<m;j++){sum+=p[j]*u[j];}t=sum/h;for(j=0;j<m;j++){w[j]=q[j]-t*u[j];}for(i=0;i<m;i++){for(j=0;j<m;j++){a[i][j]=a[i][j]-w[i]*u[j]-u[i]*p[j];}}}}}}}}}}for(i=0;i<n;i++){for(j=0;j<n;j++){if(a[i][j]<=E)a[i][j]=0;}}}void solve(double a[n][n],double s1[2],double s2[2],int m)//解方程函数{double s,t,det;t=a[m-1][m-1]*a[m-2][m-2]-a[m-1][m-2]*a[m-2][m-1];s=a[m-1][m-1]+a[m-2][m-2];det=s*s-4*t;if(det<0){det=-det;det=sqrt(det);s1[0]=s/2;s1[1]=det/2;s2[0]=s/2;s2[1]=-det/2;}else{det=sqrt(det);s1[0]=(s+det)/2;s1[1]=0;s2[0]=(s-det)/2;s2[1]=0;}}void Gauss(double a[n][n], double x[n])//定义列主元高斯消去法函数{int i,j,k,t;double sum,m[n][n],max,p;for(k=0;k<(n-1);k++){max=fabs(a[k][k]);t=k;for(i=(k+1);i<n;i++){if(max<fabs(a[i][k])){max=fabs(a[i][k]);t=i;}}if(t!=k){for(j=k;j<n;j++){p=a[k][j];a[k][j]=a[t][j];a[t][j]=p;}}for(i=(k+1);i<n;i++){m[i][k]=a[i][k]/a[k][k];for(j=(k+1);j<n;j++){a[i][j]=a[i][j]-m[i][k]*a[k][j];}}}x[n-1]=1;for(k=(n-2);k>=0;k--){sum=0;for(j=(k+1);j<n;j++){sum+=a[k][j]*x[j];}x[k]=(-sum)/a[k][k];}sum=0;for(i=0;i<n;i++){sum+=x[i]*x[i];}sum=sqrt(sum);for(i=0;i<n;i++){x[i]=x[i]/sum;}}三、处理结果四、结果讨论（1）此次编程涉及的函数众多，需要特别注意各函数之间的协调；（2）采用全局变量n节省了输入次数，简化编程；Q R矩阵，还需要另外的编辑函数来求解除Q和R，然后再求出Q R。