第八章 矩阵的特征值与特征向量的数值解法某些工程计算涉及到矩阵的特征值与特征向量的求解。

如果从原始矩阵出发，先求出特征多项式，再求特征多项式的根，在理论上是无可非议的。

但一般不用这种方法，因为了这种算法往往不稳定.常用的方法是迭代法或变换法。

本章介绍求解特征值与特征向量的一些方法。

§1 乘幂法乘幂法是通过求矩阵的特征向量来求特征值的一种迭代法，它适用于求矩阵的按模最大的特征值及对应的特征向量。

定理8·1 设矩阵A n ×n 有n 个线性无关的特征向量X i （i=1,2,…,n ），其对应的特征值λi (i =1,2,…,n)满足|λ1|>|λ2|≧…≧|λn |则对任何n 维非零初始向量Z 0,构造Z k = AZ k-1(k=1,2,11()lim()k j k k jZ Z λ→∞-= （8·1）其中(Z k )j 表示向量Z k 的第j 个分量。

证明 ： 只就λi 是实数的情况证明如下。

因为A 有n 个线性无关的特征向量X i ，（i = 1,2,…,n 用X i （i = 1,2,…,n ）线性表示，即Z 0=α1X 1 + α2X 2 +…+用A 构造向量序列{Z k }其中 21021010,,k k k Z AZ Z AZ A Z Z AZ A Z -=====， (8.2)由矩阵特征值定义知AX i =λi X i (i=1,2, …,n)，故0112211122211121k k k k k n nk k k n n nknki i i i Z A Z A X A X A X X X X X X ααααλαλαλλλααλ===+++=+++⎡⎤⎛⎫⎢⎥=+⎪⎢⎥⎝⎭⎣⎦∑ （8.3）同理有11111121k nk i k i i i Z X X λλααλ---=⎡⎤⎛⎫⎢⎥=+ ⎪⎢⎥⎝⎭⎣⎦∑ （8.4） 将（8.3）与（8.4）所得Z k 及Z k-1的第j 个分量相除，设α1≠0，并且注意到 |λi |<|λ1|(i=1,2,…,n)得11()lim()k j k k jZ Z λ→∞-= 证毕定理8·1的证明过程实际上是给出了矩阵的按模最大特征值的计算方法：1) 先任取一非零向量Z 0,一般可取Z 0=(1,1,1)T ； 2) 按（8.2）式计算Z k =AZ k-1(k=1,2,…)； 3) 当K 足够大时,即可求出11()()k j k jZ Z λ-=,为了减少λ1对于所选的第j 个分量的依赖性,还可用各个分量比的平均值来代替,即111()()nk jjk jZ Z nλ=-=∑关于对应于λ1的特征向量的计算：由（8.1）知，当k 充分大时，Z k =λ1Z k-1，又由迭代式Z k = AZ k-1，可知AZ k-1=λ1Z k-1故由特征值定义知Z k-1即为λ1对应的特征向量，或Z k =λ1Z k-1为λ1对应的特征向量。

这种求矩阵的按模最大特征值及其对应特征向量的方法称为乘幂法。

应用乘幂法计算A 的按模最大特征值λ1和对应特征向量时,由（8.3）易知11121knki ki i i Z X X λλααλ=⎡⎤⎛⎫⎢⎥=+⎪⎢⎥⎝⎭⎣⎦∑ 当|λ1|>1或|λ1|<1时,Z k将迭代向量Z k 先规范化 用max(Z)表示向量Z k 01122+…+αn X n （α1≠0）构造与（8.2）对应的向量序列。

()()()()()()()()10100110220202122020200100max 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.3）可知()()()()112101121max max max max kni i i k i k ik kk nk i i i i i X X Z A Z X k Z X A Z X X λααλλααλ==⎛⎫+⎪⎝⎭===→→∞⎛⎫⎛⎫⎪+⎪ ⎪⎝⎭⎝⎭∑∑Y （8.7）由（8.3）和（8.6）()()()()()0011001112111111211max max max max max max k k k k kn 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） 也就是说,在满足定理的条件下,规范化的向量序列Y k 仍收敛到A 的按模最大特征值对应的特征向量；而向量序列Z k 的绝对值最大的分量收敛到A 的按模最大的特征值λ1。

