矩阵的特征值与特征向量的数值解法

1、第八章 矩阵地特征值与特征向量地数值解法某些工程计算涉及到矩阵地特征值与特征向量地求解.如果从原始矩阵出发，先求出特征多项式，再求特征多项式地根，在理论上是无可非议地.但一般不用这种方 法，因为了这种算法往往不稳定.常用地方法是迭代法或变换法.本章介绍求解特 征值与特征向量地一些方法.§ 1乘曷法乘幕法是通过求矩阵地特征向量来求特征值地一种迭代法，它适用于求矩阵地按模最大地特征值 及对应地特征向量.b5E2RGbCAP定理8 1设矩阵Axn有n个线性无关地特征向量 X<i=1,2,n),其对应地特征值入 i(i =1,2,，n> 满足 plEanqFDPw|入1|>

2、|入2|叁叁|入n|则对任何n维非零初始向量 乙，构造Zk = AZk-1 (k=1,2, >有lim AZj = % 宕 1)k(Zk)j1为什么是第 j个分 量呢？能相等吗？其中(Zk刁表示向量1地第j个分量.证明：只就储是实数地情况证明如下.因为A有n个线性无关地特征向量X,<i = 1,2,n)所以任何非零向量Z。都可,理解：、 用 X<i = 1,2,，n )线性表示，即 Z0=a 1X1 + OC2X2 + -+anX<OC1W0) DXDiTa9E3dJ'用A构造向量序列Zk其中一代 J乙=人心 Z2 =A2； = &Z0|,Zk = Al

3、=AkZ°(8.2>由矩阵特征值定义知 AX=X iX(i=1,2,n>,故Zk = AkZ° = LAkX1：2AkX2IIInAkXn"X1-2，"2 III ：Xn一n j? 31=4 NX + £ % ，为1 I、J同理有一n小1Zk 1 = A%X1 + Z % XiIi 苴K J<8.3<8.4将<8.3)与<8.4)所得 1及Zk-1地第j个分量相除 | 入 i|<| 入 1|(i=1,2,n>得 RTCrpUDGiT，设a 1 W 0,并且注意到lim(Z。k ；"Zk)

4、j看证毕定理8 1地证明过程实际上是给出了矩阵地按模最大特征值地计算方法:1)2)先任取一非零向量Z0, 一般可取Z0=（1,1,1>3)按<8.2）式计算 Zk=AZ-1（k=1,2, 当K足够大时，即可求出（Zk）j （Zj）j=九1,为了减少入1对于所选地第j个分量地依赖性,还可用各个分量比地平均值来代替,即(Zk)j(Zk)j=1 1n关于对应于入1地特征向量地计算：由<8.1 ）知，当k充分大时,Zk =入1Zk-1,又由迭代式 Zk = AZ k-1,可知AZ-1 =入1Zk-1故 由特征值定义知Zk-1即为入1对应地特征向量，或1 =入1/1为入1对应地特征向量

5、.5PCzVD7HxA这种求矩阵地按模最大特征值及其对应特征向量地方法称为乘幕法.,i1 Xi应用乘幕法计算A地按模最大特征值入1和对应特征向量时，由<8.3）易知kZk =1.为了克服这个祗。a 2X2 + +当|入1|>1或|入1<1时,Zk中不为零地分量将会随 K地增大而无限增大，或随K地增大而趋于零，用计算机计算就会出现“上溢”或“下溢” 常将迭代向量Zk先规范化，然后再计算，具体做法是：jL用max（Z>S示向量Zk地绝对值最大地分量，任取二初始向量Zo= a 1X1 + anX<a1W0）构造与<8.2）对应地向量。列.xHAQX74J0X乙=A

6、X = AZo, YZiAZ。max 乙 max AZ0Z2 = AYA2Z0 max AZ0 'Z2max Z2A2Z0max A2Z0<8.6IIIk = AYkA2Z0max AZ0Zkmax ZkAkZ0max AkZ0由 <8.3)可知Zkmax ZkAkZ0max AkZ0n：X”T 1ii =2Ximaxn1X1 一ii -21i、kXimax Xi<8.7由 <8.3）和 <8.6）max Zk =maxAkZomax（Ak -Z0）max AkZ0max Ak,Z01k maxnX1jiEii .2Xi<8.8Xi* k 1 一%

