常微分方程的数值解法.ppt

《常微分方程的数值解法.ppt》由会员分享，可在线阅读，更多相关《常微分方程的数值解法.ppt（93页珍藏版）》请在三一办公上搜索。

1、Tel：86613747E-mail：授课:68学分：4,第七章 常微分方程的数值解法,问题提出 倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式,其中直径D为常数.由于体积V与相对于容器底部的任意高度H的函数关系明确,因此在容器上可以方便地标出容器刻度,而对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下表是经过测量得到部分容器高度与直径的关系.,x,1,o,根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴后,问题即为如何根据任意高度x标出容器体积V的刻度,由微元思想分析可知,H 0 0.2 0.4 0.6 0.

2、8 1.0D 0.04 0.11 0.26 0.56 1.04 1.17,其中x表示高度,直径D是高度x的函数,记为D(x),因此得到如下微分方程初值问题,第七章 常微分方程的数值解法,只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是(7.0)中的函数D(x)无解析表达式,我们根本无法求出其解析解.怎样求解不能用解析方法求解的微分方程初制值问题,就是本章要讨论的重点.,(7.0),第七章 常微分方程的数值解法,7.1 引言 包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数

3、为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数都是一次的,则称它是线性的,否则称为非线性的。,在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。譬如,这个一阶微分方程就不能用初等函数及其积分来表达它的解。,再如，方程,的解,虽然有表可查,但对于表上没有给出 的值,仍需插值方法来计算,从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主

4、要讨论一阶常微分方程初值问题,(7.1),在区间a x b上的数值解法。,可以证明,如果函数在带形区域 R=axb,-y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使,对R内任意两个 都成立,则方程(7.1)的解 在a,b上存在且唯一。,常微分方程表示方法,在微分方程中,自变量的个数只有一个,称为常微分方程.自变量的个数为两个或两个以上的微分方程叫偏微分方程。,例如,当x=0时,y=1,可得c=1特解,当x=0时,y=1,可得c=-1特解,两边积分,通解,常微分方程解法回顾,-,例：设降落伞从跳伞塔下落后，所受空气阻力与速度 成正比，并设降落伞离开跳伞

5、塔时（t=0）速度为 零，求降落伞下落速度与时间的函数关系。解：设降落伞下落速度为v(t),降落伞在空中下落时，同时受到重力P与阻力R的作用，重力大小为mg，方向与v一致；阻力大小为kv（k为比例系数），方向与v相反，从而降落伞所受外力为 F=mg-kv根据牛顿第二定律 F=ma(a 为加速度），得函数v(t)应满足的方程为,按题意，初始条件为,分离变量后得,两端积分,考虑到mg-kv0,即,或,将初始条件,代如入上式得,于是所求特解为,7.2 数值方法的基本思想 对常微分方程初值问题(7.1)式的数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 的近似值。相邻两个

6、节点的间距 称为步长，步长可以相等，也可以不等。本章总是假定h为定数，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。,对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息 计算 的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数 进行不同的离散化处理。,对于初值问题的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一

7、类是计算yi+1时只用到xi+1,xi 和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到，即前面k步的值，此类方法称为多步法；其代表是亚当斯法。,7.3 欧拉（Euler）法7.3.1 Euler公式 欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)代表通过点 的一条称之为微分方程的积分曲线。积分曲线上每一点 的切线的斜率 等于函数 在这点的值。,Euler法的求解过程是:从初始点P0(即点(x0,y0)出发,作积分曲线y=y(x)在P0点上切线(其斜

8、率为),与x=x1直线,相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x0,y0)为斜率的切线方程为 当 时,得,这样就获得了P1点的坐标。,同样,过点P1(x1,y1),作积分曲线y=y(x)的切线交直线x=x2于P2点,切线 的斜率=直线方程为,当 时,得,当 时,得,由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,Pn。对已求得点以=为斜率作直线,取,从图形上看,就获得了一条近似于曲线y=y(x)的折线。,这样,从x0逐个算出对应的数值解,通常取(常数),则Euler法的计算格式,i=0,1,n(7.2),还可

