《数值分析》计算实习报告第二题院系：机械工程及自动化学院学号：姓名：2017年11月一、题目要求试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知a ij ={sin (0.5i +0.2j ) i ≠j1.52cos (i +1.2j ) i =j(i,j =1,2, (10)说明：1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10−12。

2.打印以下内容： （1）全部源程序；（2）矩阵A 经过拟上三角化后所得的矩阵A (n−1)； （3）对矩阵A (n−1)实行QR 方法迭代结束后所得的矩阵；（4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,⋯,10)，其中R i =Re(λi )，I i = Im(λi ) 。

若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。

3.采用e 型数输出实型数，并且至少显示12位有效数字。

二、算法设计思路和方案1. 将矩阵A 拟上三角化得到矩阵A (n−1)为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n−1)，然后用QR 方法计算A (n−1)的全部特征值，而A (n−1)的特征值就是A 的特征值。

具体算法如下：记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2,,;,1,,)ij a i n j r r n ==+。

对于1,2,,2r n =-执行（1）若()(2,3,,)r ira i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。

（2）计算r d =()()1,sgn r r r r r c a d +=- (若()1,0r r r a +=，则取r r c d =)2()1,r r r r r r h c c a +=-（3）令()()()()1,2,0,,0,,,,T r r r n r r rr r rnru ac aaR ++=-∈。

（4）计算()()(1)()///r T r r r r r r rTr r r rr r r rr r T T r r r rp A u h q A u h t p u h q t u A A u u p ωω+====-=--（5）继续。

当此算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n−1)。

2.用带双步位移的QR 方法求矩阵A 全部特征值为了加速收敛，采用带双步位移的QR 方法求解A 的全部特征值，具体算法如下：（1）由步骤1得到拟上三角矩阵(1)n A -；给定精度水平0ε> 和最大迭代次数L 。

（2）记(1)(1)1=n ij n nA A a -⨯⎡⎤=⎣⎦，令k = 1,m = n 。

（3）如果(),1k m m a ε-≤，则得到A 的一个特征值()k mm a ，置m:=m-1（降阶），转（4）；否则转（5）。

（4）如果2m ≤，则得到矩阵A 的一个特征值()11k a ，或两个特征值12,s s ，转（9）；否则转（3）。

（5）如果()1,2k m m a ε--≤，则得到A 的两个特征值12,s s ，置m : = m – 2（降阶），转（4）；否则转（6）。

（6）如果k = L ，则终止计算，未得到A 的全部特征值；否则转（7）。