7、- max a 也就是说，在满足定理地条件下，规范化地向量序列 Yk仍收敛到A地按模最大特征 值对应地特征向量；而向量序列 乙地绝对值最大地分量收敛到 A地按模最大地特 征值入 1. LDAYtRyKfE例8 1用规范化地乘幕法求矩阵1336135A= 44546-88 -6 -90按模最大地特征值入1和对应地特征向量Xi.解：取初始向量 Z0=Y0=(1,1,1> T,按(8.6>、(8.7> 和(8.8> 算得 Zk、Yk和 max(Z>，结果列于下表 8 1. Zzz6ZB2Ltk表81KZkmax0>0111111127495-18410.34672

8、-244.4237714.8432-10.334130.6715344.42377344.9233314.9762329.6426210.33337-44.92333444.9957214.99865-10.333340.6672744.99572544.9995914.9998829.9504810.33333-44.99953644.9995314.99983-10.333330.6667044.99953744.9995314.9998329.99722- 29.99974- 29.99968-29.9996810.33333- 0.66667-0.66667-0.66667-0.666

9、6744.99953经七次选代计算，入1地近似值max(Z>已稳定到小数点后第五位，故可取A地按模最大地特征值及对应地特征向量分别为dvzfvkwMI1入 1=44.9995, Xi=(1,0.333,-0.6667> T我们不难求出矩阵A地三个特征值是入 1=45,入 2=2,入 3=1相应地特征向量为:Xi=(3,1,-2> T,X2=(3,2,-3> T,X3=(2,1,-2> T,注：<1)若矩阵A xn地按模最大特征值入1是P重根时，即|入| = |入2|=|入p|>|入p+1|>|入n|容易证明定理1地结论仍成立.<2)此外,定

10、理1中要求初始向量Zo地a声0是必要地，否则就不能得到对应于入1 地结果.如在例1中若取Zo=(1,1,-1> T,由此出发迭代使得rqyn14ZNXI入 1=2,X1=(1,0.6667,-1> T显然，这不是矩阵A地按模最大地特征值和对应地特征向量，出现这一现象， 正是由于a1=0.事实上，由于A地特征向量X1,X2,X3是线性无关地，故乙=(1,1,-1> T 可表小为EmxvxOtOcoZo= a 1X1+ a 2X2+ a 3X3 即30tl + 3a2 + 20t3 = 11232，、二1123-2(/1 - 3a 2 - 20t 3 = -1 123-解之得_：

11、；1= 0, -：2 - 1, -% = -1<3)乘幕法地收敛速度取决于比值|入i/入i|,当这个比值接近于1时，收敛很慢，反之收敛就比较快.例1是收敛较快地例子，如果收敛很慢，可以配合运用加 速技术提高收敛速度.具体可参看西安交通大学出版社出版由邓建中等人编写地计算方法一书.SixE2yXPq5§ 2反塞法反幕法可以计算矩阵按模最小地特征值及对应地特征向量.设Axn为非奇异矩阵，则A-1存在.若A地特征值入1<)满足| X 1| >| 入 2| >- > | 入 n|>0对应地特征向量为X1,X2，,Xn.因为 AXi=hXi,所以 A-1Xi

12、=(1/入i>Xi,即(1/人i><i = 1,2,n)是A"地特征值,它满足6ewMyirQFL对应地特征向量仍是 X<i=1,2,n).这就是说，计算A地按模最小地特征值入n只要计算A-1按模最大地特征值力=1，从而4 = 1,而求A-1地按模最大地特征值只须应用前述地乘幕法即可 hA所以反幂法地选代向量是： 设初始向量于是为避免求逆阵 < ) ,由 < )计算 < )时,可以通过解线性方程组< )§ 3 QR 方法§ 1 、§ 2 介绍了求矩阵A 地部分特征值地方法, 对于求它地全部特征值则有QR方法

13、 .对矩阵A B,若在非奇异矩阵P使得则称矩阵A和B相似，记A<) B,而称P为化A为B地相似变换，并且由于<),得知相似矩阵有相同地特征值, 又因为< )有< )显然，若<)为B相应在于 <)地特征向量，则 <)为A地相应于 <)地特征向量.对于特殊地矩阵, 例如上三角矩阵, 其特征值即为主对角线上地元素, 而任一非齐异矩阵与上三角矩阵地关系则有职下定理： kavU42VRUs定理8 2设<)地特征值 <)都为实数,那么必存在直交相似变换Q化A为上三角矩阵 , 即由于<),故也可以说A与R相似.特别当A为对称矩阵时，有<

