1、數(shù) 值 分 析 計算實習作業(yè)二學 院： 17系專 業(yè)： 精密儀器及機械姓 名： 張大軍學 號： 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經(jīng)過擬上三角化所得的矩陣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);/巧妙的使用老師上課要求上機調試練習的列主元高斯消去法求解實特征值對應的特征函數(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;