9、用数值微分、数值积分法和泰勒展开法推导Euler格式。以数值积分为例进行推导。将方程 的两端在区间 上积分得，,选择不同的计算方法计算上式的积分项,就会得到不同的计算公式。,(7.3),公式推导,用左矩形方法计算积分项,代入(7.3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式,由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。,例7.1 用欧拉法解初值问题,取步长h=0.2,计算过程保留4位小数,解:h=0.2,欧拉迭代格式,当 k=0,x1=0.2时，已知x0=0，y0=1，有 y(0.2)y1=0.21(401)0.8当 k=1,x2=0.4时

10、，已知x1=0.2，y1=0.8，有 y(0.4)y2=0.20.8(40.20.8)0.6144当 k=2,x3=0.6时，已知x2=0.4，y2=0.6144，有 y(0.6)y3=0.20.6144(4-0.40.6144)=0.4613,例7.2 用Euler法求初值问题 y=x y 2 y(0)=0 的数值解（取 h=0.1 n=5)解：f(x,y)=x-y 2;x0=y0=0;h=0.1 由Euler法的递推公式得：yn+1=yn+0.1(xn y2n)yn=0 n=0,1,2,3,4,5由上式计算所得数据,n 0 1 2 3 4 5 xn 0 0.1 0.2 0.3 0.4 0.

11、5 yn 0 0 0.01 0.02999 0.05999 0.0995,二、隐式Euler(尤拉)格式,7.3.2 梯形公式为了提高精度,对方程 的两端在区间上 积分得，改用梯形方法计算其积分项，即,(7.4),代入(7.4)式,并用近似代替式中即可得到梯形公式,(7.5),由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（7.5）比欧拉公式(7.2)的精度高一个数值方法。,(7.5),(7.5)式的右端含有未知的yi+1,它是一个关于yi+1的函数方程,这类数值方法称为隐式方法。相反地,欧拉法是关于yi+1的一个直接的计算公式，这类数值方法称为显式方法。,例7.3 用梯形公式求下面初

12、值问题的解 在x=0.01上的值y(0.01)y=y y(0)=1 解：取h=0.01,由梯形公式得,y=exe0.01=1.010050167,7.3.3 两步欧拉公式 对方程 的两端在区间上 积分得,(7.6),改用中矩形公式计算其积分项，即,代入上式,并用yi近似代替式中y(xi)即可得到两步欧拉公式,(7.7),前面介绍过的数值方法,无论是欧拉方法,还是梯形方法，它们都是单步法,其特点是在计算yi+1时只用到前一步的信息yi;可是公式(7.7)中除了yi外,还用到更前一步的信息yi-1,即调用了前两步的信息,故称其为两步欧拉公式,欧拉法的局部截断误差 衡量求解公式好坏的一个主要标准是求

13、解公式的精度,因此引入局部截断误差和阶数的概念。定义7.1 在yi准确的前提下,即 时,用数值方法计算yi+1的误差,称为该数值方法计算时yi+1的局部截断误差。对于欧拉公式，假定，则有,而将真解y(x)在xi处按二阶泰勒展开,因此有,定义7.2 数值方法的局部截断误差为,则称这种数值方法的阶数是P。步长(h1)越小,P越高,则局部截断误差越小,计算精度越高。欧拉公式的局部截断误差为,欧拉方法仅为一阶方法。两步欧拉公式比欧拉公式精度也是高一个数值方法,设,前两步准确,则两步欧拉公式,把y(xi-1)在xi处展开成Taylor级数，即,将,在xn处按二阶泰勒展开,略去余项,用近似值yn代替y(x

14、n)有,由,再将y(xi+1)在xi处进行泰勒展开,(7.8),(7.9),所以,由(7.8)和(7.9)可得两步欧拉公式的局部截断误差为,即,7.3.5 改进的欧拉公式 显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。先用欧拉公式(7.2)求出一个初步的近似值,称为预测值,它的精度不高,再用梯形公式(7.5)对它校正一次,即迭代一次,求得yi+1,称为校正值,这种预测-校正方法称为改进的欧拉公式：,（7.10）,预测,校正,可以证明,公式(7.10)的精度为二阶。这是一种一步显式格式,它可以表

