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

个人收集整理 勿做商业用途第八章 矩阵的特征值与特征向量的数值解法 某些工程计算涉及到矩阵的特征值与特征向量的求解。如果从原始矩阵出发，先求出特征多项式，再求特征多项式的根，在理论上是无可非议的。但一般不用这种方法，因为了这种算法往往不稳定。常用的方法是迭代法或变换法。本章介绍求解特征值与特征向量的一些方法。1 乘幂法乘幂法是通过求矩阵的特征向量来求特征值的一种迭代法，它适用于求矩阵的按模最大的特征值及对应的特征向量。 为什么是第j个分量呢？能相等吗？定理81 设矩阵Ann有n个线性无关的特征向量Xi（i=1，2，n）,其对应的特征值i(i =1,2，,n)满足|1|2n则对任何n维非零初始向量Z0,构造Zk = AZk1（k=1，2,） 有 (81)其中（Zk)j表示向量Zk的第j个分量。证明 : 只就i是实数的情况证明如下.理解：对比迭代因为A有n个线性无关的特征向量Xi，（i = 1，2,n）所以任何非零向量Z0都可用Xi（i = 1，2，,n)线性表示，即Z0=1X1 + 2X2 +nXn（10） 用A构造向量序列Zk其中 (8.2）由矩阵特征值定义知AXi=iXi(i=1，2, ，n)，故 (8。3）同理有 (8.4）将（8。3)与（8。4）所得Zk及Zk1的第j个分量相除，设10，并且注意到 |i|1或|1|p+1n|容易证明定理1的结论仍成立。 (2)此外,定理1中要求初始向量Z0的10是必要的，否则就不能得到对应于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 反幂法反幂法可以计算矩阵按模最小的特征值及对应的特征向量。设Ann为非奇异矩阵，则A-1存在。若A的特征值1(）满足1|2|n0对应的特征向量为X1，X2， ,Xn.因为AXi=iXi，所以A1Xi=（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的相应于（）的特征向量。 对于特殊的矩阵,例如上三角矩阵,其特征值即为主对角线上的元素，而任一非齐异矩阵与上三角矩阵的关系则有职下定理：定理82设（）的特征值（）都为实数，那么必存在直交相似变换Q化A为上三角矩阵，即由于（）,故也可以说A与R相似。特别当A为对称矩阵时，有（）这里的直交矩阵Q若能知道，即可求生物电A的特征值,但Q的求得并不那么容易,由此矩阵A的特征值也不可能直接求得。一般可由矩阵A通过直交相似变换构造矩阵列(），使其逐步逼近上三角矩阵R，从而求得矩阵A的满足精度要求的近似特征值及相应的特征向量.定理83任一（）总可分解为一个直交矩阵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 0A = 5 13 0 1 0 21） 2 1 0B = 1 2 1 0 -1 22） 7 3 2C = 3 4 1 -2 -1 33） 2 4 6D = 3 9 15 4 16 364）2 用QR方法求下列矩阵 1 1 0A = 1 1 1 0 1 11（） 3 2B = 4 52(） 12 -20 41C = 9 15 63 20 50 353（）的全部特征值勤（精确到10-2)。第九章 常微分方程的数值解法本章讨论一阶常微分方程的初值问题（）（）这类问题在工程计算中是常见的，例如,对于等截面均匀排风风道,风道内静压分布有如下规律：（）（)我们知道,只要函数（）适当光滑，理论上就可以保证初值问题（91）（92）的解（）存款额并且是唯一的。虽然求解常微分方程有各种各样的解析方法但解菥方法只能用来求解一些特殊类型的方程,大量从实际问题当中归结出来的微分方程主要靠数值解法。 所谓数值解法,就是寻求初值问题（）的解（)在一系列离散结点（）上的近似值(）相邻两个结点间(）称作步长,今后如不特别申明,总假定步长（)为定数下面就介绍几种常邮的数值解法：1 欧拉（Euler）方法初值问题（)的解,在几何上是通过点（)的一条曲线（)。欧拉法的求解过程是：先过点()作曲线的切线，该切线与直线（）相交于点（），再用（）作为曲线上点（）的纵坐标(）近似值。如图91所示。（）因为过（）点以（）为斜率的切线方程为(）当(）时得（）即取()，然后，再过（）点，以(）为斜率作直线（）当（)时得（）即取（）一般地，如果已求出（）则过此点,以（)为斜率作直线（）当(）时得取（）通过上述过程，就可逐步求出点（）所对应的数值解（）欧拉法的几何意义，是用一条折线近似代替曲线（)。欧拉（Euler）法（也叫欧拉折线法）的计算格式为（)（）欧拉法是最古老的一种数值解法，它体现了数值方法的基本思想民，但精度很低,单独用它来作计算往往不能满足确度要求。2改进的欧拉方法同一种计算格式往往可以通过多种途径构造出来，本节与下一节就会看到这一点.为了提高精度，本节以改过的欧拉方法为例,介绍构造计算格式的数值积分方法。 交方程（91)的两端从（）到（）求积分,得到（)为要通过这个积分关系式得（)的近似值,只要近似地求出积分项（）即可。选择不同的近似方法计算这个积分项会得到不同的计算格式. 例如：用矩形分式计算积分项（）代入（94）得()若用（）代替上式中的（）并交右端的值作为（）的近似值（）。这样建立起来的格式就是欧拉法的计算格式（93).由于用矩形公式求积分值很粗糙,故导出欧拉格式精度也很低，不难证明,欧拉格式（93）的截断误差为(）即()为了改造精度，我们必用梯形法计算左端积分（）（）将其中的（)分别用（)代替,则有下列计算格式（）(95）式被称为解常微分方程的梯形法则.应该注意，格式（93)与(95）有本质上的区别，欧拉格式（93）是个直接的计算公式,这类格式称作显式的,而梯形法则（95）则由于其右端含有未知的（)故被除数称作是隐式的。它实际上是关于（），可以用选代法求解（参看第五章）,不过计算量比较大.我们将综合使用这两种格式，先用欧拉格式求得一个初步的近似值（)，称为预报值,然后用(）代替形法则右端的（)再直接计算。得到校正值(),这样建立起来的预报一校正系统称作改进的欧拉格式。预报(）校正(）格式（96）的每一步需要两次调次调用函数(),它可以改写成下列形式；（） i = n ?图92描述了改进的欧拉方法.开 始读入数据 1=iX0+h=x1y0+h(x0，y0)=ypy0+hf(x1,yp)=yc(yp+yc)/2=y1印x1,y1i+1=ix1=x0y1=y0i=n？ 否 是结 果 图92欧拉法每一步只需对（）调用一次，而改进的欧拉法则不然，需对（)用两次,其计算量比欧拉法增加了一倍，付出这种代价的目的是为了提高精度。不难证明,改进的欧拉格式（96）的截断误差为（），即（）由此可见,它比欧拉格式的截断误差提高了一阶.例91 解初值部题（）取步长（）试求从（)到（）各结点上的数值解.解 我们分别用两种格式进行计算,这里欧拉式的具体形式是（）而改进的欧拉格式是（）计算结果见表91表91结点欧拉法改进欧拉法准确解00.10。20.30.40.50。60。70.80.91.011。11.1918181.2774381.3582131.4351331。5089661.5803381。6497831.7177791.7847711.0959091.1840971.2662011.343361。4164021。4859561。5525141.6164751.6781661.73786711。0954451.1832161.2649111。3416411.4142141.4832401。5491931。6124521.6733201。732051这个问题有解析解（),按这个解算出的值列在表91第四列。同准确解比较,第二列欧拉格式的结果在致只有两位有效数字，而第三列改进的欧拉格式的结果则有三位有数字.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)上的平均斜率。由此得知，只要对平均斜率提供一种算法，由(97)式便相应地得到一种计算格式。欧拉格式由于仅取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）内任一点