北航数值分析第二次大作业--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的全部特征值；
北航数值分析计算实习题目二 矩阵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分解求特征值算法
数 值 分 析（B ） 大 作 业（二）1、算法设计：①矩阵的拟上三角化：对实矩阵A 进行相似变换化为拟上三角矩阵(1)A n -，其变换矩阵采用Householder 矩阵，变换过程如下：若()a 0(2,,)r iri r n ==+ ，则r H I =； 否则，(r)(r)Tr r+1,r n,r s =(0,,0,a ,a ) ，(r)r r+1,r r 2c = -sgn(a )||s ||，()()()r r r r+11,2,u =s -c e (0,,0,,,,)r r r r r r r r nr a c a a ++=- ，2T r r r 2H =I-2u u /r u ，(1)()r r r r A H A H +=。
当2r n =-时，得(1)(2)(1)222112n n n n n n A H A H H H A H H ------=== ，令12n-2P=H H H 又r H 是对称正交矩阵，于是n-1T A =P AP 成立，因而n-1A 与 A 相似。
②矩阵的QR 分解：矩阵的QR 分解过程与拟上三角化过程相似，在这里不再重复其原理。
③求全部特征值矩阵拟上三角化后利用带双步位移的QR 方法，采用书本Page 63页具体算法实现。
为了使编程方便，并体会goto 语句使用的灵活性，程序中的跳转均使用goto Loop*实现。
④求A 的相应于实特征值的特征向量求实特征值对应的特征向量，即是求解线性方程组(λI-A)=0i i x ，(1,,)i n = 。
因此，为得到全部实特征值对应的特征向量，解线性方程组的过程要循环n 次（n 为矩阵阶数）。
线性方程组的求解采用列主元素Gauss 消去法。
#include <stdio.h>#include <math.h>#define ERR 1.0e-12 //误差限#define N10 //矩阵行列数#define L 1.0e5 //最大迭代次数double A[N][N]={0};void Init_A() //初始化矩阵{int i,j;for(i=0;i<N;i++)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));}}/*void Display_A(){int i,j;for(i=0;i<N;i++){for(j=0;j<N;j++)printf("%8.4f",A[i][j]);printf("\n");}}*/int Sgn(double a){if(a>=0)return 1;elsereturn -1;}void On_To_The_Triangle() //矩阵拟上三角化{int i,j,r,flag=0;double cr,dr,hr,tr,temp;double ur[N],pr[N],qr[N],wr[N];for(r=1;r<=N-2;r++){flag=0;for(i=r+2;i<=N;i++)if(A[i-1][r-1]!=0){flag=1;break;}if(0==flag)continue;dr=0;for(i=r+1;i<=N;i++)dr+=A[i-1][r-1]*A[i-1][r-1];dr=sqrt(dr);if(0==A[r][r-1])cr=dr;else cr=-Sgn(A[r][r-1])*dr;hr=cr*cr-cr*A[r][r-1];for(i=1;i<=r;i++)ur[i-1]=0;ur[r]=A[r][r-1]-cr;for(i=r+2;i<=N;i++)ur[i-1]=A[i-1][r-1];for(i=1;i<=N;i++){temp=0;for(j=1;j<=N;j++)temp+=A[j-1][i-1]*ur[j-1];pr[i-1]=temp/hr;}for(i=1;i<=N;i++){temp=0;for(j=1;j<=N;j++)temp+=A[i-1][j-1]*ur[j-1];qr[i-1]=temp/hr;}temp=0;for(i=1;i<=N;i++){temp+=pr[i-1]*ur[i-1];tr=temp/hr;}for(i=1;i<=N;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i<=N;i++)for(j=1;j<=N;j++)A[i-1][j-1]=A[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}}void Get_Roots(double eigenvalue[][2],int m,double ss,double tt) //求一元二次方程的根{double discriminant=ss*ss-4*tt; //if(discriminant<0){*(*(eigenvalue+m-2))=0.5*ss;*(*(eigenvalue+m-2)+1)=0.5*sqrt(-discriminant);*(*(eigenvalue+m-1))=0.5*ss;*(*(eigenvalue+m-1)+1)=-0.5*sqrt(-discriminant);}else{*(*(eigenvalue+m-2))=0.5*(ss+sqrt(discriminant));*(*(eigenvalue+m-2)+1)=0;*(*(eigenvalue+m-1))=0.5*(ss-sqrt(discriminant));*(*(eigenvalue+m-1)+1)=0;}}void Get_Mk(double mk[][N],int m,double ss,double tt) //获取Mk,用于带双步位移的QR分解{int i,j,k;for(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(k=0;k<m;k++)*(*(mk+i)+j)+=A[i][k]*A[k][j];*(*(mk+i)+j)-=ss*A[i][j];if(j==i)*(*(mk+i)+j)+=tt;}}void QR_Reslove(double mk[][N],int m) //QR分解{int i,j,r,flag=0;double cr,dr,hr,tr,temp;double ur[N],vr[N],pr[N],qr[N],wr[N];double B[N][N],C[N][N];for(i=0;i<m;i++)for(j=0;j<m;j++){B[i][j]=*(*(mk+i)+j);C[i][j]=A[i][j];}for(r=1;r<=m-1;r++){flag=0;for(i=r+1;i<=m;i++)if(B[i-1][r-1]!=0){flag=1;break;}if(0==flag)continue;dr=0;for(i=r;i<=m;i++)dr+=B[i-1][r-1]*B[i-1][r-1];dr=sqrt(dr);if(0==B[r-1][r-1])cr=dr;else cr=-Sgn(B[r-1][r-1])*dr;hr=cr*cr-cr*B[r-1][r-1];for(i=1;i<r;i++)ur[i-1]=0;ur[r-1]=B[r-1][r-1]-cr;for(i=r+1;i<=m;i++)ur[i-1]=B[i-1][r-1];for(i=1;i<=m;i++){temp=0;for(j=1;j<=m;j++)temp+=B[j-1][i-1]*ur[j-1];vr[i-1]=temp/hr;}for(i=0;i<m;i++)for(j=0;j<m;j++)B[i][j]-=ur[i]*vr[j];for(i=1;i<=m;i++){temp=0;for(j=1;j<=m;j++)temp+=C[j-1][i-1]*ur[j-1];pr[i-1]=temp/hr;}for(i=1;i<=m;i++){temp=0;for(j=1;j<=m;j++)temp+=C[i-1][j-1]*ur[j-1];qr[i-1]=temp/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++)C[i-1][j-1]=C[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}for(i=0;i<m;i++)for(j=0;j<m;j++)A[i][j]=C[i][j];}void Display_Eigenvalue(double value[][2]) //显示特征值{int i;for(i=0;i<N;i++){printf("λ%d=%8.4f",i+1,*(*(value+i)));if(*(*(value+i)+1)>0)printf("+%8.4f",*(*(value+i)+1));else if(*(*(value+i)+1)<0)printf("%8.4f",*(*(value+i)+1));printf("\n");}printf("\n");}int QR_With_Double_Step_Displacement(double eigenvalue[][2]) //带双步位移QR分解求特征值{int i,j,k=1,m=N;double s,t;double Mk[N][N];for(i=0;i<N;i++)for(j=0;j<2;j++)eigenvalue[i][j]=0;do{k++;if(m==1){eigenvalue[m-1][0]=A[m-1][m-1];m--;continue;}else if(m==2){s=A[m-2][m-2]+A[m-1][m-1];t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-1][m-2];Get_Roots(eigenvalue,m,s,t); //求一元二次方程的根m=0;continue;}else if(m==0)return 0;else if(fabs(A[m-1][m-2])<=ERR){eigenvalue[m-1][0]=A[m-1][m-1];m--;continue;}else{s=A[m-2][m-2]+A[m-1][m-1];t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-2][m-1];if(fabs(A[m-2][m-3])<=ERR){Get_Roots(eigenvalue,m,s,t); //求一元二次方程的根m-=2;continue;}else{Get_Mk(Mk,m,s,t); //获取Mk,用于带双步位移的QR分解QR_Reslove(Mk,m); //QR分解}}}while(k<L);printf("Errer\n"); //迭代过程不成功，报错}void Select_Principal_Element(int k) //Gauss消元法求解方程组之选主元{int i,max=k;double trans[N];for(i=k+1;i<N;i++)if(fabs(A[i][k])>fabs(A[max][k]))max=i;if(k==max)return;else{for(i=k;i<N;i++){trans[i]=A[k][i];A[k][i]=A[max][i];A[max][i]=trans[i];}}}void Eliminant(int k) //Gauss消元法求解方程组之消元{int i,j;double rate;for(i=k+1;i<N;i++){rate=A[i][k]/A[k][k];for(j=k;j<N;j++)A[i][j]=A[i][j]-rate*A[k][j];}}void Back_Substitution(double Eigenvector[]) //Gauss消元法求解方程组之回代{int i,j,k;double Temp_M[N][N+1];for(i=0;i<N-1;i++){for(j=i;j<N;j++)Temp_M[i][j]=A[i][j];Temp_M[i][N]=0;}Temp_M[N-1][N-1]=1;Temp_M[N-1][N]=1;for(k=N-1;k>=0;k--){*(Eigenvector+k)=Temp_M[k][N];for(i=k;i<N-1;i++)*(Eigenvector+k)=*(Eigenvector+k)-*(Eigenvector+i+1)*Temp_M[k][i+1];*(Eigenvector+k)=*(Eigenvector+k)/Temp_M[k][k];}}void Display_Eigenvector(double Eigenvector[],double eigenvalue) //显示特征值对应特征向量{int i;printf("When λ=%8.4f\nEigenvector=\n",eigenvalue);for(i=0;i<N;i++)printf("%8.4f",*(Eigenvector+i));printf("\n");}void Get_Eigenvector(double value[][2]) //利用Gauss消元法求解特征向量{int i,j;double eigenvalue;double Eigenvector[N];for(i=0;i<N;i++)if(*(*(value+i)+1)==0){Init_A(); //初始化矩阵eigenvalue=*(*(value+i));for(j=0;j<N;j++)A[j][j]-=eigenvalue;for(j=0;j<N-1;j++){Select_Principal_Element(j); //Gauss消元法求解方程组之选主元Eliminant(j); //Gauss消元法求解方程组之消元}Back_Substitution(Eigenvector); //Gauss消元法求解方程组之回代Display_Eigenvector(Eigenvector,eigenvalue); //显示特征值对应特征向量}}main(){double eigenvalue[N][2];Init_A(); //初始化矩阵On_To_The_Triangle(); //矩阵上三角化QR_With_Double_Step_Displacement(eigenvalue); //带双步位移QR分解求特征值printf("Contactme:****************\n\n");Display_Eigenvalue(eigenvalue); //显示特征值Get_Eigenvector(eigenvalue); //利用Gauss消元法求解特征向量return 0;}PS: Forfurtherdetails,pleasecontactme:****************。
数值分析大作业
第二次计算实验：SVD及其应用梁杰存2014310739航博1431.方法求矩阵A奇异值分解一个途径是求解A T A的特征值，但因为舍入误差容易丢掉小奇异值。
因此通常先将矩阵上双对角化，即构造正交阵Q和W,使得Q T AW=B(upper−bidiagnal)。
这一过程可以通过逐次Household变换或逐次Given’s变换完成，还有一种基于待定系数法思想的Lanczos算法。
由于Linpack中SVD算法需要输入上双对角矩阵，本文采用Lanczos 算法实现上双对角化。
1.1.隐式零移位QR法（implicit zero-shift QR）与传统移位QR迭代算法不同，隐式零移位QR算法不进行移位，并且第一步构造右乘Given’s变换矩阵GR(1,2)将上双对角矩阵B的(1,2)位置上的元素12b消零，而不是传统方法中引入一个非零元素21b。
但这一步可能会使原来为零的b12变为非零。
第二步左乘Given’s阵GL(1,2)使得12b为0，但可能会使为零b13变为非零。
与上述步骤类似，将b13变为0后可能会使b23非零。
如下图所示，重复上述步骤最终将恢复为上双对角矩阵，即完成一步隐式零移位QR迭代。
反复迭代，矩阵B将趋近与对角阵阵，对角元即特征值。
图1隐式QR迭代1.2.分而治之(Divide-and-conquer)分而治之算法将上双对角阵B分成有两个互相独立对角块矩阵与另一矩阵之和，即：B=B100B2+b m vvT=Q1Σ1Q1T00Q2Σ1Q2T+b m vv T =Q100Q2(Σ100Σ1+b m uuT)Q1T00Q2T所以矩阵B的特征值与矩阵D+ρu u T的特征值相同，其中D=Σ100Σ1为对角阵，又：det D+ρu u T−λI=det((D−λI)(I+ρD−λ−1u u T))由于D−λI非奇异，则det I+ρD−λ−1u u T=1+ρu T D−λ−1u=1+ρu i2d i−λ=0ni=1在每个d i与d i+1之间分布着一个特征值，可用牛顿法快速找到该特征值。
北航数值分析第二次大作业--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编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北航数值分析报告大作业第三题(fortran)
北航数值分析报告大作业第三题(fortran)“数值分析“计算实习大作业第三题——SY1415215孔维鹏一、计算说明1、将x i=0.08i，y j=0.5+0.05j分别代入方程组（A.3）得到关于t，u，v，w的的方程组，调用离散牛顿迭代子函数求出与x i，y j对应的t i，u j。
2、调用分片二次代数插值子函数在点(t i,u j)处插值得到z(x i,y j)=f(x i,y j)，得到数表(x i,y j,f(x i,y j))。
3、对于k=1,2,3,4?，分别调用最小二乘拟合子函数计算系数矩阵c rs 及误差σ，直到满足精度，即求得最小的k值及系数矩阵c rs。
4、将x i?=0.1i，y j?=0.5+0.2j分别代入方程组（A.3）得到关于t?，u?，v?，w?的的方程组，调用离散牛顿迭代子函数求出与x i?，y j?对应的t i?，u j?，调用分片二次代数插值子函数在点(t i?,u j?)处插值得到z?(x i?,y j?)=f(x i?,y j?)；调用步骤3中求得的系数矩阵c rs求得p(x i?,y j?)，打印数表(x i?,y j?,f(x i?,y j?),p(x i?,y j?))。
二、源程序（FORTRAN）PROGRAM SY1415215DIMENSIONX(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21), C(6,6) DIMENSIONX1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组，得到z=f（t，u）的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDO ENDDO!*************分片二次插值得到z=f（x，y）***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDO ENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P（x，y）**************N=11M=21WRITE (1,'(" ","K和σ分别为：")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E) WRITE (1,'(I3,2X,E20.13)') K-1,EIF(ETA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIF ENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J) ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A) DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I)) ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标，W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R************************* IF(UX(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(UY(M-1)-T/2) THENR=M-1 ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(VN) P=N IF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDO APX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*G ENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,M V(K,J)=V(K,J)+Z(I,J)*G ENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDO U(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*G APY(2)=APY(2)+(Y(I))*G*G ENDDO APY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*G DO J=1,PU(J,K)=U(J,K)+V(J,I)*G ENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDO ENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDO DO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2) IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K) ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L) ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)): X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.0790.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.3001.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.001.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769-0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454-0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.2620.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.1501.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718-0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.2250.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.3010.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+00 0.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477-0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.3630.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.2501.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739-0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.1990.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517-0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662-0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674-0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.2230.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.4940.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.4501.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.1680.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+00 0.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.5680.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.1510.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.4430.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为：0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs（按行）为：0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-01 0.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02 数表（x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00-0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+00 0.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+010.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.7001.500 -0.53E+00 -0.80E+00 0.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。
北航数值分析大作业 第二题 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为虚数，不需要求特征向量。
北航数值分析大作业第二题
数值分析第二次大作业史立峰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。
数值分析大作业QR分解
数值分析大作业QR分解题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵20010-1-24A=0-2131431?? ? ? ? ???为例进行具体的求解计算。
一、算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。
其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。
二、算法设计：1、原始实矩阵A Rn n∈拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)rij a i n j r r n ==+。
对于r=1,2,…,n -2执行(1) 若()(2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。
(2) 计算令()2()()1,1,1,sgn(0,sgn()=1)r r r r r rr r r r r c aa a ρ+++=-=，（若则取(3) 令-0=r n r u u ?? ?，()()()-1,2,1(,,...,)r r r Tn r r r r r r nr u a c a a ρ++=- (4) 计算(r)(r)(r)Tn-r r+1,r r r+2,r nr r Tn-r T n-r n-r n-r n-r r+1r 1u =(a -c ,a ,...,a )ρI H =I -2uu =H H =I -2u u A =HA H(5) 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。
2、拟上三角矩阵QR 分解的求原矩阵的全部特征值记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，对于k =1,2,…,n -1执行： (1) 分解k-1k-1k-1A =R Q选取旋转矩阵1P =R(2,1,θ)，使(1)1k-1A =PA ，从而使第一列次对角元(1)2,1=0a ；再选取旋转矩阵2P =R(3,2,θ)，使(2)(1)1A =PA ，从而使第一列次对角元(2)3,2=0a ……如此进行下去，最多经过n-1步，(n-1)A必然变为上三角矩阵k-1R ，即k-1n-121k-1-1=-1-1-1k-1n-12112n-1R =P P P A =P P P P P PQ （）k-1Q 必为正交矩阵，且满足k-1k-1k-1A =R Q(2) 计算k k-1k-1A =R Q(3) 上述两过程反复进行直到k A 变为近似舒尔矩阵，对角线元素即为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 三角化过程前进了一步。
北航数值分析报告大作业二
数值分析大作业（二）学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1 算法设计方案首先将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解，计算出RQ 。
要得出矩阵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 Λ执行 若()n r r i a r ir,,3,2)(Λ++=全为零，则令)()1(r r A A =+,转5；否则转2。
计算()∑+==nri r ir r 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 +-=令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0ΛΛ。
计算r r T r r h u A p /)(=r r rr 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(继续。
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。
北航 数值分析第二次大作业(带双步位移的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)，再求解方程组。
北航数值分析大作业第二题
“数值分析“计算实习大作业第二题——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. 问题的提出：（1）用幂法、反幂法求矩阵A =[a ij ]20×20的按摸最大和最小特征值，并求出相应的特征向量。
其中 a ij ={sin (0.5i +0.2j ) i ≠j 1.5cos (i +1.2j ) i =j要求：迭代精度达到10−12。
（2）用带双步位移的QR 法求上述的全部特征值，并求出每一个实特征值相应的特征向量。
2. 算法的描述：(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为：{ 任取非零向量u 0=(h 1(0)，⋯，h n (0))T|h r (k−1)|=max 1≤j≤n |h r (k−1)| y ⃑ k−1=u ⃑ k−1|h r (k−1)| u ⃑ k =Ay ⃑ k−1=(h 1(k )，⋯，h n (k ))T βk =sgn (h r (k−1))h r (k ) (k =1，2，⋯) 终止迭代的控制选用≤ε。
幂法的使用条件为n ×n 实矩阵A 具有n 个线性无关的特征向量x 1，x 2，⋯，x n ，其相应的特征值λ1，λ2，⋯，λn 满足不等式|λ1|>|λ2|≥|λ3|≥⋯≥|λn |或λ1=λ2=⋯=λm|λ1|>|λm+1|≥|λm+2|≥⋯≥|λn |幂法收敛速度与比值|λ2λ1|或|λm+1λ1|有关，比值越小，收敛速度越快。
(2) 反幂法反幂法用于计算n ×n 实矩阵A 按摸最小的特征值，其迭代格式为：{任取非零向量u 0∈R nηk−1=√u ⃑ k−1T u ⃑ k−1 y ⃑ k−1=u ⃑ k−1ηk−1⁄ Au ⃑ k =y ⃑ k−1 βk =y ⃑ k−1u ⃑ k (k =1，2，⋯) 每迭代一次都要求解一次线性方程组Au ⃑ k =y ⃑ k−1。
当k 足够大时，λn ≈1βk ，y ⃑ k−1可近似的作为矩阵A 的属于λn 的特征向量。
QR分解、SVD分解、最小二乘问题数值上机报告
QR分解、SVD分解、最小二乘问题数值上机报告练习6.13 随机构造一个可逆方阵，利用不同的方法给出它的QR分解，观察所得列向量的正交性、CPU时间和所谓的向后稳定性。
解：本题中所得列向量的正交性由||Q T Q−I||2来进行刻画。
向后稳定性由||A−QR||2来进行刻画。
实验步骤：1、首先分别编写出CGS、MGS、Householder、Givens 等不同方法进行矩阵的QR分解的程序。
2、在主程序中调用这些子程序。
其中在编写利用Householder矩阵进行矩阵的QR分解时，可以利用H矩阵作为秩一修正矩阵在计算复杂度方面的优势，即 H∗g=g−b−1(u T g)u,只需2n+1次乘除法运算。
实验数据：（1）取实验矩阵A为满秩的200*200的随机矩阵时得到的数据如下：（2）取实验矩阵A为满秩的400*400的随机矩阵时得到的数据如下：实验数据分析：1.从向后稳定性上来看，CGS和MGS给出的QR乘积均能很好地近似原始矩阵A，因此都是向后稳定的算法。
而Householder方法和Givens方法的向后稳定性略差。
2.从所求得列向量的正交性来看，CGS方法最差，MGS比CGS稍好，Householder方法和Givens方法的表现最好。
Householder方法给出的Q具有更好的列正交性。
3.从CPU时间上来看，Givens方法所耗时间最长，这是合理的，因为对应于一个Householder矩阵，Givens方法中要构造一系列Givens矩阵。
理论分析，对于稠密阵，Householder方法所需的乘除法次数是Givens方法的一半。
且由讲义所指出，有如下结果：N opt House≈∑2(n+1−k)(m+1−k)≈mn2−13n3nk=1N opt GS ≈∑2(k −1)m ≈mn 2n k=1.而由Householder 镜像变换法的计算复杂度略低于MGS 方法，Householder 方法所耗的CPU 时间略低于MGS 方法，故是合理的。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《数值分析A》计算实习题目二姓名学号数值分析第二次大作业一、算法设计方案首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），算法如课本P61。
使用QR分解法对矩阵A（n-1）进行QR分解，算法如课本P59,进而求出所得矩阵的Q、R、RQ矩阵。
然后对A（n-1）进行带双步位移的QR分解求矩阵的全部特征值，采用以下几步进行：第一步：判断是否a m,m-1(k)<=ε ,若不是,则进入第四步。
若是，则a m,m-1(k)为特征值，m=m-1，若m=1，则进入第二步，若m=2进入第三步，否则转第四步。
第二步：m=1，则a11(k)为特征值，转向结束步。
第三步：m=2,则可以求出A的两个特征值s1和s2，转向结束步。
第四步：判断是否a m-1,m-2(k)<=ε，若不是，进入第五步。
若是，则得到A的两个特征值s1和s2，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。
第五步：判断是否达到循环上限，若达到，则结束，否则进入第六步。
第六步：对A进行双步位移QR分解，这里的算法如课本P64，分解后转向计数步。
计数步：对循环次数进行计数，并转向第一步。
结束步：显示所求得的特征值。
最后对实特征值利用列主元高斯消元法求解其对应的特征向量，算法如课本p17.二、源程序代码#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;}}}}}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矩阵进行分解并传递参数矩阵Q QRdecomposition(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++){for(k=0;k<10;k++){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++){{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;}}for(k=j;k<10;k++){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;{c=c+A[i][j]*X[l][j];}X[l][i]=-c/A[i][i];}}}三、程序输出结果、1612091037 曾庆超20。