15、示为嵌套形式。,（7.11）,或者表示成下列平均化形式,（7.12）,7.3.6 改进欧拉法算法实现（1）计算步骤 输入,h,N 使用以下改进欧拉法公式进行计算 输出，并使 转到 直至n N 结束。,（2）改进欧拉法的流程图,(3)程序实现(见附录A A-15 改进欧拉法计算常微 分方程初值问题),例7.2 用改进欧拉法解初值问题,区间为0,1,取步长h=0.1,解:改进欧拉法的具体形式,本题的精确解为,计算见P158列表所示,例7.3 对初值问题,证明用梯形公式求得的近似解为,并证明当步长h0时,yn收敛于精确解证明:解初值问题的梯形公式为,整理成显式,反复迭代,得到,由于，有,证毕,7.4

16、 龙格-库塔（Runge-Kutta）法7.4.1 龙格-库塔(Runge-Kutta)法的基本思想 Euler公式可改写成,则yi+1的表达式与 y(xi+1)的Taylor展开式的前两项完全相同,即局部截断误差为。改进的Euler公式又可改写成,上述两组公式在形式上有一个共同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1,而且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉公式:每步计算一次f(x,y)的值,为一阶方法。改进欧拉公式需计算两次f(x,y)的值，它是二阶方法。它的局部截断误差为。,于是可考虑用函数f(x,y)在若干点上的函数值的线性

17、组合来构造近似公式，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导,又提高了计算方法精度的阶数。或者说,在 这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格库塔（Runge-Kutta）法的基本思想。,7.4.2 二阶龙格库塔法 在 上取两点xi和,以该两点处的斜率值k1和k2的加权平均(或称为线性组合)来求取平均斜率k*的近似值K，即,式中:k1为xi点处的切线斜率值，k2为 点处的切线斜率值,比照改进的欧拉法,将 视为，即可得,对

18、常微分方程初值问题(7.1)式的解 y=y(x),根据微分中值定理，存在点，使得,式中,K可看作是y=y(x)在区间 上的平均斜率。所以可得计算公式为：,（7.14）,将y(xi)在x=xi处进行二阶Taylor展开：,（7.15）,也即,（7.13）,将 在x=xi处进行一阶Taylor展开：,将以上结果代入（7.14）得：,（7.16）,对式(7.15)和(7.16)进行比较系数后可知,只要,（7.17）,成立,格式(7.14)的局部截断误差就等于,有2阶精度,式(7.17)中具有三个未知量,但只有两个方程,因而有无穷多解。若取,则p=1，这是无穷多解中的一个解，将以上所解的值代入式(7.

19、14)并改写可得,不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（7.17）有一簇形如上式的计算格式，这些格式统称为二阶龙格库塔格式。因此改进的欧拉格式是众多的二阶龙格库塔法中的一种特殊格式。,若取,则，此时二阶龙格-库塔法的计算公式为,此计算公式称为变形的二阶龙格库塔法。式中 为区间 的中点。,7.4.3 三阶龙格-库塔法,为了进一步提高精度，设除 外再增加一点,并用三个点,的斜率k1,k2,k3加权平均得出平均斜率k*的近似值，这时计算格式具有形式:,（7.18）,为了预报点 的斜率值k3,在区间 内有两个斜率值k1和k2可以用,可将k1,k2加权平均得出 上的平均斜率,从而得到 的

20、预报值,于是可得,运用Taylor展开方法选择参数,可以使格式(7.18)的局部截断误差为,即具有三阶精度，这类格式统称为三阶龙格库塔方法。下列是其中的一种，称为库塔（Kutta）公式。,（7.19）,7.4.4 四阶龙格库塔法,如果需要再提高精度，用类似上述的处理方法，只需在区间 上用四个点处的斜率加权平均作为平均斜率k*的近似值，构成一系列四阶龙格库塔公式。具有四阶精度，即局部截断误差是。由于推导复杂，这里从略，只介绍最常用的一种四阶经典龙格库塔公式。,（7.20）,7.4.5 四阶龙格库塔法算法实现(1)计算步骤 输入,h,N 使用龙格库塔公式（7.20）计算出y1 输出，并使 转到 直

