数值分析计算实习二.docx

1、数值分析计算实习二题目：试求矩阵的全部特征值，并对其中的每个实特征值求相应的特征向量，已知说明：1.要求迭代的精度水平为。2.打印以下内容：(1)采用带双步位移的QR分解法，说明算法设计方案和思路；(2)全部源程序；(3)矩阵经过拟上三角化后所得的矩阵；(4)对矩阵实行QR方法迭代结束后所得的矩阵；(5)矩阵A的全部特征值,其中，。若是实数，则令=0；(6) 的相应于实特征值的特征向量；(7)讨论：发现的现象与遇到的问题。3.采用e型数输出实型数，并且至少显示12位有效数字。算法设计方案：a) 声明二维1010数组空间，对其进行赋值得到矩阵；b) 将矩阵拟上三角化得到矩阵，具体算法见课本P61

2、-62；c) 对矩阵进行双步位移QR分解迭代，具体算法大致同课本P63-65，但有改进。改进部分为：先判断迭代矩阵右下角是否为二阶分块矩阵，若是，再求解其一对特征值。d) 根据求得的实特征值计算对应的特征向量。因为特征值已知，只需要求解以下线性方程组的一非零解即可：为了方便起见，直接利用QR分解函数分解矩阵B，则原方程组等价于然后利用直接法解出方程组即可。其中矩阵R右下角元素必为0，则使对应行未知数为1，然后回代即可。当矩阵R对应特征值的重数为2时，在矩阵R中除了右下角元素为0外，必有另一个对角元素为0，则令最后行对应未知数为0，该行对应未知数为1.当重数大于2时依次类推，即可求出对应于特定重

3、数特征值的特征向量。结果：1、 拟上三角矩阵如下。2、 迭代最终矩阵如下，迭代次数为14次。3、矩阵所有特征值如下。从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。4、对应于实特征值的特征向量如下问题发现与讨论1、 双步位移QR迭代过程中迭代矩阵始终为拟上三角矩阵。2、 双步位移QR方法与基本QR方法比较基本QR方法迭代过程中迭代矩阵始终为拟上三角矩阵。下图为基本QR方法得到的分块上三角矩阵：下图为基本QR方法得到的特征值：在基本QR迭代中发现，当矩阵收敛到分块上三角矩阵后，再次迭代它仍旧为分块上三角矩阵，然而矩阵中的元素改变了，不过块的阶数和位置不会改变。基本QR迭代收

4、敛速度为545次。可以看出虽然两者收敛到的分块上三角矩阵不同，但得到的最终特征值是一样的。也从而可初步验证所求得得特征值的正确性。收敛速度上基本QR方法明显慢于双步位移QR方法。源程序：#include stdafx.h#include stdio.h#include math.h#include stdlib.h#define N 10#define Ep 1e-12int _tmain(int argc, _TCHAR* argv) int ar0(double aNN, int s, int n); /判断矩阵多少位后是0 int sgn(double x); /0为1的符号函数 voi

5、d 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); /矩阵加法 void namul(double aN, double x, double bN, int n); /向量数乘 void naamul(double aNN, double x, double bNN, int n); /矩阵数乘 void amul(double

6、 aN, double bN, double cNN, int n); /向量乘法 double aneiji(double aN, double bN, int n); /向量内积 void aamul(double aN, double bNN, double cN, int n);/矩阵向量乘法 void aaamul(double aNN, double bNN, double cNN, int n); /矩阵矩阵乘法 void sol2(double a11, double a12, double a21, double a22, double(*p1)2, double(*p2)2

7、); /计算二阶矩阵的特征值 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 (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+) /将原矩阵拟上三角化

8、if (ar0(A,r+2,r) = 0) continue; else double 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 urN = 0; urr + 1 = Ar + 1r - cr; for (int i = r + 2; i N; i+) uri = Air; double prN; atra(A, N); aamul(ur

9、, 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); double BNN; amul(wr, ur, B, N); aaadd(A, B, A, N, -1); amul(ur, pr, B, N); aaadd(A, B, A, N,

10、 -1); /求一步操作后的矩阵Ar for (int i = 0; i N; i+) for (int j = 0; j N; j+) printf(% .11lf , Aij); if (j = N - 1) printf(n); printf(n); /打印拟上三角矩阵Ar double AqrNN = 0; aaadd(Aqr, A, Aqr, N, 1); double qqrNN = 0; for (int i = 0; i 10000; i+) double Aqr1NN = 0 ; aaadd(Aqr, Aqr1, Aqr1, N, 1); qrs(Aqr, qqr); dou

11、ble 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(Amm - 1) Ep) | (m=0) nmdm0 = Amm; nmdm1 = 0; m = m - 1; c