（7）记()(1,)k k ij m m A a i j m ⨯⎡⎤=≤≤⎣⎦，计算()()1,1()()()()1,1,11,21)() (k k m m mmk k k k m m mm m m m mk k k k k k k T k k k ks a a t a a a a M A sA tI I m M Q R M QR A Q A Q ------+=+=-=-+==是阶单位矩阵对作分解 （8）置k : = k + 1，转（3）。

（9）A 的全部特征值已计算完毕，停止计算。

其中，k M 的QR 分解与1k A +的计算用下列算法实现：记(1)()11[],[],r k ij m m r ij m m k B M b B b C A ⨯⨯====。

对于1,2,,1r m =-执行（1）若()()1,2,,r ir b i r r m =++全为零，则令11,r r r r B B C C ++==，转（5）；否则转（2）。

（2）计算()()()()1,2()sgn 0,r r r r rr r r r r r r r r r rrd c b d a c d h c c b +==-===-若则取 （3）令()()()()1,0,,0,,,,Tr r r m r rr r r r mr u b c b b R +=-∈。

（4）计算11////T r r r r T r r r r T r r r r r r r rTr r r rr r r rT T r r r r r rv B u h B B u v p C u h q C u h t p u h q t u C C u u p ωω++==-====-=--（5）继续。

此算法执行完后，得到1k m A C +=。

3.用列主元素Gauss 消去法求特征向量记λ为矩阵A 的实特征值，x 为对应的特征向量，则A x =λx ，即(A - λ)x = 0的解即为矩阵A 的特征向量。

因此对于特征向量的求解可采用列主元素Gauss 消去法。

0A I λ-=，经检查由于矩阵A 无重特征值，所以每个特征值对应的特征向量的线性空间为1，矩阵A I λ-的秩为N-1，经过列主元素Gauss 消去法消元后，矩阵有且仅有最后一行全为0。

此时设定1n x =-，再依次求出()1,2,,1k x k n n =--，就可得到对应λ的其中一个特征向量12(,,,)T n x x x x =，对应λ的全部特征向量可表示为k x 。

最后为了得到标准形式，将向量x 正交化。

具体算法如下。

1. 消元过程 对于1,2,,1k n =-执行（1）选行号k i ，使()()max k k k i k ik k i na a ≤≤=。

（2）交换()k kj a 与()(,1,,)k k i j a j k k n =+所含的数值。

（3）对于1,2,,i k k n =++计算()()(1)()()/ (1,2,,)k k ik ik kkk k k ijij ik kim a a aam a j k k n +==-=++2. 回代过程()()()111,2,,1n n k k k kj j kk j k x x a x a k n n =+=-⎛⎫=-=-- ⎪⎝⎭∑最终得到的向量12(,,,1)T n x x x =-即为A 对应于实特征值λ的特征向量。

3. 将向量正交化。

三、源程序代码#include<windows.h>#include<iostream> #include<iomanip> #include<math.h> using namespace std;#define N 10 #define E 1e-12#define MAX 10000 //最大迭代次数//创建原始矩阵Avoid create_A(double a[N][N]) { int i, j; for (i = 0; i<N; i++) for (j = 0; j<N; j++) if (i != j) a[i][j] = sin(0.5*(i + 1) + 0.2*(j + 1)); else a[i][j] = 1.52*cos((i + 1) + 1.2*(j + 1)); }//符号函数 int sgn(double a) { int z; if (a>E) z = 1; else12(,,,1)T n x x x =-z = -1;return z;}//矩阵中小于或等于E的元素置0void zhi0(double A[N][N]){int i, j;for (i = 0; i<N; i++)for (j = 0; j<N; j++)if (fabs(A[i][j]) <= E)A[i][j] = 0;}//二次多项式求deltadouble delta(double b, double c){double m;m = b*b - 4 * c;return m;}//求二次方程两个根s1，s2void s1s2(double A[N][N], int m, double lamdaR[N], double lamdaI[N]) {double s, t, x;s = A[m - 1][m - 1] + A[m][m];t = A[m - 1][m - 1] * A[m][m] - A[m][m - 1] * A[m - 1][m];x = delta(s, t);if (x>E){x = sqrt(x);lamdaR[m] = s / 2 - x / 2;lamdaR[m - 1] = s / 2 + x / 2;lamdaI[m] = 0;lamdaI[m - 1] = 0;}else{x = sqrt(fabs(x));lamdaR[m] = s / 2;lamdaR[m - 1] = s / 2;lamdaI[m] = -x / 2;lamdaI[m - 1] = x / 2;}}//矩阵的拟上三角化void nishangsjh(double A[N][N]){int i, j, k, flag;double d, c, h, t;double u[N], p[N], q[N], w[N];for (i = 0; i<N - 2; i++){flag = 0;for (j = i + 2; j<N; j++)if (fabs(A[j][i])>E){flag = 1;break;}if (flag == 0)continue;for (j = i + 1, d = 0; j<N; j++)d = d + A[j][i] * A[j][i];d = sqrt(d);c = -1 * sgn(A[i + 1][i])*d;h = c*c - c*A[i + 1][i];for (j = i + 2; j<N; j++)u[j] = A[j][i];for (j = 0; j<i + 2; j++)u[j] = 0;u[i + 1] = A[i + 1][i] - c;for (j = 0; j<N; j++){for (k = i + 1, p[j] = 0; k<N; k++)p[j] = A[k][j] * u[k] + p[j];p[j] = p[j] / h;}for (j = 0; j<N; j++){for (k = i + 1, q[j] = 0; k<N; k++)q[j] = A[j][k] * u[k] + q[j];q[j] = q[j] / h;}for (j = 0, t = 0; j<N; j++)t = t + p[j] * u[j];t = t / h;w[j] = q[j] - t*u[j];for (j = 0; j<N; j++)for (k = 0; k<N; k++)A[j][k] = A[j][k] - w[j] * u[k] - u[j] * p[k];}zhi0(A);}//对矩阵A进行一次基本QR方法迭代void QRdiv(double A[N][N],double RQ[N][N]){int i, j, k, flag;double Q[N][N], R[N][N];double d, c, h;double u[N], w[N], p[N];for (i = 0; i<N; i++)for (j = 0; j<N; j++)if (i == j)Q[i][j] = 1;elseQ[i][j] = 0;for (i = 0; i<N; i++)for (j = 0; j<N; j++)R[i][j] = A[i][j];for (i = 0; i<N - 1; i++){flag = 0;for (j = i + 1; j<N; j++)if (fabs(R[j][i])>E){flag = 1;break;}if (flag == 0)continue;for (j = i, d = 0; j<N; j++)d = d + R[j][i] * R[j][i];d = sqrt(d);c = -1 * sgn(R[i][i])*d;h = c*c - c*R[i][i];for (j = i + 1; j<N; j++)u[j] = R[j][i];u[j] = 0;u[i] = R[i][i] - c;for (j = 0; j<N; j++)for (k = 0, w[j] = 0; k<N; k++)w[j] = Q[j][k] * u[k] + w[j];for (j = 0; j<N; j++)for (k = 0; k<N; k++)Q[j][k] = Q[j][k] - w[j] * u[k] / h;for (j = 0; j<N; j++){for (k = i, p[j] = 0; k<N; k++)p[j] = R[k][j] * u[k] + p[j];p[j] = p[j] / h;}for (j = 0; j<N; j++)for (k = 0; k<N; k++)R[j][k] = R[j][k] - u[j] * p[k];}zhi0(R);zhi0(Q);for (i = 0; i<N; i++)for (j = 0; j<N; j++)for (k = 0, RQ[i][j] = 0; k<N; k++)RQ[i][j] = RQ[i][j] + R[i][k] * Q[k][j]; }//矩阵M的QR分解void QRdiv_M(double M[N][N], double A[N][N], int &m) {int i, j, k, flag;double d, c, h, tr;double u[N], w[N], p[N], q[N], v[N];for (i = 0; i<m; i++){flag = 0;for (j = i + 1; j<N; j++)if (fabs(A[j][i])>E){flag = 1;break;}if (flag == 0)continue;for (j = i, d = 0; j <= m; j++)d = d + M[j][i] * M[j][i];d = sqrt(d);c = -1 * sgn(M[i][i])*d;h = c*c - c*M[i][i];for (j = i + 1; j <= m; j++)u[j] = M[j][i];for (j = 0; j<i; j++)u[j] = 0;u[i] = M[i][i] - c;for (j = 0; j <= m; j++){for (v[j] = 0, k = i; k <= m; k++)v[j] = v[j] + M[k][j] * u[k];v[j] = v[j] / h;}for (k = i; k <= m; k++)for (j = 0; j <= m; j++)M[k][j] = M[k][j] - u[k] * v[j];for (j = 0; j <= m; j++){for (p[j] = 0, k = i; k <= m; k++)p[j] = p[j] + A[k][j] * u[k];p[j] = p[j] / h;}for (j = 0; j <= m; j++){for (q[j] = 0, k = i; k <= m; k++)q[j] = q[j] + A[j][k] * u[k];q[j] = q[j] / h;}for (tr = 0, j = i; j <= m; j++)tr = tr + p[j] * u[j];tr = tr / h;for (j = 0; j <= m; j++)w[j] = q[j] - tr*u[j];for (j = 0; j <= m; j++)for (k = 0; k <= m; k++)A[j][k] = A[j][k] - (w[j] * u[k] + u[j] * p[k]);}}//带双步位移的QR方法求特征值int tezhengzhi(double A[N][N], double lamdaR[N], double lamdaI[N]) {int L = 0, m = N - 1;int i, j, k;double s, t;double M[N][N];int flag;loop3:if (fabs(A[m][m - 1]) <= E){lamdaR[m] = A[m][m]; lamdaI[m] = 0;m--;goto loop4;}elsegoto loop5;loop4:if (m<2){if (m == 0){lamdaR[0] = A[0][0];lamdaI[0] = 0;}if (m == 1)s1s2(A, m, lamdaR, lamdaI);flag = 1;goto loop9;}elsegoto loop3;loop5:if (fabs(A[m - 1][m - 2]) <= E){s1s2(A, m, lamdaR, lamdaI);m = m - 2;goto loop4;}elsegoto loop6;loop6:if (L == MAX){flag = 0;goto loop9;}elsegoto loop7;loop7:s = A[m - 1][m - 1] + A[m][m];t = A[m - 1][m - 1] * A[m][m] - A[m][m - 1] * A[m - 1][m];for (i = 0; i <= m; i++)for (j = 0; j <= m; j++)if (i == j){for (k = 0, M[i][j] = 0; k <= m; k++)M[i][j] = A[i][k] * A[k][j] + M[i][j];M[i][j] = M[i][j] - s*A[i][j] + t;}else{for (k = 0, M[i][j] = 0; k <= m; k++)M[i][j] = A[i][k] * A[k][j] + M[i][j];M[i][j] = M[i][j] - s*A[i][j];}QRdiv_M(M, A, m);goto loop8;loop8:L++;goto loop3;loop9:return L;}//列主元的高斯消元法求解特征向量void gauss(double x[N][N], double tzR[N], int &nR, double lamdaR[N], double lamdaI[N]) {double A[N][N];double sum, ar, ch, m[N];int i, j, k, p, ik;for (i = 0, nR = 0; i<N; i++)if (fabs(lamdaI[i]) <= E)tzR[nR++] = lamdaR[i];for (p = 0; p<nR; p++){create_A(A);for (i = 0; i<N; i++)for (j = 0; j<N; j++)if (i == j)A[i][j] = A[i][j] - tzR[p];for (k = 0; k<N - 1; k++){for (ik = k, i = k + 1; i<N; i++)if (A[ik][k]<A[i][k])ik = i;for (j = k; j<N; j++){ch = A[ik][j];A[ik][j] = A[k][j];A[k][j] = ch;}for (i = k + 1; i<N; i++){m[i] = A[i][k] / A[k][k];for (j = k + 1; j<N; j++)A[i][j] = A[i][j] - m[i] * A[k][j];}}x[p][N - 1] = -1;for (k = N - 2; k >= 0; k--){for (ar = 0, j = k + 1; j<N; j++)ar = ar + A[k][j] * x[p][j];x[p][k] = -ar / A[k][k];}for (j = 0, sum = 0; j<N; j++)sum = sum + x[p][j] * x[p][j];sum = sqrt(sum);for (j = 0; j<N; j++)x[p][j] = x[p][j] / sum;}}//打印N*N阶矩阵void print_A(double a[][N]){int i, j;printf("\n");for (i = 0; i<N; i++){for (j = 0; j<int(N / 2); j++)cout << setw(22) << a[i][j];}printf("\n");for (i = 0; i<N; i++){for (j = int(N / 2); j<N; j++)cout << setw(22) << a[i][j];printf("\n");}printf("\n");}//打印矩阵所有特征值void print_tzz(double lamdaR[N], double lamdaI[N]){int i;printf("\n");for (i = 0; i<N; i++){cout << "\tλ" << setw(2) << i + 1 << " = ";cout << setw(20) << lamdaR[i];if (fabs(lamdaI[i])>E)if (lamdaI[i]>E)cout << " + " << setw(20) << lamdaI[i] << setw(2) << "i";elsecout << " - " << setw(20) << -lamdaI[i] << setw(2) << "i";printf("\n");}printf("\n");}//打印实特征值的特征向量void print_tzxl(double x[N][N], double tzR[N], int nR){int i, j;for (i = 0; i<nR; i++){cout << endl << " 实特征值：" << setw(25) << tzR[i] << endl;cout << " 特征向量：" << endl;for (j = 0; j<int(N / 2); j++)cout << setw(22) << x[i][j];printf("\n");for (j = int(N / 2); j<N; j++)cout << setw(22) << x[i][j];}}int main(){double A[N][N], RQ[N][N], x[N][N];double lamdaR[N], lamdaI[N], tzR[N];int flag, nR;system("Color f0");cout << setiosflags(ios::scientific) << setiosflags(ios::right) << setprecision(12) << endl;create_A(A);cout << " 矩阵A的拟上三角化矩阵A(n-1)：" << endl;nishangsjh(A);print_A(A);printf("\n");cout << " 对矩阵A(n-1)进行QR分解，矩阵R*Q：" << endl;QRdiv(A, RQ);print_A(RQ);printf("\n");cout << " 矩阵A(n-1)经过QR迭代得到的最终矩阵：" << endl;flag = tezhengzhi(A, lamdaR, lamdaI);print_A(A);printf("\n");if (flag == 0)cout << "未求出所有特征值!" << endl;else{cout << " 矩阵A的全部特征值：" << endl;print_tzz(lamdaR, lamdaI);printf("\n");cout << " 实特征值对应的特征向量：" << endl;gauss(x, tzR, nR, lamdaR, lamdaI);print_tzxl(x, tzR, nR);}return 0;}四、计算结果1. 拟上三角矩阵A(n−1)如下：2. 对矩阵A(n−1)进行一次基本QR方法迭代，得到如下矩阵：3. 矩阵A(n−1)经过QR迭代得到的最终矩阵如下：4. 矩阵A的全部特征值如下：5.对应于实特征值的特征向量如下：五、结果分析与讨论1. 由计算结果1和结果2可知，对拟上三角矩阵A(n−1)进行一次基本QR方法迭代，所得矩阵仍然为拟上三角阵，这与两个矩阵相似的理论推导结果相符，证明了拟上三角化算法的正确性。