北航数值分析大作业 第二题 QR分解.doc

北航2009级研究生 《数值分析B》计算实习题目 （第 二 题：QR分解） 设计文档与源程序 2009年11月16日 打 印 内 容 1 算法的设计方案 （1）运行平台 （2）算法设计 2 全部源代码 3 输出结果： 3.1 矩阵A经过拟上三角化后所得的矩阵A(n-1)； 3.2 对矩阵A(n-1)进行QR分解方法后所得的矩阵； 3.3 矩阵A的全部特征值λｉ； 3.4 A的相对应于实特征值的特征向量。 1 算法的设计方案 1.1 运行与开发平台 操作系统：Windows 7； 开发平台：VC++ 6.0； 工程类型：Win32 Console Application； 工程名：QR_EigenValue； 1.2 算法设计 设计思想：本题目要求采用带双步位移的QR分解法求矩阵A={aij}10×10的全部特征值。首先由已知公式计算出A各元素值；然后对矩阵A进行拟上三角化得到矩阵A(n-1)；然后对得出的A(n-1)进行QR分解，得到一个分块上三角矩阵以及矩阵的各个特征值λi；最后，采用高斯消去法，可以求得对应于各实特征值的特征向量。 具体算法如下：（精度eps=le-12，最大迭代次数L=50，n=10） (一)、计算矩阵A 对i,j=0到10，执行： if i=j,a[i][j]=1.5*cos*(i+1.2j);else a[i][j]=sin*(0.5*i+0.2j); (二)、矩阵拟上三角化函数 矩阵拟上三角化函数为：void Hessenberg_A(double a[][n],double a_n1[][n])。 定义矩阵ar = a，对于r=1，2，……，n-2执行以下算法： 第一步：判断a[i-1][r-1]是否全为零（i=r+2,……,n）,如果是，则令a[i-1][r-1] = a[i-1][r-1]，并执行循环；否则，跳出第一步的循环。 第二步：对i=r+1，……，n计算 d_r = sqrt(ar[i-1][r-1] * ar[i-1][r-1]) cr=--sgn(d_r) h_r = c_r*c_r - c_r*ar[r][r-1] 第三步：对 i=r+1，……，n，执行 for(i=r+1;i=n;i++) ur[i-1] = ar[i-1][r-1]; ur[r] = ur[r] - c_r; 第四步：计算 pr=A’*ur/hr qr=A*ur/hr tr=pr’*ur/hr wr=qr-tr*ur A=A-wr*ur-ur*pr’ (三)、对A(n-1)做带双步位移的QR分解 定义的带双步位移的QR分解函数为int Db_QR_Method(double Ak[n][n], double lamda[n][2], FILE *pFile) 记Ak[i][j] = A(n-1)=a_n1[i][j]，令k=0,m=n; (2) if | A[m][m-1]=eps|,则得到A的一个特征值lamda[m][0]=A[m][m]，置m=m-1,转(3);否则转(4); (3) if m=1,则得到A的一个特征值lamda[m][0]=A[1][1]，转(10);if m=0，则直接转(10)；if m1,则转(2); (4) 调用了解二次方程函数，用来求二阶子阵的两个特征值，即计算二次方程的两个根， x*x- s*x+ t=0 （*） 其中s = A[m-1][m-1]+A[m][m]; t]= A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]； (5) if m=2,则得到A的两个特征值，转(10)；否则转(7)； (6) if |A[m-1][m-2]) |=eps,则得到A的两个特征根，置m=m-2，转(3);否则转(7)； (7) if k=L,跳出程序，未能得到结果；否则转(8)； (8) 对i,j=1,…,m，计算 s= A[m-1][m-1]+A[m][m]; t= A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]; Mk[i][j] = Mk[i][j]-s*A[i][j]+t*I[i][j] Mk=Qk*Rk（对Mk作QR分解，可调用QR分解函数） Ak+1=Qk’*Ak*Qk (9) 置k=k+1，转(2)； (10) 得到所有的特征值，停止计算。 Mk的QR分解算法： 对r=1到m-1执行 (1) if |B[i][r]|=eps，转(2); (2) 计算 dr=sqrt(A[i][r]* A[i][r]) cr=--sgn(A[r+1][r]) hr = cr*(cr-A[r+1]