北航研究生数值分析大作业二(2014年)资料

1、数 值 分 析 计算实习作业二学 院： 17系专 业： 精密仪器及机械姓 名： 张大军学 号： DY1417114一、程序设计方案程序设计方案流程图如图1所示。（注：由本人独立完成，并且有几处算法很巧妙）图1.程序设计方案流程图2、 程序源代码#include <iostream.h>#include <iomanip.h>#include <math.h>#define N 10#define E 1.0e-12#define MAX 10000int main()void nishangsanjiaohua(double (*A)10);void QRf

2、enjie(double (*A)10,double (*Q)N,double (*R)N);void zhengli(double (*A)10);void subuQR(double (*A)10,double *RR,double *II);void tezhengxl(double (*a)N,double T);double A1010=0,Q1010=0,R1010=0;double B10=0,C10=0;int i,j;for(i=1;i<=10;i+)for(j=1;j<=10;j+)if(i!=j)Ai-1j-1=sin(0.5*i+0.2*j);elseAi-

3、1j-1=1.52*cos(i+1.2*j); /对实矩阵A进行拟上三角化nishangsanjiaohua(A);zhengli(A);cout<<"矩阵A经过拟上三角化所得的矩阵A(n-1):"<<endl;for(i=0;i<N;i+)for(j=0;j<N;j+)cout<<setiosflags(ios:scientific)<<setprecision(12)<<setw(21)<<Aij;cout<<'n'<<endl;/拟上三角化后进行

4、的QR分解QRfenjie(A,Q,R);zhengli(R);cout<<"矩阵A(n-1)三角化得到的Q矩阵:"<<endl;for(i=0;i<N;i+)for(j=0;j<N;j+)cout<<setiosflags(ios:scientific)<<setprecision(12)<<setw(21)<<Qij;cout<<'n'<<endl;cout<<"矩阵A(n-1)三角化得到的R矩阵:"<<

5、endl;for(i=0;i<N;i+)for(j=0;j<N;j+)cout<<setiosflags(ios:scientific)<<setprecision(12)<<setw(21)<<Rij;cout<<'n'<<endl;/求解A矩阵的全部特征值subuQR(A,B,C);zhengli(A);cout<<"矩阵A(n-1)双步位移QR迭代后RQ阵:"<<endl;for(i=0;i<N;i+)for(j=0;j<N;j+)c

6、out<<setiosflags(ios:scientific)<<setprecision(12)<<setw(21)<<Aij;cout<<'n'<<endl;cout<<"矩阵A(n-1)双步位移QR迭代后求出的所有特征值:"<<endl;for(i=0;i<N;i+) cout<<setiosflags(ios:scientific)<<setprecision(12)<<setw(21)<<Bi<

7、;<"+"<<setw(21)<<Ci<<"i"<<endl;cout<<'n'<<endl;for(i=1;i<=10;i+)for(j=1;j<=10;j+)if(i!=j)Ai-1j-1=sin(0.5*i+0.2*j);elseAi-1j-1=1.52*cos(i+1.2*j);/A相应于实特征值的特征向量cout<<"矩阵A(n-1)双步位移QR迭代后求出的所有实特征值所对应的特征向量:"<<en

8、dl;for(i=0;i<N;i+)if(Ci=0)cout<<""<<i<<"对应的特征向量"<<endl; tezhengxl(A,Bi);return 1;void zhengli(double (*A)10)int i,j;for(i=0;i<N;i+)for(j=0;j<N;j+)if(fabs(Aij)<=E)Aij=0;void nishangsanjiaohua(double (*A)10)double d,c,h,sum,t;double uN,pN,qN,wN;i