12、ontinue; else if (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; else double 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,

13、m + 1); aaadd(M, D, M, m + 1, -1); for (int i = 0; i m + 1; i+) Mii = Mii + t; for (int r = 0; r m; r+) /M的QR分解和Ar的计算 if (ar0(A, r + 1, r) = 0) continue; else double d = 0; for (int i = r; i m+1; i+) d = d + pow(Mir, 2); double dr = sqrt(d); double cr = -sgn(Mrr)*dr; double hr = pow(cr, 2) - cr*Mrr;

14、 double urN = 0 ; urr = Mrr - cr; for (int i = r + 1; i 10000) /算法失败判定 printf(算法失败n); break; for (int i = 0; i N; i+) for (int j = 0; j N; j+) printf(% .11lf , Aij); if (j = N - 1) printf(n); for (int i = 0; i N; i+) printf(%1.11e , nmdi0); printf(%1.11en, nmdi1); for (int i = 0; i Ep) continue; for

15、 (int j = i + 1; j Ep) continue; else if (fabs(nmdj0 - nmdi0) Ep) nmdj1 = 1; for (int k = 0; k Ep) continue; for (int i = 0; i N; i+) /求矩阵A-I*nmd for (int j = 0; j N; j+) if (i = j) Aij = 1.52*cos(i+1 + 1.2*(j+1) - nmdk0; else Aij = sin(0.5*(i+1) + 0.2*(j+1); double qNN; double yNN; int p1 = 0; qrs(

16、A,q); /对A-I*nmd进行QR分解得到新矩阵R Gauss(A,y,p); /齐次不满秩的高斯消去法,求非0根 printf(对应于%d实特征值的特征向量为n, k+1); for (int j = 0; j p0; j+) for (int i = 0; i N; i+) printf(%1.12en, yij); return 0; int ar0(double aNN, int s,int n) for (int i=s; i Ep) return 1; break; return 0;int sgn(double x) if (x 0) return 1; else retur

17、n -1;void atra(double aNN,int n) double x; for (int i = 0; i n; i+) for (int j = i+1; j n; j+) x = aij; aij = aji; aji = x; void aadd(double aN, double bN,double cN,int n,int co) if (co=1) for (int i = 0; i n; i+) ci = ai + bi; else if (co = -1) for (int i = 0; i n; i+) ci = ai - bi; else printf(向量加

18、减选择错误n); void aaadd(double aNN, double bNN,double cNN,int n,int co) if (co = 1) for (int i = 0; i n; i+) for (int j = 0; j n; j+) cij = aij + bij; else if (co = -1) for (int i = 0; i n; i+) for (int j = 0; j n; j+) cij = aij - bij; else printf(矩阵加减选择错误n); void namul(double aN, double x,double bN,int

19、 n) for (int i = 0; i n; i+) bi = x*ai;void naamul(double aNN, double x,double bNN,int n) for (int i = 0; i n; i+) for (int j = 0; j n; j+) bij = x*aij;void amul(double aN, double bN, double cNN,int n) for (int i = 0; i n; i+) for (int j = 0; j n; j+) cij = ai * bj;double aneiji(double aN, double bN

20、,int n) double t=0; for (int i = 0; i n; i+) t = t + ai * bi; return t;void aamul(double aN, double bNN, double cN,int n) for (int i = 0; i n; i+) ci = 0; for (int j = 0; j n; j+) ci =ci+ bij * aj; void aaamul(double aNN, double bNN, double cNN, int n) for (int i = 0; i n; i+) for (int j = 0; j n; j

21、+) cij = 0; for (int s = 0; s n; s+) cij = cij + ais * bsj; void sol2(double a11, double a12, double a21, double a22, double(*p1)2, double(*p2)2) double deta; deta = a11*a11 + a22*a22 - 2 * a11*a22 + 4 * a12*a21; if (deta 0) (*p1)0 = (a11 + a22) / 2; (*p2)0 = (a11 + a22) / 2; (*p1)1 = sqrt(0 - deta)

22、 / 2; (*p2)1 = 0 - sqrt(0 - deta) / 2; else (*p1)0 = (a11 + a22) / 2 + sqrt( deta) / 2; (*p1)1 = 0; (*p2)0 = (a11 + a22) / 2 - sqrt( deta) / 2; (*p2)1 = 0; void qrs(double aNN,double qNN) for (int i = 0; i N; i+) for (int j = 0; j N; j+) if (i = j) qij = 1; else qij = 0; for (int r = 0; r N - 1; r+) if (ar0(a, r + 1, r) = 0) continue; else double d=0, dr; double cr; double hr; for (int i = r