第八章矩阵的特征值与特征向量的数值解法.doc

个人收集整理 勿做商业用途
第八章 矩阵的特征值与特征向量的数值解法
某些工程计算涉及到矩阵的特征值与特征向量的求解。如果从原始矩阵出发，先求出特征多项式，再求特征多项式的根，在理论上是无可非议的。但一般不用这种方法，因为了这种算法往往不稳定。常用的方法是迭代法或变换法。本章介绍求解特征值与特征向量的一些方法。
§1 乘幂法
乘幂法是通过求矩阵的特征向量来求特征值的一种迭代法，它适用于求矩阵的按模最大的特征值及对应的特征向量。
为什么是第j个分量呢？能相等吗？
定理8·1 设矩阵An×n有n个线性无关的特征向量Xi（i=1，2，…，n）,其对应的特征值λi(i =1,2，…,n)满足
|λ1|〉|λ2｜≧…≧｜λn｜
则对任何n维非零初始向量Z0,构造Zk = AZk—1（k=1，2,…） 有
(8·1)
其中（Zk)j表示向量Zk的第j个分量。
证明 : 只就λi是实数的情况证明如下.
理解：对比迭代
因为A有n个线性无关的特征向量Xi，（i = 1，2,…,n）所以任何非零向量Z0都可用Xi（i = 1，2，…,n)线性表示，即Z0=α1X1 + α2X2 +…+αnXn（α1≠0）
用A构造向量序列｛Zk}其中
(8.2）
由矩阵特征值定义知AXi=λiXi(i=1，2, …，n)，故
(8。3）
同理有
(8.4）
将（8。3)与（8。4）所得Zk及Zk—1的第j个分量相除，设α1≠0，并且注意到 |λi|<|λ1｜（i=1,2,…,n)得
证毕
定理8·1的证明过程实际上是给出了矩阵的按模最大特征值的计算方法：
1) 先任取一非零向量Z0,一般可取Z0=（1,1,1)T；
2) 按（8.2）式计算Zk=AZk-1（k=1，2,…）;
3) 当K足够大时，即可求出，为了减少λ1对于所选的第j个分量的依赖性,还可用各个分量比的平均值来代替，即

关于对应于λ1的特征向量的计算：
由（8。1）知，当k充分大时,Zk =λ1Zk—1,又由迭代式Zk = AZk-1，可知AZk—1=λ1Zk-1故由特征值定义知Zk—1即为λ1对应的特征向量,或Zk =λ1Zk—1为λ1对应的特征向量。
这种求矩阵的按模最大特征值及其对应特征向量的方法称为乘幂法。
应用乘幂法计算A的按模最大特征值λ1和对应特征向量时，由（8.3)易知

无穷模！
当|λ1|>1或|λ1｜<1时，Zk中不为零的分量将会随K的增大而无限增大,或随K的增大而趋于零，用计算机计算就会出现“上溢”或“下溢”。为了克服这个缺点,常将迭代向量Zk先规范化，然后再计算,具体做法是：
用max（Z)表示向量Zk的绝对值最大的分量，任取一初始向量Z0=α1X1 + α2X2 +…+αnXn（α1≠0）构造与（8.2）对应的向量序列.
（8。6)
由（8.3)可知
（8.7）
由（8.3)和（8.6）
（8。8）
也就是说,在满足定理的条件下，规范化的向量序列Yk仍收敛到A的按模最大特征值对应的特征向量；而向量序列Zk的绝对值最大的分量收敛到A的按模最大的特征值λ1.
例8·1 用规范化的乘幂法求矩阵

