第9章常微分方程初值问题数值解法.ppt

第9章 常微分方程初值问题数值解法,9.1 引言 9.2 简单的数值方法与基本概念 9.3 龙格-库塔方法 9.4 单步法的收敛性与稳定性 9.5 线性多步法 9.6 方程组和高阶方程,9.1 引 言,科学技术中常常需要求解常微分方程的定解问题. 这类问题最简单的形式，是本章将着重考察的一阶方程的初值问题,我们知道，只有f(x, y)适当光滑譬如关于y满足利普希茨(Lipschitz)条件,理论上就可以保证初值问题的解yf(x)存在并且唯一.,虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.,所谓数值解法, 就是寻求解y(x)在一系列离散节点,上的近似值 y1,y2,yn,yn+1,. 相邻两个节点的间距hn=xn+1-xn称为步长. 今后如不特别说明，总是假定 hi=h(i=1,2,)为定数, 这时节点为xn=x0+nh(i=0,1,2,) (等距节点).,初值问题的数值解法有个基本特点，他们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进. 描述这类算法，只要给出用已知信息yn,yn-1,yn-2,计算yn+1的递推公式.,首先，要对微分方程离散化，建立求解数值解的递推公式. 一类是计算yn+1时只用到前一点的值yn，称为单步法. 另一类是用到yn+1前面 k 点的值yn,yn-1, yn-k+1，称为k步法. 其次，要研究公式的局部截断误差和阶，数值解yn与精确解y(xn)的误差估计及收敛性，还有递推公式的计算稳定性等问题.,9.2 简单的数值方法与基本概念,9.2.1 欧拉法与后退欧拉法,我们知道，在xy平面上，微分方程(1.1)式的解y=f(x)称作它的积分曲线，积分曲线上一点(x, y)的切线斜率等于函数f(x, y)的值. 如果按f(x, y)在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.,基于上述几何解释，我们从初始点P0(x0, y0)出发,先依方向场在该点的方向推进到x=x1上一点P1，然后再从P1点依方向场在该点的方向推进到 x=x2 上一点P2 , 循环前进做出一条折线P0 P1 P2.,一般地，设已做出该折线的顶点Pn,过Pn(xn, yn)依方向场的方向再推进到Pn+1(xn+1, yn+1)，显然两个顶点Pn,Pn+1的坐标有关系,这就是著名的(显式)欧拉(Euler)公式. 若初值y0已知，则依公式(2.1)可逐次逐步算出各点数值解.,即,例1 用欧拉公式求解初值问题,解 取步长h=0.1，欧拉公式的具体形式为,其中xn=nh=0.1n (n=0,1,10), 已知y0 =1, 由此式可得,依次计算下去，部分计算结果见下表.,与准确解 相比，可看出欧拉公式的计算结果精度很差.,欧拉公式具有明显的几何意义, 就是用折线近似代替方程的解曲线，因而常称公式(2.1)为欧拉折线法.,还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，,按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线.从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.,为了分析计算公式的精度，通常可用泰勒展开将y(xn+1)在xn处展开，则有,在yn=y(xn)的前提下，f(xn,yn )=f(xn,y(xn)=y(xn).于是可得欧拉法(2.1)的公式误差为,称为此方法的局部截断误差.,如果对方程(1.1)从xn到xn+1积分，得,右端积分用左矩形公式hf(xn,y(xn)近似，再以yn代替y(xn)，yn+1代替y(xn+1)也得到欧拉公式(2.1)，局部截断误差也是(2.3).,称为(隐式)后退的欧拉公式.,如果右端积分用右矩形公式hf(xn+1,y(xn+1)近似，则得到另一个公式,后退的欧拉公式与欧拉公式有着本质的区别, 后者是关于yn+1的一个直接计算公式，这类公式称作是显式的；前者公式的右端含有未知的yn+1，它实际上是关于yn+1的一个函数方程,这类方程称作是隐式的.,显式与隐式两类方法各有特点，考虑到数值稳定性等其他因素，人们有时需要选用隐式方法，但使用显式算法远比隐式方便.,隐式方程通常用迭代法求解，而迭代过程的实质是逐步显式化.,设用欧拉公式,给出迭代初值 ，用它代入(2.5)式的右端，使之转化为显式，直接计算得,然后再用 代入(2.5)式，又有,如此反复进行，得,由于f(x, y)对y满足Lipschitz条件(1.3). 由(2.6)减(2.5)得,由此可知，只要hL1，迭代法(2.6)就收敛到解.关于后退欧拉方法的公式误差，从积分公式看到它与欧拉法是相似的.,9.2.2 梯形方法,为得到比欧拉法精度高的计算公式，在等式(2.4) 右端积分用梯形求积公式近似, 并用yn代替y(xn), yn+1代替y(xn+1)，则得,称为矩形方法.,矩形方法是隐式单步法，用迭代法求解，同后退的欧拉方法一样，仍用欧拉法提供迭代初值，则矩形迭代公式为,为了分析迭代过程的收敛性, 将(2.7)与(2.8)相减, 得,于是有,使得,则当k时有 , 这说明迭代过程(2.8)是收敛的.,9.2.3 单步法的局部截断误差与阶,初值问题(1.1),(1.2)的单步法可用一般形式表示为,其中多元函数与f(x, y )有关，当含有yn+1时，方法是隐式的，若不含yn+1则为显式方法，所以显式单步法可表示为,(x, y, h)称为增量函数，例如对欧拉法(2.1)有,它的局部截断误差已由(2.3)给出, 对一般显式单步法则可如下定义.,定义1 设y(x)是初值问题(1.1),(1.2)的准确解, 称,为显式单步法(2.10)的局部截断误差.,Tn+1之所以称为局部的，是假设在xn前各步没有误差.当yn=y(xn)时，计算一步，则有,所以，局部截断误差可理解为用方法(2.10)计算一步的误差，也即公式(2.10)中用准确解y(x)代替数值解产生的公式误差. 根据定义, 显然欧拉法的局部截断误差为,即为(2.3)的结果. 这里 称为局部截断误差主项. 显然Tn+1=O(h2). 一般情形的定义如下,定义2 设y (x)是初值问题的准确解，若存在最大整数p使显式单步法(2.10)的局部截断误差满足,则称方法(2.10)具有p阶精度.,若将(2.10)展开式写成,则 称为局部截断误差主项.,以上定义对隐式单步法(2.9)也是适用的.例如，对后退欧拉法(2.5)其局部截断误差为,这里p=1是1阶方法，局部截断误差主项为,同样对梯形法(2.7)有,所以梯形方法(2.7)是二阶的. 其局部截断误差主项为,9.2.4 改进的欧拉公式,我们看到，梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式(2.9)进行实际计算时，每迭代一次，都要重新计算函数f(x, y )的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测. 为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.,具体地说，我们先用欧拉公式求得一个初步的近似值 ，称之为预测值，此预测值 的精度可能很差，再用梯形公式(2.7)将它校正一次，即按(2.8)式迭代一次，这个结果称之为校正值.,这样建立的预测校正系统通常称为改进的欧拉公式:,或表为下列平均化形式,(2.13),预测,校正,例2 用改进的欧拉法解例1中的初值问题(2.2).,解 仍取步长h=0.1,改进的欧拉公式为,部分计算结果见下表,同例1中的欧拉法的计算结果比较，明显改善了精度.,例 (两种方法的精度比较),用欧拉公式和改进的欧拉公式解下述初值问题，并与其准确解y=e-x+x进行比较.,解 取步长h=0.1，xk=kh(k=0,1,6).用两种方法进行计算对应结果及绝对误差见下表,9.3 龙格库塔方法,对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精度方法的途径.,9.3.1 显式龙格库塔法的一般形式,上节给出了显式单步法的表达式(2.10), 其局部截断误差为O(hp+1)，对欧拉法Tn+1=O(h2)，即方法为p=1阶，若用改进欧拉法(2.13)，它可表为,此时增量函数为,它比欧拉法的(xn, yn, h)=f(xn, yn), 增加了计算一个右函数f 的值，可望 p=2.若要使得到的公式阶数p更大, 就必须包含更多的f 值. 实际上从方程(1.1)等价的积分形式(2.4) ，即,若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积节点，为此可将(3.3)的右端用求积公式表示为,一般说来，点数r 越多，精度越高，上式右端相当于增量函数(xn, yn, h)，为得到便于计算的显式方法，可类似于改进欧拉法(3.1),(3.2)，将公式表示为,其中,这里ci, i, ij均为常数. (3.4)和(3.5)称为r级显式龙格-库塔(Runge-Kutta)法, 简称R-K方法.,当r=1, (xn, yn, h)=f(xn, yn)时，就是欧拉法，此时方法的阶为p=1. 当r=2时，改进欧拉法(3.1),(3.2)是其中一种，下面将证明其阶p=2. 要使公式(3.4),(3.5)具有更高的阶p，就要增加点数r. 下面我们只就r=2推导R-K方法. 并给出 r=3,4 时的常用公式，其推导方法与r=2时类似，只是计算较复杂.,9.3.2 二阶显式R-K方法,对r=2的R-K方法，由(3.4),(3.5)式可得如下计算公式,这里 c1, c2, 2, 21 均为待定常数，我们希望适当选取这些系数，使公式阶数 p 尽量高. 根据局部截断误差定义，推导出(3.6)的局部截断误差为,其中,这里yn=y(xn), yn+1=y(xn+1). 为得到Tn+1的阶p，要将上式各项在(xn, yn)处做泰勒展开，由于f(x, y )是二元函数，故要用二元泰勒展开，各项展开式为,将以上结果代入(3.7)，则有,要使公式(3.6)具有p=2阶，必须使,即,(3.9)的解是不唯一的. 可令c2=a0，则得,这样得到的公式称为二阶R-K方法.,则由此可以看出在改进的欧拉公式中相当于取(xn,yn), (xn+1,yn+1)两点处斜率的平均值，近似代替平均斜率，其精度比欧拉公式提高了.,如取a=1/2，则c1= c2=1/2, 2=21=1. 这就是改进的欧拉公式(3.1).,称为中点公式(变形的欧拉公式)，相当于数值积分的中矩形公式.也可以表示为,如取a=1，则c1=0, c2=1, 2=21=1/2. 得计算公式,对r=2的R-K公式(3.6)能否使局部误差提高到O(h4)? 为此 需把K2多展开一项，从(3.8)的 看到展开式中的项 是不能通过选择参数削掉的，实际上要使 h3 的项为零，需增加3个方程，要确定4个参数c1, c2, 2, 21，这是不可能的. 故r2的显式R-K方法的阶只能是p=2，而不能得到三阶公式.,9.3.3 三阶与四阶显式R-K方法,要得到三阶显式R-K方法，必须r=3. 此时计算(3.4), (3.5)的公式表示为,其中c1, c2, c3及2, 21, 3, 31, 32均为待定常数，公式(3.11)的局部截断误差为,只要K1, K2将按二元泰勒展开，使Tn+1O(h4)，可得待定参数满足方程,这是8个未知数6个方程的方程组，解不是唯一的. 可以得到很多公式. 满足条件(3.12)的公式(3.11)统称为三阶R-K公式. 下面只给出其中一个常见的公式.,此公式称为库塔三阶方法.,继续上述过程，经过较复杂的数学演算，可以导出各种四阶R-K公式，下列经典公式是其中常用的一个：,四阶R-K方法的每一步需要计算四次函数值f，可以证明其局部截断误差为O(h5). (例3见书p349),然而值得指出的是，龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质. 反之，如果解的光滑性差，那么，使用龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法. 实际计算时，我们应当针对问题的具体特点选择合适的算法.,9.3.4 变步长的龙格-库塔方法,单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题.,在选择步长时，需要考虑两个问题： 1. 怎样衡量和检验计算结果的精度？ 2. 如何依据所获得的精度处理步长？,我们考察四阶R-K公式(3.13) ，从节点xn出发，先以h为步长求出一个近似值，记为 ，由于公式的局部截断误差为O(h5)，故,然后将步长折半，即取为步长 ,从xn跨两步到xn+1，再求得一个近似值 ，每跨一步的局部截断误差是 ，因此有,比较(3.14)式和(3.15)式我们看到，步长折半后，误差大约减少到1/16，即有,由此易得下列事后估计式,这样，我们可以通过检查步长，折半前后两次计算结果的偏差,来判定所选的步长是否合适，具体地说，将区分以下两种情况处理：,这种通过加倍或折半处理步长的方法称为变步长方法.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的.,1.对于给定的精度,如果,我们反复将步长折半计算,直至为止,这时取最终得到的 作为结果；,2.如果为止，这时再将步长折半计算一次，就得到所要的结果.,9.4 单步法的收敛性与稳定性,9.4.1 收敛性与相容性,数值解法的基本思想是，通过某种离散化手段将微分方程转化为差分方程，如单步法,它在点xn处的解为yn，而初值问题在点xn处的精确解为y(xn)，记en=y(xn)-yn称为整体截断误差. 收敛性就是讨论当 x=xn 固定且 时en0的问题.,定义3 若一种数值方法对于固定的xn=x0+nh, 当h0时有yny(xn)，其中y(x)是(1.1),(1.2)的准确解，则称该方法是收敛的.,显然数值方法收敛是指en=y(xn)-yn0，对单步法(4.1)有下述收敛性定理：,定理1 假设单步法(4.1)具有p阶精度，且增量函数(x, y ,h)关于y满足利普希次条件,又设初值y0是准确的, 即y0=f(x0), 则其整体截断误差,证明 设以yn+1表示取yn=y(xn)用公式(4.1)求得的结果，即,则y(xn)-yn+1为局部截断误差，由于所给方法具有p阶精度，按定义2，存在定数C ，使,又由式(4.4)与(4.1)，得,利用利普希次条件(4.2)，有,从而有,即对整体截断误差en=y(xn)-yn成立下列递推关系式,据此不等式反复递推，可得,由此可以断定，如果初值是准确的，即e0=0，则(4.3)式成立. 定理证毕.,再注意到当x=x0+nhT时,最终得下列估计式,依据这一定理，判断单步法(4.1)的收敛性，归结为验证增量函数能否满足利普希次条件(4.2).,对于欧拉方法，由于其增量函数 就是f(x, y), 故当f(x, y)关于y满足利普希次条件时它是收敛的.,再考察改进的欧拉方法，其增量函数已由(3.2)式给出，假定f(x, y)关于y满足利普希次条件，即,这时有,即,设限定步长hh0(h0为定数)，上式表明关于y的利普希次常数为,因此改进的欧拉方法也是收敛的.,类似地, 不难验证其它龙格-库塔方法的收敛性.,定理1表明p1时单步法收敛, 并且当y(x)是初值问题(1.1),(1.2)的解, (4.1)具有p阶精度时, 则有展开式,所以p1的充分必要条件是 ，而 ，于是可给出如下定义：,定义4 若单步法(4.1)的增量函数满足,以上讨论表明p阶方法(4.1)当p1时与(1.1), (1.2)相容，反之相容方法至少是1阶的.,于是由定理1可知方法(4.1)收敛的充分必要条件是此方法是相容的.,则称单步法(4.1)与初值问题(1.1),(1.2)相容.,9.4.2 绝对稳定性与绝对稳定域,前面关于收敛性的讨论有个前提，必须假定数值方法本身的计算是准确的. 实际情形并不是这样，差分方程的求解还会有计算误差. 譬如由于数字舍入而引起的小扰动. 这类小扰动在传播过程中会不会恶性增长，以至于“淹没”了差分方程的“真解”呢？这就是差分方程的稳定性问题. 在实际计算时，我们希望某一步产生的扰动值，在后面的计算中能够被控制，甚至是逐步衰减的.,定义5 若一种数值方法在节点值yn上大小为的扰动，于以后各节点值ym(mn)上产生的偏差均不超过，则称该方法是稳定的.,下面以欧拉法为例考察计算稳定性.,例4 用欧拉公式求解初值问题,解 用欧拉法解方程y=-100y 得,其准确解 是一个按指数曲线衰减很快的函数.,若取步长h=0.025，则欧拉公式的具体形式为,计算结果见表, 明显计算过程不稳定, 但取h=0.005, yn+1=0.5yn, 则计算过程稳定.,对后退的欧拉公式，取h=0.025时，则计算公式为yn+1=(1/3.5)yn .计算结果见表, 这时计算过程是稳定的.,例题表明稳定性不但与方法有关，也与步长h有关，当然与方程中的f(x, y)有关. 为了只考察数值方法本身，通常只检验数值方法用于解模型方程的稳定性，模型方程为,其中为复数，这个方程分析较简单，对一般方程可以通过局部线性化化为这种形式，例如在(x, y)的邻域，可展开为,略去高阶项，再做变换即可得到的形式 u=u. 对于m个方程的方程组, 可线性化为y=Ay, 这里A为m×m雅可比矩阵(fi/yj)，若A有m个特征值1,2,m，其中i可能是复数，所以，为了使模型方程结果能推广到方程组，方程(4.8)中为复数. 为保证微分方程本身的稳定性，还应假定Re()0.,下面先研究欧拉方法的稳定性. 模型方程y=y的欧拉公式为,设在节点yn上有一扰动值n，它的传播使节点值yn+1产生大小为的扰动值n+1，假设用yn*=yn+n，按欧拉公式得出yn+1*=yn+1+n+1的计算过程不再有新的误差，则扰动值满足,可见扰动值满足原来的差分方程(4.9). 这样，如果差分方程的解是不增长的，即有,则它就是稳定的. 这一论断对于下面将要研究的其它方法同样适用.,显然，为要保证差分方程(4.9)的解是不增长的，只要选取h充分小，使,在=h的复平面上，这是以(-1,0)为圆心，1为半径的单位圆. 称为欧拉法的绝对稳定域，一般情形可由下面定义.,定义6 单步法(4.1)用于解模型方程y=y，若得到的解yn+1=E(h)yn，满足|E(h)|1，则称方法(4.1)是绝对稳定的. 在=h的平面上, 使|E(h)|1的变量围成的区域，称为绝对稳定区域，它与实轴的交称为绝对稳定区间.,对欧拉法E(h)=1+h，其绝对稳定域为|1+h|1,绝对稳定区间为-20,在例5中=-100,-2-100h0,即0h2/100=0.02为稳定区间，在例4中取h=0.025，故它是不稳定的，当取h=0.005时它是稳定的.,对二阶R-K方法，解模型方程(4.1)可得到,故,绝对稳定域由|E(h)|1得到，于是可得绝对稳定区间为-2h0，即0h2/.,类似可得三阶及四阶R-K方法的E(h)分别为,由|E(h)|1可得到相应的绝对稳定域. 当为实数时，则得绝对稳定区间，它们分别为,三阶显式R-K方法：,四阶显式R-K方法：,从以上讨论可知显式R-K方法的绝对稳定域均为有限域，都对步长h有限制. 如果h不在所给的绝对稳定区间内，方法就不稳定.,例4 分别取h=0.1及h=0.2，用经典的四阶R-K方法(3.1)计算初值问题,解 本例=-20,h分别为-2及-4，前者在绝对稳定区间内，后者则不在，用四阶R-K方法计算其误差见下表,从以上结果看到，如果步长h不满足绝对稳定条件，误差增长很快.,对隐式单步法,可以同样讨论方法的绝对稳定性,例如对后退欧拉法，用它解模型方程可得,故,由|E(h)|1，这是以(1,0)为圆心，1为半径的单位圆外部. 故方法的绝对稳定区间为-h0. 当0时，则0h，即对任何步长均为稳定的.,对隐式梯形法，它用于解模型方程(4.8)得,故,对Re()0有|E(h)|1，故绝对稳定域为=h的左半平面，绝对稳定区间为-h0，即0h时隐式梯形法均是稳定的.,9.5 线性多步法,在逐步推进的求解过程中，计算yn+1之前事实上已经求出了一系列的近似值y0,y1,yn，如果充分利用前面多步的信息来预测yn+1，则可以期望会获得较高的精度. 这就是构造所得线性多步法的基本思想.,构造多步法的主要途径基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1.1)两端积分后利用插值求积公式得到. 本节主要介绍基于泰勒展开的构造方法.,9.5.1 线性多步法的一般公式,如果计算yn+k时，除用yn+k-1的值，还要用到yn+i (i=0,1,k-2)的值，则称此方法为线性多步法. 一般的线性多步法公式可表示为,其中yn+1为y(xn+1)的近似，fn+i=f(xn+i, yn+i), 这里xn+i=xn+ih，i, i为常数， 0及0不全为零，则称(5.1)为线性k步法，计算时需先给出前面k个近似值y0,y1,yk-1，再由(5.1)逐次求出yk,yk+1,.,如果k=0，则(5.1)称为显式k步法，这时yn+k可直接由(5.1)算出；如果k0, 则(5.1)称为隐式k步法,求解时与梯形法(2.7)相同, 要用迭代法方可算出yn+k. (5.1)中系数i及i可根据方法的局部截断误差及阶确定，其定义为,定义7 设y(x)是初值问题(1.1), (1.2)的准确解，线性多步法(5.1)在xn+k上局部截断误差为,若Tn+k=O(hp+1)，则称方法(5.1)是p阶的，p1则称方法(5.1)与方程(1.1)是相容的.,由定义7，对Tn+k在xn处泰勒展开，由于,代入(5.2)得,其中,若在公式(5.1)中选择系数i及i，使它满足,由定义可知此时所构造的多步法是p阶的，且,称右端第一项为局部截断误差主项, cp+1称为误差常数.,根据相容性定义, p1，即c0=c1=0，由(5.4)得,故线性多步法(5.1)与微分方程(1.1)相容的充分必要条件是(5.6)成立.,显然，当k=1时，若1=0，则由(5.6)可求得,0=1，0=1.,此时公式(5.1)为,即为欧拉法. 从(5.4)可求得c2=1/20，故方法为1阶精度，且局部截断误差为,这和第2节给出的定义及结果是一致的.,对k=1，若10，此时方法为隐式公式，为了确定系数0,0,1，可由c0=c1=c2=0解得0=1, 0=1=1/2.于是得到公式,即为梯形公式.,由(5.4)可求得c2=-1/12，故p=2，所以梯形法是二阶方法，其局部截断误差主项是,这与第2节中讨论是一致的.,对k2的多步法公式都可利用(5.4) 确定系数i,i,并由(5.5)给出局部截断误差，下面只就若干常用的多步法导出具体公式.,9.5.2 阿当姆斯显式与隐式公式,p362自学.,9.5.3 米尔尼方法与辛普森方法,p366自学.,9.5.4 汉明方法,p367自学.,9.5.5 预测-校正方法,p368自学.,9.5.6 构造多步法公式的注记和例,p371自学.,9.6 方程组和高阶方程,9.6.1 一阶方程组,前面我们研究了单个方程y=f 的数值解，只要把y 和f 理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形.,考察一阶方程组,的初值问题，初始条件给为,若采用向量的记号，记(向量),则上述方程组的初值问题可表示为,求解这一初值问题的四阶龙格-库塔公式为(向量),式中(向量),或表示为(分量),其中(分量),这里yin是第i个因变量yi(x)在节点xn=x0+nh的近似值.,为了帮助理解这一公式的计算过程，我们再考察两个方程的特殊情形,这时四阶龙格-库塔公式具有形式,其中,这是一步法，利用节点xn上的值yn, zn，由(6.3)式顺序计算K1,L1,K2,L2,K3,L3,K4,L4，然后代入(6.2)式即可求得节点xn+1上的yn+1, zn+1.,9.6.2 化高阶方程为一阶方程组,关于高阶微分方程(或方程组)的初值问题，原则上总可以归结为一阶方程组来求解. 例如，考察下列m阶微分方程,初始条件为,只要引进新的变量,即可将m阶方程(6.4)化为如下的一阶方程组：,初始条件(6.5)则相应地化为,不难证明初始条件(6.4),(6.5)和(6.6),(6.7)是彼此等价的.,特别地，对于下列二阶方程的问题,引进新变量z=y,即可化为下列一阶方程组的初值问题:,针对这个问题应用四阶龙格-库塔公式(6.2)，有,由(6.3)式可得,如果消去K1,K2,K3,K4，则上述格式可表示为,这里,9.6.3 刚性方程组,p378自学.,课程结束! 希望同学们好好复习! 争取考个好成绩!,再见!,