北航数值分析大作业二

《北航数值分析大作业二》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、北京航咳航夭大學BEIHANGUNIVERSITY.数值分析大作业二院（系）名称机械工程及自动化学院学号ZY1107201姓名丁洪利2011年11月一、算法设计方案1、首先利用拟上三角化程序将矩阵A化为上三角矩阵在程序中用函数uppertriangularmatrix表示，其算法步骤如下记A(1)=A,并记A(r)的第r列到第n列元素为a,〉(i=1,2,...,〃;j=r,r+1,...,/?)。对于r=1,2,.../-2执行如下算法：⑴若曙(心厂+2,厂+3,…加全为零，则令A(r+1)=Ar,转(5)；否则传(

2、2)。(2)计算Vi=r^leJdr(ifa乩SO)Cr~[-drdf如/〉0)h,=c；-c"；)“(3)令“严(0,・・・,0卫乩-c,监肿斛r⑷计算pr=A(r)Tur/hrqr=A{r)ur/hrtr=pTrur/hrA(r+[)=A{r)-G)ruJ-up7(5)继续2、利用双步位移的QR分解求矩阵的全部特征值在程序中用函数twostepQR表示，其中需要调用函数MQRdecompositon(对矩阵Mk进彳亍QR分解)，算法如下Step1:调用uppertriangularmatrix函数得到的矩阵A(n~

3、[},给定精度水平£=严和迭代的最大次数L=10000；Step2:记£=肝1=「疋)],令k=l,m=10;L)」nxnStcp3:判断是否a^m_}<£,若是,贝I」4鳥一]为特征值,m=m-1,gotostcp4,否则gotostep5;Step4:女口果m=l,贝为特征值，gotostep11；若m=0,gotostep11；如果m>l,gotostep3;Step5(step6):如果m=2,则得到A的两个特征值必和s?(耳和耳为右下角二阶子阵对应的特征方程/-(%+_]+amm)A+detD=0的两个根)，

4、gotostepll;否则gotostep7;Step7:如果疋'、叶2§£，则得到A的两个特征值两个特征值5,和S'(5,和S'定义如上)，令m=m・2,gotostep4;否则gotostep8;Step&如果k>10000,输出“计算终止,未得到全部的特征值”,gotostepll;否则gotostep9;Step9：CAk(1

5、+l,gotostep3;StcplkA的所有特征值已经计算完毕，停止计算。3、利用QR分解算法來计算对双步位移中的进行QR分解，即函数MQRdecompositon的算法如卜•：记B{=Mk=,B=b-p,C]=Ako对r=1,2,...,/?-2执行如下算-J」mxmLJ」mxm法(1)若b：；)(i=r+l,r+2…沖)全为零，则令B(r+I)=Br9C(r+t)=Cr,转(5)；否则传⑵；(2)计算c=k(ifb^r<0)Cr~-dr{if/忧〉0)hr=c；~crb^⑶令叭=(0,・・・,0,砒-皿翊)T⑷

6、计算Vr=B(r)Tur/hrB(r+])=B{r)-ut.vJp,=C^Tur/hrqr=C^ur!hrtr=pTrur!hrG)r=qr-trurC(r+r}=C(r)-a)f.u,T-urPr(5)继续。经过此算法，就得到A,+1=Cmo4^利用QR分解算法对经过拟上三角化uppertriangularmatrix函数得到的矩阵川心)进行QR分解。算法设计如下iflA(1)=A，并记A(r)=[a^l,0=7。对于厂=1,2,...,77-1执行(1)若砧)(心厂+1,厂+2…加全为零，则令转(5)；否则转⑵;(

7、2)计算<=JE(4r))2Vi=rc饥SO)Cr~[-dr(ifd乩〉0)(3)令冷=(0,・・・,0,砒-c丿乩，础产⑷计算©=QrUr-©叮仇"P>.=Arur/hr4+i=~urpj经过上述运算之后，就得到0=0和上三角矩阵然后利用R和0相乘可计算R0。5、解出矩阵的特征值人之后利用选主元的Gauss消去法解线性方程组(A-刃)X=0所得的解向量兀.即为对应于特征值A,的特征向量。由于(4-V)为非满秩矩阵，故这样变换得到的矩阵％=0,所以令x10=l,这样可以得到特征向量。解线性方程组利用函数gauss进行

8、。其算法步骤如下消元过程⑴选行号L，使

9、哦

10、弋瓷略。(2)交换砧)与47(J=k,k+1,…/)0(3)对于(i=R+1,R+2,...,〃)计算卜常)%]硝+i)=硝)一加=R+1,R+2,.../)回代过程xio=°耳=一£a^Xj/=〃-1/-1,・・・,1)7=^+1经过上述过程解线性方程组(A-V)^=O所得的解向量兀即