21、至n N 结束。,(2)四阶龙格库塔算法流程图,程序实现(见附录A A-16 四阶龙格-库塔法计 算常微分方程初值问题),例7.4 取步长h=0.2，用经典格式求解初值问题,解:由四阶龙格-库塔公式可得,可同样进行其余yi的计算。本例方程的解为，数值解yi与准确解y(xi)的对照表见教材P163所示，从表中看到所求的数值解具有4位有效数字。,龙格库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的算法。,7.4.6 变步长的龙格-库塔

22、法 在微分方程的数值解中，选择适当的步长是非常重要的。单从每一步看，步长越小，截断误差就越小；但随着步长的缩小，在一定的求解区间内所要完成的步数就增加了。这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。以经典的四阶龙格-库塔法(7.20)为例。从节点xi出发，先以h为步长求出一个近似值，记为，由于局部截断误差为，故有,当h值不大时，式中的系数c可近似地看作为常数。,然后将步长折半,即以为 步长,从节点xi出发,跨两步到节点xi+1,再求得一个近似值,每跨一步的截断误差是,因此有,这样,由此可得,这表明以 作为 的近似值，其误差可用步长折半前后

23、两次计算结果的偏差,来判断所选步长是否适当,当要求的数值精度为时：,（1）如果，反复将步长折半进行计算，直至为止，并以上一次步长的计算结果作为。这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要反复判断，增加了计算工作量，但在方程的解y(x)变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。,7.5 亚当姆斯方法7.5.1 亚当姆斯格式 龙格-库塔方法是一类重要算法，但这类算法在每一步都需要先预报几个点上的斜率值，计算量比较大。考虑到计算yi+1之前已得出一系列节点上 的斜率值，能否利用这些已知值来减少计算量呢？这就是亚当姆斯（Adams）方法的设计

24、思想。,设用xi,xi+1两点的斜率值加权平均作为区间 上的平均斜率，有计算格式,（7.21）,选取参数，使格式（7.21）具有二阶精度。,将 在xi处Taylor展开,代入计算格式(7.21)化简,并假设,因此有,与y(xi+1)在xi处的Taylor展开式,相比较,需取,才使格式(7.21)具有二阶精度。这样导出的计算格式,称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式。,和四阶亚当姆斯格式。,(7.22),这里和下面均记,上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点 的斜率值来预报区间 上的平均斜率是个外推过程，效果不够理想。为了进一步改善精度，变外推为内插，即增加节点

25、xi+1的斜率值来得出 上的平均斜率。譬如考察形如,（7.23）,的隐式格式,设(7.23)右端的Taylor展开有,可见要使格式(7.23)具有二阶精度,需令,这样就可构造二阶隐式亚当姆斯格式,其实是梯形格式。类似可导出三阶隐式亚当姆斯格式,和四阶隐式亚当姆斯格式,（7.24）,7.5.2 亚当姆斯预报-校正格式 参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（7.22）和隐式（7.24）相结合，用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报-校正格式,（7.25）,预报：,校正：,这种预报-校正格式是四步法，它在计算yi+1时不但用到前一步的信息，而且要用到再前面三步的

26、信息，因此它不能自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格库塔法提供开始值。,例7.5 取步长h=0.1,用亚当姆斯预报-校正公式求解 初值问题,的数值解。,解:用四阶龙格-库塔公式求出发值,计算得：,表中的,yi和y(xi)分别为预报值、校正值和准确解(),以比较计算结果的精度。,再使用亚当姆斯预报-校正公式(7.25),见教材P166列表算得其余的计算结果,7.6 算法的稳定性及收敛性 7.6.1 稳定性 稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差

27、等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。,当在某节点上xi的yi值有大小为的扰动时，如果在其后的各节点 上的值yi产生的偏差都不大于，则称这种方法是稳定的。稳定性不仅与算法有关，而且与方程中函数f(x,y)也有关，讨论起来比较复杂。为简单起见，通常只针对模型方程,来讨论。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程都不稳定，也就很难用于其他方程的求解。,先考察显式Euler方法的稳定性。模型方程的Euler公式为,将上式反复递推后

