北航数值分析A第二题.docx

1、北航数值分析A第二题数值分析A计算实习第二题一、设计方案 1.定义矩阵。由函数A(double (*a)N)定义矩阵中的元素，方便初始化。 2.对矩阵进行拟上三角化，得到。由函数AN(double (*a)N)对进行拟上三角化，减少计算量。 3.对采用带双步位移的分解法求解特征值。在每次分解前，如果可以直接得到矩阵的一个特征值或者一对特征值，则对矩阵降阶。函数SQR(double (*a)N,double (*lc)2)代表带双步位移的分解，其中lcN0存储特征值的实部，lcN1存储特征值的虚部。函数HAH(double (*a)N,double *u,double h,int n)通过参数u

2、N、h的传递对矩阵左乘、右乘Householder矩阵。 4.采用列主元Gauss消去法求实特征值对应的特征向量。对初始化后，采用列主元素Gauss消去法求的解，消元过程中，矩阵的最后一行元素为0，令向量的最后一个元素为1，回代可求得对应于特征值的特征向量。所用函数为G(double *x,double lc)。 5.对初始化，研究拟上三角化对一般分解法求得收敛的分块上三角阵所需分解次数的影响，对比分块上三角阵一阶、二阶对角块给出的特征值与带双步位移的分解法求解结果的一致性。所用函数为QR(double (*a)N,double (*q)N,double (*r)N)，采用一阶对角块的元素、二

3、阶对角块的迹组成的向量，通过判断，确定分解的收敛性。二、源代码#include#include #define N 10 /矩阵阶数#define K 50000 /最大迭代次数#define E 1e-12 /精度void A(double (*a)N) /初始化矩阵A int i,j; for(i=0;iN;i+) for(j=0;jN;j+) aij=i!=j?sin(0.5*i+0.2*j+0.7):1.5*cos(i+1.2*j+2.2);void HAH(double (*a)N,double *u,double h,int n) /对矩阵A左乘右乘Householder矩阵 do

4、uble pN,qN,t,wN,s1,s2; int i,j; for(i=0;in;i+) s1=0.0; s2=0.0; for(j=0;jn;j+) s1+=aji*uj/h; s2+=aij*uj/h; pi=s1; qi=s2; t=0.0; for(i=0;in;i+) t+=pi*ui; t/=h; for(i=0;in;i+) wi=qi-t*ui; for(i=0;in;i+) for(j=0;jn;j+) aij-=(wi*uj+ui*pj);void AN(double (*a)N) /拟上三角化 double uN,d,c,h; int r,i,t; for(r=1;r

5、=N-2;r+) t=0; for(i=r+2;i=N;i+) if(fabs(ai-1r-1)E) t+=1; if(tN-r-1) d=0.0; for(i=r+1;i=N;i+) d+=ai-1r-1*ai-1r-1; d=sqrt(d); c=arr-1=0?d:-d; h=c*(c-arr-1); for(i=1;i=N;i+) ui-1=i1) if(fabs(an-1n-2)E) lcN-n0=an-1n-1; lcN-n1=0.0; printf(QR分解次数为%d，矩阵阶数为%d时，由最后1行倒数第2个元素为0，求得特征值n%d=%.12en,k,n,N-n+1,lcN-n0

6、); printf(此时矩阵A的%d阶顺序主子式为n,n); for(i=0;in;i+) for(j=0;j(nN/2?n:N/2);j+) printf(%20.12et,fabs(aij)N/2) for(i=0;in;i+) for(j=N/2;jn;j+) printf(%20.12et,fabs(aij)E?0:aij); printf(n); printf(n); n-=1; else /解特征值s1，s2 b=-an-2n-2-an-1n-1; c=an-2n-2*an-1n-1-an-1n-2*an-2n-1; d=b*b-4*c; if(d0) s1r=-b/2.0; s2

