chap7常微分方程的数值解法

7:ODE,1,第7章常微分方程（组）的数值解法,刘东毅天津大学理学院数学系,7:ODE,2,第7章常微分方程（组）的数值解法,主要目标：掌握常微分方程初值问题数值解法的基本理论掌握计算机上的常用算法,主要内容：初值问题计算格式的建立Runge-Kutta方法一阶常微分方程组与高阶方程的数值解法,7:ODE,3,第7章常微分方程（组）的数值解法,在科学研究和工程实践中会遇到很多微分方程，虽然从理论上可以证明其解的存在性，但其解的解析表达式往往是很难求解的，或者即使可以写出来，但也难于计算，此时，只能借助数值解来解决问题.常微分方程（组）定解问题是自然科学和工程技术领域中常见的数学模型.本章介绍求解此类问题的基本理论和数值解法。,7:ODE,4,定义7.0.1若存在常数L0,使得对一切的xa,b及y,均有,则称f(x,y)在D上关于y满足Lipschitz条件,其中L称为Lipschitz常数.,7:ODE,5,下面，我们给出常微分方程初值问题解的存在惟一性定理。,定理假设f(x,y)C(D),且关于y满足Lipschitz条件,则一阶常微分方程初值问题(7.0.1)存在唯一解.接下来，在此前提下,我们讨论上述初值问题(7.0.1)的数值解法。,7:ODE,6,然后在节点上建立逼近于原初值问题的计算格式(或差分格式),由此计算出原问题的解y(x)在节点x1,x2,.,xN处的近似值:y1,y2,.,yN，称它们为常微分方程初值问题的数值解.本章我们特别约定：y(xk)表示在节点xk的真值，yk表示相应的近似值。相邻两个节点的距离hn=xn+1-xn称为步长,通常取定步长,即节点xn=x0+nh,h0,n=0,1,N.,其基本思想是在区间a,b上引入一系列节点,7:ODE,7,7.1初值问题计算格式的建立,1.数值微分方法,在等距节点下讨论问题.利用两点数值微分公式,7.1.1计算格式的建立,将上式代入初值问题(7.0.1),有,(7.1.1),7:ODE,8,略去余项,并以数值解yn,yn+1代替y(xn)及y(xn+1),则得差分方程,上式称为Euler公式.利用此式可由初值y0出发按“步进式”方法,逐步求得数值解y1,y2,.,yN.,由于计算yn+1时,只用到它前一步的结果yn,这类公式称为单步法.又因为其关于yn+1是显式形式,故称该Euler公式为显格式.,7:ODE,9,如果利用下列数值微分公式,由类似的可导出,上述公式称为后退的Euler公式,此公式为单步法公式.又因为它关于yn+1成隐式形式,所以该公式为隐式公式，简称隐格式.,7:ODE,10,类似地，可导出,上述公式称为Euler两步法公式.这因为，当计算yn+1时,要用到yn-1与yn.显然它也是显格式.,如果利用下列三点数值微分公式,(7.1.3),7:ODE,11,设y(x)C2a,b,由Taylor公式有,由于故上式即为,略去余项,并以yn,yn+1代替y(xn)及y(xn+1),得到的差分方程正是Euler公式.,(7.1.4),2.Taylor展开法,7:ODE,12,3.数值积分方法,对，在区间xn,xn+1上积分，得,则有,对上式中的积分采用不同的数值积分公式可得到不同的差分方程.例如,对上式的积分采用左矩形数值积分公式,可得到Euler公式.,7:ODE,13,若对此式的积分采用数值积分梯形公式,则有,若略去余项,以yn,yn+1代替y(xn)及y(xn+1),得到的差分方程,7:ODE,14,上式称为梯形公式.由于它关于yn+1成隐式形式,故其为隐格式.隐格式求解比较困难,当yn已知时,要求yn+1，需解关于yn+1的非线性方程.在实际应用时,上式常与Euler公式联合使用,构成如下计算格式:,(7.1.6),7:ODE,15,隐式梯形公式的迭代格式,(7.1.7),由上式可以得到一个序列:,k=0,1,关于此序列的收敛性,有如下的定理.,7:ODE,16,定理7.1.1设f(x,y)在区域D上关于y满足Lipschitz条件,即,其中L为Lipschitz常数,当步长时,对任意的初值按格式（7.1.7）,生成的序列收敛于梯形公式（7.1.6）,的解.,7:ODE,17,为了减少计算量,可采用预测-校正格式.方法是先用Euler公式求得一个初始近似值称为预测值,再把带入梯形公式右端计算一次求得yn+1称之为校正值,即,上式称为预测-校正公式或改进的Euler公式.上式也可写成如下形式:,7:ODE,18,例7.1.1：利用Euler公式与改进的Euler公式求解初值问题（步长h=0.1）,解：由步长h=0.1，知节点设数值解为利用Euler公式得,7:ODE,19,计算结果见下表(见书P227表7.1),此初值问题的解析解为,从上表可以看出,数值解yn与解析解y(xn)比较,yn精度较差.,7:ODE,20,解此问题的改进的Euler公式为,同Euler公式比较,改进的Euler法显然精度提高了.,由于误差大小是评价计算格式优劣的重要依据,故需要给出有关误差的概念.,计算结果见下表(见书P228表7.2),7:ODE,21,7.1.2截断误差与方法的精度,定义7.1.1称误差en+1=y(xn+1)-yn+1为数值方法在点xn+1的截断误差,又称整体截断误差.,设yk=y(xk)(k=0,1,.,n)，则为数值方法在点xn+1的局部截断误差.,7:ODE,22,整体截断误差en+1是在没有引进舍入误差的情况下,纯粹因为不准确的计算格式造成的,故又称为方法误差.它不仅与x=xn+1这一步的计算有关,而且和xn,xn-1,.,x1这几步的计算都有关系.,局部截断误差是假设xn之前各数值解没有误差,仅由xn到xn+1这一步计算由计算格式引起的误差.,7:ODE,23,如Euler公式,在点xn+1的整体截断误差en+1=y(xn+1)-yn+1,局部截断误差,当前的整体截断误差=上一步的整体截断误差+上一步整体截断误差引起的导函数f(x,y)的误差+当前的局部截断误差。,7:ODE,24,定义7.1.2若某数值方法的局部截断误差为则称该方法具有P阶精度,或称其为P阶方法.,可以证明:,Euler方法的局部截断误差其具有一阶精度.梯形方法的局部截断误差其具有二阶精度.改进的Euler方法的局部截断误差具有二阶精度.,7:ODE,25,7.2Runge-Kutta方法,继续讨论前面的Taylor展开法。,设yC2a,b,由Taylor公式有,由故上式即为,略去余项,并以yn,yn+1代替y(xn)及y(xn+1),得到Euler公式.,7:ODE,26,进一步假设yCp+1a,b,由Taylor公式有,其中,由故式(7.2.1)即为,7:ODE,27,7.2Runge-Kutta方法,略去余项,并以数值解yn,yn+1代替(7.2.3)中的解析解y(xn)及y(xn+1),可得到一个差分方程，即,其中余项可写成,注：这里,7:ODE,28,在(7.2.3)中略去余项,用yn,yn+1代替y(xn)及y(xn+1),其中,称(7.2.4)式为求解常微分方程初值问题数值解Taylor的格式.,7:ODE,29,由于局部截断误差,可知它是一个p阶方法。当p=1时，上式正是Euler公式。但当p2时，需要计算f(x,y(x)的高阶导数，特别是对于复杂函数f(x,y(x)的求导，这无疑是大大增加计算量，这是它最大的缺点。因此高阶的Taylor方法是不实用的。,德国数学家C.Runge及M.W.Kutta提出了一种改进策略，得到了至今还作为高精度单步法被广泛使用龙格-库塔法(Runge-Kuttamethod)。,7:ODE,30,7.2.1Runge-Kutta方法的基本思想,Runge-Kutta方法是利用f在某些点处函数值的线性组合替代（7.2.4）步长h后面括号中的因子来构造差分方程,从而避免了高阶导数的计算,这就是Runge-Kutta方法的基本思想.,用f在某些点处函数值的线性组合替代这一部分,7:ODE,31,其一般形式为:,其中r是上式中调用f的个数,r称为级数，为待定参数,适当确定这些参数,可使之具有尽可能高的精度.如局部截断误差满足,7:ODE,32,7.2.2二阶Runge-Kutta方法,考虑r=2的情况,此时有,利用二元函数的一阶Taylor公式，即全微分公式,希望适当选择参数,使上式的局部截断,误差为,即为二阶方法.,下面将yn+1与y(xn+1)作比较,7:ODE,33,从而有,将上式代入,再由,7:ODE,34,得到,在下面要将yn+1与y(xn+1)作比较，使它们的局部截断误差满足,为此考虑y(xn+1)。,7:ODE,35,再根据y(xn+1)在点xn的一元3阶Taylor展开式,由刚才已得到的,让它们满足,7:ODE,36,即由,左式含有四个未知元三个方程,因此解不唯一.参数满足左式的一族公式统称二阶Runge-Kutta公式.,可得参数应满足下列方程组:,7:ODE,37,可见,二阶Runge-Kutta公式,每计算一步需要两次调用f的函数值.,7:ODE,38,7.2.3四阶Runge-Kutta方法,当r=4时,类似地可导出四阶Runge-Kutta公式,这种公式也有一族,其中常用地有:,标准(经典)的Runge-Kutta方法,7:ODE,39,Gill公式,Gill公式是标准的Runge-Kutta公式的改进形式,这种算法可节省存储单元,并能控制舍入误差的增长.,四阶Runge-Kutta公式,每一步计算需四次调用f的函数值,计算量较大,但其局部截断误差可达O(h5),精度较高.,7:ODE,40,解:此问题的计算公式为,计算结果如下:,显然在计算量大致相同的情况下,标准的Runge-Kutta法比改进的Euler方法精确度更高.（参见p227和p228的表）,例7.2.1用标准的四阶Rung-Kutta法解初值问题:取步长h=0.2.,7:ODE,41,7.5一阶常微分方程组与高阶方程初值问题的数值解法,7.5.1一阶常微分方程组初值问题,7:ODE,42,写成向量形式：,其中,注意:在形式上(7.5.2)与(7.0.1)一样,所以可以把求解常微分方程初值问题的各种数值方法推广到方程组上来.,7:ODE,43,利用向量值函数的微积分理论,很容易推导出一阶常微分方程组初值问题的数值解法.,如Euler公式,其中,7:ODE,44,(7.5.3)的分量形式为,或,7:ODE,45,四阶标准的Runge-Kutta公式,设则(7.5.5)的分量形式为,7:ODE,46,四阶标准的Runge-Kutta公式的分量形式,其中。,7:ODE,47,例7.5.1：试写出用中点公式解下列初值问题的计算公式：,解:令,则,再取,由向量形式的中点公式,7:ODE,48,中点公式的向量形式,上述中点公式的分量计算形式为,分量计算形式为,7:ODE,49,整理得分量计算格式,7:ODE,50,7.5.2高阶常微分方程初值问题的数值解法,高阶常微分方程初值问题的一般形式为,引入新变量,则(7.5.7)可化为如下一阶常微分方程组初值问题.,7:ODE,51,等价的一阶常微分方程组初值问题：,于是，可采用前面7.5.1所介绍的方法求此问题的数值解.,7:ODE,52,例7.5.2写出用标准四阶Runge-Kutta公式求解,的计算公式.,解:令，将此问题转化为一阶常微分方程组初值问题,再令则,7:ODE,53,求此问题数值解的标准四阶Runge-K