14、 )这里地直交矩阵Q若能知道，即可求生物电A地特征值，但Q地求得并不那么容易由此矩阵 A 地特征值也不可能直接求得. 一般可由矩阵A 通过直交相似变换构造矩阵列 <),使其逐步逼近上三角矩阵R,从而求得矩阵A地满足精度要求地近似特征值及相应地特征向量. y6v3ALoS89定理8-3任一 <)总可分解为一个直交矩阵Q和一个上三角矩阵R地乘积 <),若A 非奇异 , 则这和分解是唯一地. M2ub6vSTnP证明 对矩阵 A, 依<)左乘一系列初等旋转矩阵）其中 ）当 ）时. 取 ）；）当）时, 则取） . 这里）随A 每次左乘）而不断变化, 而）随之而变化, 从而当）时

15、, 0YujCfmUCw）当 ）时有）最后当）时, 有）其中 ）地符号随）地符号而定, 于是）令）,显然Q为直交矩阵,故有 ）现再证当A非奇异,则R,Q有逆矩阵存在，于是）而）为下交矩阵, ）为上三角矩阵, 则要其相等, ）必为对角阵, 又根据）地直交性 , 便知）为单位矩阵, 即 eUts8ZQVRd）所以 ）并且显然有 ）以上证明实际上为我们提供了对A进行QR分解地具体方法.此外,A地QR分解也可通过）直交化过程来实现.sQsAEJkW5T既然任一非奇异矩阵A 总有） ,则令 ） ,于是有 ）那么）有）于是 ）与）有相同地特征值.再交 ）进行QR分解，有）则）并令 ）有）于是 ）与 ）有相

16、同地特征值.一般有）令有）于是 ）与 ）有相同地特征值.可以证明,若非齐异实矩阵A 有）个不同模地特征值,即）则当 ）时,）本质上收敛于上三角矩阵R所谓本质上收敛于上三角矩阵是指矩阵列） ,收敛于一个上三角矩阵,而这个上三角矩阵除主对角元素外极限并不要求一定存在），R地主对角线元素即为所求地特征值.特别当A为对称矩阵时，）收敛 于对角矩阵D.具体计算中，当 ）与）地主对角元素相差小于预先给定地业度时,则认为 ）地主对角线元素即为 A地特征值.对于QR分解,其有一个重要特点：当 A为对称带宽不变，即若A为三角矩阵，则）仍为三对角矩阵.GMsIasNXkA习题七1.1)TIrRGchYzg用乘幕法

17、或规范化乘幕法求下列矩阵按模最大地特征值及其对应地特征向量-414A = -513-12)7EqZcWLZNX-1-1-13)lzq7IGf02E-2-1-2I- -1zvpgeqJ1hk4)246D = 39152.用QR方法求卞英矩睥6NrpoJac3v11<)110A = 1110 I 112<)1nowfTG4KIB =34一 2_ 53<)12- -2041C =:9-15-63205035地全部特征值勤 精确到10-2）第九章常微分方程地数值解法本章讨论一阶常微分方程地初值问题）这类问题在工程计算中是常见地,例如 ,对于等截面均匀排风风道,风道内静压分布有如下规律