例8·1 用规范化的乘幂法求矩阵13361354454688690⎡⎤⎢⎥=⎢⎥⎢⎥---⎣⎦A 按模最大的特征值λ1和对应的特征向量X 1。

解：取初始向量Z 0=Y 0=(1,1,1)T ，按(8.6)、(8.7)和(8.8)算得Z k 、Y k 和max(Z k )，结果列于下表8—1. 表经七次选代计算,λ1的近似值max(Z 7)已稳定到小数点后第五位，故可取A 的按模最大的特征值及对应的特征向量分别为λ1=44.9995，X 1=(1,0.333，-0.6667)T我们不难求出矩阵A 的三个特征值是λ1=45，λ2=2，λ3=1相应的特征向量为：X 1=(3,1,-2)T ，X 2=(3,2,-3)T ，X 3=(2,1,-2)T ，注：（1）若矩阵A n ×n 的按模最大特征值λ1是P 重根时,即|λ1|=|λ2|=…=|λp |>|λp+1|≥|λn |容易证明定理1的结论仍成立。

（2）此外，定理1中要求初始向量Z 0的α1≠0是必要的，否则就不能得到对应于λ1的结果。

如在例1中若取Z 0=(1,1,-1)T ，由此出发迭代便得λ1=2，X 1=(1,0.6667,-1)T显然,这不是矩阵A 的按模最大的特征值和对应的特征向量,出现这一现象,正是由于α1=0。

事实上，由于A 的特征向量X 1，X 2，X 3是线性无关的,故Z 0=(1,1,-1)T 可表示为Z 0=α1X 1+α2X 2+α3X 3 即1231231233321212321ααααααααα⎧++=⎪++=⎨⎪---=-⎩解之得1230,1,1ααα===-（3）乘幂法的收敛速度取决于比值 |λ1/λ1|，当这个比值接近于1时,收敛很慢,反之收敛就比较快。

例1是收敛较快的例子,如果收敛很慢,可以配合运用加速技术提高收敛速度。

具体可参看西安交通大学出版社出版由邓建中等人编写的《计算方法》一书。

§2 反幂法反幂法可以计算矩阵按模最小的特征值及对应的特征向量。

设A n ×n 为非奇异矩阵，则A -1存在。

若A 的特征值λ1（）满足|λ1|≥|λ2|≥…≥|λn |>0对应的特征向量为X 1,X 2, …,X n 。

因为AX i =λi X i ，所以A -1X i =(1/λi )X i ，即(1/λi )（i=1,2,…,n ）是A -1的特征值，它满足111nλλ≥≥对应的特征向量仍是X i （i=1,2,…,n ）。

这就是说,计算A 的按模最小的特征值λn 只要计算A -1按模最大的特征值1nλλ=，从而1n λλ=，而求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为三角矩阵,则（）仍为三对角矩阵.习题七11）2） 3）4）2．用QR1（）2（）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）的每一步需要两次调次调用函数（）,它可以改写成下列形式；（）图9—2—2欧拉法每一步只需对（）调用一次,而改进的欧拉法则不然,需对（）用两次,其计算量比欧拉法增加了一倍,付出这种代价的目的是为了提高精度.不难证明,改进的欧拉格式（9·6）的截断误差为（）,即（）由此可见,它比欧拉格式的截断误差提高了一阶.例9·1 解初值部题（）取步长（）试求从（）到（）各结点上的数值解.解我们分别用两种格式进行计算,这里欧拉式的具体形式是（）而改进的欧拉格式是（）计算结果见表9—1表9—1这个问题有解析解（）,按这个解算出的值列在表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一个点的斜率值(X i,y i)作为斜率的近似值,所以精度很低.再考察改进的欧拉格式,它可改写为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）内任一点。