第八章 矩阵的特征值与特征向量的数值解法
某些工程计算涉及到矩阵的特征值与特征向量的求解。如果从原始矩阵出发,先求出特征多项式,再求特征多项式的根，在理论上是无可非议的。但一般不用这种方法，因为了这种算法往往不稳定．常用的方法是迭代法或变换法。本章介绍求解特征值与特征向量的一些方法。
§1 乘幂法
乘幂法是通过求矩阵的特征向量来求特征值的一种迭代法,它适用于求矩阵的按模最大的特征值及对应的特征向量。
定理８·1 设矩阵Ａn ×n 有n 个线性无关的特征向量X ｉ（ｉ=1,２,…,n)，其对应的特征值λi (ｉ ＝1,2，…，n)满足
｜λ１|＞|λ2｜≧…≧|λn ｜
则对任何ｎ维非零初始向量Z ０，构造Ｚk ＝ AZ k-1
11()lim
()k j k k j
Z Z λ→∞
-=
(8·1)
其中(Ｚk )ｊ表示向量Z k 的第ｊ个分量。 证明 : 只就λｉ是实数的情况证明如下。
因为A 有n 个线性无关的特征向量X i ,(i = １,2,…,n)用X ｉ(ｉ = 1,2,…，n)线性表示,即Z 0=α1X 1 + α2Ｘ2 +用A 构造向量序列{Z k }其中 ? 21021010,
,k k k Z AZ Z AZ A Z Z AZ A Z -=====， (8.2）
由矩阵特征值定义知ＡＸi =λi X i (i=1,２, …,n),故
?
0112211122211121k k k k k n n
k k k n n n
k
n
k
i i i i Z A Z A X A X A X X X X X X ααααλαλαλλλααλ===++
+=+++????
??=+
??????
?
∑
(８.3）
同理有
1
1
11
1121k n
k i k i i i Z X X λλααλ---=?
?
????=+ ?????
?
?
∑ (8．４)
将(8.３）与(8．4)所得Ｚk 及Z k-１的第j 个分量相除，设α1≠0，并且注意到 |λi |<|λ1|(i=1,２,…，n ）得
11()lim
()k j k k j
Z Z λ→∞
-= 证毕
定理8·1的证明过程实际上是给出了矩阵的按模最大特征值的计算方法： 1) 先任取一非零向量Z 0,一般可取Z 0=（１,１，１)T ; 2) 按(8.2）式计算Z k =ＡＺk-1(ｋ＝１,2,…)； 3) 当Ｋ足够大时，即可求出
11()()k j k j
Z Z λ-=,为了减少λ1对于所选的第j 个分
量的依赖性,还可用各个分量比的平均值来代替,即
1
11()()
n
k j
j
k j
Z Z n
λ=-=∑
关于对应于λ1的特征向量的计算:
由（8.1)知，当k 充分大时,Ｚk =λ1Z k-1，又由迭代式Z k = AZ k －1,可知AZ k-1=λ1Ｚk-1故由特征值定义知Z k-１即为λ1对应的特征向量，或Ｚk =λ1Z k-1为λ1对应的特征向量。
这种求矩阵的按模最大特征值及其对应特征向量的方法称为乘幂法。
应用乘幂法计算A 的按模最大特征值λ１和对应特征向量时,由(8．3）易知
11121k
n
k
i k
i i i Z X X λλααλ=????
??=+
???????∑ 当|λ1|＞１或|λ1｜＜１时，Z k
K 的增大而趋于零,用计算机计算就会出现常将迭代向量Z k 先规范化, 用ｍa ｘ（Z)表示向量Z k 011２Ｘ2 +…+αｎＸn (α1≠0）构造与(8.２）对应的向量序列。
()()()()()
()()()
10
10011022
020
2
122020200100max max max max max max max max k k k k k k k Z AZ Z AY AZ Z AZ A Z Z A Z Z AY AZ Z A Z A Z Z A Z Z AY AZ Z A Z -?====?
?
??====??
????====
??
，Y ，Y ，Y (8．6) 由(8.３）可知
()()
()()11210
1121max max max max k
n
i i i
k i k i k k
k n k i i i i i X X Z A Z X k Z X A Z X X λααλλααλ==??
+ ???==
=→→∞????
?
+ ? ?????
∑∑Y (8．７）
由(8.3）和（8.６）
()()
()()
()001
1001112111111211max max max max max max k k k k k
n k i i i i k n k k i i i i A Z A Z A Z A Z X X k X X Z λλααλλααλλλ--=--=?? ?== ???
???? ?+ ? ?????=→→∞????
?
+ ? ?????
∑∑max (8.8） 也就是说,在满足定理的条件下，规范化的向量序列Ｙk 仍收敛到Ａ的按模最大特征值对应的特征向量；而向量序列Z k 的绝对值最大的分量收敛到A 的按模最大的特征值λ１。
例8·1 用规范化的乘幂法求矩阵
13361354454688690??
??=??
??---??
A 按模最大的特征值λ1和对应的特征向量Ｘ1。
解：取初始向量Z 0＝Y ０=(1,１,１)T ,按(8．６)、(８.7)和(８．8)算得Z ｋ、Y k 和max(Z k )，结果列于下表8—１． 表
经七次选代计算,λ1的近似值m ａx(Ｚ7)已稳定到小数点后第五位，故可取A 的按模最大的特征值及对应的特征向量分别为
λ1=４4.９995,X 1＝(1,０．３３3，-0.666７）T
我们不难求出矩阵A 的三个特征值是
λ1=45，λ2＝2，λ3=1
相应的特征向量为:
X １=（3,１,-2)T ，Ｘ2=(３,２,－3)T ，Ｘ3=(２,1，-2）T ，
注：(1)若矩阵A n ×n 的按模最大特征值λ1是Ｐ重根时，即
|λ１|=|λ２|=…=｜λp |>|λp ＋１|≥|λn |
容易证明定理1的结论仍成立。 （２）此外,定理１中要求初始向量Z ０的α1≠0是必要的,否则就不能得到对应于λ1的结果。如在例1中若取Z 0=（1,1,－１)Ｔ，由此出发迭代便得
λ1=2,Ｘ1=(1,０.666７，－1)T
显然，这不是矩阵A 的按模最大的特征值和对应的特征向量,出现这一现象,正是由于α1=０。事实上，由于A 的特征向量Ｘ1,X ２,Ｘ3是线性无关的,故Z 0=（１，1,-1)Ｔ可表示为
Z 0=α１X 1+α２Ｘ2+α３X 3 即
?1231231233321212321ααααααααα?++=?
++=??
---=-? 解之得
?1230,1,1ααα===-
(3)乘幂法的收敛速度取决于比值 ｜λ１/λ１|，当这个比值接近于1时，收敛很慢，反之收敛就比较快。例1是收敛较快的例子,如果收敛很慢，可以配合运用加速技术提高收敛速度。具体可参看西安交通大学出版社出版由邓建中等人编写的《计算方法》一书。
§2 反幂法
反幂法可以计算矩阵按模最小的特征值及对应的特征向量。 设Ａｎ×n 为非奇异矩阵，则A -1存在。若A 的特征值λ1()满足
|λ1|≥|λ2|≥…≥|λn ｜>０
对应的特征向量为X 1,X 2， …，X n 。因为AX i =λi X ｉ，所以A -1X i =(1/λｉ)X i ，即(1/λi )(ｉ＝1,２,…,n )是A -１的特征值,它满足 ?
1
1
1
n
λλ≥≥
对应的特征向量仍是X i (i=１,2,…，ｎ)。
这就是说,计算A 的按模最小的特征值λn 只要计算A －1按模最大的特征值 ?1
n
λλ=
,从而1
n λλ
=
，而求A -1的按模最大的特征值只须应用前述的乘幂法即
可。.
所以反幂法的选代向量是: 设初始向量 于是
为避免求逆阵(），由()计算()时,可以通过解线性方程组()．
§3 QR 方法
§1、§2 介绍了求矩阵A 的部分特征值的方法,对于求它的全部特征值则有QR 方法．
对矩阵A 、Ｂ,若在非奇异矩阵P 使得
则称矩阵A 和B 相似，记Ａ（)B ，而称P 为化A 为B 的相似变换,并且由于(）, 得知相似矩阵有相同的特征值,又因为 （) 有 ()
显然，若(）为B 相应在于()的特征向量,则（)为A 的相应于()的特征向量. 对于特殊的矩阵,例如上三角矩阵,其特征值即为主对角线上的元素,而任一非齐异矩阵与上三角矩阵的关系则有职下定理：
定理8·2设（)的特征值()都为实数,那么必存在直交相似变换Q 化Ａ为上三角矩阵,即
由于（）,故也可以说A与R相似.特别当A为对称矩阵时，有
()
这里的直交矩阵Q若能知道，即可求生物电A的特征值,但Q的求得并不那么容易,由此矩阵A的特征值也不可能直接求得.一般可由矩阵Ａ通过直交相似变换构造矩阵列(）,使其逐步逼近上三角矩阵R，从而求得矩阵A的满足精度要求的近似特征值及相应的特征向量.
定理8·３任一()总可分解为一个直交矩阵Q和一个上三角矩阵Ｒ的乘积(),若Ａ非奇异，则这和分解是唯一的.
证明对矩阵A,依（）左乘一系列初等旋转矩阵
(）
其中()
当（)时.取（）；（）当（)时,则取（).这里(）随Ａ每次左乘(）而不断变化,而(）随之而变化,从而当（)时,
（）
当（）时有
（）
最后当()时,有
()
其中（）的符号随(）的符号而定，于是
（）
令（），显然Q为直交矩阵,故有
(）
现再证当A非奇异,则R，Q有逆矩阵存在,于是
（)
而（)为下交矩阵,（）为上三角矩阵,则要其相等，（）必为对角阵,又根据（）的直交性,便知（）为单位矩阵,即
(）
所以(）
并且显然有()
以上证明实际上为我们提供了对A进行ＱR分解的具体方法.此外，A的QＲ分解也可通过()直交化过程来实现.
既然任一非奇异矩阵A总有（）,则令()，于是有()
那么（)有（)于是（)与(）有相同的特征值.再交（)进行QR分解，有
（）
则()
并令()
有（)
于是（)与（）有相同的特征值.一般有
（）
()
令
有（）
()
于是()与（)有相同的特征值.可以证明,若非齐异实矩阵A有()个不同模的特征值,即
（）
则当()时,（）本质上收敛于上三角矩阵R（所谓本质上收敛于上三角矩阵是指矩阵列（）,收敛于一个上三角矩阵,而这个上三角矩阵除主对角元素外极限并不要求一定存在），R的主对角线元素即为所求的特征值.特别当Ａ为对称矩阵时,()收敛于对角矩阵D．具体计算中，当（)与()的主对角元素相差小于预先给定的业度时，则认为（)的主对角线元素即为A的特征值．对于QR分解,其有一个重要特点:当Ａ为对称带宽不变,即若A为三角矩阵,则(）仍为三对角矩阵．
习题七
1
1)
2)
3）
４）
2．用QR
1()
２()
3（)
的全部特征值勤(精确到10-2）.
第九章常微分方程的数值解法
本章讨论一阶常微分方程的初值问题
（)
（)
这类问题在工程计算中是常见的,例如,对于等截面均匀排风风道，风道内静压分布有如下规律:
（)
（)
我们知道,只要函数（)适当光滑，理论上就可以保证初值问题(9·1)—（９·２)的解（)存款额并且是唯一的.虽然求解常微分方程有各种各样的解析方法但解菥方法只能用来求解一些特殊类型的方程，大量从实际问题当中归结出来的微分方程主要靠数值解法.
所谓数值解法,就是寻求初值问题（)的解()在一系列离散结点
(）
上的近似值
（)
相邻两个结点间（）称作步长，今后如不特别申明,总假定步长（）为定数
下面就介绍几种常邮的数值解法：
§1 欧拉（Euｌｅｒ）方法
初值问题（）的解,在几何上是通过点（)的一条曲线（）．欧拉法的求解过程是：先过点（）作曲线的切线，该切线与直线()相交于点（),再用（）作为曲线上点(）的纵坐标（）近似值．如图９—1所示.
(）
因为过(）点以（)为斜率的切线方程为
(）
当（）时得
（)
即取（),然后，再过()点,以（)为斜率作直线
（）
当（)时得
()
即取(）
一般地,如果已求出(）则过此点，以(）为斜率作直线
(）
当（)时得
取()
通过上述过程,就可逐步求出点（)所对应的数值解(）
欧拉法的几何意义，是用一条折线近似代替曲线（)．欧拉（Ｅuler)法(也叫欧拉折线法）的计算格式为
(）
（）
欧拉法是最古老的一种数值解法，它体现了数值方法的基本思想民，但精度很低,单独用它来作计算往往不能满足确度要求.
§2改进的欧拉方法
同一种计算格式往往可以通过多种途径构造出来,本节与下一节就会看到这一点. 为了提高精度,本节以改过的欧拉方法为例,介绍构造计算格式的数值积分方法．交方程（9·1)的两端从（）到（)求积分,得到
(）
为要通过这个积分关系式得()的近似值,只要近似地求出积分项(）即可.选择不同的近似方法计算这个积分项会得到不同的计算格式.
例如:用矩形分式计算积分项
（）
代入（9·4)得
()
若用（）代替上式中的(）并交右端的值作为(）的近似值(）.这样建立起来的格式就是欧拉法的计算格式(９·3).由于用矩形公式求积分值很粗糙,故导出欧拉格式精度也很低,不难证明,欧拉格式(9·3)的截断误差为()
即（)
为了改造精度,我们必用梯形法计算左端积分
()
（)
将其中的（)分别用（)代替，则有下列计算格式
（)
(9·5)式被称为解常微分方程的梯形法则.
应该注意,格式(９·3）与(9·5)有本质上的区别,欧拉格式(９·3)是个直接的计算公式,这类格式称作显式的,而梯形法则（９·5)则由于其右端含有未知的()故被除数称作是隐式的.它实际上是关于（），可以用选代法求解(参看第五章)，不过计算量比较大．我们将综合使用这两种格式,先用欧拉格式求得一个初步的近似值（），称为预报值,然后用（）代替形法则右端的(）再直接计算．得到校正值(）,这样建立起来的预报一校正系统称作改进的欧拉格式.
预报(）
校正（）
格式（９·6）的每一步需要两次调次调用函数（），它可以改写成下列形式；(）
欧拉法每一步只需对(）调用一次,而改进的欧拉法则不然,需对（）用两次,其计算量比欧拉法增加了一倍,付出这种代价的目的是为了提高精度.不难证明,改进的欧拉格式（9·６)的截断误差为（)，即
（)
由此可见,它比欧拉格式的截断误差提高了一阶.
例9·１解初值部题（)
取步长(）试求从(）到()各结点上的数值解.
解我们分别用两种格式进行计算,这里欧拉式的具体形式是
(）
而改进的欧拉格式是
（)
计算结果见表9—1
表9—1
第9章矩阵特征值问题的数值 方法 9.1 特征值与特征向量 9.2 Hermite矩阵特征值问题 9.3 Jacobi方法 9.4 对分法 9.5 乘幂法 9.6 反幂法 9.7 QR方法
9.1 特征值与特征向量设A是n阶矩阵，x是非零列向量. 如果有数λ存在，满足， (1) 那么，称x是矩阵A关于特征值λ的特征向量.
如果把(1)式右端写为 ，那么(1)式又可写为: x λ ()0 I A x λ-=||0 I A λ-=即1110 ()||...n n n f I A a a a λλλλλ--=-=++++记 它是关于参数λ的n 次多项式，称为矩阵A 的特 征多项式， 其中a 0=(-1)n ｜A ｜. (2)
显然，当λ是A的一个特征值时，它必然 是的根. 反之，如果λ是的根，那么齐次方程组(2)有非零解向量x，使(1)式 成立. 从而，λ是A的一个特征值. A的特征值也称为A的特征根 . ()0 fλ= ()0 fλ=
矩阵特征值和特征向量有如下主要性质： 定理9.1.1 n阶矩阵A是降秩矩阵的充分必要 条件是A有零特征值. 定理9.1.2 设矩阵A与矩阵B相似，那么它们 有相同的特征值. 定理9.1.3 n阶矩阵A与A T有相同的特征值. 定理9.1.4 设λ ≠λj是n阶矩阵A的两个互异特 i 征值，x、y分别是其相应的右特征向 量和左特征向量，那么，x T y=0 .
9.2 Hermite矩阵特征值问题?设A为n阶矩阵，其共轭转置矩阵记为A H. 如果A=A H，那么，A称为Hermite矩阵.
求矩阵特征值算法及程序简介 1.幂法 1、幂法规范化算法 (1)输入矩阵A、初始向量( 0)，误差eps； (2) k 1； (3)计算V(k)A(k 1)； (4)m k max(V(k)) ,m k1max( V ( k 1))； (5) (k)V(k)/m k； (6)如果m k m k 1eps,则显示特征值1和对应的特征向量x(1) ),终止； (7)k k 1, 转(3) 注：如上算法中的符号max(V )表示取向量V 中绝对值最大的分量。本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。 2、规范化幂法程序 Clear[a,u,x]; a=Input[" 系数矩阵A="]; u=Input[" 初始迭代向量u(0)="]; n=Length[u]; eps=Input[" 误差精度eps ="]; nmax=Input[" 迭代允许最大次数nmax="]; fmax[x_]:=Module[{m=0,m1,m2}, Do[m1=Abs[x[[k]]]; If[m1>m,m2=x[[k]];m=m1], {k,1,Length[x]}]; m2] v=a.u; m0=fmax[u]; m1=fmax[v]; t=Abs[m1-m0]//N; k=0; While[t>eps&&k