28、，可得,或,式中,要使yi有界，其充要条件为,即,由于，故有,（7.26）,可见，如欲保证算法的稳定，显式Euler格式的步长h的选取要受到式（7.26）的限制。的绝对值越大，则限制的h值就越小。,用隐式Euler格式，对模型方程 的计算公式为，可化为,由于,则恒有,故恒有,因此，隐式Euler格式是绝对稳定的（无条件稳定）的（对任何h0）。,7.6.2 收敛性 常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值yi来近似替代y(xi)，这种近似替代是否合理，还须看分割区间 的长度h越来越小时，即 时，是否成立。若成立，则称该方法是收敛的，否则称为不收敛。这里仍以Euler

29、方法为例，来分析其收敛性。Euler格式,设 表示取 时,按Euler公式的计算结果,即,Euler方法局部截断误差为：,设有常数，则,（7.27）,总体截断误差,（7.28）,又,由于f(x,y)关于y满足李普希茨条件,即,代入(7.28)上式，有,再利用式（7.27），式（7.28）,即,上式反复递推后，可得,（7.29）,设（T为常数）,因为,所以,把上式代入式（7.29），得,若不计初值误差，即，则有,（7.30）,式(7.30)说明,当时,即,所以Euler方法是收敛的，且其收敛速度为，即具有一阶收敛速度。同时还说明Euler方法的整体截断误差为，因此算法的精度为一阶。,7.7 一阶

30、常微分方程组与高阶方程 我们已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。7.7.1 一阶常微分方程组对于一阶常微分方程组的初值问题,（7.31）,可以把单个方程 中的f 和y看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。,设 为节点上的近似解，则有改进的Euler格式为,预报：,校正：,（7.32）,又，相应的四阶龙格库塔格式（经典格式）为,（7.33）,式中,（7.34）,把节点xi上的yi和zi值代入式(7.34),依次算出,然后把它们代入式(7.33),算出节点x

31、i+1上的yi+1 和zi+1值。对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理,只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初值问题。,例7.6 用改进的Euler法求解初值问题,取步长h=0.1，保留六位小数。,解:改进的Euler法公式为,预报：,校正：,将 及h=0.1代入上式,得,由初值,计算得,7.7.2 高阶方程组 高阶微分方程(或方程组)的初值问题,原则上都可以归结为一阶方程组来求解。例如,有二阶微分方程的初值问题,（7.35）,在引入新的变量 后,即化为一阶方程组初值问题:,（7.36）,式（7.36）为一个一阶方程组的初值问题，对此可用中介

32、绍的方法来求解。例如应用四阶龙格-库塔公式得,（7.37）,（7.38）,消去，上式简化为：,（7.39）,（7.40）,上述方法同样可以用来处理三阶或更高阶的微分方程（或方程组）的初值问题,例7.7 求解下列二阶微分方程的初值问题,取步长h=0.1,解:先作变换：令，代入上式，得一阶方程组,用四阶龙格-库塔方法求解,按式(7.37)及(7.38)进行计算：取步长,时,然后计算 时的y2和z2；依此类推，直到i=9时的y10和z10，即可得到其数值解。,本章小结,本章介绍了常微分方程初值问题的基本数值解法。包括单步法和多步法。单步法主要有欧拉法、改进欧拉法和龙格库塔方法。多步法是亚当姆斯法。它

33、们都是基于把一个连续的定解问题离散化为一个差分方程来求解，是一种步进式的方法。用多步法求常微分方程的数值解可获得较高的精度。实际应用时，选择合适的算法有一定的难度，既要考虑算法的简易性和计算量，又要考虑截断误差和收敛性、稳定性。,龙格-库塔法较为常用，适用于多步方法中作初值计算和函数f(x,y)较为简单的场合。四阶标准龙格库塔法精度高，程序简单，易于改变步长，比较稳定，也是一个常用的方法，但计算量较大。当函数f(x,y)较为复杂，可用显式亚当姆斯方法或亚当姆斯预测校正方法，不仅计算量较小，稳定性也比较好，但不易改变步长。一般采用龙格库塔法提供初值y1,y2,y3，然后用亚当姆斯外推公式求得预测值，再由亚当姆斯内插值求得校正值yi+1，如此求得的值近似程度好且节省计算量，是一种较好的方法。,作业习题一7.17.9,