数值分析报告.doc

北航2008级研究生《数值分析A》计算实习题2008年11月15日目录1引言12算法设计方案13特殊情况处理44计算结果45结论76参考文献77附录81引言算法背景：幂法Jacobi方法及QR方法是求矩阵特征值和特征向量的常用数值方法它们都是造构造迭代产生的矩阵序列来达到目的的。幂法计算简单特别适用于高阶稀疏矩阵但其收敛速度不能令人满意要想加快幂法的收敛速度可采用反幂法及位移技术。Jacobi方法是古典方法它收敛快精度高便于并行计算且算法稳定。用Jacobi方法求出的特征向量在较好的正交性不过它的计算量较大当阶数n增大时收敛速度减慢因此Jacobi方法适用于求低阶的对称矩阵的全部特征值和特征向量。由J.G.Francis提出的QR算法的基本思想源于LR算法即对矩阵分解获得一个相似于原矩阵的序列使其收敛到一个易于求得特征值的形式。LR算法比QR算法收敛速度快但不稳定。QR方法是60年代发展起来的被人们称为数值数学最值得注意的算法之一它是目前求任意矩阵全部特征值和特征向量最有效的方法。利用矩阵的QR分解通过逆序相乘产生对原矩阵的一系列正交相似变换使其变化为一个近似的上三角矩阵来求全部特征值。这里QR分解是指将矩阵化为一个正交矩阵Q和一个上三角矩阵左乘的形式。但是基本QR算法的收敛往往过慢因此通常采用带原点位移的QR算法。又由于A是实矩阵它可能有复特征值如果采用一般带原点位移的QR算法变换中带有复位移量则造成复数运算。故而采用带双步位移的QR算法可以减少迭代次数加速收敛避免复数运算在计算上是比较经济的。2算法设计方案2.1题目试用带双步位移的QR分解法求矩阵A=的全部特征值并对其中的每一个实特征值求相应的特征向量。要求程序输出矩阵A经过拟上三角化后所得到的矩阵；对矩阵进行QR分解方法结束后所得的矩阵；矩阵A的全部特征值及相应于实特征值的特征向量。采用e型输出数据并且至少显示12位有效数字。要求迭代的精度水平为=。已知2.2算法2.2.1用带双步位移的QR方法求实矩阵全部特征值的算法（1）使用矩阵的拟上三角化的算法把矩阵A化为拟上三角矩阵；给定精度水平和迭代最大次数L。【矩阵A的拟上三角化的算法如下：记A并记。对于r=12n-2执行（a）若全为零则令转（e）；否则转（b）。（b）计算（若则取）（c）令（d）计算（e）继续】（2）记令k=1m=n。（3）如果则得到A的一个特征值置m=n-1转（4）；否则转（5）。（4）如果m=1则得到A的一个特征值转（11）；如果m=0则直接转（11）；如果m1则转（3）。（5）求二阶子阵的两个特征值和即计算二次方程的两个根和。（6）如果m=2则得到A的两个特征值和转（11）；否则转（7）。（7）如果则得到A的两个特征值和置m=m-2转（4）；否则转（8）。（8）如果k=L则计算终止未得到A的全部特征值否则转（9）。（9）记计算（为m阶单位矩阵）（对作QR分解）【对作QR分解与的计算算法如下：记。对于r=12m-1执行（a）若全为零则令转（e）；否则转（b）。（b）计算（若则取）（c）令（d）计算（e）继续】（10）置k=k+1转（3）。（11）A的全部特征值以计算完毕停止计算。2.2.2列主元素Gauss消去法求特征值对应的特征向量的算法列主元素Gauss消去法具体算法如下：记1.消元过程对于12n-1执行（1）选行号使（2）交换与以及与所含的数值（3）对于i=k+1k+2n计算2.回代过程3特殊情况处理(1)由于拟上三角化和QR分解的算法有相同的部分故将两部分的子程序集成在一起（当m=1时为矩阵的拟上三角化程序和当m=0时为QR分解程序）(2)当输出“未求出矩阵A的全部特征值”时检查L的大小增大L的数值再计算。4计算结果4.1矩阵A经过拟上三角化后得到的矩阵（由于一行放不下十个元素故将矩阵每一行的十个元素分成两行列出）矩阵A经过拟上三角化后所得的矩阵A(n-1)：-8.82751675883e-001-9.93313649183e-002-1.10334928599e+000-7.60044358564e-0011.54910107991e-001-1.94659186287e+000-8.78243638293e-002-9.25588938718e-0016.03259944053e-0011.51886095647e-001-2.34787836242e+0002.37237010494e+0001.81929082221e+0003.23780410155e-0012.20579844032e-0012.10269266255e+0001.81613808610e-0011.27883908999e+000-