数值分析报告.doc

1、北航2008级研究生数值分析A计算实习题2008年11月15日目 录1 引言12 算法设计方案13 特殊情况处理44计算结果45结论76参考文献77附录81 引言算法背景：幂法，Jacobi方法及QR方法是求矩阵特征值和特征向量的常用数值方法，它们都是造构造迭代产生的矩阵序列来达到目的的。幂法计算简单，特别适用于高阶稀疏矩阵，但其收敛速度不能令人满意，要想加快幂法的收敛速度可采用反幂法及位移技术。Jacobi方法是古典方法，它收敛快，精度高，便于并行计算且算法稳定。用Jacobi方法求出的特征向量在较好的正交性，不过它的计算量较大，当阶数n增大时收敛速度减慢，因此Jacobi方法适用于求低阶的

2、对称矩阵的全部特征值和特征向量。由J.G.Francis提出的QR算法的基本思想源于LR算法，即对矩阵分解获得一个相似于原矩阵的序列，使其收敛到一个易于求得特征值的形式。LR算法比QR算法收敛速度快，但不稳定。QR方法是60年代发展起来的，被人们称为数值数学，最值得注意的算法之一，它是目前求任意矩阵全部特征值和特征向量最有效的方法。利用矩阵的QR分解，通过逆序相乘产生对原矩阵的一系列正交相似变换，使其变化为一个近似的上三角矩阵来求全部特征值。这里QR分解是指将矩阵化为一个正交矩阵Q和一个上三角矩阵左乘的形式。但是基本QR算法的收敛往往过慢，因此通常采用带原点位移的QR算法。又由于A是实矩阵，它

3、可能有复特征值，如果采用一般带原点位移的QR算法，变换中带有复位移量，则造成复数运算。故而采用带双步位移的QR算法，可以减少迭代次数加速收敛，避免复数运算，在计算上是比较经济的。2 算法设计方案2.1 题目 试用带双步位移的QR分解法求矩阵A=的全部特征值，并对其中的每一个实特征值求相应的特征向量。要求程序输出矩阵A经过拟上三角化后所得到的矩阵；对矩阵进行QR分解方法结束后所得的矩阵；矩阵A的全部特征值及相应于实特征值的特征向量。采用e型输出数据，并且至少显示12位有效数字。要求迭代的精度水平为=。已知 2.2 算法2.2.1用带双步位移的QR方法求实矩阵全部特征值的算法（1）使用矩阵的拟上三

4、角化的算法把矩阵A化为拟上三角矩阵；给定精度水平和迭代最大次数L。【矩阵A的拟上三角化的算法如下：记A，并记。对于r=1，2，n-2执行（a）若全为零，则令，转（e）；否则转（b）。（b）计算，（若，则取），（c）令（d）计算，（e）继续】（2）记，令k=1，m=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的两个特征

5、值和，置m=m-2，转（4）；否则转（8）。（8）如果k=L，则计算终止，未得到A的全部特征值，否则转（9）。（9）记，计算（为m阶单位矩阵）（对作QR分解）【对作QR分解与的计算算法如下：记，。对于r=1，2，m-1执行（a）若全为零，则令，转（e）；否则转（b）。（b）计算，（若，则取），（c）令（d）计算，（e）继续】（10）置k=k+1，转（3）。（11）A的全部特征值以计算完毕，停止计算。2.2.2 列主元素Gauss消去法求特征值对应的特征向量的算法列主元素Gauss消去法具体算法如下：记，1.消元过程对于1，2，n-1执行（1）选行号，使（2）交换与以及与所含的数值（3）对于i=

6、k+1,k+2, ,n计算，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.5491

7、0107991e-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-6.38057812440e-001-4.15407560380e-001-1.055685482

8、41e-0161.72827459997e+000-1.17146764279e+000-1.24383926270e+000-6.39975834174e-001-2.00283307904e+0002.92494720612e-001-6.41283006839e-0019.78399762128e-0022.55776357416e-001-5.39338381277e-017-3.16587386518e-017-1.29166953413e+000-1.11160351340e+0001.17134682410e+000-1.30735603002e+0001.80369917775

9、e-001-4.24638535837e-0017.98895523930e-0021.60881992807e-0011.53346499662e-0175.96332140618e-017-3.27306215536e-0171.56012629853e+0008.12504939751e-0014.42175683292e-001-3.58861612814e-0024.69174231367e-001-2.73659505009e-001-7.35933465775e-0021.30056273723e-016-3.09706001089e-0178.55612756587e-0170