按模最大的特征值λ1和对应的特征向量X1.
解:取初始向量Z0=Y0=（1,1,1）T，按（8。6)、（8。7)和（8.8)算得Zk、Yk和max（Zk)，结果列于下表8-1.
表8—1
K
Zk
Yk
max(Zk）
0
1
2
3
4
5
6
7
1
274
44.42377
44。92333
44。99572
44.99959
44。99953
44。99953
1
95
14。8432
14。97623
14。99865
14.99988
14.99983
14.99983
1
—184
-29.64262
—29。95048
-29。99722
—29。99974
—29。99968
—29。99968
1
1
1
1
1
1
1
1
1
0。34672
0。33413
0.33337
0.33334
0。33333
0。33333
0。33333
1
-0。67153
-0.66727
-0.66670
—0.66667
-0。66667
—0。66667
-0.66667

44。42377
44.92333
44.99572
44。99953
44.99953
44。99953

经七次选代计算,λ1的近似值max（Z7）已稳定到小数点后第五位,故可取A的按模最大的特征值及对应的特征向量分别为
λ1=44。9995,X1=(1,0。333,—0。6667)T
我们不难求出矩阵A的三个特征值是
λ1=45，λ2=2，λ3=1
相应的特征向量为：
X1=（3，1,—2）T，X2=（3，2,-3)T，X3=(2,1，—2）T，
注：(1）若矩阵An×n的按模最大特征值λ1是P重根时,即
｜λ1｜=|λ2｜=…=|λp｜>|λp+1｜≥｜λn|
容易证明定理1的结论仍成立。
(2)此外,定理1中要求初始向量Z0的α1≠0是必要的，否则就不能得到对应于λ1的结果。如在例1中若取Z0=(1,1，-1）T，由此出发迭代便得
λ1=2,X1=（1，0.6667,—1）T
显然,这不是矩阵A的按模最大的特征值和对应的特征向量，出现这一现象，正是由于α1=0。事实上，由于A的特征向量X1，X2，X3是线性无关的,故Z0=(1，1，—1)T可表示为
Z0=α1X1+α2X2+α3X3
即

解之得

（3）乘幂法的收敛速度取决于比值 ｜λ1/λ1｜，当这个比值接近于1时，收敛很慢,反之收敛就比较快。例1是收敛较快的例子,如果收敛很慢,可以配合运用加速技术提高收敛速度。具体可参看西安交通大学出版社出版由邓建中等人编写的《计算方法》一书.
§2 反幂法
反幂法可以计算矩阵按模最小的特征值及对应的特征向量。
设An×n为非奇异矩阵，则A-1存在。若A的特征值λ1(）满足
｜λ1|≥|λ2｜≥…≥|λn｜>0
对应的特征向量为X1，X2， …,Xn.因为AXi=λiXi，所以A—1Xi=（1/λi)Xi ,即(1/λi)（i=1，2，…,n）是A-1的特征值，它满足

对应的特征向量仍是Xi（i=1，2，…,n).
这就是说,计算A的按模最小的特征值λn只要计算A-1按模最大的特征值
，从而，而求A-1的按模最大的特征值只须应用前述的乘幂法即可.。

所以反幂法的选代向量是： 设初始向量

于是

为避免求逆阵（）,由（)计算(）时，可以通过解线性方程组（)。
§3 QR方法
§1、§2 介绍了求矩阵A的部分特征值的方法,对于求它的全部特征值则有QR方法.
对矩阵A、B，若在非奇异矩阵P使得

则称矩阵A和B相似,记A(）B，而称P为化A为B的相似变换,并且由于（),
得知相似矩阵有相同的特征值，又因为
（）
有
（）
显然,若（）为B相应在于（）的特征向量,则（)为A的相应于（）的特征向量。
对于特殊的矩阵,例如上三角矩阵,其特征值即为主对角线上的元素，而任一非齐异矩阵与上三角矩阵的关系则有职下定理：
定理8·2设（）的特征值（）都为实数，那么必存在直交相似变换Q化A为上三角矩阵，即

