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

第八章矩阵的特征值与特征向量的数值解法某些工程计算涉及到矩阵的特征值与特征向量的求解。如果从原始矩阵出发，先求出特征多项式，再求特征多项式的根，在理论上是无可非议的。但一般不用这种方法，因为了这种算法往往不稳定.常用的方法是迭代法或变换法。本章介绍求解特征值与特征向量的一些方法。§1乘幂法乘幂法是通过求矩阵的特征向量来求特征值的一种迭代法，它适用于求矩阵的按模最大的特征值及对应的特征向量。定理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,…)有11()lim()kjkkjZZ（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}其中21021010,,kkkZAZZAZAZZAZAZ，(8.2)由矩阵特征值定义知AXi=λiXi(i=1,2,…,n)，故0112211122211121kkkkknnkkknnnknkiiiiZAZAXAXAXXXXXX（8.3）同理有11111121knkikiiiZXX（8.4）将（8.3）与（8.4）所得Zk及Zk-1的第j个分量相除，设α1≠0，并且注意到|λi||λ1|(i=1,2,…,n)得11()lim()kjkkjZZ证毕定理8·1的证明过程实际上是给出了矩阵的按模最大特征值的计算方法：理解：对比迭代为什么是第j个分量呢？能相等吗？1)先任取一非零向量Z0,一般可取Z0=(1,1,1)T；2)按（8.2）式计算Zk=AZk-1(k=1,2,…)；3)当K足够大时,即可求出11()()kjkjZZ,为了减少λ1对于所选的第j个分量的依赖性,还可用各个分量比的平均值来代替,即111()()nkjjkjZZn关于对应于λ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）易知11121knkikiiiZXX当|λ1|1或|λ1|1时,Zk中不为零的分量将会随K的增大而无限增大,或随K的增大而趋于零，用计算机计算就会出现“上溢”或“下溢”。为了克服这个缺点,常将迭代向量Zk先规范化，然后再计算，具体做法是：用max(Z)表示向量Zk的绝对值最大的分量，任取一初始向量Z0=α1X1+α2X2+…+αnXn（α1≠0）构造与（8.2）对应的向量序列。10100110220202122020200100maxmaxmaxmaxmaxmaxmaxmaxkkkkkkkZAZZAYAZZAZAZZAZZAYAZZAZAZZAZZAYAZZAZ，Y，Y，Y（8.6）由（8.3）可知1121001121maxmaxmaxmaxkniiikikikkknkiiiiiXXZAZXkZXAZXXY（8.7）无穷模！由（8.3）和（8.6）0011001112111111211maxmaxmaxmaxmaxmaxkkkkknkiiiiknkkiiiiAZAZAZAZXXkXXZmax（8.8）也就是说,在满足定理的条件下,规范化的向量序列Yk仍收敛到A的按模最大特征值对应的特征向量；而向量序列Zk的绝对值最大的分量收敛到A的按模最大的特征值λ1。例8·1用规范化的乘幂法求矩阵13361354454688690A按模最大的特征值λ1和对应的特征向量X1。解：取初始向量Z0=Y0=(1,1,1)T，按(8.6)、(8.7)和(8.8)算得Zk、Yk和max(Zk)，结果列于下表8—1.表8—1KZkYkmax(Zk)01234567127444.4237744.9233344.9957244.9995944.9995344.9995319514.843214.9762314.9986514.9998814.9998314.999831-184-29.64262-29.95048-29.99722-29.99974-29.99968-29.999681111111110.346720.334130.333370.333340.333330.333330.333331-0.67153-0.66727-0.66670-0.66667-0.66667-0.66667-0.6666744.4237744.9233344.9957244.9995344.9995344.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即1231231233321212321解之得1230,1,1（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的特征值，它满足111n对应的特征向量仍是Xi（i=1,2,…,n）。这就是说,计算A的按模最小的特征值λn只要计算A-1按模最大的特征值1n，从而1n，而求A-1的按模最大的特征值只须应用前述的乘幂法即可。.所以反幂法的选代向量是：设初始向量于是为避免求逆阵（）,由（）计算（）时,可以通过解线性方程组（）.§3QR方法§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．用乘幂法或规范化乘幂法求下列矩阵按模最大的特征值及其对应的特征向量1）-4140A=-5130-1022）3）4）2．用QR方法求下列矩阵1（）2（）3（）的全部特征值勤（精确到10-2）.第九章常微分方程的数值解法本章讨论一阶常微分方程的初值问题（）（）这类问题在工程计算中是常见的,例如,对于等截面均匀排风风道,风道内静压分布有如下规律：（）2-10B=-12-10-1273-2C=34-1-2-13246D=391541636110A=11101132B=4512-2041C=9-15-63205035（）我们知道,只要函数（）适当光滑,理论上就可以保证初值问题（9·1）—（9·2）的解（）存款额并且是唯一的.虽然求解常微分方程有各种各样的解析方法但解菥方法只能用来求解一些特殊类型的方程,大量从实际问题当中归结出来的微分方程主要靠数值解法.所谓数值解法,就是寻求初值问题（）的解（）在一系列离散结点（）上的近似值（）相邻两个结点间（）称作步长,今后如不特别申明,总假定步长（）为定数下面就介绍几种常邮的数值解法：§1欧拉（Euler）方法初值问题（）的解,在几何上是通过点（）的一条曲线（）.欧拉法的求解过程是：先过点（）作曲线的切线,该切线与直线（）相交于点（）,再用（）作为曲线上点（）的纵坐标（）近似值.如图9—1所示.（）因为过（）点以（）为斜率的切线方程为（）当（）时得（）即取（）,然后,再过（）点,