《北航-数值分析第二次大作业(带双步位移的QR方法)》由会员分享，可在线阅读，更多相关《北航-数值分析第二次大作业(带双步位移的QR方法)（18页珍藏版）》请在金锄头文库上搜索。
1、一、 算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。总体思路：1） 初始化矩阵首先需要将需要求解的矩阵输入程序。为了防止矩阵在后面的计算中被破坏保存A。2） 对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。由于矩阵的QR分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double aN。3） 对拟上三角化后的矩阵进行QR分解
2、对拟上三角化的矩阵进行QR分解会大大减小计算量。用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。4） 对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。这里用函数intTwoStepDisplacement_QR()来实现。这是用来存储计算得到的特征值
3、的二维数组。考虑到特征值可能为复数，因此将所有特征值均当成复数处理。此函数中，QR分解部分用子函数void QR_decompositionMk()实现。这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。5） 计算特征向量得到特征值后，计算实特征值相应的特征向量。这里判断所得特征值的虚数部分是否为零。求实特征值的特征向量采用求解相应的方程组(A-I)x=0)的方法。因此先初始化矩阵Array，计算(A-I)，再求解方程组。方程组的求解采用列主元的高斯消去法，由函数intGauss_column(double aN,double b,double X)实现。由
4、于此给定矩阵的特殊性，其没有重根，所有对应于每一特征值只有一个特征向量，因此可以用高斯消去法求解此奇异的线性方程组。首先用高斯消去将矩阵(A-I)化为上三角阵，其最后一行必定全为零。因此在反代的过程中，只要令x的最后一个元素为“1”，即可得到方程组的一个基础解系，从而也就是一个特征向量。6） 输出有关结果根据题目要求，需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。由于输出结果要求都要保留12位有效数字，所以将结果输出到文件result.txt中更有利于数据的打印。程序中构造了两个输出函数专门来解决不同输出结果的问题，void print_lambda
5、(complex lambda);void print_matrix(double matN)。二、 程序源代码：#include stdafx.h#include stdlib.h#include math.h#define L 100 /定义最大迭代次数#define N 10 /定义矩阵大小#define err 1e-12 /定义误差限/定义一个结构体来表示复数typedefstruct complexdouble rpart;double ipart;FILE *pReadFile;/主函数int _tmain(intargc, _TCHAR* argv)inti,j,t;doubl
6、e yN,XN,aNN,ANN,BNN,QNN,RNN,RQNN;struct complex lambdaN;/声明要调用的函数void QuasiTriangularization(double aN);void QR_decomposition(double AN,double QN,double RN);void QR_decompositionMk(double mkN,int m, double aN);void print_lambda(complex lambda);void print_matrix(double matN);void multi_matrix(double
7、aN,double bN,double cN);intTwoStepDisplacement_QR(double aN,complex lambda);intGauss_column(double aN,double b,double X);/矩阵初始化for (i = 0; i N; i+)for (j = 0; j N; j+)if (i = j)aij = 1.5 * cos(i+1) + 1.2 * (j+1);Aij = aij;else aij = sin(0.5 * (i+1) + 0.2 * (j+1);Aij = aij;for (i = 0; i N; i+)yi = 0;
8、/对矩阵a拟上三角化QuasiTriangularization(a);/打开文件result.txtpReadFile = fopen( result.txt, w );/打印结果到文件result.txt中fprintf(pReadFile,拟上三角化之后的矩阵A%d%d:n,N,N);/printf(拟上三角化之后的矩阵A%d%d:n,N,N);print_matrix(a);/对拟上三角化后的矩阵a，进行QR分解QR_decomposition(a,Q,R);fprintf(pReadFile,Q%d%d:n,N,N);/printf(Q%d%d:n,N,N);print_matrix
9、(Q);fprintf(pReadFile,R%d%d:n,N,N);/printf(R%d%d:n,N,N);print_matrix(R);multi_matrix(R,Q,RQ);fprintf(pReadFile,RQ%d%d:n,N,N);/printf(RQ%d%d:n,N,N);print_matrix(RQ);/双步位移QR分解求全部特征值TwoStepDisplacement_QR(a,lambda);/利用校正的列主元素高斯消元法求每个实特征值对应的特征向量for (t = 0; t N; t+)if (lambdat.ipart = 0)for (i = 0; i N;
10、i+)for (j = 0; j N; j+)if (i = j)Bij = Aij - lambda1.rpart;elseBij = Aij;Gauss_column(B,y,X);fprintf(pReadFile,n矩阵与特征值=%.12e对应的特征向量为X%d：n,lambdat.rpart,N);/printf(n矩阵与特征值=%.12e对应的特征向量为X%d：n,lambdat.rpart,N);for (i = 0; i N; i+)fprintf(pReadFile,%.12e ,i,Xi);/printf(X%d: %.12e ,i,Xi);fclose( pReadFil
11、e );scanf(%d,&i);return 0;/主函数/*矩阵的拟上三角化输入n阶方阵a，将a拟上三角化无返回值*/void QuasiTriangularization(double aN)intr,i,j;double tr,hr,cr,dr,sum;double urN,prN,qrN,wrN;for (r = 0; r N-2; r+)/判断air(i=r+2,r+3,.,n-2)是否全为零sum = 0;/变量sum使用前清零for (i = r+2; i N; i+)sum = sum | air;/如果不是全部都是零，则计算if (sum)/计算drsum = 0;for
12、(i = r+1; i 0)cr = -dr;else cr = dr;/计算hrhr = cr * cr - cr * ar+1r;/取向量urfor (i = 0; i N; i+)if (i r+1)uri = 0;else if (i = r+1)uri = air - cr;elseuri = air;/计算向量qrfor (i = 0; i N; i+)sum = 0;for (j = r+1; j N; j+)sum += aij * urj;qri = sum / hr;/计算向量prfor (i = 0; i N; i+)sum = 0;for (j = r+1; j N; j+)sum += aji * urj;pri = sum / hr;/计算trsum = 0;for (i = r+1; i N; i+)sum += pri *
