数值分析第二题北航 大作业

1、数值分析计算实习题目二数值分析第二题目 录数值分析第二题11 引言：21.1 矩阵的拟上三角化21.2 矩阵的特征值求解31.3 矩阵的特征向量求解52 算法的程序实现62.1 主程序62.2 子程序的实现83计算结果93.2矩阵Q、R 以及乘积RQ103.3 各实特征值及其相对应的特征向量114 实验结论12附录 源代码131 引言1.1 矩阵的拟上三角化为了减少计算量，对矩阵利用Householder矩阵进行相似变换，把化为上三角矩阵A(n-1)。对拟上三角化，得到拟上三角矩阵，具体算法如下：记，并记的第r列至第n列的元素为。对于执行1. 若全为零，则令,转5；否则转2。2. 计算3. 令

2、。4. 计算 5. 继续。1.2 矩阵的特征值求解使用带双步位移的QR方法计算矩阵的全部特征值，也是的全部特征值，具体算法如下：1. 给定精度水平和迭代最大次数。2. 记，令。3. 如果，则得到的一个特征值，置（降阶），转4；否则转5。4. 如果，则得到的一个特征值，转11；如果，则转3。5. 求2阶子阵的两个特征值和，即计算二次方程的两个根和。6. 如果，则得到的两个特征值和，转11；否则转7。7. 如果，则得到的两个特征值和，置（降阶），转4；否则转88. 如果，则计算终止，未得到的全部特征值；否则转9。9. 记，计算10. 置，转3。11. 的全部特征值已计算完毕，停止计算。其中，的分解

3、与的计算用下列算法实现：记。对于执行1. 若全为零，则令,转5；否则转2。2. 计算3. 令。4. 计算5. 继续。此算法执行完后，就得到。1.3 矩阵的特征向量求解用列主元素Gauss消去法计算矩阵对应于实特征值的特征向量，具体算法如下：记1. 消元过程对于执行(1) 选行号，使。(2) 交换与所含的数值。(3) 对于计算2. 回代过程最终得到的向量的即为对应于实特征值的特征向量。2 算法的程序实现2.1 主程序#include <iostream>#include"matrix-tool.h"#include"simi-trian.h"#

4、include"QR_decompose.h"#include"Gauss_by_max.h"#define threhold (1e-12)using namespace std; int n = 10;extern double *a; extern double *Mk; extern double *a_lambda; extern double *b_lambda; double *spec_u;/高斯求解方程组的右端常数列向量 double *G_b; int main() simi_Initial();cout<<"处

5、理之前的A矩阵为："<<endl; Print_Matrix(a,n);cout<<"拟上三角化后的A矩阵为："<<endl;pseudo_triangle(a,n); BI_STEP_QR(a,n);/*特征向量求取代码段*/初始化spec_u10spec_u = new double n;G_b = new double n; char i;simi_Initial();/为G_b 赋值for(i=0; i<n; i+) (*(G_b+i) = 0; for(i=0; i<n; i+) if( (*(b_lamb

