微机原理第7章

1、常微分方程初值问题的数值解法常微分方程初值问题的数值解法7.1 7.1 欧拉法及其改进欧拉法及其改进7.2 7.2 龙格龙格- -库塔方法库塔方法7.3 7.3 线性多步法线性多步法第七章第七章7.0 7.0 概述概述7.4 7.4 数值算例数值算例本章着重讨论一阶常微分方程初值问题本章着重讨论一阶常微分方程初值问题 00yxyyxfdxdy）（），（bxa的数值解法。的数值解法。7.0 7.0 概述概述 常微分方程初值问题的数值解是求上述初常微分方程初值问题的数值解是求上述初值问题的解值问题的解y y( (x x) )在区间在区间 a a，b b 中的点中的点列列），（nihxxiii211

2、上的近似值上的近似值 . .以下设以下设iyih不变，记为不变，记为h-h-步长步长 。)(nxy 假设解假设解y y( (x x) )在区间在区间 a a, ,b b 上是存在而且唯一的，上是存在而且唯一的，并且具有充分的光滑度，因此，要求并且具有充分的光滑度，因此，要求f f( (x,yx,y) )也充也充分光滑。初值问题的解析解分光滑。初值问题的解析解( (理论解）用理论解）用 表示表示, ,数值解法的精确解用数值解法的精确解用 表示。表示。ny常微分方程数值解法一般分为：常微分方程数值解法一般分为： （1）一步法：在计算 时，只用到 , 和 ，即前一步的值。 nx1nyny1nx（2）

3、多步法：计算 时，除用到 , 和 以外，还要用 和 ,即前 k步的值。 1nypny1nxnxnypnx)0;2 , 1(kkp（3）显式格式与隐式格式。7.1 7.1 欧拉法与梯形法欧拉法与梯形法设节点为设节点为 ，得欧拉方法计算公得欧拉方法计算公式为：式为： ），（3 , 2100nnhxxn），）（，（32101nyxhfyynnnn一、欧拉一、欧拉(Euler)(Euler)法法下面通过几种常用的方法来推导该公式。下面通过几种常用的方法来推导该公式。1 1、泰勒展开法、泰勒展开法假设在假设在 附近把附近把y y（x x）做做TaylorTaylor展开，有：展开，有：nx)( 2)(

4、)()(2nnnnxyxhyxyhxyh取取h h的线性部分的线性部分, ,并用并用 表示表示 的近似值的近似值, ,得得ny)(nxy), 2 , 1 , 0)(,(1nyxhfyynnnn2 2、数值积分法、数值积分法 从 到 +h对等式 y(t )=f(t,y(t) 进行积分得到 nxnx1)(,()()(nnxxnndttytfxyhxy再利用左矩形公式，得再利用左矩形公式，得)(,()()(nnnnxyxhfxyhxy从而得到从而得到EulerEuler公式。公式。),()(nnnyxfxy由由)(,()(1)()(1nnnnnxxxyxfxyxxyynn3 3、数值微分法、数值微分

5、法4 4、几何方法、几何方法过点过点( (x xn n, ,y yn n) )作以作以f f( (x xn n , ,y yn n) )为斜率的直线方程：为斜率的直线方程：)(,(nnnnxxyxfyy 将将x x= =x xn+n+1 1处该直线上的函数值做为处该直线上的函数值做为y y( (x xn n+1+1) )的近似的近似值，则有值，则有EulerEuler公式。这实质上是在每个小区间上利公式。这实质上是在每个小区间上利用折线来代替曲线的结果，故用折线来代替曲线的结果，故EulerEuler法法又称又称EulerEuler折线折线法法。二、梯形法二、梯形法在式在式 中中, ,将积分用

6、将积分用梯形公式来代替，则有梯形公式来代替，则有 1)(1nnxxnndttytfxyxy）（，（）),()(,( 12)(,()(,(2)()(13111nnnnnnnnxxyfhxyxfxyxfhxyxy从而得到梯形公式：从而得到梯形公式：, 2 , 1 , 0),(),(2111nyxfyxfhyynnnnnn 梯形方法关于梯形方法关于y yn n+1+1是隐式的，而是隐式的，而EulerEuler方法是显方法是显式的。一般情形下不容易从上式解出式的。一般情形下不容易从上式解出y yn n+1+1，因而可将，因而可将上式与上式与EulerEuler公式联合使用，即公式联合使用，即), 2

7、 , 1 , 0;, 2 , 1 , 0)(,(),(2),()(11)1(1)0(1 nkyxfyxfhyyyxhfyyknnnnnknnnnn使用上式时，先用第一式算出使用上式时，先用第一式算出x xn n+1+1处处y yn n+1+1的初始近似的初始近似)0(1ny再用第二式反复迭代，得到数列再用第二式反复迭代，得到数列0)1(1kkny用用)(1)1(1knknyy来控制迭代次数，这里来控制迭代次数，这里 为允许误差。把满足误差要求为允许误差。把满足误差要求的的)1(1kny可以证明可以证明, ,当当f f( (x,yx,y) )满足满足LipschitzLipschitz条件条件,

8、 ,即：即：12,LhLyf且( (L L为为LipschitzLipschitz常数常数) )时时, ,上述数列收敛。上述数列收敛。作为作为y y( (x xn n+1+1) )的近似值的近似值y yn n+1+1. .类似地可以得出类似地可以得出y yn+2n+2,y,yn+3n+3,证明：证明：由由),(),(2111nnnnnnyxfyxfhyy和和),(),(2)(11 11knnnnnknyxfyxfhyy有：有：),(),(2111)(111)1(1nnknnnknyxfyxfhyy)(,(211)(11nknnyyxyfh)(211)(11)(1之间与介于nknnknyyyyh

9、L反复使用不等式有：反复使用不等式有：1)1(1nknyy1)(121nknyyhL0211)0(11nnkyyhL ）（ 实用中实用中, ,在在h h 取得较小时取得较小时, ,用梯形公式计算用梯形公式计算, ,第第二式只迭代一次就结束二式只迭代一次就结束, ,得到得到EulerEuler预估预估- -校正格式校正格式: : ), 2 , 1 , 0)(,(),(2),()0(111)0(1nyxfyxfhyyyxhfyynnnnnnnnnn第一式称为第一式称为预估公式预估公式，第二式称为，第二式称为校正公式校正公式。 三、三、EulerEuler预估预估- -校正格式校正格式四、方法的误差

10、估计、收敛性和稳定性四、方法的误差估计、收敛性和稳定性 定义1： 为 某一数值方法在xn处的整体截断误差(不考虑舍入误差的影响)。 定义2：对单步法，在 的假设下， 称为在 处的局部截断误差。 nnnyxy）（nx）（nnxyy 11nnnyxyR）（nxRemark1: Euler法的局部截断误差为法的局部截断误差为:,)(2)(1211 nnnnnnnxxyhyxyRRemark2: 梯形方法的局部截断误差为:,)(12)(1311 nnnnnnnxxyhyxyR用泰勒展开法推导用泰勒展开法推导EulerEuler预估校正预估校正格式的局部截断误差格式的局部截断误差改写改写EulerEul

11、er预估校正公式为：预估校正公式为： ),(),(2121121211kyhxhfkyxhfkkkyynnnnnn)( )(,(),(1nnnnnxhyxyxhfyxhfk在在)(nnxyy 的假定下，的假定下，而而),(12kyhxhfknn)()(,()(,() )(,()(,(211hOxyxfykxyxfxhxyxfhkxyhxhfnnnnnnnn)()(,()(,()(,()(,(32hOxyxfyxyxfxyxfxhxyxhfnnnnnnnn)()( )( 32hOxyhxhynn)()( 21)( )()( )( )( 21,32321121hOxyhxhyyhOxyhxhyxh

12、yyyykknnnnnnnnn有代入将而而)()( ! 2)(321hOxyhxyhxyxynnnn）（）因此有因此有)(311hOyxyRnnn）（故故EulerEuler预估校正方法为的局部截断误差阶为预估校正方法为的局部截断误差阶为O O( (h h3 3) )。 定义定义3 3：若一个方法的局部截断误差为：若一个方法的局部截断误差为 ，则称该方法为则称该方法为p p阶方法阶方法, ,或称该方法具有或称该方法具有p p阶精度阶精度。 ）1(phO截断误差截断误差RemarkRemark：EulerEuler方法是一阶方法，梯形法和方法是一阶方法，梯形法和EulerEuler预估校正法是二

13、阶方法。预估校正法是二阶方法。整体截断误差与局部截断误差的关系整体截断误差与局部截断误差的关系且局部截断误差有界：且局部截断误差有界：则则EulerEuler法的整体截断误差法的整体截断误差 n n满足估计式：满足估计式： 其中其中L L为李普希兹常数，为李普希兹常数，b-ab-a为求解区间长度，为求解区间长度，定理定理：如果：如果f f( (x x，y y) )满足李普希兹（满足李普希兹（LipschitzLipschitz）条件）条件2121),(),(yyLyxfyxf), 2 , 1(21|22nMhRn) 1(2)(20)(LabLabneLhMe。)(max2xyMbxa 证明略证

14、明略。收敛性与稳定性收敛性与稳定性收敛性定义收敛性定义：如果某一数值方法对于任意固定的：如果某一数值方法对于任意固定的x xn n= =x x0 0+ +nhnh，当，当h h0(0(同时同时n n ) )时有时有y yn n y y( (x xn n) )，则称该方法则称该方法收敛收敛。定义定义 用一个数值方法，求解微分方程初值问题时，对给定的步长h0，若在计算 时引入误差 （也称扰动），但由此引起计算后面的 时的误差按绝对值均不增加，则称这个数值方法是稳定的稳定的。 nyn), 2 , 1(kyknRemarkRemark：该定理表明，整体截断误差比局部截断：该定理表明，整体截断误差比局部

15、截断误差低一阶。对其它方法，也有类似的结论。误差低一阶。对其它方法，也有类似的结论。稳定性定义稳定性定义 稳定性稳定性RemarkRemark：由于稳定性问题比较复杂，通常的：由于稳定性问题比较复杂，通常的做法是将满足李普希兹条件的微分方程模型做法是将满足李普希兹条件的微分方程模型化。设化。设 f f y y= = = =常数，此时微分方程为线性常数，此时微分方程为线性方程方程 y y = = y y。为保证微分方程的稳定性，。为保证微分方程的稳定性，假定假定 00。讨论某方法的稳定性，就是讨论。讨论某方法的稳定性，就是讨论该方法对模型方程的稳定性。该方法对模型方程的稳定性。 稳定性结论稳定性

16、结论EulerEuler法的稳定性条件是：法的稳定性条件是：梯形法梯形法是绝对稳定的。是绝对稳定的。EulerEuler预估校正格式的稳定性条件是：预估校正格式的稳定性条件是：对对非线性方程非线性方程，应视，应视，此时，此时 将是变化的。将是变化的。的变化将引起的变化将引起 h h的变化，的变化，属于绝对稳定区域，则认为对属于绝对稳定区域，则认为对20 h1)(2112hhyfyfyfhh如果步长如果步长h h固定，固定，此时，若此时，若此方程而言，方法是稳定的。此方程而言，方法是稳定的。7.2 7.2 泰勒展开法与龙格泰勒展开法与龙格- -库库塔（塔（RungeKuttaRungeKutta

17、）方法）方法问题问题：利用泰勒展开法推导高阶单步的：利用泰勒展开法推导高阶单步的求解常微分方程初值问题的数值方法。求解常微分方程初值问题的数值方法。从提高截断误差阶的阶数入手。从提高截断误差阶的阶数入手。 假定初值问题的解假定初值问题的解y y（x x）及函数及函数f f（x x，y y）是充分光滑的，则是充分光滑的，则: :)()(!)( ! 2)( 1)(21pnppnnnnhOxyphxyhxhyxyxy）（）（)()(,(!)(,( ! 2)(,()(1) 1(2pnnppnnnnnhOxyxfphxyxfhxyxhfxy当当n n 充分小时，略去余项充分小时，略去余项 ，则有，则有p

18、 p阶计算公式阶计算公式 )(1phO一、一、Taylor Taylor 方法方法 0) 1(21)(,(!),( 2),(yxyxfphyxfhyxhfyynnppnnnnnn其中，其中， ),(nnnyxfy ),(),(),( nnynnnnxnyxfyxfyxfy),(2 22nnyxyyyxyxxnyxfffffffffy 上式称为上式称为p p阶阶TaylorTaylor方法方法。特别地，当。特别地，当p p1 1时，就是时，就是EulerEuler公式。当公式。当p p2 2时，得二时，得二阶阶TaylorTaylor方法：方法：),(2212),(2nnyxyxnnnnnnnf

19、 ffhyxhfyyhyhyy 当当TaylorTaylor方法的阶数方法的阶数p p取的较大时，需计取的较大时，需计算算f f( (x x，y y) )的高阶导数值，计算量较大。特别的高阶导数值，计算量较大。特别当当f f( (x x，y y) )较复杂时，较复杂时，y y(x)(x)的高阶导数会很复的高阶导数会很复杂。因此杂。因此TaylorTaylor方法很少单独使用，但可以方法很少单独使用，但可以用它来启发思路。用它来启发思路。二、二、RungeRungeKutta Kutta 方法方法 基本思想基本思想：用不同点的函数值作线性组合，构：用不同点的函数值作线性组合，构造近似公式，把近似

20、公式和解的造近似公式，把近似公式和解的TaylorTaylor展开比较，展开比较，使前面的若干项吻合，从而使近似公式达到一定使前面的若干项吻合，从而使近似公式达到一定的阶数。一般的显式的阶数。一般的显式R-KR-K方法，可以写成方法，可以写成 11,22112321313312122122111NNNNNnNnNnnnnnnNNnnkkkyhxhfkkkyhxhfkkyhxhfkyxhfkkckckcyy，（），（），（），（其中，其中， 为常数，选取这些常数的原则是，为常数，选取这些常数的原则是，要求第一式的右端在要求第一式的右端在 处泰勒展开后，按处泰勒展开后，按h 的的幂次重新整理，得到

21、幂次重新整理，得到 iijiC，），nnyx(3322113121hhhyynn！与微分方程的解的与微分方程的解的TaylorTaylor展开式展开式 3213121)(hfhfhfxyxynnnnn！）（有尽可能多的项重合，即要求有尽可能多的项重合，即要求 ，nnnfff 321 上述公式叫做上述公式叫做N N级的级的Runge-KuttaRunge-Kutta方法，其局方法，其局部截断误差为部截断误差为。）1(NhO其中其中，nnnfff 表示表示),(),(),()(nnnnnxyxyyxfxy 显然，显然，EulerEuler法是一级一阶法是一级一阶R-KR-K方法。方法。 下面以二级

22、下面以二级R-KR-K公式为例，来说明公式为例，来说明R-KR-K方法的推导方法的推导过程。过程。 复习泰勒公式复习泰勒公式)()(! 11)()(*11*1*nnnxxpfxxpfpfpf 211122211222111121111 21)xx)(p(f)xx)(xx)(p(f)xx)(xx)(p(f)xx)(xx)(p(f)xx)(xx)(p(f)xx)(p(f!*nn*nn*nn*n*nn*n*nn*n*），（），（12122122111kyhxhfkyxhfkkckcyynnnnnn对对 要求适当选取系数系数 ， ， 和 ，使当 时，上式的局部截断误差为1C2C221）（nnxyy ）

23、3(hO将将 在在 处展开，有处展开，有 K2),(yxnn),(),(),(),(121212122nnynnxnnnnyxfkyxhfyxfhkyhxhfk)(),(),(2),(! 2132121121222hOyxfkyxfhkyxfhnnyynnxynnxx将 代入 有， 21,kk1ny)(),()(21),()(),(21),()(),()()()(423221232122322222122222122111hOyxfxyhcyxfxyhcyxfhcyxfxyhcyxfhcxyhcxyhcxykckcyynnyynnnxynnnxxnnynnnxnnnnn)(),()(),()(

24、2),(! 21),()(),()(3222212122222212hOyxfxyhyxfxyhyxfhyxfxyhyxhfxyhnnyynnnxynnnxxnnynnnxn)(2121)()()(43),(222122122222),(2122221hOhfcffcfchffcfchxyccxynnnnyxyyxyxxyxyxnn而而y y( (x xn n+1+1) )在在x xn n处的处的TaylorTaylor展式为：展式为：)()(2! 3121)()()()(! 31)(21)()()(43),(22),(4321hOhfffffffffhfffhxyxyhOhxyhxyhxyx

25、yxynnnnyxyxyyyxyxxyxyxnnnnnnn )()(61)2161()31()2161()21()21()()1 ()(43),(222122122222),(212222111hOhfffffcffcfchffcfchxyccyxynnnnyxyxyyyxyxxyxyxnnnnR故故 212112122221cccc由于二级方法由于二级方法 ，方程也可写为，方程也可写为 02c221222121211cccc要使要使R Rn n尽量小，应首先让尽量小，应首先让h h，h h2 2项的系数为零，即项的系数为零，即 上述方程所确定的解都能使二级上述方程所确定的解都能使二级 R-K

26、R-K方方法成为一个二阶方法。法成为一个二阶方法。此时此时R Rn n为为)()(61)2)(8161(43),(22hOhfffffffffcnnyxyxyyyxyxxnR为使为使R Rn n最小，应令最小，应令081612c从而有从而有局部截断误差为：局部截断误差为：)()(64),(3hOffffhnnyxyxynR32,32,41,4321212cc故有二阶故有二阶R-KR-K方法为：方法为：)32,32(),()3(41121211kyhxhfkyxhfkkkyynnnnnn该方法又称为二阶该方法又称为二阶HeunHeun（休恩）方法。（休恩）方法。EulerEuler预估预估- -

27、校正格式校正格式 ),(),()(21121211kyhxhfkyxhfkkkyynnnnnn若取若取1,2121221则ccRemark1Remark1：我们可以构造无穷多个二级我们可以构造无穷多个二级R-KR-K方方法，这些方法的截断误差均为法，这些方法的截断误差均为O O( (h h3 3),),即都是二即都是二阶方法。其中二阶阶方法。其中二阶HeunHeun方法是截断误差项数最方法是截断误差项数最少，且允许少，且允许f f 任意变化的情况下截断误差最小任意变化的情况下截断误差最小的二阶方法。的二阶方法。Remark2Remark2：二级二级R-KR-K方法不可能达到三阶方法不可能达到三

28、阶Remark3Remark3：同样可构造其他阶的同样可构造其他阶的R-KR-K方法，它方法，它们都有无穷多组解，且三级们都有无穷多组解，且三级R-KR-K方法阶数不超方法阶数不超过过3 3，四级，四级R-KR-K方法阶数不超过方法阶数不超过4 4。Remark4Remark4：更高阶的方法由于计算量较大，一更高阶的方法由于计算量较大，一般不再采用。般不再采用。标准（经典）四级四阶标准（经典）四级四阶R-KR-K公式公式 ),()2,2()2,2(),()22(61342312143211kyhxhfkkyhxhfkkyhxhfkyxhfkkkkkyynnnnnnnnnn关于关于R-KR-K方

29、法计算量的讨论方法计算量的讨论 二阶二阶R-KR-K方法需计算两个函数值，四阶方法需计算两个函数值，四阶R-KR-K方法需方法需计算四个函数值，但精度要比二阶方法高出两阶。因计算四个函数值，但精度要比二阶方法高出两阶。因此，要达到同样的精度，用低阶方法需步长取得比较此，要达到同样的精度，用低阶方法需步长取得比较小，但若用高阶方法则可以将步长取得大一些，从而小，但若用高阶方法则可以将步长取得大一些，从而降低计算量。降低计算量。四阶经典四阶经典R RK K方法的稳定性条件是方法的稳定性条件是1)(241)(61)(211432hhhh关于关于R-KR-K方法稳定性的讨论方法稳定性的讨论二阶二阶R

30、RK K方法的稳定性条件是方法的稳定性条件是 1)(2112hh 线性多步法的线性多步法的基本思想基本思想：如果充分利用：如果充分利用前面多步的信息来预测前面多步的信息来预测y yn n+ +k k，则可以期望获，则可以期望获得较高的精度。得较高的精度。7.3 7.3 线性多步法线性多步法 前面的前面的R RK K方法是增加一些非节点处的函方法是增加一些非节点处的函数值的计算来提高单步法的精度的，这样使计数值的计算来提高单步法的精度的，这样使计算量增加了许多。本节介绍多步法，是在不过算量增加了许多。本节介绍多步法，是在不过分增加计算量的情况下取得较高的计算精度。分增加计算量的情况下取得较高的计

31、算精度。 线性多步法公式的线性多步法公式的构造构造一般用两种方法，一般用两种方法，即即TaylorTaylor展开法与数值积分法。展开法与数值积分法。线性多步法的一般形式线性多步法的一般形式)(110111101rnrnnnrnrnnnffffhyyyy式中式中jiknknknrkyxff,), 1 , 0 , 1)(,(都为实数，且都为实数，且0rr。当。当 1 1 0 0时上式为隐式方法，当时上式为隐式方法，当 1 10 0时，上式为显示方法。由于求时，上式为显示方法。由于求y yn n+1+1用到前面用到前面y yn n, ,y yn n- -1 1,y yn n-r-r等等r r1 1

32、个值，且关于个值，且关于y yn n- -j j和和f fn n- -j j( (j j=0,1,2,=0,1,2,r r) )都是线性的，因此称上式为线性都是线性的，因此称上式为线性r r1 1步方法。步方法。 一、用数值积分方法构造线性多步法用数值积分方法构造线性多步法将将 方程两端从方程两端从 积分得积分得 )(,(xyxfy 1nknxx到001)()(,()()(1yxydxxyxfxyxynknxxknn（1）对对 取等距插值节点取等距插值节点 ，)(,()(xyxfxFknnnnxxxx11,。)(,(,),(,(11nnknknxyxfxyxf对应对应的函数值为的函数值为 如果

33、如果k k取取不同的值，以及不同的值，以及F F( (x x) )取不同的插值多项式近似，由取不同的插值多项式近似，由上式就可以推导出不同的线性多步公式。上式就可以推导出不同的线性多步公式。其插值余项为：其插值余项为： 1. 1. 阿达姆斯阿达姆斯(Adams)(Adams)外插公式外插公式 在在(1)(1)式中取式中取k k=0=0，并选择，并选择x xn n, ,x xn-n-1 1, , x xn n-2-2, ,x xn n-3-3作为作为插值节点，作函数插值节点，作函数F F( (x x) )的三次插值多项式：的三次插值多项式： 30303)()(iinijjjninjnxFxxxx

34、xLnnnnnnxxxxxxxxxxFxR3321)4(3)()()()(! 41)(把把F F( (x x)=)=L L3 3( (x x)+ )+ R R3 3( (x x) )代入（代入（1 1）式，有）式，有 11)()()()(331nnnnxxxxnndxxRdxxLxyxy略去上式右端第三项，得略去上式右端第三项，得), 5 , 4 , 3()()()(131ndxxLxyxynnxxnn对于上式积分部分用变量代换对于上式积分部分用变量代换x=xx=xn n+th+th，并注意到，并注意到hxxxxxxnnnnnn32211则则)(9)(37)(59)(5524)2)(2(3)(

35、) 3)(1(2)() 3)(2(2)() 3)(2)(1(! 3)()(3213211031nnnnnnnnxxxFxFxFxFhhdttttxFtttxFtttxFtttxFdxxLnn从而得到线性四步从而得到线性四步AdamsAdams显式公式：显式公式：), 5 , 4 , 3(),(9),(37),(59),(55243322111nyxfyxfyxfyxfhyynnnnnnnnnn其局部截断误差就是数值积分的误差其局部截断误差就是数值积分的误差111)()()()(! 41)(32)4(3nnnnnnxxnnnxxndxxxxxxxxxFdxxRR因因( (x x- -x xn n

36、)()(x x- -x xn n-1-1)()(x x- -x xn n-2-2)()(x x- -x xn n-3-3) )在在 x xn n, ,x xn n+1+1 上不上不变号，并设变号，并设F F(4)(4)( (x x) )在在 x xn n, ,x xn n+1+1 上连续，利用积分中上连续，利用积分中值定理，存在值定理，存在 n n x xn n, ,x xn n+1+1 ，使得，使得)(720251)(720251)()()()(! 41)5(5)4(532)4(11nnxxnnnnyhFhdxxxxxxxxxFRnnnn 因为插值多项式因为插值多项式L L3 3( (x x

37、) )是在是在 x xn n3 3, ,x xn n 上作出上作出的，而积分区间为的，而积分区间为 x xn n, ,x xn n+1+1 ，故上式称为，故上式称为AdamsAdams外外插公式。插公式。2. 2. 阿达姆斯阿达姆斯(Adams)(Adams)内插公式内插公式 若在若在(1)(1)式中取式中取k k=2=2，并选择，并选择x xn n1 1, ,x xn n, , x xn n-1-1, ,x xn n-2-2作为插值节点，作函数作为插值节点，作函数F F( (x x) )的三次插值多项式。的三次插值多项式。类似于上面的外插公式，有类似于上面的外插公式，有)(720195199

38、24)5(52111nnnnnnnnyhRffffhyy该公式也称为该公式也称为AdamsAdams内插公式，为三步隐式方法。内插公式，为三步隐式方法。3. 3. 阿达姆斯阿达姆斯(Adams)(Adams)预估预估- -校正公式校正公式 由于由于AdamsAdams内插公式是隐式方法，故用它做计内插公式是隐式方法，故用它做计算需使用迭代法。通常把算需使用迭代法。通常把AdamsAdams外插公式与内插公外插公式与内插公式结合起来使用，先由前者提供初值，再由后者式结合起来使用，先由前者提供初值，再由后者进行修正，即进行修正，即), 5 , 4 , 3;, 2 , 1 , 0(519),(924

39、93759552421)(11)1(1321)0(1nkfffyxfhyyffffhyynnnknnnknnnnnnn), 5 , 4 , 3(519,(92493759552421)0(111321)0(1nfffyxfhyyffffhyynnnnnnnnnnnnn当183Lyfh在求解区域内成立时，迭代收敛。在求解区域内成立时，迭代收敛。 若上式中的第二式只迭代一次，便得到若上式中的第二式只迭代一次，便得到AdamsAdams预预估估- -校正格式。校正格式。 Taylor Taylor展开法更具一般性，不仅可以构造用数展开法更具一般性，不仅可以构造用数值积分法得出的数值方法，而且还可导出

40、积分法得值积分法得出的数值方法，而且还可导出积分法得不到的方法。它比积分法更加灵活。下面仅举一例不到的方法。它比积分法更加灵活。下面仅举一例说明如何用这种方法构造线性多步法。说明如何用这种方法构造线性多步法。二、用二、用TaylorTaylor展开法构造线性多步公式展开法构造线性多步公式 首先以首先以x xn n-1-1, ,x xn n, ,x xn n+1+1为节点，构造形如为节点，构造形如)(110111101nnnnnnfffhyyy的公式。假设上式右边的公式。假设上式右边) 1, 1()()(,(), 1()(nnnixyxyxffnnixyyiiiiii将函数将函数)(),(),(

41、111nnnxyxyxy在在x x= =x xn n处展开，有处展开，有: : 5) 5(4) 4(321)(! 51)(! 41)(! 31)(! 21)(! 11)()(hxyhxyhxyhxyhxyxyxynnnnnnn 5) 6(4) 5(3) 4(21)(! 51)(! 41)(! 31)(! 21)(! 11)()(hxyhxyhxyhxyhxyxyxynnnnnnn 5) 6(4) 5(3) 4(21)(! 51)(! 41)(! 31)(! 21)(! 11)()(hxyhxyhxyhxyhxyxyxynnnnnnn代入给定公式并按代入给定公式并按h的幂次整理得到下式：的幂次整

42、理得到下式： )()!41!41!51()()!31!31!41()()2121!31()()!21()()()()()5(5111)4(4111311121111011101nnnnnnnxyhxyhxyhxyhxyhxyy将上式与将上式与 5)5(4)4(321)(!51)(!41)(!31)(!21)(! 11)()(hxyhxyhxyhxyhxyxyxynnnnnnn比较，选择系数比较，选择系数 i i( (i i=0,1)=0,1)和和 i i( (i i=-1,0,1)=-1,0,1)使两式中使两式中关于关于h h的同次幂的系数有尽可能多的项相等。故有：的同次幂的系数有尽可能多的项

43、相等。故有：2416161241:61212161:2121:1:1:11141113111210111100hhhhh求解上述方程组，得出求解上述方程组，得出 0 0， 1 1， 1 1， 0 0 ， 1 1。所。所得到的算式的局部截断误差为得到的算式的局部截断误差为O O( (h h5 5) )。)4(31111nnnnnfffhyyReamrkReamrk：我们也可以只要求前面几个方程组成立，：我们也可以只要求前面几个方程组成立，如要求前面如要求前面4 4个方程组成立时，所得算式的局部截断个方程组成立时，所得算式的局部截断误差为误差为O O( (h h4 4) )。如令。如令 0 00

44、0，带入上式的前，带入上式的前4 4个方程，个方程，解得解得 1 1=1=1， 1 1= = 1 1=1/3=1/3 ， 0 0=4/3,=4/3,于是得到计算公式于是得到计算公式为：为： 此时上式中第此时上式中第5 5式也恰巧成立。可以得到上式式也恰巧成立。可以得到上式得截断误差为：得截断误差为：)()(901)()()2412411201(12016)5(56)5(5111hOxyhhOxyhRnnn称上式为称上式为辛浦生（辛浦生（SimpsonSimpson）公式）公式，它可由数值积，它可由数值积分方法而得到。分方法而得到。我们也可以用类似的方法构造其它的线性多步法，我们也可以用类似的方

45、法构造其它的线性多步法，如前面的如前面的AdamsAdams公式等。公式等。三、出发值的计算三、出发值的计算 使用线性使用线性k k步法求解初值问题时，需要知道步法求解初值问题时，需要知道k k个个出发值出发值y y0 0, ,y y1 1,y yk k-1-1才能进行计算。然而初值问题才能进行计算。然而初值问题只提供一个只提供一个y yn n, ,还有还有k k1 1个出发值，需要通过别的方个出发值，需要通过别的方法计算出来。常用的方法是一步方法。由于初值对法计算出来。常用的方法是一步方法。由于初值对于确定微分方程的解有重要作用，因而在求解数值于确定微分方程的解有重要作用，因而在求解数值解时

46、，对出发值的精度也必须有相应的要求。为了解时，对出发值的精度也必须有相应的要求。为了保证多步方法的精确度，用于计算出发值的一步方保证多步方法的精确度，用于计算出发值的一步方法的阶数至少不低于多步方法的阶。法的阶数至少不低于多步方法的阶。 理论上讲，可用理论上讲，可用TaylorTaylor展开法和展开法和Runge-Runge-KuttaKutta方法，计算出发值。但由于方法，计算出发值。但由于TaylorTaylor展开展开法要计算高阶导数值，故最常用的方法还是选法要计算高阶导数值，故最常用的方法还是选择与多步法同阶的择与多步法同阶的Runge-KuttaRunge-Kutta方法。一旦出方

47、法。一旦出发值计算出来，线性多步法的计算量（特别是发值计算出来，线性多步法的计算量（特别是显式公式）就会很小，因为每次只须计算一次显式公式）就会很小，因为每次只须计算一次f f值。值。RemarkRemark：有关线性多步法的整体截断误差、：有关线性多步法的整体截断误差、收敛性及数值稳定性的讨论可参考有关文献。收敛性及数值稳定性的讨论可参考有关文献。7.47.4数值算例数值算例求解常微分方程初值问题求解常微分方程初值问题1)0(2yyxydxdy在在0,10,1上的解，取步长上的解，取步长h h0.10.1。计算结果如下图所示：计算结果如下图所示：第第7 7章例题章例题例例1 1：用欧拉预：用

48、欧拉预校方法求解初值问题（要求取步长校方法求解初值问题（要求取步长h=0.2h=0.2，计算，计算y(1.2)y(1.2)和和y(1.4)y(1.4)的近似值，小数点后保的近似值，小数点后保留留5 5位小数）：位小数）：1) 1 (0sin2yxyyy解解:欧拉预欧拉预校格式为校格式为:), 2 , 1 , 0)(,(),(2),()0(111)0(1nyxfyxfhyyyxhfyynnnnnnnnnn由由y(1)=y0=1计算得计算得:715488.0)2.1(63171.01)0(1yyy52611.0)4 .1 (47696.02)0(2yyy)sinsin( 1 . 0)sin( 2 . 0