北航数值分析大作业二(纯原创,高分版)

一、算法设计方案：按题目要求，本程序运用带双步位移的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)，再求解方程组。

应用数理统计作业二学号：姓名：电话：二〇一四年十二月对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日左右结束。

《数值分析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。

数值分析第二题1 算法设计方案要想得出该题的答案首先要将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解。

要得出矩阵A 的全部特征值先对A 进行QR 的双步位移得出特征值。

采用列主元的高斯消元法求解特征向量。

1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构，因此在进行QR 分解前对矩阵A 进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。

具体算法如下所示，记A A =)1(，并记)(r A 的第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)( ++=全为零，则令)()1(r r A A =+,转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. 令()n Tr 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 -=ωTrr T r r r r p u u A A --=+ω)()1(5. 继续。

1.2 A 的QR 分解具体算法如下所示，记)1(1-=n A A ，并记[]nn r ij r a A ⨯=)(，令I Q =1 对于1,,2,1-=n r 执行1.若()n r r i a r ir,,3,1)( ++=全为零，则令r r Q Q =+1r r A A =+1,转5；否则转2。

2.计算()∑==nri r irr a d 2)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,2r r r r r r a c c h -=3.令()n Tr nrr r r r r r r r R a a c a u ∈-=+)()(,2)(,,,,,0,,0 。

北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。

我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。

我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。

在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。

通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。

第二次大作业是关于数值积分的方法。

数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。

在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。

我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。

在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。

通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。

第三次大作业是关于常微分方程的数值解法。

常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。

在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。

我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。

在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。

通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。

总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。

通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。

虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。

《数值分析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为虚数，不需要求特征向量。

数值分析第二题 梁进明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 -。

数理统计大作业（二）全国各省发展程度的聚类分析及判别分析指导教师院系名称材料科学与工程院学号学生姓名2015 年 12 月21 日目录全国各省发展程度的聚类分析及判别分析 (1)摘要： (1)引言 (1)1实验方案 (2)1.1数据统计 (2)1.2聚类分析 (3)1.3判别分析 (4)2结果分析与讨论 (5)2.1聚类分析结果 (5)2.2聚类分析结果分析： (8)2.3判别分析结果 (9)2.4 Fisher判别结果分析： (11)参考文献： (16)全国各省发展程度的聚类分析及判别分析摘要：利用SPSS软件对全国31个省、直辖市、自治区(浙江、安徽、甘肃除外)的主要经济指标进行多种聚类分析，分析选择最佳聚类类数，并对浙江、湖南、甘肃进行类型判别分析。

通过这两个方法对全国各省进行发展分类。

本文选取了7项社会发展指标作为决定发展程度的影响因素，其中经济因素为主要因素，同时评估城镇化率和人口素质因素。

各项数据均来自2014年国家统计年鉴。

分析结果表明：北京市和上海市和天津市为同一类；江苏省和山东省和广东省为同一类型；河北、湖北、河南、湖南、四川、辽宁为同一类；其余的为另一类。

关键词：聚类分析、判别分析、发展引言聚类分析是根据研究对象的特征对研究对象进行分类的多元统计分析技术的总称。

它直接比较各事物之间的性质，将性质相近的归为一类，将性质差别较大的归入不同的类。

系统聚类分析又称集群分析，是聚类分析中应用最广的一种方法，它根据样本的多指标（变量）、多个观察数据，定量地确定样品、指标之间存在的相似性或亲疏关系，并据此连结这些样品或指标，归成大小类群，构成分类树状图或冰柱图。

判别分析是根据多种因素(指标)对事物的影响来实现对事物的分类，从而对事物进行判别分类的统计方法。

判别分析适用于已经掌握了历史上分类的每一个类别的若干样品，希望根据这些历史的经验（样品），总结出分类的规律性（判别函数）来指导未来的分类。

2020级数值分析第二次作业(非线性方程求根)参考答案和评分标准班级学号姓名一(20分)用二分法求方程3()10f x x x =--=在区间[1.0,1.5]内的一个实根，且要求有3位有效数字。

试完成：(1)估计需要二分的次数；(8分)解：容易知道方程在[1.0,1.5]有且仅有一个实根。

记此实根为*x ，根据二分法误差估计公式有12()1*22k k k b a x x ++--≤=要使得近似解有3位有效数字，只需要有22111022k -+≤⨯从而可得6k ≥，即满足精度要求的二分次数为6次。

(2)将计算过程中数据填入表1.(中间过程填写到小数点后面3位)(12分，每个k x 得2分，其它空不计分)表1题1计算过程kk a k b kx 0 1.0 1.5 1.251 1.25 1.5 1.3752 1.25 1.375 1.3163 1.313 1.375 1.3494 1.313 1.344 1.3285 1.313 1.328 1.32061.3201.3281.324二.(10分)为了计算方程()3sin 2120f x x x =--=的根，某同学将()0f x =改写为14sin 23x x =+，并建立迭代公式114sin 23k k x x +=+。

请问此迭代公式在R 上是否全局收敛的吗？说明理由。

证明：(1)对任意的x R ∈，有11113()4sin 2,333x x R ϕ⎡⎤=+∈⊆⎢⎣⎦；(2)对任意的x R ∈，有22'()cos 2133x x ϕ=≤<;从而可知，迭代格式在R 上全局收敛。

三.(20分)设有方程3()10f x x x =--=，试回答下列问题：(1)确定方程3()10f x x x =--=实根的数目；(4分)解：由2'()31f x x =-可知函数3()1f x x x =--的单调递增区间是,33⎛⎛⎫-∞-⋃+∞ ⎪ ⎪⎝⎭⎝⎭，单调递减区间是,33⎛- ⎝⎭。

北京航空航天大学数值分析历年部分试题整理
(2001
2002
2003
2006
2008
2009)
2010年3月16日
北航数值分析2001年试题
北航数值分析2002年试题
北航数值分析2003年试题
北航数值分析2006年试题
北航数值分析2008年试题（版本一）
北航数值分析2008年试题（版本二）
北航数值分析2008年试题答案
一.选择题
1.C
2.A
3.D
4.D
5.A
6.D
二.填空题
1.错误！未找到引用源。

2.错误！未找到引用源。

3.错误！未找到引用源。

4.错误！未找到引用源。

5.0
6.0
7.错误！未找到引用源。

北航数值分析2009年试题
感谢手工抄写试卷的zhp、jhw、zby同学，以上试卷答案为本人和同学讨论所得，不保证正确性，请同学们参考使用，祝大家取得好成绩(^_^)
2010年3月16日。

一、算法设计方案：按题目要求，本程序运用带双步位移的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)，再求解方程组。

