北航数值分析大作业第二题

1、数值分析第二次大作业史立峰SY1505327一、 方案（1）利用循环结构将(i,j=1,2,10)进行赋值，得到需要变换的矩阵A；（2）然后，对矩阵A利用Householder矩阵进行相似变换，把A化为上三角矩阵A(n-1)。对A拟上三角化，得到拟上三角矩阵A(n-1)，具体算法如下：记A(1)=A，并记A(r)的第r列至第n列的元素为。对于执行1. 若全为零，则令A(r+1) =A(r),转5；否则转2。2. 计算3. 令。4. 计算 5. 继续。（3）使用带双步位移的QR方法计算矩阵A(n-1)的全部特征值，也是A的全部特征值，具体算法如下：1. 给定精度水平和迭代最大次数。2. 记，令。

2、3. 如果，则得到的一个特征值，置（降阶），转4；否则转5。4. 如果，则得到的一个特征值，转11；如果，则转3。5. 求2阶子阵的两个特征值和，即计算二次方程的两个根和。6. 如果，则得到的两个特征值和，转11；否则转7。7. 如果，则得到的两个特征值和，置（降阶），转4；否则转88. 如果，则计算终止，未得到的全部特征值；否则转9。9. 记，计算10. 置，转3。11. 的全部特征值已计算完毕，停止计算。其中，的分解与的计算用下列算法实现：记。对于执行1. 若全为零，则令,转5；否则转2。2. 计算3. 令。4. 计算5. 继续。此算法执行完后，就得到。（4）计算Q,R，一边求R*Q矩阵。

3、记对于执行1.若全为零，则令转5；否则转 2.2.计算3.令 。4.计算5.继续当此算法执行完毕后，就得到正交矩阵和上三角矩阵6.然后计算出矩阵（5）用列主元素Gauss消去法计算矩阵对应于实特征值的特征向量，具体算法如下：记1. 消元过程对于执行(1) 选行号，使。(2) 交换与所含的数值。(3) 对于计算2. 回代过程最终得到的向量的即为对应于实特征值的特征向量。二、源程序#include<iostream.h>#include<math.h>#include<iomanip.h>#define n 10#define E 1.0e-12void Hou

4、seholder(double ann);/拟上三角化函数double sgn(double a);/符号函数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);/解方程函数void Gauss(double ann, double xn);/定义列主元高斯消去

5、法函数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);/调用拟上三角化函数cout<<endl<<endl<<"对矩阵A拟上三角化前七列的结果:"<<endl;fo

6、r(i=0;i<n;i+)/k=0;/为了显示方便，每行显示三个元素，使用k来实现/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、达三个元素cout<<endl;cout<<"对矩阵A拟上三角化后三列的结果:"<<endl;for(i=0;i<n;i+)/k=0;/为了显示方便，每行显示三个元素，使用k来实现/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)/判断每行是否到达三个元素/cout<<endl;cout<<endl;cout<<endl;QRfenjie(a,q,r);QR(a,L);cout<<endl<<"对矩阵A进行QR分解后所得矩阵前七列的结果:"<<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分解后所得矩阵三列的结果:"<<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<<"实特征值对应的特征向量分别是："<<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<<""<<"对应的特征向量是:"<<endl;for(j=0;j<n;j+)cout<<xj<<endl;else continue;void Householder(double ann)/拟上三角化函数 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; /以上是拟上三角化的结果 double sgn(double a)/符号函数if(a>0)return 1;else if(a<0)return -1;elser

18、eturn 0;/以上是符号函数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);/调用解方程函数，将求解结果存到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); /调用矩阵相乘函数构造矩阵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)/解方程函数double s,t,det;t=am-1m-1*am-2m-2-am-1m-2*am-2m-1;s=am-1m-1+am-2