北航數(shù)值分析大作業(yè)第二題

1、數(shù)值分析第二次大作業(yè)史立峰SY1505327一、 方案（1）利用循環(huán)結(jié)構(gòu)將(i,j=1,2,10)進行賦值，得到需要變換的矩陣A；（2）然后，對矩陣A利用Householder矩陣進行相似變換，把A化為上三角矩陣A(n-1)。對A擬上三角化，得到擬上三角矩陣A(n-1)，具體算法如下：記A(1)=A，并記A(r)的第r列至第n列的元素為。對于執(zhí)行1. 若全為零，則令A(yù)(r+1) =A(r),轉(zhuǎn)5；否則轉(zhuǎn)2。2. 計算3. 令。4. 計算 5. 繼續(xù)。（3）使用帶雙步位移的QR方法計算矩陣A(n-1)的全部特征值，也是A的全部特征值，具體算法如下：1. 給定精度水平和迭代最大次數(shù)。2. 記，令。

2、3. 如果，則得到的一個特征值，置（降階），轉(zhuǎn)4；否則轉(zhuǎn)5。4. 如果，則得到的一個特征值，轉(zhuǎn)11；如果，則轉(zhuǎn)3。5. 求2階子陣的兩個特征值和，即計算二次方程的兩個根和。6. 如果，則得到的兩個特征值和，轉(zhuǎn)11；否則轉(zhuǎn)7。7. 如果，則得到的兩個特征值和，置（降階），轉(zhuǎn)4；否則轉(zhuǎn)88. 如果，則計算終止，未得到的全部特征值；否則轉(zhuǎn)9。9. 記，計算10. 置，轉(zhuǎn)3。11. 的全部特征值已計算完畢，停止計算。其中，的分解與的計算用下列算法實現(xiàn)：記。對于執(zhí)行1. 若全為零，則令,轉(zhuǎn)5；否則轉(zhuǎn)2。2. 計算3. 令。4. 計算5. 繼續(xù)。此算法執(zhí)行完后，就得到。（4）計算Q,R，一邊求R*Q矩陣。

3、記對于執(zhí)行1.若全為零，則令轉(zhuǎn)5；否則轉(zhuǎn) 2.2.計算3.令 。4.計算5.繼續(xù)當(dāng)此算法執(zhí)行完畢后，就得到正交矩陣和上三角矩陣6.然后計算出矩陣（5）用列主元素Gauss消去法計算矩陣對應(yīng)于實特征值的特征向量，具體算法如下：記1. 消元過程對于執(zhí)行(1) 選行號，使。(2) 交換與所含的數(shù)值。(3) 對于計算2. 回代過程最終得到的向量的即為對應(yīng)于實特征值的特征向量。二、源程序#include<iostream.h>#include<math.h>#include<iomanip.h>#define n 10#define E 1.0e-12void Hou

4、seholder(double ann);/擬上三角化函數(shù)double sgn(double a);/符號函數(shù)void QRfenjie(double ann,double Qnn,double Rnn);void QR(double ann,double Ln2);/帶雙步位移的QR分解void MxM(double Mnn,double Ann,double Bnn,int m);/矩陣相乘void solve(double ann,double s12,double s22,int m);/解方程函數(shù)void Gauss(double ann, double xn);/定義列主元高斯消去