数理统计第二次大作业材料行业股票的聚类分析与判别分析2015年12月26日材料行业股票的聚类分析与判别分析摘要1 引言2 数据采集及标准化处理2.1 数据采集本文选取的数据来自大智慧软件的股票基本资料分析数据，从材料行业的股票中选取了30支股票2015年1月至9月的7项财务指标作为分类的自变量，分别是每股收益（单位：元）、净资产收益率（单位：%）、每股经营现金流（单位：元）、主营业务收入同比增长率（单位：%）、净利润同比增长率（单位：%）、流通股本（单位：万股）、每股净资产（单位：元）。

各变量的符号说明见表2.1，整理后的数据如表2.2。

表2.1 各变量的符号说明自变量符号每股收益（单位：元）X1净资产收益率（单位：%）X2每股经营现金流（单位：元）X3主营业务收入同比增长率（单位：%）X4净利润同比增长率（单位：%）X5流通股本（单位：万股）X6每股净资产（单位：元）X7表2.2 30支股票的财务指标股票代码X1 X2 X3 X4 X5 X6 X7 武钢股份600005-0.0990-2.81-0.0237-35.21-200.231009377.98 3.4444宝钢股份6000190.1400 1.980.9351-14.90-55.011642427.88 6.9197山东钢铁600022-0.11650.060.0938-20.5421.76643629.58 1.8734北方稀土6001110.0830 3.640.652218.33-24.02221920.48 2.2856杭钢股份600126-0.4900-13.190.4184-36.59-8191.0283893.88 3.4497抚顺特钢6003990.219310.080.1703-14.26714.18112962.28 1.4667盛和资源6003920.0247 1.84-0.2141-5.96-19.3739150.00 1.2796宁夏建材6004490.04000.510.3795-22.15-92.3447818.108.7321宝钛股份600456-0.2090-2.53-0.3313-14.81-6070.2043026.578.1497山东药玻6005290.4404 5.26 1.2013 6.5016.7825738.018.5230国睿科技6005620.410011.53-0.2949 3.3018.9416817.86 3.6765海螺水泥600585 1.15169.05 1.1960-13.06-25.33399970.2612.9100华建集团6006290.224012.75-0.57877.90-6.4034799.98 1.8421福耀玻璃6006600.790014.250.9015 3.6017.27200298.63 6.2419宁波富邦600768-0.2200-35.02-0.5129 3.1217.8813374.720.5188马钢股份600808-0.3344-11.710.3939-21.85-689.22596775.12 2.6854亚泰集团6008810.02000.600.1400-23.63-68.16189473.21 4.5127博闻科技6008830.503516.71-0.1010-10.992612.8023608.80 3.0126新疆众和6008880.0523 1.04-0.910662.64162.0464122.59 5.0385西部黄金6010690.0969 3.940.115115.5125.5712600.00 2.4965中国铝业601600-0.0700-2.920.2066-9.0882.79958052.19 2.3811明泰铝业6016770.2688 4.66-1.09040.8227.8640770.247.4850金隅股份6019920.1989 3.390.3310-10.05-39.01311140.26 6.7772松发股份6032680.35007.00-0.3195-4.43-9.622200.00 6.0244方大集团0000550.0950 5.66-0.480939.2920.6742017.94 1.6961铜陵有色0006300.0200 1.220.6132 3.23-30.74956045.21 1.5443鞍钢股份000898-0.1230-1.870.7067-27.32-196.21614893.17 6.4932中钢国际0009280.572714.45-0.4048-14.33410.2441286.57 4.2449中材科技0020800.684610.27 1.219547.69282.1740000.00 6.8936中南重工0024450.1100 4.300.340518.8445.0950155.00 2.70302.2 数据的标准化处理由于不同的变量之间存在着较大的数量级的差别，因此要对数据变量进行标准化处理。