由于（）,故也可以说A与R相似。特别当A为对称矩阵时，有
（）
这里的直交矩阵Q若能知道，即可求生物电A的特征值,但Q的求得并不那么容易,由此矩阵A的特征值也不可能直接求得。一般可由矩阵A通过直交相似变换构造矩阵列(），使其逐步逼近上三角矩阵R，从而求得矩阵A的满足精度要求的近似特征值及相应的特征向量.
定理8·3任一（）总可分解为一个直交矩阵Q和一个上三角矩阵R的乘积（）,若A非奇异,则这和分解是唯一的.
证明 对矩阵A，依（）左乘一系列初等旋转矩阵

(）
其中（）

当（）时.取（）;（）当(）时,则取（）.这里(）随A每次左乘（）而不断变化，而（)随之而变化,从而当（）时,
（)
当（）时有
（)
最后当（）时,有
（)
其中()的符号随()的符号而定,于是
(）
令（），显然Q为直交矩阵，故有
(）
现再证当A非奇异,则R,Q有逆矩阵存在，于是
（）
而（）为下交矩阵，（）为上三角矩阵,则要其相等，（）必为对角阵，又根据（）的直交性，便知（）为单位矩阵,即
（）
所以()
并且显然有()
以上证明实际上为我们提供了对A进行QR分解的具体方法。此外，A的QR分解也可通过(）直交化过程来实现。
既然任一非奇异矩阵A总有（),则令（），于是有（）

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

习题七

1． 用乘幂法或规范化乘幂法求下列矩阵按模最大的特征值及其对应的特征向量

—4 14 0
A = —5 13 0
—1 0 2
1）

2 —1 0
B = —1 2 —1
0 -1 2
2）

7 3 —2
C = 3 4 —1
-2 -1 3
3）

2 4 6
D = 3 9 15
4 16 36

4）

2． 用QR方法求下列矩阵

1 1 0
A = 1 1 1
0 1 1

1（）

3 2
B =
4 5

2(）

12 -20 41
C = 9 —15 —63
20 50 35

3（）

的全部特征值勤（精确到10-2)。

第九章 常微分方程的数值解法
本章讨论一阶常微分方程的初值问题
（）
（）
这类问题在工程计算中是常见的，例如,对于等截面均匀排风风道,风道内静压分布有如下规律：
（）
（)
我们知道,只要函数（）适当光滑，理论上就可以保证初值问题（9·1）—（9·2）的解（）存款额并且是唯一的。虽然求解常微分方程有各种各样的解析方法但解菥方法只能用来求解一些特殊类型的方程,大量从实际问题当中归结出来的微分方程主要靠数值解法。
所谓数值解法,就是寻求初值问题（）的解（)在一系列离散结点
（）
上的近似值
(）
相邻两个结点间(）称作步长,今后如不特别申明,总假定步长（)为定数
下面就介绍几种常邮的数值解法：
§1 欧拉（Euler）方法

初值问题（)的解,在几何上是通过点（)的一条曲线（)。欧拉法的求解过程是：先过点()作曲线的切线，该切线与直线（）相交于点（），再用（）作为曲线上点（）的纵坐标(）近似值。如图9—1所示。

（）

因为过（）点以（）为斜率的切线方程为
(）
当(）时得
（）
即取()，然后，再过（）点，以(）为斜率作直线
（）
当（)时得
（）
即取（）
一般地，如果已求出（）则过此点,以（)为斜率作直线
（）
当(）时得
取（）
通过上述过程，就可逐步求出点（）所对应的数值解（）
欧拉法的几何意义，是用一条折线近似代替曲线（)。欧拉（Euler）法（也叫欧拉折线法）的计算格式为
（)
（）
欧拉法是最古老的一种数值解法，它体现了数值方法的基本思想民，但精度很低,单独用它来作计算往往不能满足确度要求。
§2改进的欧拉方法
同一种计算格式往往可以通过多种途径构造出来，本节与下一节就会看到这一点.
为了提高精度，本节以改过的欧拉方法为例,介绍构造计算格式的数值积分方法。
交方程（9·1)的两端从（）到（）求积分,得到
（)
为要通过这个积分关系式得（)的近似值,只要近似地求出积分项（）即可。选择不同的近似方法计算这个积分项会得到不同的计算格式.
例如：用矩形分式计算积分项
（）
代入（9·4）得
()
若用（）代替上式中的（）并交右端的值作为（）的近似值（）。这样建立起来的格式就是欧拉法的计算格式（9·3).由于用矩形公式求积分值很粗糙,故导出欧拉格式精度也很低，不难证明,欧拉格式（9·3）的截断误差为(）
即()
为了改造精度，我们必用梯形法计算左端积分
（）
（）
将其中的（)分别用（)代替,则有下列计算格式
（）
(9·5）式被称为解常微分方程的梯形法则.
应该注意，格式（9·3)与(9·5）有本质上的区别，欧拉格式（9·3）是个直接的计算公式,这类格式称作显式的,而梯形法则（9·5）则由于其右端含有未知的（)故被除数称作是隐式的。它实际上是关于（），可以用选代法求解（参看第五章）,不过计算量比较大.我们将综合使用这两种格式，先用欧拉格式求得一个初步的近似值（)，称为预报值,然后用(）代替形法则右端的（)再直接计算。得到校正值(),这样建立起来的预报一校正系统称作改进的欧拉格式。
预报(）
校正(）
格式（9·6）的每一步需要两次调次调用函数(),它可以改写成下列形式；
（）
i = n ?
图9—2描述了改进的欧拉方法.
开 始

读入数据

1=〉i

X0+h=〉x1
y0+h(x0，y0)=>yp
y0+hf(x1,yp)=〉yc
(yp+yc)/2=>y1

印x1,y1

i+1=〉i
x1=〉x0
y1=>y0

i=n？
否

是
结 果

图9—2

欧拉法每一步只需对（）调用一次，而改进的欧拉法则不然，需对（)用两次,其计算量比欧拉法增加了一倍，付出这种代价的目的是为了提高精度。不难证明,改进的欧拉格式（9·6）的截断误差为（），即
（）
由此可见,它比欧拉格式的截断误差提高了一阶.
例9·1 解初值部题（）
取步长（）试求从（)到（）各结点上的数值解.
解 我们分别用两种格式进行计算,这里欧拉式的具体形式是
（）
而改进的欧拉格式是
（）
计算结果见表9—1
表9—1

