数值分析计算实习二教程.docxVIP

题目：? 试求矩阵A=[aij]10×10的全部特征值，并对其中的每个实特征值求相应的特征向量，已知 aij=sin0.5i+0.2j i≠j1.52cosi+1.2j i=j (i,j=1,2,…,10) 说明：? 1.要求迭代的精度水平为ε=10-12。 2.打印以下内容：? (1)采用带双步位移的QR分解法，说明算法设计方案和思路； (2)全部源程序； (3)矩阵A经过拟上三角化后所得的矩阵A(n-1)； (4)对矩阵A(n-1)实行QR方法迭代结束后所得的矩阵； (5)矩阵A的全部特征值λi=(Ri,Ii)(i=1,2,?,10),其中Ri=Re(λi)， Ii=Im(λi) 。若λi是实数，则令Ii=0； (6)A的相应于实特征值的特征向量； (7)讨论：发现的现象与遇到的问题。 3.采用e型数输出实型数，并且至少显示12位有效数字。 算法设计方案： 声明二维10×10数组空间，对其进行赋值得到矩阵A； 将矩阵A拟上三角化得到矩阵A(n-1)，具体算法见课本P61-62； 对矩阵A(n-1)进行双步位移QR分解迭代，具体算法大致同课本P63-65，但有改进。改进部分为：先判断迭代矩阵右下角是否为二阶分块矩阵，若是，再求解其一对特征值。 根据求得的实特征值计算对应的特征向量。因为特征值已知，只需要求解以下线性方程组的一非零解即可： Bx=0 B=A-λI 为了方便起见，直接利用QR分解函数分解矩阵B，则原方程组等价于 Rx=0 然后利用直接法解出方程组即可。其中矩阵R右下角元素必为0，则使对应行未知数为1，然后回代即可。当矩阵R对应特征值的重数为2时，在矩阵R中除了右下角元素为0外，必有另一个对角元素为0，则令最后行对应未知数为0，该行对应未知数为1.当重数大于2时依次类推，即可求出对应于特定重数特征值的特征向量。 结果： 拟上三角矩阵A(n-1)如下。 迭代最终矩阵如下，迭代次数为14次。 3、矩阵所有特征值如下。 从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。 4、对应于实特征值的特征向量如下 问题发现与讨论 双步位移QR迭代过程中迭代矩阵Ak+1始终为拟上三角矩阵。 双步位移QR方法与基本QR方法比较 基本QR方法迭代过程中迭代矩阵Ak+1始终为拟上三角矩阵。 下图为基本QR方法得到的分块上三角矩阵： 下图为基本QR方法得到的特征值： 在基本QR迭代中发现，当矩阵收敛到分块上三角矩阵后，再次迭代它仍旧为分块上三角矩阵，然而矩阵中的元素改变了，不过块的阶数和位置不会改变。 基本QR迭代收敛速度为545次。 可以看出虽然两者收敛到的分块上三角矩阵不同，但得到的最终特征值是一样的。也从而可初步验证所求得得特征值的正确性。收敛速度上基本QR方法明显慢于双步位移QR方法。 源程序： #include stdafx.h #include stdio.h #include math.h #include stdlib.h #define N 10 #define Ep 1e-12 int _tmain(int argc, _TCHAR* argv[]) { int ar0(double a[N][N], int s, int n); //判断矩阵多少位后是0 int sgn(double x); //0为1的符号函数 void atra(double a[N][N], int n); //矩阵转置 void aadd(double a[N], double b[N], double c[N], int n, int co); //向量加法 void aaadd(double a[N][N], double b[N][N], double c[N][N], int n, int co); //矩阵加法 void namul(double a[N], double x, double b[N], int n); //向量数乘 void naamul(double a[N][N], double x, double b[N][N], int n); //矩阵数乘 void amul(double a[N], double b[N], double c[N][N], int n); //向量乘法 double aneiji(double a[N], double b[N], int n); //向量内积 void aamul(double a[N], double b[N][N], double c[N], int n);//矩阵向量乘法 void aaamul(double a[N][N],