数值分析第二次大作业史立峰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。

“数值分析“计算实习大作业第二题——SY1415215孔维鹏一、计算说明本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。

在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L；2、记，令k=1，m=n；3、如果，则可直接计算出最后1或2个特征值，转8，否则转4；4、如果，则可得一个特征值，置m=m-1；转3，否则转5；5、如果，则可得两个特征值，置m=m-2；转3，否则转6；6、记，计算7、k=k+1,转38、A的全部特征值已经求出，停止计算。

二、计算源程序（FORTRAN）PROGRAM SY1415215_2PARAMETER (N=10)DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部REAL(8) A,A1,A2,C,Q,R,CME=1E-12 !精度水平L=1000 !迭代最大次数OPEN(1,FILE='数值分析大作业第二题计算结果.TXT')DO I=1,NDO J=1,NIF(I==J) THENA(I,J)=*COS(I+*J)ELSEA(I,J)=SIN*I+*J)ENDIFENDDOENDDOA1=AWRITE(*,"('矩阵A为：')")WRITE(1,"('矩阵A为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"(' ')")WRITE(1,"(' ')")ENDDO!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALL HESSENBERG(A,N)WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")WRITE(1,"('拟上三角化后矩阵A（n-1）为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵A2=ACALL QRD(A2,N,Q,R)WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") Q(I,J)WRITE(1,"(2X,,2X,\)") Q(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOWRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") R(I,J)WRITE(1,"(2X,,2X,\)") R(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用带双步位移的QR方法求解矩阵A（n-1）的特征值K=1M=NDO WHILE(K<=L)IF(M<=2) THENIF(M==1) THENC(1,M)=A(M,M)ELSE IF(M==2) THENCALL CALCUS(A,N,M,C)ENDIFEXITELSE IF(ABS(A(M,M-1))<E) THENC(1,M)=A(M,M)M=M-1ELSE IF(ABS(A(M-1,M-2))<E) THENCALL CALCUS(A,N,M,C)M=M-2ELSECALL CALM(A,M,N)ENDIFK=K+1ENDDOWRITE(*,"('矩阵A的全部特征值为：')")WRITE(1,"('矩阵A的全部特征值为：')")DO J=1,NWRITE(*,",'+',,'i')") C(1,J),C(2,J)WRITE(1,",'+',,'i')") C(1,J),C(2,J)ENDDO!使用反幂法求解A的相应于实特征值的特征向量J=1DO I=1,NIF(C(2,I)==0)THENCR(J)=C(1,I)J=J+1ENDIFENDDOJC=J-1WRITE(*,"('矩阵A的实特征值为：')")WRITE(1,"('矩阵A的实特征值为：')")DO I=1,JCWRITE(*,"") CR(I)WRITE(1,"") CR(I)ENDDODO II=1,JCDO I=1,NDO J=1,NIF(I==J) THENA(I,J)=A1(I,J)-CR(II)ELSEA(I,J)=A1(I,J)ENDIFENDDOENDDOCALL INPOVERMETHOD(A,N,CM)WRITE(*,"('与实特征值',,'对应的特征向量为：')") CR(II) WRITE(1,"('与实特征值',,'对应的特征向量为：')") CR(II) DO I=1,NWRITE(*,"(2X,,2X,\)") CM(I)WRITE(1,"(2X,,2X,\)") CM(I)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOCLOSE(1)END!***************拟上三角化子函数*************************！SUBROUTINE HESSENBERG(A,N)DIMENSION A(N,N),P(N),Q(N),W(N),U(N),AT(N,N)REAL(8) A,P,Q,W,U,ATREAL(8) S0,S1,S2,S3,S4,TDO L=1,N-2JUDGE=0DO I=L+2,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIFENDDOIF(JUDGE==0) THENA=ACYCLEELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L+1,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L+1,L)==0)THENCR=DRELSECR=-SGN(A(L+1,L))*DRENDIF!计算HRHR=CR**2-CR*A(L+1,L)!给u赋值IF(I<L+1) THENU(I)=0ELSE IF(I==L+1) THEN U(I)=A(I,L)-CRELSE IF(I>L+1) THEN U(I)=A(I,L)ENDIFENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS1=0DO J=1,NS1=S1+AT(I,J)*U(J)ENDDOP(I)=S1/HRENDDO!计算QDO I=1,NS2=0DO J=1,NS2=S2+A(I,J)*U(J)ENDDOQ(I)=S2/HRENDDO!计算TS3=0DO I=1,NS3=S3+P(I)*U(I)ENDDOT=S3/HR!计算WDO I=1,NW(I)=Q(I)-T*U(I)!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J)ENDDOENDDOENDIFENDDORETURNEND!***************符号函数子程序*****************！FUNCTION SGN(X)REAL(8) XIF(X>0) THENSGN=1ELSE IF(X<0) THENSGN=-1ELSE IF(X==0) THENSGN=0ENDIFEND!*********计算二阶子阵特征值s1，s2子函数*******！SUBROUTINE CALCUS(X,N,M,Y)DIMENSION X(N,N),Y(2,N)REAL(8) A,B,C,D,X,YA=1C=X(M-1,M-1)*X(M,M)-X(M-1,M)*X(M,M-1)B=-(X(M-1,M-1)+X(M,M))D=B**2-4*CIF(D>=0) THENY(1,M)=(-B-SQRT(D))/2Y(1,M-1)=(-B+SQRT(D))/2ELSEIF(D<0) THENY(1,M)=-B/2Y(1,M-1)=-B/2Y(2,M)=-SQRT(-D)/2Y(2,M-1)=-SQRT(-D)/2ENDIFRETURNEND!*********计算Mk，Ak+1子函数************！SUBROUTINE CALM(A,M,N)DIMENSION A(N,N),MK(M,M),X(M,M),QK(M,M),RK(M,M),S1(M,M),S2(M,M),QKT(M,M) REAL(8) A,MK,X,QK,RK,QKTREAL(8) S0,S1,S2DO I=1,MDO J=1,MIF(I==J) THENX(I,J)=1ELSEX(I,J)=0ENDIFENDDOENDDOS=A(M-1,M-1)+A(M,M)T=A(M-1,M-1)*A(M,M)-A(M,M-1)*A(M-1,M)DO I=1,MDO J=1,MS0=0DO K=1,MS0=S0+A(I,K)*A(K,J)ENDDOMK(I,J)=S0-S*A(I,J)+T*X(I,J)ENDDOENDDO!对Mk做QR分解CALL QRD(MK,M,QK,RK)DO I=1,MDO J=1,MQKT(I,J)=QK(J,I)ENDDOENDDODO I=1,MDO J=1,MS1(I,J)=0DO K=1,MS1(I,J)=S1(I,J)+QKT(I,K)*A(K,J)ENDDOENDDOENDDOA=S1DO I=1,MDO J=1,MS2(I,J)=0DO K=1,MS2(I,J)=S2(I,J)+A(I,K)*QK(K,J)ENDDOENDDOENDDOA=S2RETURNEND!************QR分解子程序***************！SUBROUTINE QRD(A,N,Q,R)DIMENSION A(N,N),AT(N,N),Q(N,N),U(N),W(N),P(N),R(N,N) REAL(8) A,AT,Q,U,W,P,RREAL(8) DR,S0,CR,HR,S1,S2DO I=1,NDO J=1,NIF(I==J) THENQ(I,J)=1ELSEQ(I,J)=0ENDIFENDDOENDDODO L=1,N-1JUDGE=0DO I=L+1,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIF!A(I,L)中有一个不为零,判断条件为真，跳出循环转ENDDOIF(JUDGE==0) THENQ=QA=ACYCLE!A(I,L)全为零，结束本循环，进入下一个ELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L,L)==0)THENCR=DRELSECR=-SGN(A(L,L))*DRENDIF!计算HRHR=CR**2-CR*A(L,L)!给u赋值DO I=1,NIF(I<L) THENU(I)=0ELSE IF(I==L) THENU(I)=A(I,L)-CRELSE IF(I>L) THENU(I)=A(I,L)ENDIFENDDO!计算WDO I=1,NS1=0DO J=1,NS1=S1+Q(I,J)*U(J)ENDDOW(I)=S1ENDDO!计算Q（r+1）DO I=1,NDO J=1,NQ(I,J)=Q(I,J)-W(I)*U(J)/HRENDDOENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS2=0DO J=1,NS2=S2+AT(I,J)*U(J)ENDDOP(I)=S2/HRENDDO!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-U(I)*P(J)ENDDOENDDOENDIFENDDOQ=QR=ARETURNEND!*************运用反幂法求解矩阵A实特征值的特征向量***********！SUBROUTINE INPOVERMETHOD(A,N,Y)DIMENSION A(N,N),U(N),Y(N),U1(N,N),L1(N,N)REAL(8) E,Z,Z1,Z2,S1,S2,BREAL(8) A,U,Y,U1,L1U(I)=1ENDDO!任取非零向量UCALL DETA(A,N,U1,L1)Z2=EIZ1=E=1K=1DO WHILE (E>1E-12)S1=0DO I=1,NS1=S1+U(I)**2ENDDOB=SQRT(S1) !1DO I=1,NY(I)=U(I)/BENDDO!2CALL DOOLITTLE(U1,L1,Y,N,U) !3利用DOOLITTLE分解法法求解Au=yS2=0DO I=1,NS2=S2+Y(I)*U(I)ENDDOZ1=Z2Z2=S2 !4E=ABS(1/Z2-1/Z1)/ABS(1/Z2) !判断是否满足精度K=K+1ENDDORETURNENDSUBROUTINE DOOLITTLE(U,L,B1,N,X)DIMENSION B(N),U(N,N),X(N),Y(N),B1(N)REAL(8) L(N,N)REAL(8) B,U,X,Y,B1REAL(8) S1,S2,S3,S4B=B1Y(1)=B(1)S3=0DO M=1,I-1S3=S3+L(I,M)*Y(M)ENDDOY(I)=B(I)-S3ENDDOX(N)=Y(N)/U(N,N)DO I=N-1,1,-1S4=0DO M=I+1,NS4=S4+U(I,M)*X(M)ENDDOX(I)=(Y(I)-S4)/U(I,I)ENDDORETURNENDSUBROUTINE DETA(A1,N,U,L) DIMENSION A(N,N),U(N,N),A1(N,N) REAL(8) L(N,N)REAL(8) X,S1,S2REAL(8) A,U,A1X=1A=A1!对矩阵A进行Doolittle分解DO K=1,NDO J=K,NS1=0DO M=1,K-1S1=S1+L(K,M)*U(M,J)ENDDOU(K,J)=A(K,J)-S1A(K,J)=U(K,J)ENDDOIF (K==N) THENEXITELSEDO I=K+1,NS2=0DO M=1,K-1S2=S2+L(I,M)*U(M,K)ENDDOL(I,K)=(A(I,K)-S2)/U(K,K)A(I,K)=L(I,K)ENDDOENDIFENDDORETURNEND三、计算结果矩阵A为：+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +01 +00 +00 +00 +00 +00 +00+00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +01拟上三角化后矩阵A（n-1）为：+00 +01 +00 +00 +01 +00 +00 +00+01 +01 +01 +00 +00 +01 +00 +01 +00 +00+01 +01 +01 +00 +01 +00 +00 +00 +00+01 +01 +01 +01 +00 +00 +00 +00+01 +00 +00 +00 +00 +00+00 +01 +00 +00 +00 +00+00 +00 +00+00 +00 +00 +00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得Q为：+00 +00 +00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00 +00 +00+00 +00+00 +00 +00 +00 +00+00 +00 +00 +00+00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得R为：+01 +01 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +01 +01 +00 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +01 +01 +00 +00 +00+00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00+00 +00 +00 +00 +00 +00+00 +00 +00 +00矩阵A的全部特征值为：+01+ +00i+01++00i+01++00i+01+ +00i+01+ +00i+00++00i+00++00i+00+ +00i+00+ +00i+ +00i矩阵A的实特征值为：+01+01+01+00+00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00 +00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值+01对应的特征向量为：+00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00。