6、da+i) = 0 ) Gauss_root( a, G_b, n, (*(a_lambda+i),spec_u);cout<<"属于"<<(*(a_lambda+i)<<" i "<<(*(b_lambda+i)<<"的特征向量为： "<<endl;Print_Vector(spec_u,n); system("pause"); return 0;2.2 子程序的实现2.2.1 “matrix-tool.h”void Matrix_M_Vec

7、tor(double *A,double *b,double *y,int n);函数实现 N*N的矩阵乘以N*1的向量，得到N*1向量Yvoid Matrix_Converted_M_Vector(double *A,double *b,double *y,int n);函数实现转置后N*N矩阵乘以N*1向量,得到N*1向量Yvoid Column_M_Row(double *a,double *b,double *A,int n); 函数实现列 乘以 行 得到一个矩阵Avoid Row_M_Column(double *a,double *b,double *y,int n); 函数实现行

8、 乘以 列 得到一个数值yvoid Print_Matrix(double *a,int n); 函数实现将矩阵a显示出来void Print_Vector(double *b,int n); 函数实现将向量b显示出来2.2.2 “Gauss_by_max.h”void Gauss_root(double *a,double *b,int n,double lambda,double *u)这个函数用于使用高斯列主元消元法实现方程组的求解。M_switch(double *a,double *b,int n,int k,int flag)这个函数用于实现按照列主元的绝对值大小重新排列方程组2.

9、2.3 “QR_decompose.h”void BI_STEP_QR(double *a,int n)这个函数用于实现矩阵的双步位移QR分解。2.2.4 “simi-trian.h”void pseudo_triangle(double *A,int n); 这个函数主要实现矩阵的拟上三角化void simi_Initial();这个函数主要用于按照题目条件实现矩阵的初始化3计算结果3.1 拟三角化后的矩阵如图1所示：图1拟三角化后的矩阵A3.2矩阵Q、R 以及乘积RQ图2 QR分解后得到的RQ矩阵图3 QR分解后得到的Q矩阵图4 QR分解后得到的R矩阵3.3 各实特征值及其相对应的特征向量

10、得到的特征向量如下图5所示：图5 计算得到A的特征值得到的A的各特征值对应的特征向量为：图6 A的各特征值对应的特征向量图7 A的各特征值对应的特征向量4 实验结论本次作业运用带双步位移的QR分解法求出了给定矩阵的全部特征值和对应于实特征值的特征向量，通过本人调试后，达到了题目的设计要求。本次调试过程中，还是发现了比较多的问题。1. 当一个数据被多个程序段进行操作的时候，需要在具体程序中，对被操作的数据另行创建副本，以避免这种操作改变了原始数据而影响了其他的程序。我在本次程序设计过程中，就出现这种问题：在对矩阵A进行QR分解求取特征值的程序和后来求取特征向量的程序中，都没有创建原始变量的副本，

11、导致程序之间互相干扰。2. 在调试过程中，由于我是使用指针来操作数组的，所以整个程序对二维数组的维数非常敏感。例如我在调试过程中，在完成矩阵A的QR分解过程中，原来是通过语句来*MkC+i*m+j= *a+i*m+j实现的，最后发现程序在m=10时，程序运行正确，但是当m=8后，程序就出现问题了。最后才知道这句赋值语句出现了问题，应该是*MkC+i*m+j= *a+i*n+j。3. 在调试过程中，我发现我写的函数间的耦合性非常大，所使用的全局变量也特别多。导致整个程序调试起来非常不方便。以后可以通过减少全局变量的个数等其他方式来降低程序的耦合性。附录 源代码1. “matrix-tool.c”

12、#include"matrix-tool.h"#include<math.h>#include<iostream>#include<iomanip>using namespace std;void Matrix_M_Vector(double *A ,double *b,double *y,int n)int i,j;double sum = 0;for(i=0; i<n; i+)for(j=0; j<n; j+)sum += (*(A+i*n+j) )*(*(b+j);(*(y+i) = sum ; sum = 0;void

13、 Matrix_Converted_M_Vector(double *A,double *b,double *y,int n) int i,j;double sum = 0;for(i=0; i<n; i+)for(j=0; j<n; j+)sum += (*(A+j*n+i) )*(*(b+j);(*(y+i) = sum ; sum = 0;void Column_M_Row(double *a,double *b,double *A,int n) /列 乘以 行 得到一个矩阵Aint i,j;for(i=0; i<n; i+)for(j=0; j<n; j+) (

14、*(A+i*n+j) ) = (*(a+i)*(*(b+j);double Row_M_Column(double *a,double *b,int n) /行 乘以 列 得到一个数值yint i;double sum = 0;for(i=0; i<n; i+) sum += (*(a+i)*(*(b+i);return sum;void Print_Matrix(double *a,int n)int i,j;for(i=0; i<n; i+)for(j=0; j<n; j+)cout<<setiosflags(ios:scientific)<<se

15、tprecision(12)<<""<<i<<""<<""<<j<<"= "<<(*(a+i*n+j)<<" "cout<<endl;cout<<endl;void Print_Vector(double *b,int n) int i; for(i=0; i<n; i+) cout<<setw(15)<<setiosflags(ios:scien

16、tific)<<setprecision(6)<<(*(b+i)<<" "<<endl;2. “Gauss_by_max.c”#include"Gauss_by_max.h"#include"matrix-tool.h"#include<iostream>#include"math.h"using namespace std;/extern double *a_lambda;/extern double *b_lambda;void M_switch(do

17、uble *a,double *b,int n,int k,int flag)int i;double temp;for(i=0; i<n; i+) temp = (*(a+k*n+i); (*(a+k*n+i) = (*(a+flag*n+i);(*(a+flag*n+i) = temp; temp = (*(b+k); (*(b+k) = (*(b+flag);(*(b+flag) = temp;return ; void Gauss_root(double *a,double *b,int n,double lambda,double *u) double *G_a;double

18、*G_b;double G_temp;double G_sum;int G_flag = 0;G_a = new doublen*n;G_b = new doublen;int k,i,j;for(i=0; i<n; i+)for(j=0; j<n; j+)if(i=j) (*(G_a+i*n+j) = (*(a+i*n+j)-lambda;else (*(G_a+i*n+j) = (*(a+i*n+j); (*(G_b+i) = (*(b+i);for(k=0; k<(n-1); k+)G_temp = fabs(*(G_a+k*n+k); G_flag = k;for(i

19、=(k+1); i<n; i+)if( G_temp<fabs(*(G_a+i*n+k) )G_flag = i;G_temp = (*(G_a+i*n+k);M_switch(G_a,G_b,n,k,G_flag); for(k=0; k<(n-1); k+) for(i=(k+1); i<n; i+)G_temp = (*(G_a+i*n+k)/(*(G_a+k*n+k);for(j=(k); j<n; j+)(*(G_a+i*n+j) -= G_temp*(*(G_a+k*n+j);(*(G_b+i) -= G_temp*(*(G_b+k); for(i=0

20、; i<n; i+) (*(u+i) = 0; (*(u+n-1) = 1;for(k=(n-2); k>=0; k-)G_sum = 0;for(j=(k+1); j<n; j+) G_sum += (*(G_a+k*n+j) * (*(u+j);(*(u+k) = ( (*(G_b+k) - G_sum )/(*(G_a+k*n+k);G_sum = 0;return ; 3. “QR_decompose.c”#include"QR_decompose.h"#include"matrix-tool.h"#include"s

21、imi-trian.h"#include<iostream>#include"math.h"#define threhold (1e-12)using namespace std; int m;extern double *a; double *a1; double *Mk;double *Mk_B;double *Mk_C;double *Mk_pr;double *Mk_qr; double *Mk_wr; double *Mk_ur;double *Mk_vr;double Mk_tr = 0;double Mk_hr = 0;double Mk

22、_dr = 0;double Mk_cr = 0;char complex_flag = 0;char lambda_flag = 0;double *a_lambda;double *b_lambda;double lambda2 = 0 ;void root_lambda(double *a,int n,int m,double *lambda) double det = 0; double delta = 0; double tra = 0; double temp = 0;/求取A右下角的2*2行列式的值 det = ( *(a+(m-2)*n+(m-2) )*( *(a+(m-1)*

23、n+(m-1) ); det -= ( *(a+(m-1)*n+(m-2) )*( *(a+(m-2)*n+(m-1) );/求取A右下角的2*2矩阵的迹 tra = ( *(a+(m-2)*n+(m-2) )+( *(a+(m-1)*n+(m-1) );/求取一元二次方程的delta值 temp = tra*tra-4*det; if(temp>=0) delta = sqrt(tra*tra-4*det); lambda0 = (tra-delta)/2; lambda1 = (tra+delta)/2; if(m=2) cout<<"一个特征根是："

24、<<(tra-delta)/2<<endl; cout<<"另一个特征根是："<<(tra+delta)/2<<endl; (*(a_lambda+lambda_flag) ) = lambda0; (*(b_lambda+lambda_flag) ) = 0; lambda_flag +;(*(a_lambda+lambda_flag) ) = lambda1; (*(b_lambda+lambda_flag) ) = 0; lambda_flag +; else complex_flag = 1; delta

25、= sqrt(-tra*tra+4*det); lambda0 = tra; lambda1 = delta;if(m=2) cout<<"一个复特征根是："<<tra/2<<"-i"<<delta/2<<endl; cout<<"另一个复特征根是："<<tra/2<<"+i"<<delta/2<<endl; (*(a_lambda+lambda_flag) ) = tra/2; (*(b_lam

26、bda+lambda_flag) ) = delta/2; lambda_flag +;(*(a_lambda+lambda_flag) ) = tra/2; (*(b_lambda+lambda_flag) ) = -delta/2; lambda_flag +; return ;void step9(double *a,int n,int m)int i,j,r;double s,t;double a_pow_2 = 0;double temp = 0;Mk = new double m*m;/初始化 cmm矩阵Mk_C = new double m*m;for(i=0; i<m;

27、i+)for(j=0; j<m; j+)(*(Mk_C+i*m+j) = (*(a+i*n+j);/ 求取ss = ( *(a+(m-2)*n+(m-2) )+( *(a+(m-1)*n+(m-1) );/求取tt = ( *(a+(m-2)*n+(m-2) )*( *(a+(m-1)*n+(m-1) );t -= ( *(a+(m-2)*n+(m-1) )*( *(a+(m-1)*n+(m-2) ) );/求取M矩阵 同时求取Mk_B矩阵for(i=0; i<m; i+)for(j=0; j<m; j+)for(r=0; r<m; r+) a_pow_2 += ( *

28、(a+i*n+r) )*( *(a+r*n+j) ); /求取Ak的平方中的元素if(i=j) temp = 1; /求取单位矩阵的元素值else temp = 0;(*(Mk+i*m+j) ) = a_pow_2 - s*( *(a+i*n+j) ) + t*temp;temp = 0;a_pow_2 = 0;/初始化 bmm矩阵 Mk_B = new double m*m;for(i=0; i<m; i+)for(j=0; j<m; j+) (*(Mk_B+i*m+j) ) = (*(Mk+i*m+j) ) ; return;/实现 pr,qr,tr,wr和A(k+1)矩阵的求

29、取void Mk_QR_step3(double *cr,double *Mk,int n,int m,int r) int i,j; /实现pr的求取 Matrix_Converted_M_Vector(cr,Mk_ur,Mk_pr,m); for(i=0; i<m;i+) (*(Mk_pr+i) /= Mk_hr; /实现qr的求取 Matrix_M_Vector(cr,Mk_ur,Mk_qr,m); for(i=0; i<m;i+) (*(Mk_qr+i) /= Mk_hr; /实现tr的求取Mk_tr = Row_M_Column(Mk_pr,Mk_ur,m); Mk_tr

30、 /= Mk_hr; /实现wr的求取 for(i=0; i<m; i+) (*(Mk_wr+i) = (*(Mk_qr+i) - Mk_tr*(*(Mk_ur+i); /实现C(R+1)的求取; C(r+1) = C（r)- wr * urT - ur*prT for(i=0; i<m; i+) for(j=0; j<m; j+)(*(cr+i*m+j) -= ( (*(Mk_wr+i)*(*(Mk_ur+j)+(*(Mk_ur+i)*(*(Mk_pr+j) ); /进行Mk_dr,Mk_cr,Mk_hr,Mk_ur的求取void Mk_QR_step2(double *c

31、r,double *br,int n,int m,int r)int i, j;double sum = 0;/ 求取 Mk_dr for(j=(r); j<m; j+) sum += pow( (*(br+j*m+r) ,2); Mk_dr = sqrt(sum); sum = 0;/ 实现Mk_cr = -sgn(ar+1r）*Mk_drif ( (*(br+(r)*m+r) <= 0) Mk_cr = Mk_dr; else Mk_cr = -Mk_dr; / 实现 Mk_hr 值的求取 Mk_hr = pow(Mk_cr,2) - Mk_cr*( *(br+(r)*m+r)

32、 ); / 实现 Mk_ur 值的求取 for(i=0; i<m; i+) if (i<r) (*(Mk_ur+i) = 0 ;else if ( i=r ) (*(Mk_ur+i) = (*(br+i*m+r) - Mk_cr;else (*(Mk_ur+i) = (*(br+i*m+r); /实现Mr_vr的求取 Matrix_Converted_M_Vector(br,Mk_ur,Mk_vr,m); for(i=0; i<m;i+) (*(Mk_vr+i) /= Mk_hr; /实现B(r+1)矩阵的求解 for(i=0; i<m; i+) for(j=0; j&

33、lt;m; j+)(*(br+i*m+j) -= (*(Mk_ur+i)*(*(Mk_vr+j); Mk_QR_step3(cr,Mk,n,m,r);void Mk_QR(double *A,double *Mk,int n,int m)int r,i,j; Mk_pr = new double m;Mk_qr = new double m;Mk_wr = new double m;Mk_ur = new double m;Mk_vr = new double m;char flag = 0; for(r=0; r<(m-1); r+) mark1: for(i=(r+1); i<

34、m; i+) if( fabs( (*(Mk+i*m+r) ) = 0 ) flag += 1; if( flag != (m-r-1) ) Mk_QR_step2(Mk_C,Mk_B,n,m,r); flag = 0; else r+;goto mark1; void BI_STEP_QR(double *a,int n)int L = 100; /预置迭代次数为1000次 int k ,i,j;a_lambda = new double n;b_lambda = new double n; k = 1;m = n; /初始化 a1nn矩阵 a1 = new double n*n; for(

35、i=0; i<n; i+)for(j=0; j<n; j+)(*(a1+i*n+j) = (*(a+i*n+j); for(j=0; j<n; j+)(*(a_lambda+j) = 0;(*(b_lambda+j) = 0;stage3:if( fabs(*(a1+(m-1)*n+m-2) < threhold )(*(a_lambda+lambda_flag) ) = (*(a1+(m-1)*n+(m-1);(*(b_lambda+lambda_flag) ) = 0;lambda_flag +;m -;goto stage4;else goto stage5;st

36、age4: if(m=1)(*(lambda+n-1) ) = (*(a1);(*(a_lambda+lambda_flag) ) = (*a1 );(*(b_lambda+lambda_flag) ) = 0;lambda_flag +;goto stage12; else goto stage3;stage5: root_lambda(a1, n, m, lambda);if(m=2)goto stage12;else goto stage7;stage7:if( fabs(*(a1+(m-2)*n+m-3) < threhold )if(complex_flag=1)(*(a_la

37、mbda+lambda_flag) ) = lambda0/2;(*(b_lambda+lambda_flag) ) = lambda1/2; lambda_flag +;(*(a_lambda+lambda_flag) ) = lambda0/2; (*(b_lambda+lambda_flag) ) = -lambda1/2; lambda_flag +;complex_flag = 0;else(*(a_lambda+lambda_flag) ) = lambda0;(*(b_lambda+lambda_flag) ) = 0; lambda_flag +;(*(a_lambda+lam

38、bda_flag) ) = lambda1;(*(b_lambda+lambda_flag) ) = 0; lambda_flag +;m -= 2;goto stage4;else goto stage8;stage8:if(k=L)goto stage11;else goto stage9;stage9: step9(a1,n,m); Mk_QR(Mk_C,Mk_B,n,m); for(i=0; i<m; i+) for(j=0; j<m; j+)(*(a1+i*n+j) = (*(Mk_C+i*m+j); /用于显示A(n-1)的QR分解后的结果， if(m=n) /显示得到

39、的乘积矩阵RQ cout<<"QR分解后得到的乘积矩阵RQ为："<<endl; Print_Matrix(a1,m); / 显示得到的Q矩阵 cout<<"QR分解后得到的Q矩阵为："<<endl; Print_Matrix(Mk_B,m); /显示得到的R矩阵 cout<<"QR分解后得到的R矩阵为："<<endl; Print_Matrix(Mk_C,m); k+;goto stage3;stage11:cout<<"未能得到A的全部特征

40、值"<<endl;stage12:cout<<"A的全部特征值已计算完毕，停止计算"<<endl;cout<<"得到的十个特征值为："<<endl;for(j=0; j<n; j+)cout<<"lambda"<<j<<""<<(*(a_lambda+j)<<" i "<<(*(b_lambda+j);cout<<endl;4. “sim

41、i-trian.c”#include<iostream>#include"simi-trian.h"#include"matrix-tool.h"#include"math.h"using namespace std; extern int n; double *a ;double *tri_pr;double *tri_qr; double *tri_wr; double *tri_ur;double tri_tr = 0;double tri_hr = 0;double tri_dr = 0;double tri_cr = 0; /实现对A=a1010矩阵的初始化和赋值 void simi_Initial () char i,j; a = new double n*n; for(i=0; i<n; i+) for(j=0; j<n; j+)if( i!=j )(*(a+i*n+j) = sin(0.5*(i+1)+0.2*(j+1);else(*(a+