10、.00000000000e+000-7.70777375519e-001-1.58305142574e+000-3.04284317680e-0012.52871244603e-001-6.70992540145e-0012.54461992908e-0011.61021672477e-016-2.21157183737e-016-3.92548490401e-0170.00000000000e+0000.00000000000e+000-7.46345345694e-001-2.70836515702e-002-9.48652189368e-0011.19587108150e-0011.92

11、926561795e-0021.36855018620e-0167.15151319080e-017-8.67261531120e-0170.00000000000e+0000.00000000000e+000-1.07207471836e-016-7.70180137436e-001-4.69762399062e-0014.98825946801e-0011.13769160378e-001-2.78085130072e-017-6.70863078836e-0178.47612717394e-0170.00000000000e+0000.00000000000e+000-4.8797742

12、8748e-0171.85482928622e-0167.01316709211e-0011.58218068848e-0013.86259461423e-001-2.12460444005e-017-1.70797975893e-016-2.79894228707e-0170.00000000000e+0000.00000000000e+0004.15375832524e-0171.22262169189e-0160.00000000000e+0004.84380760278e-0013.99277799518e-0014.2 对矩阵进行QR分解后所得的矩阵（由于一行放不下十个元素，故将矩阵

13、每一行的十个元素分成两行列出）对拟三角化后的矩阵A(n-1)进行QR分解方法结束后的矩阵：1.53147700891e+000-2.31151347770e+000-5.13352162117e-002-1.56228767543e-0013.95270218975e-001-4.60307123730e-001-4.92665761043e-0011.54615717231e-0014.58107069729e-0029.01038975186e-002-3.19420475749e+000-5.63497701019e-001-2.71483239643e-001-6.4285865942

14、9e-0011.03913640971e+0002.87072983195e+000-7.00811288884e-001-8.03485266055e-001-7.59721306153e-0012.13060816678e+0003.55229096910e-0161.78505809730e+000-1.74399102792e+0007.97374725643e-001-6.25624360652e-001-1.11600242269e+0002.59559452545e-0013.42921412680e-0015.10113644468e-001-9.75562845054e-00

15、1-2.38895193340e-016-2.87790678694e-016-4.06561715572e-001-5.06868432472e-0011.82855102452e+0005.18941002166e-001-3.59510952371e-001-1.13845216101e-001-4.51225605836e-0025.59740463210e-0011.13201004533e-016-9.76111965681e-0183.31049963415e-0161.48676607216e+0003.71340891995e-001-1.69622855201e-0011.

16、74920002091e-001-1.73910145562e-001-2.43180817397e-0011.67045385533e-001-1.95011078163e-017-7.21884669449e-017-8.32462282641e-017-4.58793017954e-017-6.55122789790e-001-9.39665990713e-0012.83922343339e-001-5.65703210957e-002-8.96845791791e-002-3.88998224402e-0011.55950784591e-0161.55024068510e-0161.3

17、9811287854e-016-4.23319394901e-017-2.36039665905e-016-4.57356256434e-001-1.25312729615e+000-2.30116600863e-0012.99531681862e-001-4.26239601973e-0014.17370606398e-0174.99418388088e-0174.85013749449e-017-1.25802674480e-017-1.44155149365e-016-6.10890682971e-017-4.66929122722e-0018.08204282850e-001-1.55

18、343840859e-0011.63069004449e-0012.43149869404e-0174.06334732673e-0172.38568246195e-017-1.93516231771e-017-7.96194310258e-017-3.33232735828e-017-1.35991797860e-016-1.59398041636e-0019.12617646163e-0023.11904439743e-003-2.01582627298e-017-5.01986969674e-017-1.50289426382e-0181.07390932683e-017-8.37186

19、191687e-017-3.67783217967e-0172.76682192903e-017-5.82867087928e-016-2.41769845651e-0017.01517104511e-0014.5 矩阵A的全部特征值矩阵A的全部特征值:6.36062787575e-001+0.00000000000e+000i5.65048899350e-002+0.00000000000e+000i9.35588907819e-001+0.00000000000e+000i-9.80530956290e-001+1.13948912743e-001i-9.80530956290e-001-

20、1.13948912743e-001i-1.48403982226e+000+0.00000000000e+000i1.57754855711e+000+0.00000000000e+000i-2.32349621021e+000+8.93040517720e-001i-2.32349621021e+000-8.93040517720e-001i3.38303961744e+000+0.00000000000e+000i4.6 矩阵A的相应于实特征值的特征向量（由于一行放不下十个元素，故将矩阵每一行的十个元素分成两行列出）特征值对应的特征向量(按行排列)(不是实特征值的,不求对应的特征向量):