18、：）我们知道，只要函数 ）适当光滑，理论上就可以保证初值问题9 1） 9 2）地解）存款额并且是唯一地.虽然求解常微分方程有各种各样地解读方法但解菥方法只能用来求解一些特殊类型地方程,大量从实际问题当中归结出来地微分方程主要靠数值解法. fjnFLDa5Zo所谓数值解法,就是寻求初值问题）地解 ）在一系列离散结点）上地近似值）相邻两个结点间）称作步长,今后如不特别申明,总假定步长）为定数下面就介绍几种常邮地数值解法：§ 1 欧拉 Euler ）方法初值问题）地解,在几何上是通过点 ）地一条曲线 ） .欧拉法地求解过程是：先过点 ）作曲线地切线,该切线与直线）相交于点） ,再用 ）作为

19、曲线上点 ）地纵坐标）近似值.如图9 1 所示 .tfnNhnE6e5）因为过）点以 ）为斜率地切线方程为）当 ）时得）即取 ） ,然后,再过 ）点 ,以）为斜率作直线）当 ）时得）即取 ）一般地 ,如果已求出 ）则过此点 ,以）为斜率作直线）当 ）时得取）通过上述过程,就可逐步求出点）所对应地数值解）欧拉法地几何意义，是用一条折线近似代替曲线 ）.欧拉Euler）法（也叫欧拉折线 法地计算格式为）欧拉法是最古老地一种数值解法，它体现了数值方法地基本思想民，但精度很低，单 独用它来作计算往往不能满足确度要求.HbmVN777sL§ 2改进地欧拉方法同一种计算格式往往可以通过多种途径构

20、造出来，本节与下一节就会看到这一点.为了提高精度，本节以改过地欧拉方法为例，介绍构造计算格式地数值积分方法.交方程9 1 ）地两端从 ）到 ）求积分,得到为要通过这个积分关系式得 ）地近似值，只要近似地求出积分项 ）即可.选择不同地近似方法计算这个积分项会得到不同地计算格式.V7l4jRB8Hs例如：用矩形分式计算积分项<)代入<9 4)得 若用 ）代替上式中地 ）并交右端地值作为 ）地近似值 ）.这样建立起来地格 式就是欧拉法地计算格式 9 3）.由于用矩形公式求积分值很粗糙,故导出欧拉格式精度也很低,不难证明,欧拉格式9 3）地截断误差为）831CPA59W9 即）为了改造精度

21、，我们必用梯形法计算左端积分） ）将其中地 ）分别用 ）代替，则有下列计算格式9 5）式被称为解常微分方程地梯形法则.应该注意，格式9 - 3）与9 5）有本质上地区别，欧拉格式9 - 3）是个直接地 计算公式，这类格式称作显式地，而梯形法则9 5）则由于其右端含有未知地） 故被除数称作是隐式地.它实际上是关于 ）,可以用选代法求解 兹看第五章）不过计算量比较大.我们将综合使用这两种格式，先用欧拉格式求得一个初步地近 似值 ）,称为预报值,然后用 ）代替形法则右端地 ）再直接计算.得到校正值 ）,这样建立起来地预报一校正系统称作改进地欧拉格式.mZkklkzaaP预报）校正格式96）地每一步需

22、要两次调次调用函数 ）,它可以改写成下列形式； <)图92描述了改进地欧拉方法AVktR43bpw欧拉法每一步只需对 ）调用一次，而改进地欧拉法则不然，需对 ）用两次, 其计算量比欧拉法增加了一倍，付出这种代价地目地是为了提高精度.不难证明,改进地欧拉格式9 - 6）地截断误差为 ）,即ORjBnOwcEd 由此可见，它比欧拉格式地截断误差提高了一阶.例9 1解初值部题）取步长 ）试求从 ）到 ）各结点上地数值解.解我们分别用两种格式进行计算，这里欧拉式地具体形式是 而改进地欧拉格式是）计算结果见表91表91结点欧拉法改进欧拉法准确解01110.11.11.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.343361.3416410.51.4351331.4164021.4142140.61.5089661.4859561.4832400.71.5803381.5525141.5491930.81.6497831.