7、r=-b/2.0; s1i=sqrt(-d)/2.0; s2i=-sqrt(-d)/2.0; else s1r=(-b+sqrt(d)/2.0; s2r=(-b-sqrt(d)/2.0; s1i=0.0; s2i=0.0; if(n=2|fabs(an-2n-3)E) lcN-n0=s1r; lcN-n1=s1i; lcN-n+10=s2r; lcN-n+11=s2i; if(n=2) if(fabs(s1i)E) printf(QR分解次数为%d，矩阵阶数为%d时，直接求得特征值n%d=%.12en%d=%.12en,k,n,N-n+1,lcN-n0,N-n+2,lcN-n+10); els

8、e printf(QR分解次数为%d，矩阵阶数为%d时，直接求得特征值n%d=%.12e+%.12ein%d=%.12e+%.12ein,k,n,N-n+1,lcN-n0,lcN-n1,N-n+2,lcN-n+10,lcN-n+11); else if(fabs(s1i)E) printf(QR分解次数为%d，矩阵阶数为%d时，由倒数第2行倒数第3个元素为0，求得特征值n%d=%.12en%d=%.12en,k,n,N-n+1,lcN-n0,N-n+2,lcN-n+10); else printf(QR分解次数为%d，矩阵阶数为%d时，由倒数第2行倒数第3个元素为0，求得特征值n%d=%.12

9、e+%.12ein%d=%.12e+%.12ein,k,n,N-n+1,lcN-n0,lcN-n1,N-n+2,lcN-n+10,lcN-n+11); printf(此时矩阵A的%d阶顺序主子式为n,n); for(i=0;in;i+) for(j=0;j(nN/2?n:N/2);j+) printf(%20.12et,fabs(aij)N/2) for(i=0;in;i+) for(j=N/2;jn;j+) printf(%20.12et,fabs(aij)E?0:aij); printf(n); printf(n); n-=2; else if(kK) s=an-2n-2+an-1n-1;

10、 t=an-2n-2*an-1n-1-an-1n-2*an-2n-1; for(i=0;iN;i+) for(j=0;jN;j+) sum=0.0; for(r=0;rN;r+) sum+=air*arj; mij=sum-s*aij+t*(i=j?1:0); for(r=1;rn;r+) sum=0.0; for(i=r+1;i=n;i+) if(fabs(mi-1r-1)E) sum+=1; if(sumn-r) d=0.0; for(i=r;i=n;i+) d+=mi-1r-1*mi-1r-1; d=sqrt(d); c=mr-1r-1=0?d:-d; h=c*(c-mr-1r-1);

11、for(i=1;i=n;i+) ui-1=ir?0:mi-1r-1; ur-1=mr-1r-1-c; for(i=0;in;i+) sum=0.0; for(j=0;jn;j+) sum+=mji*uj; vi=sum/h; for(i=0;in;i+) for(j=0;jn;j+) mij=mij-ui*vj; HAH(a,u,h,n); k+=1; else printf(未求出全部特征值n); break; if(n=1) lcN-n0=a00; lcN-n1=0.0; printf(QR分解次数为%d，矩阵阶数为%d时，直接求得特征值n%d=%.12en,k,n,N-n+1,lcN-n

12、0); printf(此时矩阵A的1阶顺序主子式为n%20.12enn,a00); void G(double *x,double lc) /列主元素gauss消去法求特征向量 double aNN,t,m; int i,j,k,ik; A(a); for(i=0;iN;i+) for(j=0;jN;j+) aij-=i=j?lc:0; for(k=1;kN;k+) t=0.0; for(i=k;i=fabs(t) t=ai-1k-1; ik=i; for(j=0;jN;j+) t=ak-1j; ak-1j=aik-1j; aik-1j=t; for(i=k+1;i=N;i+) m=ai-1k

13、-1/ak-1k-1; for(j=k+1;j0;k-) t=0.0; for(j=k+1;j=N;j+) t+=ak-1j-1*xj-1; xk-1=-t/ak-1k-1; void QR(double (*a)N,double (*q)N,double (*r)N) double d,c,h,t,s,uN,wN,pN; int i,j,k; for(i=0;iN;i+) for(j=0;jN;j+) qij=i=j?1:0; rij=aij; for(k=1;kE) d=0.0; for(i=k;i=N;i+) d+=ai-1k-1*ai-1k-1; d=sqrt(d); c=ak-1k-

14、1=0?d:-d; h=c*(c-ak-1k-1); for(i=1;i=N;i+) ui-1=ik?0:ai-1k-1; uk-1=ak-1k-1-c; for(i=0;iN;i+) t=0.0; s=0.0; for(j=0;jN;j+) t+=qij*uj; s+=rji*uj; wi=t; pi=s/h; for(i=0;iN;i+) for(j=0;jN;j+) qij-=wi*uj/h; rij-=ui*pj; HAH(a,u,h,N); void main() double aNN,lcN2,xN,t,qNN,rNN,bNN; /lc存储复特征值 int i,j,k,m=0; A

15、(a); /定义矩阵A AN(a); /对矩阵A进行拟上三角化 printf(对矩阵A进行拟上三角化n); for(i=0;iN;i+) for(j=0;jN/2;j+) printf(%20.12et,fabs(aij)E?0:aij); printf(n); printf(n); for(i=0;iN;i+) for(j=N/2;jN;j+) printf(%20.12et,fabs(aij)E?0:aij); printf(n); printf(n); SQR(a,lc); /对拟上三角矩阵A进行带双步位移的QR分解，求解特征值 for(i=0;iN;i+) if(fabs(lci1)E

16、) m+=1; printf(所有%d个实特征值及其对应的特征向量n,m); for(i=0,j=1;iN;i+) if(fabs(lci1)E) printf(%d=%.12en,j,lci0); printf( X%d=,j); G(x,lci0); /gauss消去法求特征向量 for(k=0;kN;k+) printf(t%20.12e,xk); if(k=N/2-1|k=N-1) printf(n); j+; printf(所有%d对复特征值n,(N-m)/2); for(i=0,j=m+1;iE) printf(%d,%d=%.12e%.12ein,j,j+1,lci0,lci1)

17、; i+; j+=2; A(a); /初始化矩阵A AN(a); /对矩阵A进行拟上三角化 for(k=0,t=1;kE) for(i=0;iN;i+) for(j=0;jN;j+) bij=aij; QR(a,q,r); /对拟上三角阵A进行第k次QR分解，求Q、R及乘积矩阵RQ else printf(nQR分解次数为%dn,k); break; t=0.0; for(i=0;iN;i+) if(fabs(ai+1i)E) t+=(bii-aii)*(bii-aii); else t+=(bii+bi+1i+1-aii-ai+1i+1)*(bii+bi+1i+1-aii-ai+1i+1);

18、 i+=1; t=sqrt(t); /以矩阵B-A的一阶对角块和二阶对角块的迹构成的向量的2-范数判断QR分解法收敛 printf(第%d次QR分解中，R=n,k); /输出第k次分解时的R for(i=0;iN;i+) for(j=0;jN/2;j+) printf(%20.12et,fabs(rij)E?0:rij); printf(n); printf(n); for(i=0;iN;i+) for(j=N/2;jN;j+) printf(%20.12et,fabs(rij)E?0:rij); printf(n); printf(n); printf(第%d次QR分解中，Q=n,k); /

19、输出第k次分解时的Q for(i=0;iN;i+) for(j=0;jN/2;j+) printf(%20.12et,fabs(qij)E?0:qij); printf(n); printf(n); for(i=0;iN;i+) for(j=N/2;jN;j+) printf(%20.12et,fabs(qij)E?0:qij); printf(n); printf(n); printf(第%d次QR分解中，RQ=n,k); /输出第k次分解时的RQ for(i=0;iN;i+) for(j=0;jN/2;j+) printf(%20.12et,fabs(aij)E?0:aij); print

20、f(n); printf(n); for(i=0;iN;i+) for(j=N/2;jN;j+) printf(%20.12et,fabs(aij)E?0:aij); printf(n); printf(n);三、计算结果 拟上三角矩阵 带双步位移的分解法求特征值过程 所有特征值及实特征值对应的特征向量 对采用分解法求、及，下面只给出最后一次迭代的结果四、计算结果分析 1.对矩阵进行拟上三角化可以减少分解的迭代次数。对矩阵直接进行分解，得到收敛的分块上三角阵所需次数为866次；对拟上三角阵进行分解，得到收敛的分块上三角阵所需次数为489次。 2. 分解法得到收敛的分块上三角阵，由一阶、二阶对角块求得的特征值与带双步位移的分解法求得的特征值一致。 3.带双步位移的分解法通过降阶大大减少了QR分解的次数。采用带双步位移的分解法求的特征值，只需14次QR分解。