常微分方程初值问题数值解法

1、第第9章章 常微分方程初值问题数值解法常微分方程初值问题数值解法9.1 引言引言9.2 简单的数值方法简单的数值方法9.1引言引言 包含自变量、未知函数及未知函数的导数或微包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中分的方程称为微分方程。在微分方程中, 自变量的自变量的个数只有一个个数只有一个, 称为称为常微分方程常微分方程。自变量的个数为。自变量的个数为两个或两个以上的微分方程叫两个或两个以上的微分方程叫偏微分方程偏微分方程。微分方。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数程的阶数。如果

2、未知函数y及其各阶导数及其各阶导数都是一次的都是一次的,则称它是线性的则称它是线性的,否则称为非线性的。否则称为非线性的。 )(,nyyy 2 在高等数学中，对于常微分方程的求解，给出在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如了一些典型方程求解析解的基本方法，如可分离变可分离变量法量法、常系数齐次线性方程常系数齐次线性方程的解法、常系数非齐次的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。有限的，大多数的常微分方程是不可能给出解析解。 譬如譬如 22yxy

3、 这个一阶微分方程就不能用初等函数及其积分这个一阶微分方程就不能用初等函数及其积分来表达它的解。来表达它的解。 引言引言3再如，方程再如，方程 1)0(yyy的解的解 , ,虽然有表可查虽然有表可查,但对于表上没有给但对于表上没有给出出 的值的值,仍需插值方法来计算仍需插值方法来计算xey xe引言引言4从实际问题当中归纳出来的微分方程，通常主要依从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程靠数值解法来解决。本章主要讨论一阶常微分方程初值问题初值问题 00)(),(yxyyxfy ( 1 ) 在区间在区间a x ba x b上的数值解法上的数值解法。

4、 可以证明可以证明, ,如果函数在带形区域如果函数在带形区域 R=axb,R=axb,-y y 内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件条件，即存在常数，即存在常数L(L(它与它与x,yx,y无关无关) )使使 引言引言5引言引言6对对R内任意两个内任意两个 都成立都成立, ,则方程则方程( 1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 2121),(),(yyLyxfyxf9. 2 简单的数值方法简单的数值方法 对常微分方程初值问题对常微分方程初值问题( (1)式的数值解法，就是式的数值

5、解法，就是要算出精确解要算出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值 。相邻两个节点的间距相邻两个节点的间距 称为步长，步长可称为步长，步长可以相等，也可以不等。本章总是假定以相等，也可以不等。本章总是假定h为定数，称为定数，称为为定步长定步长，这时节点可表示为，这时节点可表示为bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi, 2 , 1 ,07数值方法的基本思想数值方法的基本思想8数值解法需要把连续性的问题加以离散化，从数值解法需要把连续性的问题加以离散化，从而

6、求出离散节点的数值解。而求出离散节点的数值解。 对常微分方程数值解法的基本出发点就是离散对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用化。其数值解法有两个基本特点，它们都采用“步步进式进式”，即求解过程顺着节点排列的次序一步一步，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息地向前推进，描述这类算法，要求给出用已知信息 计算计算 的递推公式。建立这的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、类递推公式的基本方法是在这些节点上用数值积分、021,yyyyiii1iy中的导数中的导数 进行不同的离散化处理进行不

7、同的离散化处理。对于初值问题对于初值问题y00)(),(yxyyxfy数值方法的基本思想数值方法的基本思想900)(),(yxyyxfy数值微分、泰勒展开等离散化方法，对初值问题数值微分、泰勒展开等离散化方法，对初值问题的数值解法，首先要解决的问题就是的数值解法，首先要解决的问题就是如何对微分方如何对微分方程进行离散化程进行离散化，建立求数值解的，建立求数值解的递推公式递推公式。递推公。递推公式通常有两类，一类是计算式通常有两类，一类是计算yi+1时只用到时只用到xi+1, xi 和和yi，即前一步的值，因此有了初值以后就可以逐步往下即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称

8、为计算，此类方法称为单步法单步法；其代表是；其代表是龙格龙格库塔库塔法。另一类是计算法。另一类是计算y yi+1i+1时，除用到时，除用到x xi+1i+1,x,xi i和和y yi i以外，以外，还要用到还要用到 ，即前面，即前面k步的值，此类步的值，此类方法称为方法称为多步法多步法；其代表是；其代表是亚当斯法亚当斯法。 ), 2 , 1( ,kpyxpipi数值方法的基本思想数值方法的基本思想101 Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简单的）方法是解初值问题的最简单的数值方法。初值问题数值方法。初值问题的解的解y=y(x)代表通过点代表通过点 的一条称之为的一

9、条称之为微分方微分方程的积分曲线程的积分曲线。积分曲线上每一点。积分曲线上每一点 的切线的切线的斜率的斜率 等于函数等于函数 在这点的值。在这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf11Euler公式公式12通常取通常取 ( (常数常数),),hhxxiii1将方程将方程 的两端在区间的两端在区间 上积分得，上积分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy(1) 选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 ,

10、就会得到不同的计算公式。就会得到不同的计算公式。 1)(,iixxdxxyxfEuler公式公式13用用左矩形方法左矩形方法计算积分项计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii 代入代入(1)式式,并用并用yi近似代替式中近似代替式中y(xi)即可得到向前欧拉即可得到向前欧拉（Euler）公式）公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度很低，所以欧拉由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。）公式当然很粗糙。 Euler公式公式14Euler公式公式15Euler法的计算格式法的计算格式 )(),(001xyyyxh

11、fyyiiii i=0,1,n ( 2 ) 例例1 用欧拉法解初值问题用欧拉法解初值问题 1)0()6 . 00(2yxxyyy取步长取步长h=0.2 ,计算过程保留计算过程保留4位小数位小数 解解: h=0.2, 欧拉迭代格式欧拉迭代格式 2),(xyyyxf21),(iiiiiiiiyhxhyyyxhfyy)2 , 1 , 0()4(2 . 0iyxyiii例题例题16例题例题17当当 i=0, x1=0.2时，已知时，已知x0=0，y0=1，有，有 y(0.2) y1=0.21(401)0.8当当 i=1, x2=0.4时，已知时，已知x1 =0.2， y1 =0.8，有，有 y(0.4

12、) y2 =0.20.8(40.20.8)0.6144当当 i=2, x3 =0.6时，已知时，已知x2 =0.4， y2 =0.6144，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 2 梯形公式梯形公式 为了提高精度为了提高精度, ,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得， 改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii( 4 ) 代入代入(4)(4)式式,

13、 ,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 梯形公式梯形公式18 (5)式的右端含有未知的式的右端含有未知的yi+1,它是一个关于它是一个关于yi+1的函数方程的函数方程,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地。相反地,欧拉法是关于欧拉法是关于yi+1的一个直接的计算公式，的一个直接的计算公式， 这类数这类数值方法称为显式方法。值方法称为显式方法。 梯形公式梯形公式19 由于数值积分的梯形公式比矩形公式的精度高，由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（5）比欧拉公式）比欧拉公式(2)的精度高一个数值的精度高一个数值方法。

14、方法。 ),(),(2111iiiiiiyxfyxfhyy(5 ) 3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间上的两端在区间上 积分得积分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 6 ) 改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式,并用并用yi近似代替式中近似代替式中y(xi)即可得到两步欧即可得到两步欧拉公式拉公式 两步欧拉公式两步欧拉公式20 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方法无论是欧拉方法,还还是梯形方法，

15、它们都是单步法是梯形方法，它们都是单步法,其特点是在计算其特点是在计算yi+1时只用到前一步的信息时只用到前一步的信息yi;可是公式可是公式(7)中除了中除了yi外外,还用到更前一步的信息还用到更前一步的信息yi-1,即调用了前两步的信息即调用了前两步的信息,故称其为两步欧拉公式故称其为两步欧拉公式 两步欧拉公式两步欧拉公式21),(211iiiiyxhfyy (7 ) 4 欧拉法的局部截断误差欧拉法的局部截断误差 衡量求解公式好坏的一个主要标准是求解公式的衡量求解公式好坏的一个主要标准是求解公式的精度精度, 因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。 定义定义1 在

16、在yi准确的前提下准确的前提下, 即即 时时, 用数值用数值方法计算方法计算yi+1的误差的误差 , 称为该数值方称为该数值方法计算时法计算时yi+1的局部截断误差。的局部截断误差。 对于欧拉公式，假定对于欧拉公式，假定 ，则有，则有)(iixyy 11)(iiiyxyR)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy欧拉法的局部截断误差欧拉法的局部截断误差22欧拉法的局部截断误差欧拉法的局部截断误差23而将真解而将真解y(x)在在xi处按二阶泰勒展开处按二阶泰勒展开 ),()(! 2)()()(121 iiiiixxyhxyhxyxy因此有因此有 )(!2)(21

17、1yhyxyii 定义定义2 数值方法的局部截断误差为数值方法的局部截断误差为 ,则称这种则称这种数值方法的阶数是数值方法的阶数是P。步长。步长(h1) 越小越小,P越高越高, 则局则局部截断误差越小部截断误差越小,计算精度越高。欧拉公式的局部截计算精度越高。欧拉公式的局部截断误差为断误差为 , 欧拉方法仅为一阶方法。欧拉方法仅为一阶方法。)(1phO)(2hO两步欧拉公式比欧拉公式精度也是高一个数值方法两步欧拉公式比欧拉公式精度也是高一个数值方法,设设 , 前两步准确前两步准确,则两步欧拉则两步欧拉公式公式 )(iixyy )(11iixyy)(,(2)(11iiiixyxhfxyy把把y(

18、xi-1)在在xi处展开成处展开成Taylor级数，即级数，即 )()(! 2)()()()(32111hOxxxyxxxyxyxyiiiiiiii 欧拉法的局部截断误差欧拉法的局部截断误差24)()(! 2)()(32hOxyhxyhxyiii ),(2)(11iiiiyxhfxyy由由 ),(2)()(! 2)()(32iiiiiyxhfhOxyhxyhxy )()(! 2)()(32hOxyhxyhxyiii 再将再将y(xi+1)在在xi+1处进行泰勒展开处进行泰勒展开)()(! 2)()()(321hOxyhxyhxyxyiiii (8)(9)欧拉法的局部截断误差欧拉法的局部截断误差

19、25欧拉法的局部截断误差欧拉法的局部截断误差26所以所以, 由由(8)和和(9)可得两步欧拉公式的局部截断误差可得两步欧拉公式的局部截断误差为为, 即即 )()(311hOyxyii5 改进的欧拉公式改进的欧拉公式 显式欧拉公式计算工作量小显式欧拉公式计算工作量小,但精度低。梯形公但精度低。梯形公式虽提高了精度式虽提高了精度,但为隐式公式但为隐式公式,需用迭代法求解需用迭代法求解,计算计算工作量大。综合欧拉公式和梯形公式便可得到改进工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。的欧拉公式。 先用欧拉公式先用欧拉公式(2)求出一个初步的近似值求出一个初步的近似值 ,称称为预测值为预测值

20、, 它的精度不高它的精度不高, 再用梯形公式再用梯形公式(5)对它对它1iy改进的欧拉公式改进的欧拉公式27改进的欧拉公式改进的欧拉公式28),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy预测预测 校正校正 （10） 1,2 , 1 , 0ni 可以证明可以证明,公式公式(10)的精度为二阶。这是一种一的精度为二阶。这是一种一步显式格式步显式格式,它可以表示为嵌套形式它可以表示为嵌套形式。校正一次校正一次,即迭代一次即迭代一次,求得求得yi+1,称为校正值称为校正值, 这种预这种预测测-校正方法称为改进的欧拉公式：校正方法称为改进的欧拉公式：),(,(),(211

21、iiiiiiiiyxhfyxfyxfhyy（11）或者表示成下列平均化形式或者表示成下列平均化形式 )(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy（12） 1,2 , 1 , 0ni改进的欧拉公式改进的欧拉公式29例例2 用改进欧拉法解初值问题用改进欧拉法解初值问题 1)0(2yyxyy区间为区间为 0,1 ,取步长取步长h=0.1 解解: 改进欧拉法的具体形式改进欧拉法的具体形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本题的精确解为本题的精确解为 ,计算结果略。计算结果略。 xxy21)(例题例题309. 3 龙

22、格龙格-库塔方法库塔方法1 龙格龙格-库塔库塔(Runge-Kutta)法法的基本思想的基本思想 Euler公式可改写成公式可改写成 ),(111iiiiyxfKhKyy则则yi+1的表达式的表达式y(xi+1)与的与的Taylor展开式的前两项完展开式的前两项完全相同全相同,即局部截断误差为即局部截断误差为 。 改进的改进的Euler公式又可改写成公式又可改写成 )(2hO31 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点:都是用都是用f(x,y)在某些点上值的线性组合得出在某些点上值的线性组合得出y(xi+1)的近似值的近似值yi+1,而且增加计算的次数而且增加计算的次

23、数f(x,y)的次数的次数,可提高截断误可提高截断误差的阶。如欧拉公式差的阶。如欧拉公式:每步计算一次每步计算一次f(x,y)的值的值,为一为一阶方法。改进欧拉公式需计算两次阶方法。改进欧拉公式需计算两次f(x,y)的值，它是的值，它是二阶方法。它的局部截断误差为二阶方法。它的局部截断误差为 。)(3hO龙格龙格-库塔法的基本思想库塔法的基本思想32),(),()(21121211hKyxfKyxfKKKhyyiiiiii 于是可考虑用函数于是可考虑用函数f(x,y)在若干点上的函数值的在若干点上的函数值的线性组合来构造线性组合来构造近似公式近似公式，构造时要求近似公式在，构造时要求近似公式在(xi,yi)处的处的Taylor展开式与解展开式与解y(x)在在xi处的处的Taylor展开式展开式的前面几项重合，从而使近似公式达到所需要的阶数。的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导既避免求偏导,又提高了计算方法精度的阶数。或者又提高了计算方法精度的阶数。或者说说,在在 这一步内多预报几个点的斜率值，然后这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高精