9、nt r,i,j,k;for(r=0;r<=N-3;r+)k=0;/(1)做判断for(i=r+2;i<N;i+)if(Air=0)k+;if(k!=N-r-2)/(2)计算sum=0;for(i=r+1;i<N;i+)sum=sum+Air*Air;d=sqrt(sum);if(Ar+1r=0)c=d;elsec=(-1)*fabs(Ar+1r)/Ar+1r*d;h=c*c-c*Ar+1r;/(3)给u赋值for(i=0;i<N;i+)if(i<=r)ui=0;elseif(i=r+1)ui=Air-c;elseui=Air;/(4)计算for(i=0;i<

10、;N;i+)ui=ui/h;for(i=0;i<N;i+)sum=0;for(j=r+1;j<N;j+)sum=sum+Aji*uj;pi=sum;sum=0;for(j=r+1;j<N;j+)sum=sum+Aij*uj;qi=sum;sum=0;for(i=r+1;i<N;i+)sum=sum+pi*ui;t=sum;for(i=0;i<N;i+)wi=qi-t*ui*h;for(i=0;i<N;i+)for(j=0;j<N;j+)Aij=Aij-wi*uj*h-ui*h*pj;void QRfenjie(double (*A)10,double

11、(*Q)N,double (*R)N) int i,j,k,m;double d,c,h;double uN,wN,pN;for(i=0;i<N;i+)for(j=0;j<N;j+) if (i=j) Qij=1; else Qij=0;for(i=0;i<N;i+)for(j=0;j<N;j+) Rij=Aij;for(i=0;i<N-1;i+)for(j=i+1;j<N;j+) if(Rji<=E) m=m+1;if(m=(N-1-i) continue;elsefor(j=i,d=0;j<N;j+) d=d+Rji*Rji;d=sqrt(d

12、); c=-1*fabs(Rii)/Rii*d; h=c*c-c*Rii; for(j=i+1;j<N;j+) uj=Rji; for(j=0;j<i;j+) uj=0;ui=Rii-c;for(j=0;j<N;j+) for(k=0,wj=0;k<N;k+)wj=Qjk*uk+wj;for(j=0;j<N;j+) for(k=0;k<N;k+) Qjk=Qjk-wj*uk/h;for(j=0;j<N;j+) for(k=i,pj=0;k<N;k+)pj=Rkj*uk+pj; pj=pj/h;for(j=0;j<N;j+) for(k=0;

13、k<N;k+) Rjk=Rjk-uj*pk;/矩阵的QR分解 double kaifang(double b,double c) double m;m=b*b-4*c;return m;/使用双步位移QR法求实矩阵A的全部特征值void subuQR(double (*A)10,double *RR,double *II) int m=N-1,BU=3,i,j; int L=1; int k=0; double s,t,x; double MNN,BNN; int f=0; double d,c,h; double uN,wN,pN; double QNN,RNN;while(BU!=1

14、1) /编程精妙之处*if(BU=3)if(fabs(Amm-1)<=E) RRm=Amm;IIm=0; m=m-1; BU=4;elseBU=5;if(BU=4)if(m=0) RRm=Amm; IIm=0; BU=11;elseBU=3;if(BU=5)s=Am-1m-1+Amm; t=Am-1m-1*Amm-Amm-1*Am-1m;x=kaifang(s,t);if(x>=0)x=sqrt(x);RRm=(s-x)/2;RRm-1=(s+x)/2;IIm=0;IIm-1=0;elsex=sqrt(-x);RRm=s/2;RRm-1=s/2;IIm=x/2;IIm-1=-x/2

15、;BU=6; if(BU=6)if(m=1)BU=11;elseBU=7;if(BU=7) if(fabs(Am-1m-2)<=E) m=m-2;BU=4; else BU=8;if(BU=8)if(L=MAX)BU=11;elseBU=9;if(BU=9)for(i=0;i<=m;i+) for(j=0;j<=m;j+) if(i=j) for(k=0,Mij=0;k<=m;k+) Mij=Aik*Akj+Mij; Mij=Mij-s*Aij+t; else for(k=0,Mij=0;k<=m;k+) Mij=Aik*Akj+Mij; Mij=Mij-s*Ai

16、j; /以下是M的QR分解 for(i=0;i<=m;i+) for(j=0;j<=m;j+) if (i=j) Qij=1; else Qij=0; for(i=0;i<=m;i+) for(j=0;j<=m;j+) Rij=Mij; for(i=0;i<m;i+) for(j=i+1;j<=m;j+) if(Rji<=E) f=f+1; if(f=(m-i) continue; for(j=i,d=0;j<=m;j+) d=d+Rji*Rji; d=sqrt(d); c=-1*fabs(Rii)/Rii*d; h=c*c-c*Rii; for

17、(j=i+1;j<=m;j+) uj=Rji; for(j=0;j<i;j+) uj=0; ui=Rii-c; for(j=0;j<=m;j+) for(k=0,wj=0;k<=m;k+) wj=Qjk*uk+wj; for(j=0;j<=m;j+) for(k=0;k<=m;k+) Qjk=Qjk-wj*uk/h; for(j=0;j<=m;j+) for(k=i,pj=0;k<=m;k+) pj=Rkj*uk+pj;pj=pj/h; for(j=0;j<=m;j+) for(k=0;k<=m;k+) Rjk=Rjk-uj*pk;

18、for(j=0;j<=m;j+) for(k=0;k<=m;k+) Mjk=Qjk; for(i=0;i<=m;i+) for(j=0;j<=m;j+) for(k=0,Bij=0;k<=m;k+) Bij=Mki*Akj+Bij; for(i=0;i<=m;i+) for(j=0;j<=m;j+) for(k=0,Aij=0;k<=m;k+) Aij=Bik*Mkj+Aij;BU=10;if(BU=10)L=L+1;BU=3;void tezhengxl(double (*a)N,double T)void qstzxl(double (*a)N);double MNN;int i,j;for(i=0;i<N;i+)for(j=0;j<N;j+)if(i=j)Mij=aij-T;elseMij=aij;qstzxl(M);/巧妙的使用老师上课要求上机调试练习的列主元高斯消去法求解实特征值对应的特征函数void qstzxl(double (*a)N)double bN=0; double HNN=0,lN=0;double XN;double B;double sum;int i,j,m,k,z;for(k=0;k<N-1;k+)for(j=k;