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

目录 TOC \o 1-3 \h \z \u HYPERLINK \l _Toc404616009 题目： PAGEREF _Toc404616009 \h 2 HYPERLINK \l _Toc404616010 算法设计思路和方案 PAGEREF _Toc404616010 \h 2 HYPERLINK \l _Toc404616011 关于第一步矩阵A的拟上三角化 PAGEREF _Toc404616011 \h 2 HYPERLINK \l _Toc404616012 关于对矩阵An-1进行带双步位移的QR分解迭代 PAGEREF _Toc404616012 \h 3 HYPERLINK \l _Toc404616013 关于求从属于矩阵A的实特征值λi的特征向量 PAGEREF _Toc404616013 \h 3 HYPERLINK \l _Toc404616014 计算结果 PAGEREF _Toc404616014 \h 4 HYPERLINK \l _Toc404616015 发现的现象与问题: PAGEREF _Toc404616015 \h 7 HYPERLINK \l _Toc404616016 探究带双步位移的QR分解比一般QR分解节省的计算量 PAGEREF _Toc404616016 \h 7 HYPERLINK \l _Toc404616017 探究拟上三角化对QR分解迭代收敛速度的影响 PAGEREF _Toc404616017 \h 7 HYPERLINK \l _Toc404616018 关于直接单步QR分解计算发现的问题 PAGEREF _Toc404616018 \h 9 HYPERLINK \l _Toc404616019 源程序 PAGEREF _Toc404616019 \h 10  数值分析上机实习作业二 王强 SY1413315 题目： 试求矩阵A=a aij 说明： 在所用的算法中，凡是要给出精度水平的ε，都取 ε=10 打印以下内容： 采用带双步位移的QR分解法，说明算法设计方案和思路。 全部源程序。 矩阵A经过拟上三角化后的矩阵An-1 对矩阵An-1实行QR 矩阵A的全部实特征值λ A相应的实特征值的特征向量。 发现的现象与遇到的问题。 采用e型数输出实型数，并且至少显示12位有效数字。 算法设计思路和方案 该问题的求解起始主要分为三个步骤，第一步是对A拟上三角化得矩阵An-1；第二步是对矩阵An-1进行带双步位移的 关于第一步矩阵A的拟上三角化 对于此步书上已经有明确的思路和实现方法了，只是在编程的过程中注意对下式在计算的时候应注意避免矩阵乘矩阵Ar+1=HrArHr如果 Ar+1=HrArHr=I-2νrνrTAr=Ar-2ν 关于对矩阵An-1进行带双步位移的QR分解迭代 尽管书上给出了算法实现的11个步骤，然而在思路很混乱，在判断某步QR迭代后能否对得到的矩阵Ak+1求部分特征值以达到降阶方面，按照书上的判断方法，几乎不可避免的得使用goto语句，这使得程序流程变得混乱。经过仔细分析后，可以发现能否求解Ak+1的部分特征值然后实现矩阵Ak+1降阶的关键在于判断一下拟上三角阵a11a12??a1ma21a22a23?a2m0??am-1m-2am-1m-1am-1m0?0 amm-1amm是否存在下列形式的子块Bk+1=am-1m-1am-1m0ammCk+1=am-2m-2am-2m-1am-2m0am-1m-1am-1m0amm-1amm而其本质就是 关于求从属于矩阵A的实特征值λi的特征向量 在求出A的全部特征值后，求从属于矩阵A的实特征值λi的特征向量的时候，一个思路是对矩阵做平移：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≥1是特征值λi的重数，而R和矩阵B有相同的秩，因而方程Rx=0有p个线性无关的解，即为矩阵A的p重特征值λi的p个特征向量，而求解方程Rx=0的p个线性无关的解可以用如下算法实现：1. 因为矩阵R的秩为10-p，所以一定存在i1 计算结果 进行拟上三角A(n-1)为： Table SEQ Table \* ARABIC 1矩阵An-1的前