《数值分析计算实习二》由会员分享，可在线阅读，更多相关《数值分析计算实习二（14页珍藏版）》请在金锄头文库上搜索。
1、题目：试求矩阵A=aij1010的全部特征值，并对其中的每个实特征值求相应的特征向量，已知aij=sin0.5i+0.2j ij1.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的相应于实特征值的特征
2、向量；(7)讨论：发现的现象与遇到的问题。3.采用e型数输出实型数，并且至少显示12位有效数字。算法设计方案：a) 声明二维1010数组空间，对其进行赋值得到矩阵A；b) 将矩阵A拟上三角化得到矩阵A(n-1)，具体算法见课本P61-62；c) 对矩阵A(n-1)进行双步位移QR分解迭代，具体算法大致同课本P63-65，但有改进。改进部分为：先判断迭代矩阵右下角是否为二阶分块矩阵，若是，再求解其一对特征值。d) 根据求得的实特征值计算对应的特征向量。因为特征值已知，只需要求解以下线性方程组的一非零解即可：Bx=0B=A-I为了方便起见，直接利用QR分解函数分解矩阵B，则原方程组等价于Rx=0然
3、后利用直接法解出方程组即可。其中矩阵R右下角元素必为0，则使对应行未知数为1，然后回代即可。当矩阵R对应特征值的重数为2时，在矩阵R中除了右下角元素为0外，必有另一个对角元素为0，则令最后行对应未知数为0，该行对应未知数为1.当重数大于2时依次类推，即可求出对应于特定重数特征值的特征向量。结果：1、 拟上三角矩阵A(n-1)如下。2、 迭代最终矩阵如下，迭代次数为14次。3、矩阵所有特征值如下。从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。4、对应于实特征值的特征向量如下 问题发现与讨论1、 双步位移QR迭代过程中迭代矩阵Ak+1始终为拟上三角矩阵。2、 双步位移QR方
4、法与基本QR方法比较基本QR方法迭代过程中迭代矩阵Ak+1始终为拟上三角矩阵。下图为基本QR方法得到的分块上三角矩阵：下图为基本QR方法得到的特征值：在基本QR迭代中发现，当矩阵收敛到分块上三角矩阵后，再次迭代它仍旧为分块上三角矩阵，然而矩阵中的元素改变了，不过块的阶数和位置不会改变。基本QR迭代收敛速度为545次。可以看出虽然两者收敛到的分块上三角矩阵不同，但得到的最终特征值是一样的。也从而可初步验证所求得得特征值的正确性。收敛速度上基本QR方法明显慢于双步位移QR方法。源程序：#include stdafx.h#include stdio.h#include math.h#include
5、stdlib.h#define N 10#define Ep 1e-12int _tmain(int argc, _TCHAR* argv) int ar0(double aNN, int s, int n);/判断矩阵多少位后是0int sgn(double x);/0为1的符号函数void atra(double aNN, int n);/矩阵转置void aadd(double aN, double bN, double cN, int n, int co);/向量加法void aaadd(double aNN, double bNN, double cNN, int n, int co
6、);/矩阵加法void namul(double aN, double x, double bN, int n);/向量数乘void naamul(double aNN, double x, double bNN, int n);/矩阵数乘void amul(double aN, double bN, double cNN, int n);/向量乘法doubleaneiji(double aN, double bN, int n);/向量内积void aamul(double aN, double bNN, double cN, int n);/矩阵向量乘法void aaamul(double
7、 aNN, double bNN, double cNN, int n);/矩阵矩阵乘法void sol2(double a11, double a12, double a21, double a22, double(*p1)2, double(*p2)2);/计算二阶矩阵的特征值void qrs(double aNN, double qNN);void Gauss(double aNN,double yNN,int *p);double ANN;double nmdN2 = 0;for (int i = 0; i N; i+)/对原矩阵赋值for (int j = 0; j N; j+)if
8、 (i = j)Aij = 1.52*cos(i+1 + 1.2*(j+1);else Aij = sin(0.5*(i+1) + 0.2*(j+1);for (int r = 0; r N - 2; r+)/将原矩阵拟上三角化if (ar0(A,r+2,r) = 0) continue;elsedouble d = 0;for (int i = r + 1; i N; i+)d = d + pow(Air,2);double dr = sqrt(d);double cr = -sgn(Ar + 1r)*dr;double hr = pow(cr, 2) - cr*Ar + 1r;double
9、 urN = 0;urr + 1 = Ar + 1r - cr;for (int i = r + 2; i N; i+)uri = Air;double prN;atra(A, N);aamul(ur, A, pr, N);namul(pr, 1 / hr, pr, N);double qrN;atra(A, N);aamul(ur, A, qr, N);namul(qr, 1 / hr, qr, N);double tr;tr = aneiji(pr, ur, N) / hr;double wrN;namul(ur, -tr, wr, N);aadd(qr, wr, wr, N,1);dou
10、ble BNN;amul(wr, ur, B, N);aaadd(A, B, A, N, -1);amul(ur, pr, B, N);aaadd(A, B, A, N, -1);/求一步操作后的矩阵Arfor (int i = 0; i N; i+)for (int j = 0; j N; j+)printf(% .11lf , Aij);if (j = N - 1)printf(n);printf(n);/打印拟上三角矩阵Ardouble AqrNN = 0;aaadd(Aqr, A, Aqr, N, 1);double qqrNN = 0;for (int i = 0; i 10000;
11、 i+)double Aqr1NN = 0 ;aaadd(Aqr, Aqr1, Aqr1, N, 1);qrs(Aqr, qqr);double BNN = 0;aaamul(Aqr, qqr, B, N);for (int s = 0; s N; s+)for (int k = 0; k N; k+)Aqrsk = Bsk;int s=0;for (int t = 0; t N; t+)if (fabs(Aqrt+1t) Ep)s = s + 1;if (s = N-2)break;for (int i = 0; i N; i+)for (int j = 0; j -1)if (fabs(A
12、mm - 1) Ep) | (m=0)nmdm0 = Amm;nmdm1 = 0;m = m - 1;continue;elseif (fabs(Am - 1m - 2) Ep) | (m=1)sol2(Am - 1m - 1, Am - 1m, Amm - 1, Amm,(nmd+m),(nmd+(m - 1);/求解二阶矩阵的特征值m = m - 2;continue;elsedouble s = Am - 1m - 1 + Amm;double t = Am - 1m - 1 * Amm - Amm - 1 * Am - 1m;double MNN = 0;double DNN = 0;aaamul(A, A, M, m + 1);naamul(A, s, D, m + 1);aaadd(M, D, M, m + 1, -1);for (int i = 0; i m + 1; i+)Mi
