数值分析大作业QR分解
北航数值分析计算实习题目二 矩阵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 λλ---++=的两个根。
数值分析大作业
第二次计算实验：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。
数值分析7.2矩阵的正交分解与求矩阵全部特征值的QR方法
但 x y ,则存在householder阵
2
2
UU T
H I2 U 2
2
使Hx y，其中U x y。 W
x
x y
y
证：若设W U ，则有 W 1，因此
U
2
I
H 2
2
I 2WW T
(x y) x y 2
( xT
yT
I )
UU T 2 U2
2
Hx
x
2
( x2 y) x y 2
2
k
)
,
H
k
(k 2
)
,
,
H
k
(k n
)
A(k 1)
1(
k
1)
,
(k 2
1)
,
,
(k n
1)
a(2) 11 0
a(2) 1k
Hk
A(k )
Hk
0
a(k) kk
0
0
a(k) nk
a(2) 1n
a(k kn
)
a(k) nn
H
(
k1
k
)
,
H
k
(k 2
)
,
,
H
k
(k n
)
a1(12) 0 0 0
迭代格式
Ak Qk Rk Ak 1 RkQk
(k 1, 2, ).
将A A1化成相似的上三角阵（或分块上三角阵），
从而求出矩阵A的全部特征值与特征向量。
由A A1 Q1R1 ,即Q11 A R1。 于是A2 R1Q1 Q1 AQ1 ,即A2与A相似。
同理可得，Ak A (k 2, 3, )。 故它们有相同的特征值。
北航数值分析大作业二
数值分析第二题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 。
数值分析大作业之2
数值分析大作业二一、算法设计方案(一)算法流程1.将要求解的矩阵进行Householder变换，进行拟上三角化，并输出拟上三角化的结果；2.将拟上三角化后的矩阵进行带双步位移的QR分解(最大迭代次数1000次)，逐个求出特征值，输出QR法结束后得到的Q、R、RQ阵，输出求出的特征值（用实部和虚部表示）3.对于求出来的实特征值，使用带原点平移的反幂法求出对应的特征向量，输出这些特征向量。
（二）程序设计流程1. 定义精度和最大迭代次数；2. 使用容器vector进行编程（方便增减元素），使用传引用调用（提高执行效率）；3. 将各个步骤用到的数学算法封装成函数，方便调用。
具体需要的函数如下：double VectMod(vector< double > &p)：求向量的模vector< double > NumbMultiVect(vector< double > &vect, double a)：数乘向量double VectMultiVect(vector< double > &y, vector< double > &u)：求两个向量的积vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector<double > &u)：矩阵的转置乘向量vector< double > ArraMultiVect(vector< vector< double > > &B, vector< double > &u)：矩阵乘向量vector< vector< double > > VectMultiCovVect(vector< double > &a, vector< double >&b)：向量乘向量转置得到矩阵void ArraSubs(vector< vector< double > > &C, vector< vector < double > > D) ：两个矩阵相减Vector< double > VectSubs(vector< double > &a, vector< double > &b)：两个向量相减vector<vector<double>> GetM(vector<vector<double>> &A)：求解矩阵M （用于双部位移QR迭代用）double GetMode(vector < vector < double > > &B, const int r)：求解矩阵B的r列向量的模double GetModeH(vector < vector < double > > &B, const int r)：求解矩阵B的第r列向量的模，用于拟对角化vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)：一个实数乘矩阵bool IsBirZeroH(vector< vector< double > > &B, const int r)：判断B[i][r]对角线下是否为零void GausElim(vector< vector< double > > a)：列主元高斯消元法求齐次方程解向量void Stop(vector< vector< double > > &Ar)：停止，结束程序void SolutS1S2(complex< double > &s1, complex< double > &s2, vector< vector<double > > &A)：求解二阶子阵的特征值s1，s2；void Save2(complex< double > &s1, complex< double > &s2)：保存两个特征值void Save1(complex< double > &s)：保存一个特征值void JudgemBelow2(vector< vector< double > > &A, vector< vector< double > > Abk)：对于m == 1 及m == 0 的处理void Hessenberg(vector< vector< double > > &A)：矩阵拟上三角化void QRMethod(vector< vector< double > > A)：矩阵QR分解void CalculatAk(vector< vector < double > > &Ak)：带双步位移QR迭代法二、源程序#include "stdafx.h"#include <vector>#include <iostream>#include <math.h>#include <complex>#include <fstream>using namespace std;const double epsion = 1e-12;const int L = 1000;int m,n;int k = 1;vector< vector< double > > I;vector< complex< double > > Lambda;///////////////////////////以下为自定义的算法流程中用到的函数double VectMod(vector< double > &p) //求向量的模{double value = 0.0;vector< double >::size_type i,j;j = p.size();for (i=1; i<j; i++){value += p[i] * p[i];}value = sqrt(value);return value;}vector< double > NumbMultiVect(vector< double > &vect, double a) //数乘向量{int j = vect.size();vector< double > b(j, 0);for (int i=1; i<j; i++){b[i] = a * vect[i];}return b;}double VectMultiVect(vector< double > &y, vector< double > &u)//两个向量相乘{vector< double >::size_type a = y.size();double value = 0;for (vector< double >::size_type i=1; i<a; i++){value += y[i] * u[i];}return value;}vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector< double > &u)//矩阵的转置乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[j][i] * u[j];}}return vec;}}vector< double > ArraMultiVect(vector< vector< double > > &B,vector< double > &u) //矩阵乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[i][j] * u[j];}}return vec;}}vector< vector<double> > GetM(vector< vector <double> > &A){int a = A.size();double s = A[a-2][a-2] + A[a-1][a-1];double t = A[a-2][a-2] * A[a-1][a-1] - A[a-1][a-2] * A[a-2][a-1];vector<vector<double>> D(a, vector< double >(a, 0));for (int i=1; i<a; i++){{double sum = 0;for (int k=1; k<a; k++){sum += A[i][k] * A[k][j];}D[i][j] = sum - s * A[i][j] + t *(i==j ? 1.0 : 0);}}return D;}double GetMode(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}double GetModeH(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r+1; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)//数乘//向量{int b = D.size();vector< vector< double > > U(b, vector< double >(b, 0));for (int i=1; i<b; i++){{U[i][j] = a *D[i][j];}}return U;}bool IsBirZero(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+1; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}bool IsBirZeroH(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+2; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}vector< vector< double > > VectMultiCovVect(vector< double > &a,vector< double > &b)//向量乘向量的转置得到矩阵{int s1 = a.size();int s2 = b.size();if (s1 != s2){cerr << "Vectors not match in size ! ";}else{vector< vector< double > > U(s1, vector< double >(s1, 0));for (int i=1; i<s1; i++){for (int j=1; j<s1; j++){U[i][j] = a[i] * b[j];}}return U;}}void ArraSubs(vector< vector< double > > &C,vector< vector < double > > D){int a = C.size();int b = D.size();int c = C[0].size();int d = D[0].size();if (a!=b || c!=d){cerr << "Vectors not match in size !";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){C[i][j] -= D[i][j];}}}}vector< double > VectSubs(vector< double > &a,vector< double > &b){int s1 = a.size();int s2 = b.size();if(s1 != s2){cerr << "Vectors not match in size !";}else{vector< double > value(s1,0);for (int i=1;i<s1;i++){value[i] = a[i] - b[i];}return value;}}void GausElim(vector< vector< double > > a) //高斯消元{int sz = a.size();vector< int > fx, ufx,ufxp, record;vector< double > x(sz, 0);vector< vector< double > > ret;//*****消元过程********for (int k = 1; k < sz-1; k++){double max = a[k][k];int p=0;for (int i=k+1; i<sz; i++){if (abs(a[i][k]) > abs(max)){p =i;max = a[i][k];}}//选出主元行if (abs(max) >= epsion){if (p != 0){for (int j=k; j<sz; j++){double temp = a[k][j];a[k][j] = a[p][j];a[p][j] = temp;}//交换主元行}for (int i=k+1; i<sz; i++){double tt = a[i][k] / a[k][k];for (int j=k+1; j<sz; j++){a[i][j] = a[i][j] - tt * a[k][j];}}}// end of if (abs(max) >= epsion)}// end of for (int k = 1; k < sz; k++)for (int i=1; i<sz; i++){int p=0;for (int j=1; j<=i;j++){if(abs(a[j][i]) >= 10*epsion)//不为零{p=j;}}record.push_back(p);}for (int i=1; i<sz; i++){int p=0;for (int j=sz-2; j>=0; j--){if (record[j] == i){p=j + 1;}}if (p != 0){ufxp.push_back(p);ufx.push_back(i);}}for (int i=1; i<sz; i++){int p = 0;for (int j=0; j < ufxp.size(); j++){if (ufxp[j] == i){p = 1;}}if (p == 0){fx.push_back(i);}}//end of for (int i=1; i<sz; i++)int c = fx.size();//***************************************************** for (int i=0; i<c; i++){for (int j=0; j<c; j++){if (i == j){x[fx.at(j)] = 1.0 + epsion;}else{x[fx.at(j)] = 0;}}int b = ufxp.size();for (int s=b-1; s>=0; s--){double temp=0;for(int t=ufxp.at(s)+1; t<sz; t++){temp += x[t] * a[ufx.at(s)][t];}x[ufxp.at(s)] = (0 - temp) / a[ufx.at(s)][ufxp.at(s)];}ret.push_back(x);}//end of for (int i=0; i<c; i++)//******************************************************* int sz1 = ret.size();int sz2 = ret[0].size();ofstream result2("CharactVector .txt", ios::app);result2.precision(12);for (int i=0; i<sz1; i++){for (int j=1;j<sz2; j++){result2 << scientific << ret[i][j] << endl;}result2 << endl << endl;}result2 << "下一个特征值的特征向量！" << endl;}void Stop(vector< vector< double > > &Ar){int a = Lambda.size();for (int i=0; i<a; i++){if (abs(Lambda[i].imag())<=epsion){vector< vector< double > > An=Ar;ArraSubs(An,NumbMultiArra(I,Lambda[i].real()));GausElim(An);}}ofstream result("result.txt");result.precision(12);vector<complex< double > >::iterator i,j;j = Lambda.end();for (i = Lambda.begin(); i<j; i++){result << scientific <<(*i) << endl;}exit(0);}void SolutS1S2(complex< double > &s1, complex< double > &s2,vector< vector< double > > &A){double s = A[m-1][m-1] + A[m][m];double t = A[m-1][m-1] * A[m][m] - A[m][m-1] * A[m-1][m];double det = s *s - 4.0 * t;if (det > 0){s1.imag(0);s1.real ((s + sqrt(det)) / 2);s2.imag (0);s2.real((s - sqrt(det)) / 2);}else{s1.imag(sqrt(0 - det) / 2);s1.real(s / 2);s2.imag(0 - sqrt(0 - det) / 2);s2.real(s / 2);}}void Save2(complex< double > &s1, complex< double > &s2)//存储特征值s1，s2 {Lambda.push_back(s1);Lambda.push_back(s2);}void Save1(complex< double > &s){Lambda.push_back(s);}void JudgemBelow2(vector< vector< double > > &A,vector< vector< double > > Abk){if (m == 1){complex< double > lbd;lbd.imag(0);lbd.real(A[1][1]);Save1(lbd);Stop(Abk);}else if (m == 0){Stop(Abk);}}//拟上三角化void Hessenberg(vector< vector< double > > &A){int a = A.size();vector< double > u(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;for (int r=1; r<a-2; r++){if (!IsBirZeroH(A, r)){d = GetModeH(A,r);c = (A[r+1][r] > 0) ?(0-d) : d;h = c * c - c * A[r+1][r];for (int j=1; j<a; j++){if (j < r+1){u[j] = 0;}else if (j == r+1){u[j] = A[j][r] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)p = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);q = NumbMultiVect(ArraMultiVect(A, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(A, VectMultiCovVect(w, u));ArraSubs(A, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("NISHANGDANJIAO.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}}void QRMethod(vector< vector< double > > A){int a = A.size();vector< vector< double > > Q(a,vector < double > (a,0));for (int i=1; i<a; i++){Q[i][i] = 1.0;}vector< vector< double > > RQ(A);//vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);vector< double > l(a,0);double d, c, h, t;for (int r=1; r<a-1; r++){if (!IsBirZero(A, r)){d = GetMode(A,r);c = (A[r][r] > 0) ?(0-d) : d;h = c * c - c * A[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = A[j][j] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)l = ArraMultiVect(Q, u);ArraSubs(Q, VectMultiCovVect(l,NumbMultiVect(u, 1/h)));v = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);ArraSubs(A, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(RQ, u), 1 / h);q = NumbMultiVect(ArraMultiVect(RQ, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(RQ, VectMultiCovVect(w, u));ArraSubs(RQ, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("QR.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << Q[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << RQ[i][j] << " ";}result << endl;}}void CalculatAk(vector< vector < double > > &Ak){int a = Ak.size();vector< vector < double > > B(a,vector < double > (a,0));vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;B = GetM(Ak);for (int r=1; r<a-1; r++){if (!IsBirZero(B, r)){d = GetMode(B,r);c = (B[r][r] > 0) ?(0-d) : d;h = c * c - c * B[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = B[j][j] - c;}else{u[j] = B[j][r];}}// end of for (int j=1; j<=m; j++)v = NumbMultiVect(ConveArraMultiVect(B, u), 1 / h);ArraSubs(B, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(Ak, u), 1 / h);q = NumbMultiVect(ArraMultiVect(Ak, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(Ak, VectMultiCovVect(w, u));ArraSubs(Ak, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)}int _tmain(int argc, _TCHAR* argv[]){vector< vector< double > > A(11, vector< double >(11,0));vector< vector< double > > Abk;I=A;for (int i=1; i<11; i++){vector< double > temp;for (int j=1; j<11; j++){if(i != j){A[i][j] = sin(0.5 * i + 0.2 * j);I[i][j] = 0;}else{A[i][j] = 1.5 * cos(i + 1.2 *j);I[i][j] = 1.0;}}}Abk = A;n = A.size() - 1;m = n;//初始化问题Hessenberg(A);QRMethod(A);while (1){if (abs(A[m][m-1]) <= epsion){complex< double > lbdm;lbdm.imag(0);lbdm.real(A[m][m]);Save1(lbdm);m -= 1;//对A进行降维处理！！！！A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();}JudgemBelow2(A, Abk);}else{complex< double > va1, va2;SolutS1S2(va1, va2, A);if (m == 2){Save2(va1,va2);Stop(Abk);}//end of if (m == 2)if ( abs(A[m-1][m-2]) <= epsion){Save2(va1,va2);m = m - 2;//矩阵降维A.pop_back();A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();A[i].pop_back();}JudgemBelow2(A, Abk);}else{if (k == L){cerr << "Stop without solution";exit(-1);}else{CalculatAk(A);k += 1;}}}// end of if (abs(A[m][m-1]) >= epsion)}}三、实验结果（1）A经过拟上三角化后得到的矩阵-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000-7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000-8.782436382928e-002 -9.255889387184e-001 6.032599440534e-0011.518860956469e-001-2.347878362416e+000 2.372370104937e+000 1.819290822208e+0003.237804101546e-001 2.205798440320e-001 2.102692662546e+0001.816138086098e-001 1.278839089990e+000 -6.380578124405e-001-4.154075603804e-0011.0556********e-016 1.728274599968e+000 -1.171467642785e+000-1.243839262699e+000 -6.399758341743e-001 -2.002833079037e+0002.924947206124e-001 -6.412830068395e-001 9.783997621285e-0022.557763574160e-001-5.393383812774e-017 -3.165873865181e-017 -1.291669534130e+000-1.111603513396e+000 1.171346824096e+000 -1.307356030021e+0001.803699177750e-001 -4.246385358369e-001 7.988955239304e-0021.608819928069e-0011.533464996622e-017 5.963321406181e-017 0.000000000000e+0001.560126298527e+000 8.125049397515e-001 4.421756832923e-001-3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-0021.300562737229e-016 -3.097060010889e-017 0.000000000000e+0000.000000000000e+000 -7.707773755194e-001 -1.583051425742e+000-3.042843176799e-001 2.528712446035e-001 -6.709925401449e-0012.544619929082e-0011.610216724767e-016 -2.211571837369e-016 0.000000000000e+0000.000000000000e+000 6.483933100712e-017 -7.463453456938e-001-2.708365157019e-002 -9.486521893682e-001 1.195871081495e-0011.929265617952e-0021.368550186199e-016 7.151513190805e-017 0.000000000000e+0000.000000000000e+000 -2.522454384246e-017 1.072074718358e-016-7.701801374364e-001 -4.697623990618e-001 4.988259468008e-0011.137691603776e-001-2.780851300718e-017 -6.708630788363e-017 0.000000000000e+0000.000000000000e+000 -3.282796821065e-017 4.879774287475e-0171.854829286220e-016 7.013167092107e-001 1.582180688475e-0013.862594614233e-001-2.124604440055e-017 -1.707979758930e-016 0.000000000000e+0000.000000000000e+000 1.013496750890e-016 -4.153758325241e-0171.222621691887e-016 -5.551115123126e-017 4.843807602783e-0013.992777995177e-001（2）Q-3.519262579534e-001 4.427591982245e-001 -6.955982513607e-0016.486200753651e-002 3.709718861896e-001 1.855847143605e-001-1.628942319628e-002 -1.181053169648e-001 -5.255375383724e-002 -5.486582943568e-002-9.360277287361e-001 -1.664679186545e-001 2.615299548560e-001 -2.438671728934e-002 -1.394774360893e-001 -6.977585391241e-0026.124472142963e-003 4.440505443139e-002 1.975907909728e-0022.062836970533e-002-4.208697095111e-017 -8.810520554692e-001 -3.989762796959e-0013.720308728479e-002 2.127794064090e-001 1.064463557221e-001-9.343171079758e-003 -6.774200464527e-002 -3.014340698675e-002 -3.146955080444e-002-2.150178169911e-017 4.009681353156e-017 -5.371806806439e-001 -1.234945854205e-001 -7.063151608719e-001 -3.533456368496e-0013.101438948264e-002 2.248676491598e-001 1.000601783527e-0011.044622748702e-0016.113458775639e-018 -3.721194648197e-0177.068175586628e-0189.892235468621e-001 -1.239411731211e-001 -6.200358589825e-0025.442272839461e-003 3.945881637235e-002 1.755813350011e-0021.833059462907e-0025.184948268593e-017 -4.198303559531e-017 3.255071298412e-017-1.527665297025e-017 5.323610690264e-001 -6.733900344896e-0015.910581205868e-002 4.285425323867e-001 1.906901343193e-0011.990794495295e-0016.419444583601e-017 4.121668945107e-017 3.096964582901e-017-5.050360323651e-017 -7.078737686239e-017 -6.0597********e-001 -9.165783032818e-002 -6.645586508974e-001 -2.957110877580e-001 -3.0872********e-0015.455993559780e-017 -9.724896332186e-017 3.956603694870e-0171.913857001710e-018 -4.316583456405e-0172.797376665557e-0179.933396625117e-001 -9.690440311939e-002 -4.311990584470e-002-4.501694411183e-002-1.108640877071e-017 4.655237056115e-017 -1.0722********e-017 -9.470236121745e-018 4.277993227986e-017 8.866601870176e-017 -2.516725033065e-016 5.410088006061e-001 -5.817540838226e-001 -6.0734********e-001-8.470152033092e-018 9.650816729410e-017 -1.429362211392e-017 -2.789222052040e-017 -3.690793770141e-017 -8.090765462521e-017 -1.964050295697e-016 -6.772274683437e-017 -7.221591336735e-0016.917269588876e-001（3）R2.508342744917e+000 -2.185646885493e+000 -1.314609070786e+000-3.558787493836e-002 -2.609857850388e-001 -1.283121847090e+000 -1.390878610606e-001 -8.712897972161e-001 3.849367902971e-0013.353802899665e-0012.100627753398e-016 -1.961603277854e+000 2.407523727633e-0017.054714572823e-001 5.957204318279e-001 5.526978774676e-001-3.268209924413e-001 -5.769498668364e-002 2.871129330189e-001 -8.895128754189e-002-3.300197935770e-016 -1.872873410421e-016 2.404534601993e+0001.706758096328e+000 -4.239566704091e-001 3.405332305815e+000-1.050017655852e-001 1.462257102734e+000 -6.684487469283e-001 -4.027*********e-0013.0773********e-017 1.746386351950e-017 0.000000000000e+0001.577122080722e+000 6.399535133956e-001 3.468127872427e-001-5.701786649768e-002 4.014788054433e-001 -2.222476176311e-001 -6.317059236442e-0021.760039865880e-016 9.988285339980e-017 0.000000000000e+0000.000000000000e+000 -1.447846997770e+000 -1.415724007744e+000-2.806139044665e-001 -2.817910521892e-001 -4.611434881851e-0021.996629079956e-0018.804885435596e-017 4.996802050991e-017 0.000000000000e+0000.000000000000e+000 8.831633340975e-017 1.231641451542e+0001.619701003419e-001 1.962638275504e-001 5.350035621760e-001-1.509273424767e-001-7.728357669400e-018 -4.385868928757e-018 0.000000000000e+0000.000000000000e+000 -7.751835246841e-018 -1.279231078565e-017-7.753441914209e-001 -3.464514508821e-001 4.312226803504e-0011.234643696237e-001-5.603391361152e-017 -3.179943413314e-017 0.000000000000e+0000.000000000000e+000 -5.620413613517e-017 -9.274974944455e-0170.000000000000e+000 1.296312940612e+000 -4.288053318338e-0012.737334158165e-001-2.493361499851e-017 -1.414991023727e-017 0.000000000000e+0000.000000000000e+000 -2.500935953597e-017 -4.127119443934e-0170.000000000000e+000 0.000000000000e+000 -6.707396440648e-001-4.842320121884e-001-2.603055667460e-017 -1.477242832192e-017 0.000000000000e+0000.000000000000e+000 -2.610963355436e-017 -4.308689959101e-0170.000000000000e+000 0.000000000000e+000 -1.110223024625e-0167.168323926323e-002（4）RQ1.163074414164e+0002.632670934508e+000 -1.772796003272e+000-8.668899138521e-002 3.300503471047e-001 1.455162371214e+000 -9.730650448593e-001 -4.873031174655e-001 -7.756411630489e-001 -3.249201979113e-0011.836115060851e+000 1.144286420080e-001 -9.880381403133e-0015.589725694767e-001 4.694190067101e-002 -2.978478237007e-0011.617130577649e-002 6.936977702522e-001 1.367670571405e-0011.419099231519e-0025.598200524418e-016 -2.118520153533e+000 -1.876189745783e+000-5.407071940597e-001 1.171538359721e+000 -2.550323020223e+0001.691577936540e+000 1.229951613262e+000 1.387947777212e+0008.667502917242e-001-3.179059303049e-017 -5.261714527977e-017 -8.471995127808e-0014.382910468318e-001 -1.008632199185e+000 -7.959374261495e-0014.769258865577e-001 4.072683083890e-001 4.096390493527e-0013.363378940862e-001-3.514195999269e-016 3.391949386044e-017 -2.160938214545e-016 -1.432244342447e+000 -5.742284908055e-001 1.213151477723e+000 -3.457508625575e-001 -4.749853573124e-001 -3.176158274191e-001 -4.294507015032e-002-3.704735750656e-017 -1.0998********e-016 -1.953241363386e-0178.269089741494e-017 6.556779598004e-001 -9.275250974463e-0012.529079844053e-001 6.905949216976e-001 -2.374430675823e-002-2.429781119781e-001-6.420051822287e-017 3.865817713597e-017 -3.806718328740e-0172.129734126613e-017 7.853*********e-017 4.698400884876e-001-2.730776009527e-001 7.821296259798e-001 -9.580964936399e-0027.846239841323e-0021.478496020372e-016 -8.383952267131e-017 9.577540382744e-017-7.120338911078e-017 -1.220876056602e-016 -1.854471832043e-0161.287679058937e+000 -3.576058900348e-001 -4.116725408807e-0033.914268216423e-0014.477158378610e-017 -6.204017427568e-017 3.360322998507e-017-1.133882337819e-017 -2.759056827456e-017 -9.217582575819e-0172.214639994444e-016 -3.628760503545e-001 7.398980975354e-0017.241608309576e-0023.408891482791e-017 2.353495494255e-017 1.932283926688e-017-3.456910153081e-017 -2.015177337156e-017 -8.084422691634e-017 -5.839476980893e-017 5.551115123126e-017 -5.176670596524e-0024.958522909877e-002（5）特征值(6.360627875745e-001,0.000000000000e+000)(5.650488993501e-002,0.000000000000e+000)(9.355889078188e-001,0.000000000000e+000)(-9.805309562902e-001,1.139489127430e-001)(-9.805309562902e-001,-1.139489127430e-001)(1.577548557113e+000,0.000000000000e+000)(-1.484039822259e+000,0.000000000000e+000)(-2.323496210212e+000,8.930405177200e-001)(-2.323496210212e+000,-8.930405177200e-001)(3.383039617436e+000,0.000000000000e+000)（6）实特征值的特征向量4.745031936539e+0003.157868541750e+0001.729946912420e+001-1.980049331455e+000-3.187521973524e+0017.794009603201e+000-1.004255685845e+0011.670757770480e+0011.310524273054e+0011.000000000001e+000下一个实特征值对应的特征向量：-5.105003830625e+000-4.886319842360e+0009.505161576151e+000-6.788331788214e-001-9.604334110499e+000-3.0457********e+0001.574873885602e+001-7.395037126277e+000-7.109654943661e+0001.000000000001e+000下一个实特征值对应的特征向量：2.792418944529e+0001.598236841512e+000-5.207507040911e-001-1.667886451702e+000-1.225708535859e+0017.241214790777e+000-5.398214501433e+0002.841008912974e+001-1.216518754416e+0011.000000000001e+000下一个实特征值对应的特征向量：6.217350824581e-001-1.115111815226e-001-2.483773580804e+000-1.306860840421e+000-3.815605442533e+0008.117305509395e+000-1.239170883675e+000-6.800309586197e-0012.691900144840e+0001.000000000001e+000下一个实特征值对应的特征向量：-1.730784582112e+0012.408192534967e+0014.133852365119e-001-8.572044074538e+0009.287334657928e-002-7.832726042776e-002-6.374274025716e-001-3.403204760832e-001。
北航数值分析大作业 第二题 QR分解
《数值分析B》课计算实习第一题设计文档与源程序姓名：杨彦杰学号：SY10171341 算法的设计方案（1）运行平台操作系统：Windows XP；开发平台：VC6.0++；工程类型：文档视图类；工程名：Numanalysis；（2）开发描述首先新建类CMetrix，该类完成矩阵之间的相关运算，包括相乘、加减等，以主程序方便调用；题目的解算过程在视图类CNumanalysisView中实现，解算结果在视图界面中显示；（3）运行流程（4）运行界面2、全部源代码（1）类CMetrixMetrix.h文件：class CMetrix{public:double** MetrixMultiplyConst(double**A,int nRow,int nCol,double nConst);//矩阵乘常数double** MetrixMultiplyMetrix(double**A,double**mA,int nRow,int nCol);//矩阵相乘double** MetrixSubtractMetrix(double **A, double **subA, int nRow,int nCol);//矩阵减矩阵double VectorMultiplyVector(double*V,double*mulV,int nV);//向量点积double** VectorMultiplyVectortoMetrix(double*V,double*VT,int nV);//向量相乘为矩阵double* VectorSubtractVector(double*V,double*subV,int nV);//向量相减double* VectorMultiplyConst(double *V, int nV, double nConst);//向量乘常数double LengthofVector(double *V,int nV);//求向量的长度double* MetrixMultiplyVector(double**A,int nRow,int nCol,double*V,int nV);//矩阵与向量相乘double** AtoAT(double **A,int Row,int Col);//矩阵转置运算void FreeMem();CMetrix(int nRow,int nCol);uCMetrix();virtual ~CMetrix();double* vector; //过渡向量double** B; //过渡矩阵};Metrix.cpp文件：CMetrix::CMetrix(int nRow, int nCol){B = new double*[nRow];for (int i = 0;i < nCol;i++){B[i] = new double[nCol];}vector = new double[nRow];}CMetrix::~CMetrix(){delete vector;B = NULL;delete B;}double** CMetrix::AtoAT(double **A, int nRow, int nCol){for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[col][row] = A[row][col];}}return B;}double* CMetrix::MetrixMultiplyVector(double **A, int nRow, int nCol, double *V, int nV) {if (nCol != nV){AfxMessageBox("矩阵列数和向量维数不等,不能相乘!");return 0;}double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){sum += A[row][col]*V[col];}vector[row] = sum;sum = 0.0;}return vector;}double CMetrix::LengthofVector(double *V, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*V[col];}return length;}double* CMetrix::VectorMultiplyConst(double *V, int nV, double nConst){for (int col = 0;col < nV;col++){vector[col] = V[col]*nConst;}return vector;}double* CMetrix::VectorSubtractVector(double *V, double *subV, int nV){for (int col = 0;col < nV;col++){vector[col] = V[col]-subV[col];}return vector;}double** CMetrix::VectorMultiplyVectortoMetrix(double*V, double *VT, int nV){for (int row = 0;row < nV;row++){for (int col = 0;col < nV;col++){B[row][col] = V[row]*VT[col];}}return B;}double CMetrix::VectorMultiplyVector(double *V, double *mulV, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*mulV[col];}return length;}double** CMetrix::MetrixSubtractMetrix(double **A, double **subA, int nRow, int nCol) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]-subA[row][col];}}return B;}double** CMetrix::MetrixMultiplyMetrix(double **A, double **mA, int nRow, int nCol) {double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){for(int n = 0;n < nCol;n++){sum += A[row][n]*mA[n][col];}B[row][col] = sum;sum = 0.0;}}return B;}double** CMetrix::MetrixMultiplyConst(double **A, int nRow, int nCol, double nConst) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]*nConst;}}return B;}（2）类CNumanalysisViewNumanalysisview.hclass CNumanalysisView : public CEditView{…………public:double Sign(double x);void DisplayVector(double*V,int nV); // 显示向量数据void DisplayMetrix(double **A,int Row,int Col); //显示矩阵void DisplayText(CString str); //显示文本protected://{{AFX_MSG(CNumanalysisView)afx_msg void OnQRanalyze(); //运行主函数…………};Numanalysisview.cppvoid CNumanalysisView::OnQRanalyze(){//开辟空间int nRow = 10;int nCol = 10;CString str;CMetrix Metrix(nRow,nCol);double tempa = 0.0;double *V = new double[nCol]; //分配10*10矩阵空间double *ur = new double[nCol];double *pr = new double[nCol];double *qr = new double[nCol];double *wr = new double[nCol];double *tempV = new double[nCol];double **Ar = new double*[nRow];double **C = new double*[nRow];double **Cr = new double*[nRow];double **tempA = new double*[nRow];double **A = new double*[nRow];double **R = new double*[nRow];for (int col = 0;col < nRow;col++){A[col] = new double[nCol];Ar[col] = new double[nCol];C[col] = new double[nCol];Cr[col] = new double[nCol];tempA[col] = new double[nCol];R[col] = new double[nCol];}//矩阵A求解for (int i = 0;i < nRow;i++){for (int j = 0;j < nCol;j++){if(i == j)A[i][j] = 1.5*cos((i+1.0)+1.2*(j+1.0));elseA[i][j] = sin(0.5*(i+1.0)+0.2*(j+1.0));}}//--------------------拟上三角化-------------------------// double dr = 0.0,cr = 0.0,hr = 0.0,tr = 0.0;for (int r = 0;r < nCol - 2;r++){dr = 0.0;for (i = r+1;i < nCol;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+2;i < nCol;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r+1][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r+1][r])*dr;hr = cr*cr - cr*A[r+1][r]; //hrstr.Format("dr = %.6e, cr = %.6e, hr = %.6e",dr,cr,hr);for (int row = 0;row < nRow;row++) //ur{if (row < r+1)ur[row] = 0.0;else if (row == r+1)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,nRow,nCol,ur,nCol); //pr memcpy(pr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(pr,nCol,1.0/hr);memcpy(pr,tempV,nCol*8);tempV = Metrix.MetrixMultiplyVector(A,nRow,nCol,ur,nCol); //qr memcpy(qr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(qr,nCol,1.0/hr);memcpy(qr,tempV,nCol*8);tr = Metrix.VectorMultiplyVector(pr,ur,nCol)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,nCol,tr); //wr memcpy(wr,tempV,nCol*8);tempV = Metrix.VectorSubtractVector(qr,wr,nCol);memcpy(wr,tempV,nCol*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,nCol); //Arfor (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)A[row][col] = tempA[row][col];}tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}}DisplayText("矩阵A拟上三角化后所得的矩阵为：");DisplayMetrix(A,nRow,nCol);for (int row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)R[row][col] = A[row][col];}// -------------------------------------------------////--------------------带双步位移的QR分解-------------------------// int m = nCol;struct EigenVal //定义特征值结构，实数和虚数{double Realnum;double Imagnum;};EigenVal *eigenvalue = new EigenVal[m];EigenVal tmpEigen1,tmpEigen2;double b = 0.0,c = 0.0,delta = 0.0,s = 0.0,t = 0.0;double *vr = new double[m];for (int k = 1;k < 100; k++){//m代表矩阵阶数，判断式中直接用，运算中需要-1while (m > 1 && fabs(A[m-1][m-2]) <= IPSLEN)//第三步和第四步{eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;m = m - 1;}if (m == 1){eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;DisplayText("已求出A的全部特征值：");break;}b = -(A[m-2][m-2]+A[m-1][m-1]); //第五步求一元二次方程式的根s1,s2c = A[m-2][m-2]*A[m-1][m-1]-A[m-2][m-1]*A[m-1][m-2];delta =b*b - 4*c;if (delta >= 0.0){tmpEigen1.Realnum = (-b-sqrt(delta))/2;tmpEigen1.Imagnum = 0.0;tmpEigen2.Realnum = (-b+sqrt(delta))/2;tmpEigen2.Imagnum = 0.0;}else{tmpEigen1.Realnum = -b/2;tmpEigen1.Imagnum = -sqrt(fabs(delta))/2 ;tmpEigen2.Realnum = -b/2;tmpEigen2.Imagnum = sqrt(fabs(delta))/2;}if (m == 2) //第六步 m=2时结束运算{eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;DisplayText("已求出A的全部特征值：");break;}else //第七步 m > 1{if (fabs(A[m-2][m-3]) <= IPSLEN){eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;m = m - 2;continue;}}for (int row = 0;row < m;row++) //Mk求之前需要把A付给C{for (int col = 0;col < m;col++)C[row][col] = A[row][col];}double **I = new double*[m]; //第九步求Mk和Mk的QR分解for (int i = 0;i < m;i++) //求单位矩阵I,分配m*m矩阵空间{I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}s = A[m-2][m-2]+A[m-1][m-1];t = A[m-2][m-2]*A[m-1][m-1] - A[m-2][m-1]*A[m-2][m-1];tempA = Metrix.MetrixMultiplyMetrix(A,A,m,m);//A*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixMultiplyConst(A,m,m,s);//s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(Ar,A,m,m);//A*A-s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col]; }tempA = Metrix.MetrixMultiplyConst(I,m,m,-1*t);//-t*Ifor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m);//A*A - s*A + r*I for (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}delete I;//Mk的QR分解for (int r = 0;r < m - 1;r++){dr = 0.0;for (i = r;i < m;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+1;i < m;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r][r])*dr;hr = cr*cr - cr*A[r][r]; //hrfor (int row = 0;row < m;row++) //ur{if (row < r)ur[row] = 0.0;else if (row == r)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,m,m); //Btfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,m,m,ur,m); //Bt*ur memcpy(vr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(vr,m,1.0/hr); //vr = Bt*ur/hr memcpy(vr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(ur,vr,m);//Ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m); //Br-ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}tempA = Metrix.AtoAT(C,m,m); //Ctfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempV = Metrix.MetrixMultiplyVector(Cr,m,m,ur,m); //pr memcpy(pr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(pr,m,1.0/hr);memcpy(pr,tempV,m*8);tempV = Metrix.MetrixMultiplyVector(C,m,m,ur,m); //qr memcpy(qr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(qr,m,1.0/hr);memcpy(qr,tempV,m*8);tr = Metrix.VectorMultiplyVector(pr,ur,m)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,m,tr); //wr memcpy(wr,tempV,m*8);tempV = Metrix.VectorSubtractVector(qr,wr,m);memcpy(wr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,m);//Cr+1for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)C[row][col] = tempA[row][col]; }tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++){C[row][col] = tempA[row][col];if (fabs(C[row][col]) < IPSLEN){C[row][col] = 0.0;}}}}str.Format("矩阵A%d QR分解结束后所得到的矩阵为：",m);//计算结果输出DisplayText(str);DisplayMetrix(A,m,m);for (row = 0;row < m;row++) //Mk的QR分解后需要把C付给A{for (col = 0;col < m;col++)A[row][col] = C[row][col];}str.Format("迭代完成后的矩阵A%d = ",k);DisplayText(str);DisplayMetrix(A,m,m);}DisplayText("矩阵A的全体特征值如下: ");for (i = 0;i<nCol;i++){str.Format("%.6e + j%.6e",eigenvalue[i].Realnum,eigenvalue[i].Imagnum);DisplayText(str);}// -------------------------------------------------//求实特征值的特征向量，在拟上三角矩阵基础上直接求解即可////(A-egiI)X = 0.0;m = nRow;for (row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)A[row][col] = R[row][col];}double **I = new double*[m]; //求单位矩阵I,分配m*m矩阵空间double sum = 0.0;for (i = 0;i < m;i++){I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}for (i = 0;i < nRow;i++){if (eigenvalue[i].Imagnum != 0.0){str.Format("特征值%.6e+j%.6e为虚数，不需要求特征向量。
数值分析大作业QR分解
数值分析大作业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 的近似特征值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目:设计求取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 的近似特征值。
3、 带位移的QR 方法求实矩阵A Rn n⨯∈全部特征值一般情况下，QR 分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR 算法所得矩阵序列k A 的收敛。
0A =A （A 是原始矩阵的拟上三角矩阵）(1) 分解k-1-1k-1k-1A -u I=R k Q ，其中-1u k 即位移量，一般取A 的某一特征值的近似值；实际计算通常取为k-1A 的右下角元素(k-1)n,n a ，或取为右下角22⨯矩阵特征之中接近(k-1)n,n a 者。
(2) 计算k k-1k-1-1A =+R u I k Q 。
(3) 重复过程（1）（2）直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。
4、 反幂法求实矩阵A Rn n⨯∈的特征向量记通过QR 分解得到A 的某一特征值的近似值为p ，反幂法步骤如下： (1) 三角分解A -pI =LR 。
(2) 对k=1，2，…，做①当k=1时，令Tu =(1 11) 当k 1≠时，解-1Lu=z k②解Ry =u k③求-1y k 绝对值最大的元素k m ④令k k k z =y /m⑤当k k-1m -m 或k k-1z -z 小于规定的误差界时，停止计算，k z 即为所求的特征向量，k p +1/m 即为对应特征值的更为精确的取值。
三、 程序设计：1、 对生成的矩阵A 进行拟上三角化利用hohd 函数对A 进行householder 变换，得到A 的拟上三角矩阵。
2、 对拟上三角化后的矩阵进行QR 分解和带位移的QR 分解，求矩阵的全部特征值在绝对误差界为-41102⨯的条件下，利用qrfenjie 函数对拟上三角矩阵进行QR 分解，利用dp_qrfenjie 函数进行带位移的QR 分解，比较两者的收敛速度。
3、 反幂法求更精确的特征值和特征向量利用vect 函数和（2）得到的特征值的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为-61102⨯。
4、 输出相关结果。
四、 源程序：1、 hohd 函数function[A]=hohd(a) n=length(a); for i=1:n-2b=a(i+1:n,i); if b(1)>=0c=-sqrt(sum(b.^2)); elsec=sqrt(sum(b.^2)); endrho=sqrt(2*c*(c-b(1))); u1=(b-c*eye(n-i,1))/rho; H1=eye(n-i)-2*u1*u1'; H=eye(n-i);H(i+1:n,i+1:n)=H1; a=H*a*(H'); endA=a2、qrfenjie函数function A=qrfenjie(a)n=length(a);A=a;for i=1:100theta(n-1)=0;c(n-1)=0;s(n-1)=0;Q=1;for j=1:n-1theta(j)=atan(A(j+1,j)/A(j,j));c(j)=cos(theta(j));s(j)=sin(theta(j));P=eye(n);P(j,j)=c(j);P(j+1,j)=-s(j);P(j,j+1)=s(j);P(j+1,j+1)=c(j);invP=eye(n);invP(j,j)=c(j);invP(j+1,j)=s(j);invP(j,j+1)=-s(j);invP(j+1,j+1)=c(j);Q=Q*invP;R=P*A;A=R;endA=R*Q;tor=max(abs(diag(A)-diag(a)));if tor>5*10^-5a=A;elsebreak;endendi,Ak=A,lanmda=diag(A)'3、dp_qrfenjie函数function A=dp_qrfenjie(a)n=length(a);A=a;A=a-a(n,n)*eye(n);theta(n-1)=0;c(n-1)=0;s(n-1)=0;Q=1;for j=1:n-1theta(j)=atan(A(j+1,j)/A(j,j));c(j)=cos(theta(j));s(j)=sin(theta(j));P=eye(n);P(j,j)=c(j);P(j+1,j)=-s(j);P(j,j+1)=s(j);P(j+1,j+1)=c(j);invP=eye(n);invP(j,j)=c(j);invP(j+1,j)=s(j);invP(j,j+1)=-s(j);invP(j+1,j+1)=c(j);Q=Q*invP;R=P*A;A=R;endA=R*Q+a(n,n)*eye(n);tor=max(abs(diag(A)-diag(a)));if tor>5*10^-5a=A;elsebreak;endendi,Ak=A,lanmda=diag(A)'4、vect函数function z=vect(a0,p)n=length(a0);Z=zeros(n,n+2);for x=1:na=a0-p(x)*eye(n);r=zeros(n,n);r(1,1:n)=a(1,1:n);l(2:n,1)=a(2:n,1)/a(1,1);for i=2:nfor j=i:nM=0;for k=1:i-1M=l(i,k)*r(k,j)+M;endr(i,j)=a(i,j)-M;endif i<nfor j=i+1:nN=0;for k=1:i-1N=l(j,k)*r(k,i)+N;endl(j,i)=(a(j,i)-N)/r(i,i);endelsebreak;endendR=r;L=l;u=ones(n,1);mo=0;invr=inv(R);invl=inv(L);for i=1:100y=invr*u;p1=max(y);p2=max(-y);if p1>p2m=p1;elsem=-p2;endz=y/m;tor=abs(m-mo);if tor<5*10^-7break;elsemo=m;endu=invl*z;endZ(x,3:n+2)=z';Z(x,1)=i;Z(x,2)=p(x)+1/m;endZ五、运行结果：1、执行命令：>> clear,clc,a=[2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1];hohd(a);得到原始实数矩阵的拟上三角矩阵A。
A =2.0000 -1.0000 0.0000 0.0000-1.0000 1.0000 5.0000 -0.00000.0000 5.0000 -2.2000 -0.40000.0000 -0.0000 -0.4000 2.20002、执行命令：>> clear,clc,A=[2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2];qrfenjie(A);得到迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。
i =38Ak =-5.9068 0.0315 -0.0000 0.00000.0315 4.8972 -0.0000 -0.00000.0000 -0.0000 2.2138 0.00070.0000 0.0000 0.0007 1.7958lanmda =-5.9068 4.8972 2.2138 1.79583、执行命令：>> clear,clc,A=[2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2];dp_qrfenjie(A);得到用带位移的QR分解法所用的迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。
i =8Ak =-5.9068 -0.0056 -0.0000 -0.0000-0.0056 4.8973 0.0000 0.0000-0.0000 0.0000 1.7958 0.00000.0000 0.0000 0.0000 2.2138lanmda =-5.9068 4.8973 1.7958 2.21384、执行命令：>> clear,clc,a=[2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1],p=[-5.9068 4.8972 2.2138 1.7958];vect(a,p);得到Z矩阵，Z矩阵每一行的第一个数是迭代次数，第二个数是特征值，后面的数为特征向量。