北航研究生数值分析作业第二题北航研究生数值分析作业第二题：一、算法设计方案1.按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为a_init()；2.对矩阵A进行householder变换，使其拟上三角化：对应的函数为householder()；3.输出拟上三角化后的A：对应的函数为aout(int)；4.对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代次数L=500），逐个求出其特征值，对应的函数为eigen_a()；中间包含两个子程序：calc_mk()和qr_analyze()，分别用来计算矩阵M k和对M k进行QR 分解并得到A k+1;5.输出QR分解过程完毕后的A及求得的特征向量：对应的函数为aout()和eigenvalout()；6.对于在第三步中求得的每个实特征值，使用带原点平移的反幂法求出其对应的特征向量，对应的函数为eigenvec()；其中包含一个解方程(A-μI)=y k-1的程序段。

这部分也用迭代完成，仍然将最大迭代次数L设置为500；7.输出矩阵A的特征向量，结束计算：对应的函数为eigenvecout()。

算法编译环境：vlsual c++6.0二、源程序如下：#include#include#define N 10 //矩阵阶数；#define EPSL 1.0e-12 //迭代的精度水平；#define L 500 //迭代最大次数；#define OUTPUTMODE 1 //输出格式：0--输出至屏幕，1--输出至文件double a[N][N], a2[N][N], eigen[N][N]; //声明矩阵A；double sa_re[N] = {0}, sa_im[N] = {0}; //声明矩阵的特征值数组；double u_init[N] = {2,1,2,1,2,1,2,1,2,1}; //定义反幂法中使用的初始向量u；//主程序开始；int main(){FILE *p;void a_init();void householder();void equal_zero(double matrix[N][N], int);void eigenvec();int eigen_a();void aout(int);void eigenvalout(int);void eigenvecout(int);if(OUTPUTMODE){p = fopen("Result.txt", "w+");fprintf(p, "计算结果：\n");fclose(p);}a_init(); //对矩阵A进行初始化；householder(); //对矩阵A进行拟上三角化；equal_zero(a, N); //对矩阵A的元素进行归零处理，消除误差；aout(OUTPUTMODE); //输出A；if(eigen_a()) printf("迭代超过最大次数，特征值求解结果可能不正确。

