北航数值分析解析大作业二（纯原创）.docx

《北航数值分析解析大作业二（纯原创）.docx》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、目录题目：2算法设计思路和方案2关于第一步矩阵A的拟上三角化2关于对矩阵An-1进行带双步位移的QR分解迭代3关于求从属于矩阵A的实特征值λi的特征向量3计算结果4发现的现象与问题:7探究带双步位移的QR分解比一般QR分解节省的计算量7探究拟上三角化对QR分解迭代收敛速度的影响7关于直接单步QR分解计算发现的问题9源程序10数值分析上机实习作业二王强SY1413315题目：试求矩阵A=aij10×10的全部特征值，并对其中的每个实特征值求相应的特征向量，已知：aij=sin0.5i+0.2ji≠j&1.52cosi+1.2ji=j(i,j=1,2,⋯,10)说明：1.在所用的算

2、法中，凡是要给出精度水平的ε，都取ε=10-12。2.打印以下内容：(1)采用带双步位移的QR分解法，说明算法设计方案和思路。(2)全部源程序。(3)矩阵A经过拟上三角化后的矩阵An-1。(4)对矩阵An-1实行QR方法迭代结束后所得的矩阵；(5)矩阵A的全部实特征值λi=Ri,Iii=1,2,,⋯,10,其中Ri=Reλi,Ii=Imλi。若λi是实数，则令Ii=0;(6)A相应的实特征值的特征向量。(7)发现的现象与遇到的问题。3.采用e型数输出实型数，并且至少显示12位有效数字。算法设计思路和方案该问题的求解起始主要分为三个步骤，第一步是对A拟上三角化得矩阵An-1；第二

3、步是对矩阵An-1进行带双步位移的QR分解迭代，并求出全部特征值；第三部是求出所有实特征值对应的特征向量。关于第一步矩阵A的拟上三角化对于此步书上已经有明确的思路和实现方法了，只是在编程的过程中注意对下式在计算的时候应注意避免矩阵乘矩阵Ar+1=HrArHr如果设airr,(i=r+1,r+2,⋯,n不全为零，则有：sr=0,⋯,0,ar+1,rr,⋯,anrrTcr=-sgnar+1,rrsr2，若ar+1,rr=0，择取cr=sr2。ur=sr-cerνr=urur2Hr=I-2νrνrT为了避免矩阵相乘采用如下步骤计算Ar+1=HrArHr=I-2νrνrTA'r=A'r-

4、2νrνrTA'r其中A'r=Ar-2ArνrνrT1.计算中间向量xr=Arνr；2.从i=1,2,⋯,n,j=r+1,⋯,n计算aij'r=aijr-2νrixrj,A'r=aij'r；3.计算中间向量xr=νrTA'r；4.从j=1,2,⋯,n,i=r+1,⋯,n计算aijr+1=aij'r-2νrixrj,Ar+1=aijr+1；关于对矩阵AN-1进行带双步位移的QR分解迭代尽管书上给出了算法实现的11个步骤，然而在思路很混乱，在判断某步QR迭代后能否对得到的矩阵Ak+1求部分特征值以达到降阶方面，按照书上的判断方法，几乎不可避免的得使用goto语句，这使得程序流程变得混

5、乱。经过仔细分析后，可以发现能否求解Ak+1的部分特征值然后实现矩阵Ak+1降阶的关键在于判断一下拟上三角阵a11a12⋯⋯a1ma21a22a23⋯a2m0⋮⋮am-1m-2am-1m-1am-1m0⋯0amm-1amm是否存在下列形式的子块Bk+1=am-1m-1am-1m0ammCk+1=am-2m-2am-2m-1am-2m0am-1m-1am-1m0amm-1amm而其本质就是1.令k=1,m=10,A1=An-1=aijm×m以及最大迭代步数L；2.若m≤0，则结束计算，已求出A的全部特征值，判断amm-1k<10-12或am-1m-2k<10-12或m≤2是否成立，

6、成立则转3，否则转4；3.若amm-1k<10-12，则得一个特征值amm，m=m-1，降阶；若am-1m-2k<10-12，则计算矩阵：am-1m-1kam-1mkamm-1kammk的特征值得矩阵A的两个特征值，m=m-2，降阶，转2.；4.若k≤L，成立则令s=am-1m-1k+ammkt=ammkam-1m-1k-amm-1kam-1mkMk=Ak2-sAk+tIMk=QkRkAk+1=QkTAkQkk=k+1，转2，否则结束计算，为计算出矩阵A的全部特征值；关于求从属于矩阵A的实特征值ΛI的特征向量在求出A的全部特征值后，求从属于矩阵A的实特征值λi的特征向量的时候，

7、一个思路是对矩阵做平移：B=A-(λi-δ)I后对矩阵B使用反幂法，其中δ为适当选择的一个小量，但这样做计算量大，没有充分利用λi已知了这个条件，λi已知时，求矩阵A的从属于λi的特征向量等价于求以下线性方程组的一个非零解：Bx=0B=A-λiI为求该方程的一个非零解，可以先对矩阵B进行QR分解，然后在方程组两边同乘Q，可得同解方程组：Rx=0其中R=r11r12⋯r1100r22⋯r210⋮⋱⋱⋮0⋯0r1010由于在本例中λi是矩阵A的特征值，故矩阵B的秩r=10-p，其中p