北航研究生數(shù)值分析大作業(yè)二(2014年)資料

1、數(shù) 值 分 析 計(jì)算實(shí)習(xí)作業(yè)二學(xué) 院： 17系專(zhuān) 業(yè)： 精密儀器及機(jī)械姓 名： 張大軍學(xué) 號(hào)： DY1417114一、程序設(shè)計(jì)方案程序設(shè)計(jì)方案流程圖如圖1所示。（注：由本人獨(dú)立完成，并且有幾處算法很巧妙）圖1.程序設(shè)計(jì)方案流程圖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); /對(duì)實(shí)矩陣A進(jìn)行擬上三角化nishangsanjiaohua(A);zhengli(A);cout<<"矩陣A經(jīng)過(guò)擬上三角化所得的矩陣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;/擬上三角化后進(jìn)行

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相應(yīng)于實(shí)特征值的特征向量cout<<"矩陣A(n-1)雙步位移QR迭代后求出的所有實(shí)特征值所對(duì)應(yīng)的特征向量:"<<en

8、dl;for(i=0;i<N;i+)if(Ci=0)cout<<""<<i<<"對(duì)應(yīng)的特征向量"<<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)計(jì)算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)計(jì)算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法求實(shí)矩陣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);/巧妙的使用老師上課要求上機(jī)調(diào)試練習(xí)的列主元高斯消去法求解實(shí)特征值對(duì)應(yīng)的特征函數(shù)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;