一、 算法设计方案：

按题目要求，本程序运用带双步位移的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)，再求解方程组。

方程组的求解采用列主元的高斯消去法，由函数intGauss_column(double a[][N],double b[],double X[])实现。由于此给定矩阵的特殊性，其没有重根，所有对应于每一特征值只有一个特征向量，因此可以用高斯消去法求解此奇异的线性方程组。首先用高斯消去将矩阵(A-λI)化为上三角阵，其最后一行必定全为零。因此在反代的过程中，只要令x[]的最后一个元素为“1”，即可得到方程组的一个基础解系，从而也就是一个特征向量。

6）输出有关结果

根据题目要求，需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。由于输出结果要求都要保留12位有效数字，所以将结果输出到文件result.txt中更有利于数据的打印。程序中构造了两个输出函数专门来解决不同输出结果的问题，void print_lambda(complex lambda[]);void print_matrix(double mat[][N])。

二、 程序源代码：

#include \#include \#include \

#define L 100 //定义最大迭代次数 #define N 10 //定义矩阵大小 #define err 1e-12 //定义误差限 //定义一个结构体来表示复数 typedefstruct complex{ double rpart; double ipart; };

FILE * pReadFile;

//主函数

int _tmain(intargc, _TCHAR* argv[]) { inti,j,t; double y[N],X[N],a[N][N],A[N][N],B[N][N],Q[N][N],R[N][N],RQ[N][N]; struct complex lambda[N]; //声明要调用的函数 void QuasiTriangularization(double a[][N]); void QR_decomposition(double A[][N],double Q[][N],double R[][N]); void QR_decompositionMk(double mk[][N],int m, double a[][N]); void print_lambda(complex lambda[]); void print_matrix(double mat[][N]); void multi_matrix(double a[][N],double b[][N],double c[][N]); intTwoStepDisplacement_QR(double a[][N],complex lambda[]); intGauss_column(double a[][N],double b[],double X[]); //矩阵初始化 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)); A[i][j] = a[i][j]; } else{ a[i][j] = sin(0.5 * (i+1) + 0.2 * (j+1)); A[i][j] = a[i][j]; } } } for (i = 0; i < N; i++){ y[i] = 0;

}

//对矩阵a[][]拟上三角化 QuasiTriangularization(a);

//打开文件result.txt

pReadFile = fopen( \//打印结果到文件result.txt中

fprintf(pReadFile,\拟上三角化之后的矩阵A[%d][%d]:\\n\//printf(\拟上三角化之后的矩阵A[%d][%d]:\\n\print_matrix(a);

//对拟上三角化后的矩阵a[][]，进行QR分解 QR_decomposition(a,Q,R);

fprintf(pReadFile,\//printf(\print_matrix(Q);

fprintf(pReadFile,\//printf(\print_matrix(R);

multi_matrix(R,Q,RQ);

fprintf(pReadFile,\//printf(\print_matrix(RQ);

//双步位移QR分解求全部特征值 TwoStepDisplacement_QR(a,lambda);

//利用校正的列主元素高斯消元法求每个实特征值对应的特征向量 for (t = 0; t < N; t++){ if (lambda[t].ipart == 0){ for (i = 0; i < N; i++){ for (j = 0; j < N; j++){ if (i == j) B[i][j] = A[i][j] - lambda[1].rpart; else B[i][j] = A[i][j]; } } Gauss_column(B,y,X);

fprintf(pReadFile,\矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\\n\ //printf(\矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\\n\ for (i = 0; i < N; i++){ fprintf(pReadFile,\ //printf(\ } } }

fclose( pReadFile ); scanf(\return 0;

}//主函数 /************************************* 矩阵的拟上三角化 输入n阶方阵a[][]，将a[][]拟上三角化 无返回值 **************************************/ void QuasiTriangularization(double a[][N]){ intr,i,j; double tr,hr,cr,dr,sum; double ur[N],pr[N],qr[N],wr[N]; for (r = 0; r < N-2; r++){

//判断a[i][r](i=r+2,r+3,...,n-2)是否全为零 sum = 0;//变量sum使用前清零 for (i = r+2; i < N; i++){ sum = sum || a[i][r]; }

//如果不是全部都是零，则计算 if (sum){

//计算dr sum = 0;

for (i = r+1; i < N; i++){ sum += a[i][r] * a[i][r]; }

dr = sqrt(sum);

//计算cr

if (a[r+1][r] > 0) cr = -dr; else cr = dr; //计算hr

hr = cr * cr - cr * a[r+1][r]; //取向量ur[]

for (i = 0; i < N; i++){ if (i < r+1) ur[i] = 0; else if (i == r+1) ur[i] = a[i][r] - cr; else ur[i] = a[i][r]; }

//计算向量qr[] for (i = 0; i < N; i++){ sum = 0; for (j = r+1; j < N; j++) sum += a[i][j] * ur[j]; qr[i] = sum / hr;

}

//计算向量pr[] for (i = 0; i < N; i++){ sum = 0; for (j = r+1; j < N; j++) sum += a[j][i] * ur[j]; pr[i] = sum / hr; }

//计算tr sum = 0;

for (i = r+1; i < N; i++) sum += pr[i] * ur[i]; tr = sum / hr;

//计算wr[]

for (i = 0; i < N; i++){ if (i < r+1) wr[i] = qr[i]; else wr[i] = qr[i] - tr * ur[i]; }

//计算新产生的矩阵a[][] for (i = 0; i < N; i++) for (j = 0; j < N; j++) a[i][j] = a[i][j] - wr[i] * ur[j] - ur[i] * pr[j]; } } } /************************************* 矩阵的QR分解 将A[][]分解为Q[][]*R[][] 无返回值 **************************************/

void QR_decomposition(double A[][N],double Q[][N],double R[][N]){ intr,i,j; double tr,hr,cr,dr,sum; double ur[N],pr[N],wr[N];

//取矩阵R1[][]为A[][] for (i = 0; i < N; i++) for (j = 0; j < N; j++) R[i][j] = A[i][j]; //取矩阵Q1[][]为单位矩阵 for (i = 0; i < N; i++) for (j = 0; j < N; j++){ if (i == j) Q[i][j] = 1; else Q[i][j] = 0; }