《北航数值分析报告大作业第二题.pdf》由会员分享，可在线阅读，更多相关《北航数值分析报告大作业第二题.pdf（25页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、实用文档数值分析第二次大作业数值分析第二次大作业史立峰史立峰SY1505327SY1505327实用文档一、一、方案方案（1）利用循环结构将换的矩阵 A；（2）然后，对矩阵 A 利用 Householder 矩阵进行相似变换，把 A 化为上三角矩阵 A对 A 拟上三角化，得到拟上三角矩阵A(n-1)(n-1)sin(0.5i0.2 j)(i j)aij1.5cos(i1.2 j)(i j)(i,j=1,2,10)进行赋值，得到需要变。，具体算法如下：(r)i 1,2,n;j r,r 1,naij记 A(1)=A，并记 A(r)的第 r 列至第 n 列的元素为对于r 1,2,n 2执行1.若(r
2、)i r 2,r 3,nair。全为零，则令 A(r+1)=A(r),转 5；否则转 2。2.计算drir1an(r)2ir)(r)cr sgnar(r1,rdr若ar1,r 0,则取cr dr)hr cr2crar(r1,r3.令)(r)(r)ur 0,0,ar(r1,rcr,ar2,r,anrTRn。4.计算pr A(r)Tur/hrqr A(r)ur/hrTtr prur/hrr qrtrurTTA(r1)A(r)rururpr5.继续。（3）使用带双步位移的 QR 方法计算矩阵 A算法如下：1.给定精度水平 0和迭代最大次数L。(1)2.记A(1)A(n1)aijnn，令k 1,m n
3、。(n-1)的全部特征值，也是 A 的全部特征值，具体实用文档(k)3.如果am，转4；否则,m1，则得到A的一个特征值amm，置m:m1（降阶）(k)转 5。4.如果m 1，则得到的一个特征值a11，转 11；如果m 1，则转 3。5.求 2 阶子阵(k)(k)am1,m1am1,mDk(k)(k)aammm,m1(k)的两个特征值s1和s2，即计算二次方程(k)(k)2(am1.m1 amm)det Dk 0的两个根s1和s2。6.如果m 2，则得到A的两个特征值s1和s2，转 11；否则转 7。(k)7.如果am，转 4；1,m2，则得到A的两个特征值s1和s2，置m:m2（降阶）否则转
4、 88.如果k L，则计算终止，未得到A的全部特征值；否则转9。(k)9.记A(k)aijmm(1 i,j m)，计算k)(k)sa(am1,m1mmk)(k)(k)(k)ta(m1,m1ammam,m1am1,mMkAk2sAktI(I是m阶单位矩阵)MkQkRk(对Mk作QR分解)TAk1QkAkQk10.置k:k 1，转 3。11.A的全部特征值已计算完毕，停止计算。其中，Mk的QR分解与Ak1的计算用下列算法实现：(1)(r)记B1 Mkbijmm,Brbijmm,C1 Ak。对于r 1,2,m 1执行1.若bir(r)i r 1,r 2,m全为零，则令Br1 Br,Cr1 Cr,转
5、5；否则转 2。2.计算实用文档drbirm(r)2ir(r)dr若ar(r)1,r 0,则取cr drcr sgnbrr(r)hr cr2 crbrr3.令ur 0,0,brrcr,br1,r,bmr4.计算(r)(r)(r)TRm。vr BrTur/hrTBr1 Brurvrpr CrTur/hrqr Crur/hrTtrprur/hrr qrtrurTTCr1 Crrururpr5.继续。此算法执行完后，就得到Ak1 Cm。（4）计算 Q,R，一边求 R*Q 矩阵。(r)记A1 A，并记Ar aijnn,令Q1I（n阶单位阵）2,n 1执行对于r 1,1.若(r)i r 2,r 3,na
6、ir全为零，则令Qr1QrAr1Ar转 5；否则转 2.2.计算drbirm(r)2ir(r)dr若ar(r)1,r 0,则取cr drcr sgnbrr(r)hr cr2 crbrr实用文档3.令ur 0,0,brrcr,br1,r,bmr4.计算(r)(r)(r)TRm。rQrurQr1T urrQrThrprArurhrTAr1Arurpr5.继续当此算法执行完毕后，就得到正交矩阵QQn和上三角矩阵R6.然后计算出RQ矩阵（5）用列主元素 Gauss 消去法计算矩阵A对应于实特征值的特征向量，具体算法如下：记A(1)AI(为A的实特征值)1.消元过程对于k 1,2,n 1执行(1)选行号
7、ik，使ai(kkk)maxakin(k)ik。(k)(2)交换akj与ai(kkj)(j k,k 1,n)所含的数值。(3)对于i k 1,k 2,n计算(k)(k)mik aik/akk(k1)(k)(k)aij aij mikaki(j k 1,k 2,n)2.回代过程xn1xkjk1an(k)kj(k)xjakkk n 1,n 2,1T最终得到的向量(x1,x2,xn1)的即为A对应于实特征值的特征向量。实用文档二、源程序二、源程序#include#include#include#define n 10#define E 1.0e-12void Householder(double a
8、nn);/拟上三角化函数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);/定义列主元高斯消去法函数void main()in
9、t i,j,k;double ann,qrnn,qnn,rnn,Ln2,xn;for(i=0;in;i+)/录入要进行变换的矩 a，并对 q 初始赋值。/coutendl;Householder(a);/调用拟上三角化函数for(j=0;jn;j+)if(i=j)aij=1.5*cos(i+1+1.2*(j+1);elseaij=sin(0.5*(i+1)+0.2*(j+1);实用文档coutendlendl对矩阵 A 拟上三角化前七列的结果:endl;for(i=0;in;i+);cout对矩阵 A 拟上三角化后三列的结果:endl;for(i=0;in;i+);/k+;/if(k%3=0)
10、/判断每行是否到达三个元素/k=0;/为了显示方便，每行显示三个元素，使用k 来实现/cout矩阵 A 的第i+1行元素为：endl;for(j=7;jn;j+)if(fabs(aij)E)aij=0;/k+;coutendl;/coutendl;/k=0;/为了显示方便，每行显示三个元素，使用k 来实现/cout矩阵 A 的第i+1行元素为：endl;for(j=0;j7;j+)if(fabs(aij)E)aij=0;coutsetprecision(12)setiosflags(ios:scientific)aij/if(k%3=0)/判断每行是否到达三个元素coutsetprecisio
11、n(12)setiosflags(ios:scientific)aij实用文档/coutendl;coutendl;coutendl;QRfenjie(a,q,r);QR(a,L);coutendl对矩阵 A 进行 QR 分解后所得矩阵前七列的结果:endl;for(i=0;in;i+);cout对矩阵 A 进行 QR 分解后所得矩阵三列的结果:endl;for(i=0;in;i+);coutendl;for(j=7;jn;j+)if(fabs(aij)E)aij=0;coutendl;for(j=0;j7;j+)if(fabs(aij)E)aij=0;coutsetprecision(12)
12、setiosflags(ios:scientific)aijcoutsetprecision(12)setiosflags(ios:scientific)aij实用文档/*cout对矩阵 A 进行 QR 分解后所得矩阵:endl;for(i=0;in;i+);*/for(i=0;in;i+)coutendlR*Q 相乘的前七列:endl;for(i=0;in;i+)/k=0;/coutR*Q 的第i+1行元素为：endl;for(j=0;j7;j+)for(j=0;jn;j+)for(k=0,qrij=0;kn;k+)qrij=qrij+rik*qkj;coutendl;k+;if(k%3=0
13、)coutendl;k=0;cout矩阵 A 的第i+1行元素为：endl;for(j=0;jn;j+)if(fabs(aij)E)aij=0;coutsetprecision(12)setiosflags(ios:scientific)aij实用文档;if(fabs(qrij)E)aij=0;coutsetprecision(12)setiosflags(ios:scientific)qrij/k+;/if(k%3=0)/coutendl;coutendl;coutendlR*Q 相乘的后三列:endl;for(i=0;in;i+);coutendl矩阵 A 的全部特征值为endl;for(
14、i=0;in;i+)if(Li1=0)couti+1=Li0endl;coutendl;for(j=7;jn;j+)if(fabs(qrij)E)aij=0;coutsetprecision(12)setiosflags(ios:scientific)qrijelsecouti+1=Li0+Li1iendl;实用文档coutendl实特征值对应的特征向量分别是：endl;for(k=0;kn;k+)if(Lk1=0)/判断特征值是否为实特征值else continue;for(i=0;in;i+)/重新录入矩阵 AGauss(a,x);coutk+1=Lk0for(j=0;jn;j+)cout
15、xjendl;对应的特征向量是:endl;for(j=0;jn;j+)if(i=j)aij=1.5*cos(i+1+1.2*(j+1)-Lk0;elseaij=sin(0.5*(i+1)+0.2*(j+1);void Householder(double ann)/拟上三角化函数 int i,j,k;int m=0;实用文档double d,c,h,t;double un,pn,qn,wn;for(i=0;in-2;i+)for(j=i+2;jn;j+)if(aji=E)m=m+1;if(m=(n-2-i)continue;for(j=i+1,d=0;jn;j+)d=d+aji*aji;d=s
16、qrt(d);c=-1*sgn(ai+1i)*d;h=c*c-c*ai+1i;for(j=i+2;jn;j+)for(j=0;ji+2;j+)ui+1=ai+1i-c;for(j=0;jn;j+)for(k=i+1,pj=0;kn;k+)uj=0;uj=aji;实用文档pj=pj/h;pj=akj*uk+pj;for(j=0;jn;j+)for(j=0,t=0;jn;j+)t=t/h;for(j=0;jn;j+)for(j=0;jn;j+)for(k=0;kn;k+)ajk=ajk-wj*uk-uj*pk;wj=qj-t*uj;t=t+pj*uj;for(k=i+1,qj=0;k0)retur
17、n 1;else if(a0)return-1;elsereturn 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;in;i+)for(i=0;in;i+)for(i=0;in-1;i+)for(j=0;jn;j+)Rij=aij;for(j=0;jn;j+)if(i=j)Qij=1;else Qij=0;for(j=i,d=0;jn;j+)d=d+Rji*Rji;实用文档d=sqrt(d);c=-1*sgn(R
18、ii)*d;h=c*c-c*Rii;for(j=i+1;jn;j+)uj=Rji;for(j=0;ji;j+)uj=0;ui=Rii-c;for(j=0;jn;j+)for(j=0;jn;j+)for(j=0;jn;j+)for(j=0;jn;j+)for(k=0;kn;k+)Rjk=Rjk-uj*pk;for(k=i,pj=0;kn;k+)pj=Rkj*uk+pj;for(k=0;kn;k+)Qjk=Qjk-wj*uk/h;for(k=0,wj=0;kn;k+)wj=Qjk*uk+wj;pj=pj/h;/以上是 矩阵的 QR 分解void MxM(double Mnn,double Ann,
19、double Bnn,int m)/矩阵相乘double Cnn,Dnn;实用文档int i,j,g;for(i=0;im;i+)for(j=0;jm;j+)Cij=Aij;Dij=Bij;Mij=0;for(i=0;im;i+)for(j=0;jm;j+)for(g=0;gm;g+)Mij=Mij+Cig*Dgj;void QR(double ann,double Ln2)/带双步位移的 QR 方法double Mnn,s,t,s12,s22,un,vn,pn,qn,wn,sum,h,c,d;int i,j,m,r,k;for(i=0;in;i+)m=n;for(k=1;k=100;k+)i
20、f(fabs(am-1m-2)=E)Lm-10=am-1m-1;for(j=0;j2;j+)Lij=0;m=m-1;if(m=1)L00=a00;实用文档 break;elsesolve(a,s1,s2,m);/调用解方程函数，将求解结果存到s1 和 s2 中if(m=2)elseif(fabs(am-2m-3)=E)Lm-20=s10;Lm-21=s11;Lm-10=s20;Lm-11=s21;m=m-2;if(m=1)L00=a00;L00=s10;L01=s11;L10=s20;L11=s21;break;else if(m=0)break;elsecontinue;实用文档break;
21、else if(m=0)break;elsecontinue;elseif(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;im;i+)for(j=0;jm;j+)Mij=Mij-s*aij;for(i=0;im;i+)Mii=Mii+t;for(r=0;r(m-1);r+)for(i=(r+1);im;i+)if(Mir=0)continue;实用文档elsesum=0;for(i=r;im;i+)s
22、um+=pow(Mir,2);d=sqrt(sum);if(Mrr!=0)c=-sgn(Mrr)*d;else c=d;h=c*(c-Mrr);for(i=0;im;i+)for(i=0;im;i+)for(i=0;im;i+)for(j=0;jm;j+)Mij=Mij-ui*vj;sum=0;for(j=r;jm;j+)sum+=Mji*uj;vi=sum/h;if(ir)ui=Mir;实用文档for(i=0;im;i+)sum=0;for(j=r;jm;j+)sum+=aji*uj;pi=sum/h;sum=0;for(j=r;jm;j+)sum+=aij*uj;qi=sum/h;sum=
23、0;for(j=r;jm;j+)sum+=pj*uj;t=sum/h;for(j=0;jm;j+)wj=qj-t*uj;for(i=0;im;i+)for(j=0;jm;j+)aij=aij-wi*uj-ui*pj;实用文档for(i=0;in;i+)for(j=0;jn;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-2m-2;det=s*s-4*t;if(det0)de
24、t=-det;det=sqrt(det);s10=s/2;s11=det/2;s20=s/2;elsedet=sqrt(det);s21=-det/2;s10=(s+det)/2;s11=0;s20=(s-det)/2;s21=0;实用文档void Gauss(double ann,double xn)/定义列主元高斯消去法函数int i,j,k,t;double sum,mnn,max,p;for(k=0;k(n-1);k+)max=fabs(akk);t=k;for(i=(k+1);in;i+)if(t!=k)for(j=k;jn;j+)for(i=(k+1);in;i+)mik=aik/
25、akk;p=akj;akj=atj;atj=p;if(maxfabs(aik)max=fabs(aik);t=i;实用文档for(j=(k+1);j=0;k-)sum=0;for(i=0;in;i+)sum=sqrt(sum);for(i=0;in;i+)xi=xi/sum;sum+=xi*xi;sum=0;for(j=(k+1);jn;j+)sum+=akj*xj;xk=(-sum)/akk;三、处理结果三、处理结果实用文档实用文档四、结果讨论四、结果讨论（1）此次编程涉及的函数众多，需要特别注意各函数之间的协调；（2）采用全局变量 n 节省了输入次数，简化编程；（3）由于需要求出RQ矩阵，还需要另外的编辑函数来求解除Q和R，然后再求出RQ。这就造成编程负担。（4）同样的程序，使用程序语言不同，得到的结果不同。最简单的就是把 C+程序语言，改换成 C 语言，结果就不相同。而且C 语言的到的结果与课本提供的参考结果相同，C+则有些偏差。（和另一位同学对比，改换运行的到的）