数值分析一、单项选择题（共20分,每小题2分）1-110=11=12=，则Lagranage 二次插值多项式为（ ） A.2(121)(144)(100)(144)(100)(121)()101112(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ B ．2(121)(144)(100)(144)(100)(121)()111012(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ C ．2(121)(144)(100)(144)(100)(121)()121110(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ D ．2(121)(144)(100)(144)(100)(121)()101211(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ 1-210=11=12=，用Lagranage为（ ）精确到小数点后4位。

A.9.7227 B ．11.7227 C ．10.7227 D ．13.72271-3、已知(1 2 3 4)TX =，则向量X 的21, , Xx x ∞的值分别是：（ ）B. -9，212,7C. 4,5,6D. 9,4,71-4、设 2121A --⎛⎫= ⎪⎝⎭，则21,, , F A A A x ∞的值分别为（ ）A. 4B. -9，C.4,5,6D. 9,4,71-5、设节点00 (=0,1,2,...,n), (0),k x x kh k x x th t =+=+>则Newton 向前插值公式为（ ）A. 100010()()!k k nn k j f N x th f t j k -==∆+=+-∑∏ B. 110()()!k k nn n n n k j f N x th f t j k -==∆+=+-∑∏ C. 100010()()!k k nn k j f N x th f t j k -==∇+=+-∑∏ D. 110()()!k k nn n n n k j f N x th f t j k -==∇+=+-∑∏1-6、方程组⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++47401815622189622315694962424321432143214321x x x x x x x x x x x x x x x x 进行直接三角分解法得到的L 矩阵为（ ）A. 1211213321B.1613216241C.16332202102 D.12147165511-7、对方程组的系数矩阵123412312341346262414535x x x x x x x x x x x x x x ++-=⎧⎪++=-⎨++-=⎪--+=-⎩进行Crout 分解法得到的U 矩阵为（ ）A.1111363111261131-- B. 1111363111569171--C.111136611151091371-- D. 111166311223611121--1-8、1、已知642()1f x x x x =+-+，2, 2 (0,1,2,...)k x kh h k =+==，则[2,6,10,14,18,22,26,30]f =（ ）A ．5！B ．4！C ．0D ．11-9、1、已知64()f x x x =+，2, 2 (0,1,2,...)k x kh h k =+==，则[2,4,6,8,10,12,14]f =（ ）A ．5！B ．4！C ．0D ．11-10、复合Cotes 求积公式, 复合梯形求积公式和复合Simpson 求积公式的收敛阶分别为（ ） A ．5，1，3 B ．4，2 ，6 C ．6，2，4 D ．以上都不对1-11、对线性方程组1231231232211221x x x x x x x x x +-=⎧⎪++=⎨++=⎪⎩，若用Jocabi 迭代法和G-S 迭代法求解，则（ ）A.Jocabi 迭代法收敛和G-S 迭代法发散B. Jocabi 迭代法和G-S 迭代法均发散C. Jocabi 迭代法和G-S 迭代法均收敛D. Jocabi 迭代法发散和G-S 迭代法收敛1-12、对线性方程组1231213918 293x x x x x x x --=⎧⎪-+=⎨-+=⎪⎩，若用Jocabi 迭代法和G-S 迭代法求解（ ），则 B.Jocabi 迭代法收敛和G-S 迭代法发散 A. Jocabi 迭代法和G-S 迭代法均发散C. Jocabi 迭代法和G-S 迭代法均收敛D. Jocabi 迭代法发散和G-S 迭代法收敛1-13、设线性方程组为1231213918 293x x x x x x x --=⎧⎪-+=⎨-+=⎪⎩，则Jocabi 迭代格式和G-S 迭代格式分别为（ ），则(Ⅰ) 2311(1)()()1(1)()2(1)()311799917881899k k k k k k k x x x x x x x +++⎧=++⎪⎪⎪=+⎨⎪⎪=+⎪⎩(Ⅱ) 2311(1)()()1(1)(1)2(1)(1)311799917881899k k k k k k k x x x x x x x +++++⎧=++⎪⎪⎪=+⎨⎪⎪=+⎪⎩A.(Ⅰ)和(Ⅱ)B. (Ⅱ)和(Ⅰ)C.(Ⅰ)和(Ⅰ)D. (Ⅱ)和(Ⅱ)1-14、已知*x 是()f x 的 (2)m m ≥重根，则求重根的修正Newton 公式为（ ）1(). ()k k k k f x A x x mf x +=-' 10(). ()k k k f x B x x mf x +=-'111(). ()()()k k k k k k k f x C x x x x f x f x +--=--- 111()(). ()()k k k k k k k f x f x D x x f x x x -+--=--1-15、若记(),()k k k k y f x z f y ==，则对迭代格式1()k k x f x -=使用Aitken 加速后得到的新迭代迭代格式为（ ）21(()). (())2()k k k k k k k f x x A x x f f x f x x +-=--+21(()). ()(())2()k k k k k k k f x x B x f x f f x f x x +-=--+21(). 2k k k k k k k z y C x z z y x +-=--+ 21((())()). (())(())2()k k k k k k k f f x f x D x f f x f f x f x x +-=--+1-16、将积分区间[a,b]n 等分，分点为kh a x k +=，k=0,1,2,3,4....n,其中nab h -=，则复合梯形公式为( )A. ])()(4)([211∑-=++n k k b f x f a f hB.])()(2)([211∑-=++n k k b f x f a f hC.)]()(4)(2)([6102111b f x f x f a f hn k k n k k +++∑∑-=+-=D.)]()(4)(2)([6112110b f x f x f a f hn k k n k k +++∑∑-=+-=二、填空题（共20分,每空2分）2-1、根据数值方法的稳定性与算法设计原则在连加运算中要防止 ，在减法运算中要避免 ，在除法运算中要避免，在乘法运算中要避免 。

数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号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 λλ---++=的两个根。