21、-1.07037468397e-001-7.12345583256e-002-3.90237917138e-0014.46655513756e-0027.19034744275e-001-1.75815688440e-0012.26537961544e-001-3.76886150502e-001-2.95625408505e-001-2.25577972559e-0022.09017518582e-0012.00063796685e-001-3.89176061817e-0012.77939118707e-0023.93236547906e-0011.24703810952e-001-6.4

22、4810939558e-0013.02779853108e-0012.91095263314e-001-4.09436555812e-002-8.05651548269e-002-4.61113468803e-0021.50243791976e-0024.81208346184e-0023.53633892069e-001-2.08919077811e-0011.55745966395e-001-8.19670427271e-0013.50982512819e-001-2.88513852783e-002不是实特征值,不求对应的特征向量不是实特征值,不求对应的特征向量5.60118116800

23、e-001-7.79341508770e-001-1.33780114857e-0022.77409287891e-001-3.00557588301e-0032.53483408960e-0032.06284849088e-0021.10134829107e-0021.22486166525e-002-3.23620930409e-002-6.24962559917e-0021.12090045150e-0022.49666705179e-0011.31364486150e-0013.83541256120e-001-8.15944310353e-0011.24560352067e-0016

24、.83561054721e-002-2.70587401777e-001-1.00519108146e-001不是实特征值,不求对应的特征向量不是实特征值,不求对应的特征向量1.05221573436e-0012.18364197327e-0014.73017877776e-0012.60887478870e-0013.05805625140e-0012.58379783221e-001-8.73379363951e-002-4.05454033853e-001-5.09013113720e-001-2.40925044562e-0015 结论QR方法是目前求任意矩阵全部特征值和特征向量最有效

25、的方法。利用矩阵的QR分解，通过逆序相乘产生对原矩阵的一系列正交相似变换，使其变化为一个近似的上三角矩阵（拟上三角矩阵）来求全部特征值，算出特征值后对实特征值利用列主元素Gauss消去法求实特征值对应的特征向量。利用QR分解将矩阵化为一个正交矩阵Q和一个上三角矩阵R左乘的形式。从QR 算法的构造过程可以看到算法的主要计算量出现在QR分解上，如果直接对矩阵A用QR方法求全部特征值，那麽涉及的计算量是很大的，因此应该先对A作预处理。应用中常先对 做正交相似变换将其化为上Hessenberg矩阵H，然后再对H采用QR方法，可以大大减少计算量，这里Hessenberg矩阵也称为拟三角矩阵，它的非零元素

26、比三角矩阵多了一条次对角线。采用此程序计算题目所给矩阵的特征值时，迭代次数k=20，可见采用带双步位移的QR算法，可以减少迭代次数加速收敛，避免复数运算，在计算上是比较经济的。6 参考文献【1】 冯康等编：数值计算方法，国防工业出版社，1978年。【2】 颜庆津编：数值分析（第三版），北京航空航天大学出版社，2006年。【3】 谭浩强编：C程序设计（第二版），清华大学出版社，2004年。7 附录附录A 全部源程序（包括主程序和每个子程序的功能注释）/带双步位移的QR算法求矩阵特征值和特征向量程序/声明：本程序为北航2008级研究生数值分析A计算实习题第一题。# include # includ

27、e void input(int m,double n10)/数据输入子程序（将aij的数据输入二维数组）int i,j;for(i=1;im+1;i+)for(j=1;j0)b=1;elseb=-1;return b;void mul(int n1,int n2,int n3, double a10,double m110,double m210)/定义矩阵乘法的子程序int i,j,l;for(i=0;in1;i+)for(j=0;jn3;j+)aij=0;for(l=0;ln2;l+)aij=aij+m1il*m2lj;void equ(double A10,double B10)/定义

28、矩阵赋值的子程序int i,j;for(i=0;i10;i+)for(j=0;j10;j+)Aij=Bij;void sub(int n1,int n2,double m,double a10,double m110,double m210)/定义矩阵减法的子程序int i,j;for(i=0;in1;i+)for(j=0;jn2;j+)aij=0;aij=m1ij-m*m2ij;void add(int n1,int n2,double m,double a10,double m110,double m210)/定义矩阵加法的子程序int i,j;for(i=0;in1;i+)for(j=0;jn2;j+)aij=0;aij=m1ij+m*m2ij;void function(double x2,double D2)/求解二次方程X*X+a*X+b=0的子程序do