1、7.3.2 二阶龙格-库塔方法,原来改进欧拉公式k2为xi+1点的近似平均值, 现对k2重新设置取,内任一点,使用xi和xi+l两点斜率值k1,k2加权平均作为, 即,其中,为待定系数,k1仍为x1上的斜率值，即,问题是如何预测xi+l处的斜率k2.,类似于改进的欧拉方法，先利用欧拉公式得到斜率值,所以计算公式如下:,其中,为调节参数, 使其具有2阶精度,考虑截断误差,对,做二阶泰勒展开,再使用yi+l, 通过f(x,y)形式, 得到斜率值,对,作为二元函数1阶泰勒展开,其中,，余项是关于,的高阶无穷小.,对节点,即,，此时,对,利用复合函数求导：对,其中,有,可对式关于x求导：,将和带入整理
2、得：,要使,和,为一族系数与l有关,当 时,时,即为改进欧拉公式 ;,当,该公式为变形欧拉公式, 即,特殊地,在区间,上取两点，预测相应斜率值,对,作加权得平均斜率近似值,对截断误差公式，要达到,，则,的系数须为0, 从而建立方程组求系数.,总结二阶龙格-库塔构造方程:,7.3.3 高阶龙格-库塔公式：,为进一步提高精度（如三阶精度），在xi, xi+l上除xi和xi+l,外,可再增加一点,以xi, xi+l, xi+m三点处的斜率值k1,k2,k3加权平均,计算公式：,其中k1,k2仍为前面定义，即,需要求K3, 可在,上利用二阶龙格-库塔公式,得到预测值后，再利用f(x,y)得到斜率值,整
3、个计算公式为：,需求,使其具有3阶精度.,可利用截断误差公式求解，具体略, 可得条件:,该方程组得到的是一族解， 常用的三阶龙格库塔公式有：,同样地，精度还可以提高到4阶， 具体略。,在计算yi+1之前，已经求出y0、y1、yi. 是否可以用前面已算好 的yi，得到高精度的yi+1的值？,7.4线性多步法,该方法为多步法，（前面方法为单步法, 只涉及到yi、yi+1）,线性r步计算公式可写成：,yi+1=yi+hf(xi,yi),yi+1=yi+,f(xi,y(xi)+f(xi+1,y(xi+1)，,其中,例如欧拉公式：,其中,梯形公式：,式中，j, j为常数。 当-1=0时为显式公式, -1
4、0时为隐式公式.,Ri+1=y(xi+1)-,常微分方程数值解常用方法： 利用数值积分方法（另外还有微分中值定理）,y(xi+1)=y(xi)+,对,提高精度（原来只使用一点、两点插值）.,7.4.1阿当姆斯内插公式（Admas）,使用xi,xi+1作插值节点，节点少，精度不高， 可增加节点插值（但节点数不能太多）,误差余项:, 采用更好的节点插值方法构造，,需增加节点, 若先取xi,xi+1内部点，函数值未知（可通过龙格-库塔方法 计算，但相对复杂）可选取xi,xi+1外部的已知节点.,阿当姆斯内插公式： 选取节点为xi-2, xi-1, xi, xi+1来计算y(xi+1)。,构造三次插值
5、多项式代替,，近似计算,被积函数表示为,(xi-2,xi+1),（利用积分第一中值定理: xxi,xi+1内,不变号）,多项式函数可求出积分值（可通过换元法） 前两项可作为近似计算公式，最后一项作为误差余项。,具体化简为：,内插公式为:,是关于yi+1的隐式公式，截断误差为4阶.,需计算xi,xi-1,xi-2三点函数值，初值点y0,y1,y2称为三步法。 （可通过欧拉方法或龙格-库塔方法计算）,7.4.2 阿当姆斯外推公式,将隐式转为显式公式，即将节点xi+1转为xi-3,取xi-3，xi-2，xi-1，xi为插值节点，求y(xi+1).,采用插值多项式，类似地得到,其误差为:,外推公式：,
6、该公式是显式公式，截断误差为4阶.,需要节点值,即初值需,，可通过欧拉,方法或龙格库塔方法求得，称为四步法。,比较内插和外推公式：,内插公式：隐式（计算复杂），但误差相对较小（余项系数较小），且计算公式（三步法）中的4个系数相对较小，产生的误差影响小。,外推公式：显式（计算相对简单），但误差相对较大（余项系数较大），且计算（四步法）公式中的4个系数相对较大，误差影响大。,可将两公式联合使用，构成阿当姆斯校正公式： 即先通过外推公式求得,预测值,代入到内插公式中：,再将,7.5 一阶方程组和高阶方程7.5.1 一阶方程组,原来一阶方程初值问题,将y, f和 看成向量, 可应用到一阶方程组，使用下
7、列标记：,其中u和v是关于x的两个函数,和,是关于, 的两个函数，即已知,求函数u和v.,该问题可转为,为向量.,这里要求每个向量相等，可利用前面方法求u和v的数值解.,7.5.2 高阶方程,有些高阶方程可将其转化为一阶方程组：,例如二阶方程初值问题：,可引入新变量,，并将其转为一阶方程组：,再使用一阶方程组方法求解。,7.5.3 二阶线性常微分方程通解,可利用差分方程求解，另外还需增加边界条件,方程形式:,其中,为关于x的已知函数.,（1）,将a, b作n等分, 即步长为,各节点为,数值解即为求出这些节点xi上y(xi)的近似值.,首先对方程（1）在节点xi上离散化:,对xi用一阶差商和二阶
8、差商分别近似该点一阶导数和二阶导数值:,令,所以（1）式即为差分方程：,具体求解方程：,(2),所以构成方程组，有n-1个方程, 求n+1个未知数.,需补充两个方程，可由不同的边界条件来补充，如,（1）,（2）,这样可通过求解方程组方法求出,该方程由为三对角方程组，所以可使用追赶法快速求解。,第八章 矩阵的特征值及特征向量计算,对n阶方阵A=,其特征值为方程,的根, 可写成:,其中ci为系数，与矩阵A有关.,称为矩阵A的特征方程，有n个根（包括重根、复数根）,当,是A的特征值时，相应方程组,的非零解x称为矩阵A关于,的特征向量.,对特征值和特征向量的计算，可看成是方程求根和线性方程组求解问题（
9、理论上可以求出）.但根据向量和矩阵及线性代数中的知识，可以使用一些别的数值方法，如对模最大或最小特征值的求解;一般特征值的统一求解方法.,8.1问题提出,称为特征值（特征根）.,8.2 模最大特征值的求解幂法,矩阵A的特征值中，模最大的特征值称为主特征值，其模长 称为谱半径（实数中，模长为绝对值，对于复数,模长为,）,思路：,设n 阶方阵A有n个线性无关的特征向量 (前提条件),所对应的特征值,并按模的大小排列, 即,分以下两种情况讨论:,（1）,任取初始向量v0,矩阵A有n个线性无关的特征向量,任何一个n维向量都可以由它们线性表示，即,(1),现从v0出发作一系列迭代：,可得一向量序列，,并
10、利用（1）可得:,分别为,(线性代数知:不同特征值的特征向量线性无关,相同特征值的特征向量不一定线性无关),k=0,1,2,同理有：,i=2,3,n,i=2,3,n,(2),若（2）式迭代收敛，则,近似线性相关，,即为模最大的特征值.,其系数,为向量,其中,表示向量,的第j个分量,（1）它是迭代法，收敛速度主要取决于,当,越小，收敛越快。,（2）迭代过程中，主要计算,所以称为幂法。,另外防止计量出的新向量,中的分量值太大或过小，可使用,（3）初始向量对迭代过程有影响，因此通常采用先算下去再说的方法，及时调整初值进行计算。,（4）另一种情况：,类似地, 有,该方法的特点：,(即每次迭代时将Vk的
11、最大分量化为1).,规范化的幂法求解，这里略,(i=3，4，n),(i=3，4，n),同理：,将上面三式带入,可知：,上式说明：,三向量大体上线性相关。,令,则有,时 , p, q趋向于某固定值,，可利用,求出,特殊地，当,时，,可近似看成,设矩阵A为非奇异矩阵 （前提条件, 即零不是A的特征值）,并设A的特征值有,因为A为非奇异矩阵, 所以,存在（逆矩阵）,且由,可得,即说明矩阵,的特征值为,（j=1,2,n）, 并有,所以可对矩阵,，利用幂法计算模最大特征值，对其求倒数，,8.3 模最小特征的求解反幂法,即得A的模最小特征值（反幂法）,求矩阵特征值和特征向量的一般方法 对特征值重根，复根，
12、特征值可能为0的情况都适用。,若是奇异矩阵，则有零特征值,任取一个不为A特征值的数,，则有,因为矩阵,是非奇异矩阵，只要得到,的特征值和特征向量就可方便求出A的特征值和特征向量。,因此所求特征值,与原来的特征值相差一个 .,不失一般性，可设A为非奇异矩阵，对非奇异矩阵A，,可以将其分解成一个正交矩阵Q (即A-1 = AT，有AAT = ATA = E)和一个上三角矩阵R的乘积, 即,A1 = A = Q1R1，(可得R1 = Q-1A1),交换该乘积顺序得,8.4 QR分解法,思路：对矩阵A,A2 = R1Q1 = Q1-1A1Q1,A1与A2有相同的特征值.,继续迭代下去，可得一迭代序列Ak，,当Q为正交矩阵时， 可以证明：当k时，矩阵序列Ak本质上收敛于R*,R主对角线上的元素可看成是A的全部特征值. 特别地，当A为对称矩阵时，Ak收敛于对角阵.,另外，矩阵A的QR分解是唯一的 (可以利用施密特正交化方法计算，具体略).,即 A2与A1是相似矩阵,(即指R*的主对角元素趋于常数，其他元素不一定，证明略),A3 = Q2-1A2Q2 ,Ak = QkRk，Ak+1 = RkQk,如A2 = Q2R2 （利用QR分解），,