结点
欧拉法
改进欧拉法
准确解
0
0.1
0。2
0.3
0.4
0.5
0。6
0。7
0.8
0.9
1.0
1
1。1
1.191818
1.277438
1.358213
1.435133
1。508966
1.580338
1。649783
1.717779
1.78477
1
1.095909
1.184097
1.266201
1.34336
1。416402
1。485956
1。552514
1.616475
1.678166
1.737867
1
1。095445
1.183216
1.264911
1。341641
1.414214
1.483240
1。549193
1。612452
1.673320
1。732051

这个问题有解析解（),按这个解算出的值列在表9—1第四列。同准确解比较,第二列欧拉格式的结果在致只有两位有效数字，而第三列改进的欧拉格式的结果则有三位有数字.
§3 龙格一库塔（Runge-—Kutta）方法
考察差商
（）
由微分中值定理,存在()使得
y（Xi+1）
于是，利用所给方程（）得到

y（Xi+1）= y（xi)+hf(Xi+Øh，y(xi+Øh)）
这里的f（Xi+Øh,y（xi+Øh））称作区间（Xi,Xi+1)上的平均斜率。由此得知，只要对平均斜率提供一种算法，由(9·7)式便相应地得到一种计算格式。
欧拉格式由于仅取Xi一个点的斜率值（Xi，yi)作为斜率的近似值，所以精度很低。
再考察改进的欧拉格式,它可改写为

yi+1 = yi +h/2 （Ki+K2)
其中 K1 =f（Xi, yi),K2 = f（Xi+1,yi+hK1）,
可见改进的欧拉格式可这样来理解,它用Xi与Xi+1两个点的斜率值K1和K2取算术平均作为平均斜率的近似值，而Xi+1处的斜率值勤K2,则通过已知信息yi来预报.
这个处理过程启发我们,若设法在(Xi，Xi+1)内多计算几个点的斜率值，然后将它们加权平均斜率的近似值,则有可能构造出更高精度的计算格式。这就是龙格一库塔方法的基本思想。
我们首先推广改进的欧拉格式。考察区间（Xi，Xi+1）内任一点