5、法函數(shù)void main()int i,j,k;double ann,qrnn,qnn,rnn,Ln2,xn;for(i=0;i<n;i+)/錄入要進行變換的矩a，并對q初始賦值。for(j=0;j<n;j+)if(i=j)aij=1.5*cos(i+1+1.2*(j+1);elseaij=sin(0.5*(i+1)+0.2*(j+1); / cout<<endl;Householder(a);/調(diào)用擬上三角化函數(shù)cout<<endl<<endl<<"對矩陣A擬上三角化前七列的結(jié)果:"<<endl;fo

6、r(i=0;i<n;i+)/k=0;/為了顯示方便，每行顯示三個元素，使用k來實現(xiàn)/cout<<"矩陣A的第"<<i+1<<"行元素為："<<endl;for(j=0;j<7;j+)if (fabs(aij)<E)aij=0;cout<<setprecision(12)<<setiosflags(ios:scientific)<<aij<<""/k+;/cout<<endl;/if(k%3=0)/判斷每行是否到

7、達(dá)三個元素cout<<endl;cout<<"對矩陣A擬上三角化后三列的結(jié)果:"<<endl;for(i=0;i<n;i+)/k=0;/為了顯示方便，每行顯示三個元素，使用k來實現(xiàn)/cout<<"矩陣A的第"<<i+1<<"行元素為："<<endl;for(j=7;j<n;j+)if (fabs(aij)<E)aij=0;cout<<setprecision(12)<<setiosflags(ios:scien

8、tific)<<aij<<""/k+;/if(k%3=0)/判斷每行是否到達(dá)三個元素/cout<<endl;cout<<endl;cout<<endl;QRfenjie(a,q,r);QR(a,L);cout<<endl<<"對矩陣A進行QR分解后所得矩陣前七列的結(jié)果:"<<endl;for(i=0;i<n;i+)for(j=0;j<7;j+)if (fabs(aij)<E)aij=0;cout<<setprecision(12)

9、<<setiosflags(ios:scientific)<<aij<<""cout<<endl;cout<<"對矩陣A進行QR分解后所得矩陣三列的結(jié)果:"<<endl;for(i=0;i<n;i+)for(j=7;j<n;j+)if (fabs(aij)<E)aij=0;cout<<setprecision(12)<<setiosflags(ios:scientific)<<aij<<""cout

10、<<endl;/*cout<<"對矩陣A進行QR分解后所得矩陣:"<<endl;for(i=0;i<n;i+)k=0;cout<<"矩陣A的第"<<i+1<<"行元素為："<<endl;for(j=0;j<n;j+)if (fabs(aij)<E)aij=0;cout<<setprecision(12)<<setiosflags(ios:scientific)<<aij<<"&

11、quot;k+;if(k%3=0)cout<<endl;cout<<endl;*/for(i=0;i<n;i+)for(j=0;j<n;j+) for(k=0,qrij=0;k<n;k+) qrij=qrij+rik*qkj;cout<<endl<<"R*Q相乘的前七列:"<<endl;for(i=0;i<n;i+)/k=0;/cout<<"R*Q的第"<<i+1<<"行元素為："<<endl;for(j

12、=0;j<7;j+)if (fabs(qrij)<E)aij=0;cout<<setprecision(12)<<setiosflags(ios:scientific)<<qrij<<""/k+;/if(k%3=0)/cout<<endl;cout<<endl;cout<<endl<<"R*Q相乘的后三列:"<<endl;for(i=0;i<n;i+)for(j=7;j<n;j+)if (fabs(qrij)<E)ai

13、j=0;cout<<setprecision(12)<<setiosflags(ios:scientific)<<qrij<<""cout<<endl;cout<<endl<<"矩陣A的全部特征值為"<<endl;for(i=0;i<n;i+)if(Li1=0)cout<<""<<i+1<<"="<<Li0<<endl;elsecout<<&q

14、uot;"<<i+1<<"="<<Li0<<"+"<<Li1<<"i"<<endl;cout<<endl<<"實特征值對應(yīng)的特征向量分別是："<<endl;for(k=0;k<n;k+)if(Lk1=0)/判斷特征值是否為實特征值for(i=0;i<n;i+)/重新錄入矩陣Afor(j=0;j<n;j+)if(i=j)aij=1.5*cos(i+1+1.2*(j+1)-

15、Lk0;elseaij=sin(0.5*(i+1)+0.2*(j+1);Gauss(a,x);cout<<""<<k+1<<"="<<Lk0<<""<<"對應(yīng)的特征向量是:"<<endl;for(j=0;j<n;j+)cout<<xj<<endl;else continue;void Householder(double ann)/擬上三角化函數(shù) int i,j,k;int m=0;double d,c

16、,h,t;double un,pn,qn,wn;for(i=0;i<n-2;i+) for(j=i+2;j<n;j+) if(aji<=E) m=m+1; if(m=(n-2-i) continue; for(j=i+1,d=0;j<n;j+) d=d+aji*aji; d=sqrt(d);c=-1*sgn(ai+1i)*d;h=c*c-c*ai+1i;for(j=i+2;j<n;j+)uj=aji;for(j=0;j<i+2;j+)uj=0;ui+1=ai+1i-c;for(j=0;j<n;j+) for(k=i+1,pj=0;k<n;k+)

17、pj=akj*uk+pj;pj=pj/h;for(j=0;j<n;j+) for(k=i+1,qj=0;k<n;k+) qj=ajk*uk+qj;qj=qj/h;for(j=0,t=0;j<n;j+)t=t+pj*uj;t=t/h;for(j=0;j<n;j+) wj=qj-t*uj;for(j=0;j<n;j+)for(k=0;k<n;k+) ajk=ajk-wj*uk-uj*pk; /以上是擬上三角化的結(jié)果 double sgn(double a)/符號函數(shù)if(a>0)return 1;else if(a<0)return -1;elser

18、eturn 0;/以上是符號函數(shù)void QRfenjie(double ann,double Qnn,double Rnn)/矩陣的QR分解 int i,j,k;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,d=0;j<n;j+) d=d+Rji*Rji;d=sqrt(d);c=-1*sgn(Rii)*d

19、;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;k<n;k+) Rjk=Rjk-uj*pk;/以

20、上是 矩陣的QR分解 void MxM(double Mnn,double Ann,double Bnn,int m)/矩陣相乘double Cnn,Dnn;int i,j,g;for(i=0;i<m;i+)for(j=0;j<m;j+)Cij=Aij;Dij=Bij;Mij=0;for(i=0;i<m;i+)for(j=0;j<m;j+)for(g=0;g<m;g+)Mij=Mij+Cig*Dgj;void QR(double ann,double Ln2)/帶雙步位移的QR方法double Mnn,s,t,s12,s22,un,vn,pn,qn,wn,sum,

21、h,c,d;int i,j,m,r,k;for(i=0;i<n;i+)for(j=0;j<2;j+)Lij=0;m=n;for(k=1;k<=100;k+)if(fabs(am-1m-2)<=E)Lm-10=am-1m-1; m=m-1; if(m=1)L00=a00; break;else if(m=0)break;else continue;elsesolve(a,s1,s2,m);/調(diào)用解方程函數(shù)，將求解結(jié)果存到s1和s2中if(m=2)L00=s10;L01=s11;L10=s20;L11=s21;break;elseif(fabs(am-2m-3)<=E

22、)Lm-20=s10;Lm-21=s11;Lm-10=s20;Lm-11=s21;m=m-2;if(m=1)L00=a00;break;else if(m=0)break;else continue;else if(k=100)cout<<"未能得到全部特征值"<<endl;break;elset=am-1m-1*am-2m-2-am-1m-2*am-2m-1;s=am-1m-1+am-2m-2;MxM(M,a,a,m); /調(diào)用矩陣相乘函數(shù)構(gòu)造矩陣Mfor(i=0;i<m;i+)for(j=0;j<m;j+)Mij=Mij-s*aij;

23、for(i=0;i<m;i+)Mii=Mii+t;for(r=0;r<(m-1);r+)for(i=(r+1);i<m;i+)if(Mir=0)continue;elsesum=0;for(i=r;i<m;i+)sum+=pow(Mir,2);d=sqrt(sum);if(Mrr!=0)c=-sgn(Mrr)*d;else c=d;h=c*(c-Mrr);for(i=0;i<m;i+)if(i<r)ui=0;if(i=r)ui=Mrr-c;if(i>r)ui=Mir;for(i=0;i<m;i+)sum=0;for(j=r;j<m;j+)s

24、um+=Mji*uj;vi=sum/h;for(i=0;i<m;i+)for(j=0;j<m;j+)Mij=Mij-ui*vj;for(i=0;i<m;i+)sum=0;for(j=r;j<m;j+)sum+=aji*uj;pi=sum/h;sum=0;for(j=r;j<m;j+)sum+=aij*uj;qi=sum/h;sum=0;for(j=r;j<m;j+)sum+=pj*uj;t=sum/h;for(j=0;j<m;j+)wj=qj-t*uj;for(i=0;i<m;i+)for(j=0;j<m;j+)aij=aij-wi*uj-ui*pj;for(i=0;i<n;i+)for(j=0;j<n;j+)if(aij<=E)aij=0;void solve(double ann,double s12,double s22,int m)/解方程函數(shù)double s,t,det;t=am-1m-1*am-2m-2-am-1m-2*am-2m-1;s=am-1m-1+am-2