数值分析报告.doc

北航 2008 级研究生 数值分析 A 计算实 习题 2008 年 11 月 15 日 目目 录录 1 引言引言 1 2 算法设计方案算法设计方案 1 3 特殊情况处理特殊情况处理 4 4 计算计算结果结果 4 5 结论结论 7 6 参考文献参考文献 7 7 附录附录 8 1 1 引言引言 算法背景 算法背景 幂法 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 的全部特征值 并对其中的每一个实特 10 10 ij a 征值求相应的特征向量 要求程序输出矩阵 A 经过拟上三角化后所得到的矩阵 对矩阵 1 n A 进行 QR 分解方法结束后所得的矩阵 矩阵 A 的全部特征值及相应于实特 1 n A 1 2 10 iii R Ii 征值的特征向量 采用 e 型输出数据 并且至少显示 12 位有效数字 要求迭代的精度水平为 已知 12 10 ij a sin 0 50 2 1 2 10 1 5cos 1 2 ij ij i j ij ij 2 2 算法算法 2 2 1 用带双步位移的用带双步位移的 QR 方法求实矩阵全部特征值的算法方法求实矩阵全部特征值的算法 1 使用矩阵的拟上三角化的算法把矩阵 A化为拟上三角矩阵 给定精度水平 n n R 1 n A 和迭代最大次数 L 0 矩阵 A 的拟上三角化的算法如下 记A 并记 对于 r 1 2 n 2 执行 1 A rr ij n n Aa a 若全为零 则令 转 e 否则转 b 2 3 r i r airrn 1 rr AA b 计算 若 则取 2 1 n r rir i r da 1 r rrrr csign ad 1 0 r rr a rr cd 2 1 r rrrrr hcc a c 令 1 2 0 0 rrrTn rrrrrrn r uac aaR d 计算 r T rrr pAuh r rrr qA uh T rrrr tp uh rrrr qt u 1 rrTT rrrr AAuu p e 继续 2 记 令 k 1 m n 1 1 1 n ij n n AAa 3 如果 则得到 A 的一个特征值 置 m n 1 转 4 否则转 5 1 k m m a k m m a 4 如果 m 1 则得到 A 的一个特征值 转 11 如果 m 0 则直接转 11 如果 1 1 k a m1 则转 3 5 求二阶子阵的两个特征值和 即计算二次方程 1 11 1 kk mmmm k kk m mm m aa D aa 1 s 2 s 的两个根和 2 1 1 det0 kk mmm mk aaD 1 s 2 s 6 如果 m 2 则得到 A 的两个特征值和 转 11 否则转 7 1 s 2 s 7 如果 则得到 A 的两个特征值和 置 m m 2 转 4 否则转 8 1 2 k mm a 1 s 2 s 8 如果 k L 则计算终止 未得到 A 的全部特征值 否则转 9 9 记 计算 1 k kij m m Aai jm 1 1 kk mmm m saa 1 1 11 kkkk mmm mm mmm taaaa 为 m 阶单位矩阵 2 kkk MAsAtI I 对作 QR 分解 kkk MQ R k M 1 T kkkk AQ A Q 对作 QR 分解与的计算算法如下 k M 1k A 记 对于 r 1 2 m 1 执行 1 1kij m m BMb r rij m m Bb 1k CA a 若全为零 则令 转 e 否则转 b 1 2 r i r birrm 11 rrrr BB CC b 计算 若 则取 2 m r rir i r db r rr rr csign bd 0 r r r b rr cd 2 r rrrr r hcc b c 令 1 0 0 rrrTm rr rrrrm r ubc bbR d 计算 T rrrr vB uh 1 T rrrr BBu v T rrrr pC uh rrrr qC uh T rrrr tp uh rrrr qt u 1 TT rrrrrr CCuu p e 继续 10 置 k k 1 转 3 11 A 的全部特征值以计算完毕 停止计算 2 2 2 列主元素列主元素 Gauss 消去法求特征值对应的特征向量的算法消去法求特征值对应的特征向量的算法 列主元素 Gauss 消去法具体算法如下 记 1 1 2 ijij aa i jn 1 1 2 ii bb in 1 消元过程 对于 1 2 n 1 执行 1 选行号 使 k i max k kk i kik k i n aa 2 交换与以及与所含的数值 k kj a 1 k k i j ajk kn k k b k k i b 3 对于 i k 1 k 2 n 计算 kk ikikkk maa 1 1 2 kkk ijijikkj aam ajkkn 1 kkk iiikk bbm b 2 回代过程 nn nnnn xba 1 1 2 1 n kkk kkkjjkk j k xbaxaknn 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 001 1 54910107991e 001 1 94659186287e 000 8 78243638293e 002 9 25588938718e 001 6 03259944053e 0011 51886095647e 001 2 34787836242e 000 2 37237010494e 000 1 81929082221e 000 3 23780410155e 0012 20579844032e 001 2 10269266255e 000 1 81613808610e 0011 27883908999e 000 6 38057812440e 001 4 001 1 05568548241e 016 1 72827459997e 000 1 000 1 24383926270e 000 6 39975834174e 001 2 00283307904e 000 2 92494720612e 001 6 41283006839e 001 9 78399762128e 0022 55776357416e 001 5 39338381277e 017 3 16587386518e 017 1 29166953413e 000 1 11160351340e 000 1 000 1 30735603002e 000 1 80369917775e 001 4 24638535837e 001 7 98895523930e 0021 60881992807e 001 1 53346499662e 0175 96332140618e 017 3 27306215536e 017 1 56012629853e 000 8 12504939751e 001 4 42175683292e 001 3 58861612814e 002 4 69174231367e 001 2 73659505009e 001 7 35933465775e 002 1 30056273723e 016 3 09706001089e 017 8 55612756587e 0170 00000000000e 000 7 70777375519e 001 1 58305142574e 000 3 04284317680e 001 2 52871244603e 001 6 70992540145e 001 2 54461992908e 001 1 61021672477e 016 2 21157183737e 016 3 92548490401e 017 0 00000000000e 000 0 00000000000e 000 7 46345345694e 001 2 70836515702e 002 9 48652189368e 001 1 19587108150e 0011 92926561795e 002 1 36855018620e 0167 017 8 67261531120e 017 0 00000000000e 000 0 00000000000e 000 1 07207471836e 016 7 70180137436e 001 4 69762399062e 001 4 98825946801e 0011 001 2 78085130072e 017 6 70863078836e 017 8 47612717394e 0170 00000000000e 000 0 00000000000e 000 4 87977428748e 017 1 85482928622e 0167 01316709211e 0011 58218068848e 0013 86259461423e 001 2 12460444005e 017 1 70797975893e 016 2 79894228707e 017 0 00000000000e 000 0 00000000000e 000 4 0171 22262169189e 0160 00000000000e 000 4 84380760278e 0013 99277799518e 001 4 2 对矩阵进行对矩阵进行 QR 分解后所得的矩阵分解后所得的矩阵 由于一行放不下十个元素 故将矩阵每一行的十个元素分成两行列出 对拟三角化后的矩阵 A n 1 进行 QR 分解方法结束后的矩阵 1 53147700891e 000 2 31151347770e 000 5 002 1 56228767543e 001 3 95270218975e 001 4 60307123730e 001 4 92665761043e 001 1 54615717231e 0014 58107069729e 0029 01038975186e 002 3 19420475749e 000 5 63497701019e 001 2 71483239643e 001 6 42858659429e 001 1 03913640971e 000 2 87072983195e 000 7 00811288884e 001 8 03485266055e 001 7 59721306153e 001 2 000 3 55229096910e 0161 78505809730e 000 1 74399102792e 000 7 97374725643e 001 6 25624360652e 001 1 11600242269e 000 2 59559452545e 0013 42921412680e 0015 10113644468e 001 9 75562845054e 001 2 38895193340e 016 2 87790678694e 016 4 06561715572e 001 5 06868432472e 001 1 82855102452e 000 5 001 3 59510952371e 001 1 001 4 51225605836e 002 5 59740463210e 001 1 016 9 76111965681e 018 3 31049963415e 0161 48676607216e 000 3 71340891995e 001 1 69622855201e 001 1 74920002091e 001 1 73910145562e 001 2 43180817397e 001 1 67045385533e 001 1 95011078163e 017 7 21884669449e 017 8 32462282641e 017 4 58793017954e 017 6 55122789790e 001 9 39665990713e 001 2 83922343339e 001 5 65703210957e 002 8 96845791791e 002 3 88998224402e 001 1 55950784591e 0161 55024068510e 0161 39811287854e 016 4 23319394901e 017 2 36039665905e 016 4 57356256434e 001 1 25312729615e 000 2 30116600863e 001 2 99531681862e 001 4 26239601973e 001 4 0174 99418388088e 0174 85013749449e 017 1 25802674480e 017 1 44155149365e 016 6 10890682971e 017 4 66929122722e 001 8 08204282850e 001 1 55343840859e 001 1 63069004449e 001 2 43149869404e 0174 06334732673e 0172 38568246195e 017 1 93516231771e 017 7 96194310258e 017 3 33232735828e 017 1 35991797860e 016 1 59398041636e 001 9 12617646163e 0023 11904439743e 003 2 01582627298e 017 5 01986969674e 017 1 50289426382e 018 1 07390932683e 017 8 37186191687e 017 3 67783217967e 017 2 76682192903e 017 5 82867087928e 016 2 41769845651e 001 7 01517104511e 001 4 5 矩阵矩阵 A 的全部特征值的全部特征值 矩阵 A 的全部特征值 6 36062787575e 001 0 00000000000e 000i 5 65048899350e 002 0 00000000000e 000i 9 35588907819e 001 0 00000000000e 000i 9 80530956290e 001 1 001i 9 80530956290e 001 1 001i 1 48403982226e 000 0 00000000000e 000i 1 57754855711e 000 0 00000000000e 000i 2 32349621021e 000 8 93040517720e 001i 2 32349621021e 000 8 93040517720e 001i 3 38303961744e 000 0 00000000000e 000i 4 6 矩阵矩阵 A 的相应于实特征值的特征向量的相应于实特征值的特征向量 由于一行放不下十个元素 故将矩阵每一行的十个元素分成两行列出 特征值对应的特征向量 按行排列 不是实特征值的 不求对应的特征向量 1 07037468397e 001 7 12345583256e 002 3 90237917138e 001 4 46655513756e 0027 19034744275e 001 1 75815688440e 001 2 26537961544e 001 3 76886150502e 001 2 95625408505e 001 2 25577972559e 002 2 09017518582e 0012 00063796685e 001 3 89176061817e 001 2 77939118707e 0023 93236547906e 001 1 24703810952e 001 6 44810939558e 001 3 02779853108e 0012 91095263314e 001 4 09436555812e 002 8 05651548269e 002 4 61113468803e 002 1 50243791976e 0024 81208346184e 0023 53633892069e 001 2 08919077811e 001 1 55745966395e 001 8 19670427271e 001 3 50982512819e 001 2 88513852783e 002 不是实特征值 不求对应的特征向量 不是实特征值 不求对应的特征向量 5 60118116800e 001 7 79341508770e 001 1 33780114857e 002 2 77409287891e 001 3 00557588301e 003 2 53483408960e 0032 06284849088e 0021 10134829107e 0021 22486166525e 002 3 23620930409e 002 6 24962559917e 002 1 12090045150e 0022 49666705179e 0011 31364486150e 0013 83541256120e 001 8 001 1 24560352067e 0016 83561054721e 002 2 70587401777e 001 1 00519108146e 001 不是实特征值 不求对应的特征向量 不是实特征值 不求对应的特征向量 1 05221573436e 0012 0014 73017877776e 0012 60887478870e 0013 05805625140e 001 2 58379783221e 001 8 73379363951e 002 4 05454033853e 001 5 09013113720e 001 2 40925044562e 001 5 结论结论 QR 方法是目前求任意矩阵全部特征值和特征向量最有效的方法 利用矩阵的 QR 分解 通过逆序相乘产生对原矩阵的一系列正交相似变换 使其变化为一个近似的上三角矩阵 拟 上三角矩阵 来求全部特征值 算出特征值后对实特征值利用列主元素 Gauss 消去法求实特 征值对应的特征向量 利用 QR 分解将矩阵化为一个正交矩阵 Q 和一个上三角矩阵 R 左乘的 形式 从 QR 算法的构造过程可以看到算法的主要计算量出现在 QR 分解上 如果直接对矩 阵 A 用 QR 方法求全部特征值 那麽涉及的计算量是很大的 因此应该先对 A 作预处理 应 用中常先对 做正交相似变换将其化为上 Hessenberg 矩阵 H 然后再对 H 采用 QR 方法 可 以大大减少计算量 这里 Hessenberg 矩阵也称为拟三角矩阵 它的非零元素比三角矩阵多 了一条次对角线 采用此程序计算题目所给矩阵的特征值时 迭代次数 k 20 可见采用带双步位移的 QR 算法 可以减少迭代次数加速收敛 避免复数运算 在计算上是比较经济的 6 参考文献参考文献 1 冯康等编 数值计算方法 国防工业出版社 1978 年 2 颜庆津编 数值分析 第三版 北京航空航天大学出版社 2006 年 3 谭浩强编 C 程序设计 第二版 清华大学出版社 2004 年 7 附录附录 附录 A 全部源程序 包括主程序和每个子程序的功能注释 带双步位移的 QR 算法求矩阵特征值和特征向量程序 声明 本程序为北航 2008 级研究生 数值分析 A 计算实习题第一题 include include void input int m double n 10 数据输入子程序 将 aij 的数据输入二维数组 int i j for i 1 i m 1 i for j 1 j0 b 1 else b 1 return b void mul int n1 int n2 int n3 double a 10 double m1 10 double m2 10 定义矩阵乘法的子程序 int i j l for i 0 i n1 i for j 0 j n3 j a i j 0 for l 0 l n2 l a i j a i j m1 i l m2 l j void equ double A 10 double B 10 定义矩阵赋值的子程序 int i j for i 0 i 10 i for j 0 j 10 j A i j B i j void sub int n1 int n2 double m double a 10 double m1 10 double m2 10 定义矩阵减法的子程 序 int i j for i 0 i n1 i for j 0 j n2 j a i j 0 a i j m1 i j m m2 i j void add int n1 int n2 double m double a 10 double m1 10 double m2 10 定义矩阵加法的子程 序 int i j for i 0 i n1 i for j 0 j 0 x 0 0 a sqrt a a 4 b 2 x 0 1 0 x 1 0 a sqrt a a 4 b 2 x 1 1 0 else x 0 0 a 2 x 0 1 sqrt 4 b a a 2 x 1 0 a 2 x 1 1 sqrt 4 b a a 2 void qr int n int m double B 10 double C 10 矩阵的拟上三角化和 qr 分解的集成子程序 当 m 1 时为拟上三角化 m 0 时为 qr 分解 double Br 10 10 BrT 10 10 Br1 10 10 double Cr 10 10 CrT 10 10 Cr1 10 10 double vr 10 10 vrT 10 10 ur 10 10 urT 10 10 uv 10 10 double pr 10 10 prT 10 10 qr 10 10 tr 10 10 wr 10 10 double wu 10 10 up 10 10 double dr cr hr int i j r l equ Br B equ Cr C for r 1 r n m r Br Cr for i r 1 m i n 1 i l 0 if Br i 1 r 1 0 l 1 continue else break if l 1 continue else dr 0 for i r m i n 1 i dr dr Br i 1 r 1 Br i 1 r 1 dr sqrt dr cr sign Br r m 1 r 1 dr hr cr cr cr Br r m 1 r 1 for i 1 i n 1 i if i r m ur i 1 0 0 else if i r m ur i 1 0 Br i 1 r 1 cr else ur i 1 0 Br i 1 r 1 if m 0 m 0 时为 qr 分解 for i 1 i n 1 i for j 1 j n 1 j BrT i 1 j 1 Br j 1 i 1 mul n n 1 vr BrT ur for i 1 i n 1 i vr i 1 0 vr i 1 0 hr for i 1 i n 1 i vrT 0 i 1 vr i 1 0 mul n 1 n uv ur vrT sub n n 1 Br1 Br uv equ Br Br1 for i 1 i n 1 i for j 1 j n 1 j CrT i 1 j 1 Cr j 1 i 1 mul n n 1 pr CrT ur for i 1 i n 1 i pr i 1 0 pr i 1 0 hr mul n n 1 qr Cr ur for i 1 i n 1 i qr i 1 0 qr i 1 0 hr for i 1 i n 1 i prT 0 i 1 pr i 1 0 mul 1 n 1 tr prT ur tr 0 0 tr 0 0 hr sub n 1 tr 0 0 wr qr ur for i 1 i n 1 i urT 0 i 1 ur i 1 0 mul n 1 n wu wr urT mul n 1 n up ur prT sub n n 1 Cr1 Cr wu sub n n 1 Cr Cr1 up if m 1 m 1 时为拟上三角化 equ Br Cr equ C Cr void gauss int n double x 10 10 double A 10 10 double b 10 10 列主元素 gauss 消去法子程序 int i j k int ik double aa double AK 10 10 bK 10 10 m 10 10 equ bK b equ AK A for k 1 k n k aa AK k 1 k 1 ik k for i k 1 iaa aa AK i 1 k 1 ik i else continue for j 1 j n 1 j aa AK k 1 j 1 AK k 1 j 1 AK ik 1 j 1 AK ik 1 j 1 aa aa bK k 1 0 bK k 1 0 bK ik 1 0 bK ik 1 0 aa for i k 1 i n 1 i m i 1 k 1 AK i 1 k 1 AK k 1 k 1 for j k j0 k aa 0 for j k 1 j n 1 j aa aa AK k 1 j 1 x j 1 0 x k 1 0 bK k 1 0 aa AK k 1 k 1 void main 主函数 double ev 10 2 ss 2 2 double x 10 10 xx 10 10 b 10 10 double A 10 10 A1 10 10 A2 2 2 double M 10 10 M1 10 10 AA 10 10 I 10 10 double eps 1 0e 12 x1 定义精度水平 1 0e 12 double s t int i j k r 1 int n 10 m n l 1 L 40 定义迭代最大次数 L 40 char c1 c2 i FILE fp for i 1 i n 1 i 定义单位矩阵 for j 1 j n 1 j if i 1 j 1 I i 1 j 1 1 else I i 1 j 1 0 input m A 调用数据输入子程序 equ A1 A qr n 1 A1 A1 进行拟上三角化 if fp fopen data txt w NULL printf cannot write this file n fprintf fp 矩阵 A 经过拟上三角化后所得的矩阵 A n 1 按行排列 n 拟上三角化后矩阵的输 出 for i