第三章 常微分方程的差分方法

1、School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 引言引言 包含自变量、未知函数及未知函数的导数或微分的方包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中程称为微分方程。在微分方程中, 自变量的个数只有一个自变量的个数只有一个, 称为常微分方程。自变量的个数为两个或两个以上的微分称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如

2、果未知函数数的阶数称为微分方程的阶数。如果未知函数y及其各阶导及其各阶导数数都是一次的都是一次的,则称它是线性的则称它是线性的,否则称为非线性的。否则称为非线性的。 )(,nyyy 第三章第三章 常微分方程的差分方法常微分方程的差分方法School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 在高等数学中，对于常微分方程的求解，给出了一些典在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次型方程求解析解的基本方法，如可分离变量法、常系数齐次

3、线性方程的解法、常系数非齐次线性方程的解法等。但能求线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。能给出解析解。 譬如譬如 22yxy 这个一阶微分方程就不能用初等函数及其积分来表达这个一阶微分方程就不能用初等函数及其积分来表达它的解。它的解。 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 从实际问题当中归纳出来的微分方程，通常主要依靠数值解从实际问题当中

4、归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题法来解决。本章主要讨论一阶常微分方程初值问题 00)(),(yxyyxfy ( 3.1 ) 在区间在区间a x b上的数值解法上的数值解法。 可以证明可以证明,如果函数在带形区域如果函数在带形区域 R=axb,-y内连续，内连续，且关于且关于y满足李普希兹满足李普希兹(Lipschitz)条件，即存在常数条件，即存在常数L(它与它与x,y无关无关)使使 2121),(),(yyLyxfyxf对对R内任意两个内任意两个 都成立都成立,则方程则方程( 3.1 )的解的解在在 a, b 上存在且唯一。上存在且唯一。 2

5、1, yy)(xyy School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 7.2 数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题(7.1)式的数值解法，就是要算出精确式的数值解法，就是要算出精确解解y(x)在区间在区间 a,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值 相邻两个节点的间距相邻两个节点的间距 称为步长，步长可以相等，称为步长，步长可以相等，也可以不等。本章总是假定也可以不等。本章总是假定h为定数，称为定

6、步长，这时节点为定数，称为定步长，这时节点可表示为可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。的数值解。 bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi, 2 , 1 ,0School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 对常微分方程数值解法的基本出发点就是离散化。其数对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用值

7、解法有两个基本特点，它们都采用“步进式步进式”，即求解过，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息要求给出用已知信息 计算计算 的递推的递推公式。建立这类递推公式的基本方法是在这些节点上用数值公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数中的导数 进行不同的离散化处理进行不同的离散化处理。 021,yyyyiii1iyy00)(),(yxyyxfySchool of Automation En

8、gineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 对于初值问题对于初值问题的数值解法，首先要解决的问题就是如何对微分方程进行离的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算类是计算yi+1时只用到时只用到xi+1, xi 和和yi，即前一步的值，因此有了初，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格是龙格库塔法。另一类是

9、计算库塔法。另一类是计算yi+1时，除用到时，除用到xi+1,xi和和yi以外，以外，还要用到还要用到 ，即前面，即前面k步的值，此类方法步的值，此类方法称为多步法；其代表是亚当斯法。称为多步法；其代表是亚当斯法。 00)(),(yxyyxfy), 2 , 1( ,kpyxpipiSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1 欧拉（欧拉（Euler）法）法3.1.1 Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简单的数值方法。）方法是解初值问题的最简

10、单的数值方法。初值问题初值问题的解的解y=y(x)代表通过点代表通过点 的一条称之为微分方程的积的一条称之为微分方程的积分曲线。积分曲线上每一点分曲线。积分曲线上每一点 的切线的斜率的切线的斜率 等于等于函数函数 在这点的值。在这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxfSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 Euler法的

11、求解过程是法的求解过程是:从初始点从初始点P0(即点即点(x0,y0)出发出发,作积分曲线作积分曲线y=y(x)在在P0点上切线点上切线 (其斜率其斜率为为 ),与与x=x1直线直线10PP),()(000yxfxy相交于相交于P1点点(即点即点(x1,y1),得到得到y1作为作为y(x1)的近似值的近似值,如上图所示。如上图所示。过点过点(x0,y0),以以f(x0,y0)为斜率的切线方程为为斜率的切线方程为 当当 时时,得得 )(,(0000 xxyxfyy1xx )(,(010001xxyxfyy这样就获得了这样就获得了P1点的坐标。点的坐标。 School of Automation

12、Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 x0 x1School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同样同样, 过点过点P1(x1,y1),作积分曲线作积分曲线y=y(x)的切线的切线交直线交直线x=x2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为21PP)(1xy),(11y

13、xf)(,(1111xxyxfyy)(,(121112xxyxfyy当当 时时,得得 2xx School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 当当 时时,得得 Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此获得了由此获得了P2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系列的就可获得一系列的点点:P1,P1,Pn。对已求得点。对已求得点以以 = 为斜率作直线为斜率作直线 ),(nnnyxP)(nxy),

14、(nnyxf)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)的折线的折线 。 Pi+ 1 Pn y= y(x) P1 Pi Pn Pi+ 1 P0 x0 x1 xi xi+ 1 xn Pi P1 这样这样,从从x0 逐个算出逐个算出对应的数值解对应的数值解 nxxx,21nyyy,21nPPPP321Schoo

15、l of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 通常取通常取 (常数常数),则则Euler法的计算格式法的计算格式 hhxxiii1)(),(001xyyyxhfyyiiii i=0,1,n ( 3.2 ) 还可用数值微分、数值积分法和泰勒展开法推导还可用数值微分、数值积分法和泰勒展开法推导Euler格式。格式。以数值积分为例进行推导。以数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy

16、11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 ,就就会得到不同的计算公式。会得到不同的计算公式。 1)(,iixxdxxyxf(3.3)School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 用左矩形方法计算积分项用左矩形方法计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(3.3)式式,并用并用yi近似代替式中近似代替式中y(xi)即可得到向前欧

17、拉即可得到向前欧拉（Euler）公式）公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度很低，所以欧拉（由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。公式当然很粗糙。 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例3.1 用欧拉法解初值问题用欧拉法解初值问题 1)0()6 . 00(2yxxyyy取步长取步长h=0.2 ,计算过程保留计算过程保留4位小数位小数 解解: h=0.2, 欧拉迭代格式欧拉迭代格式 2),(xyyyxf21),(

18、iiiiiiiiyhxhyyyxhfyy)2 , 1 , 0()4(2 . 0iyxyiii当当 k=0, x1=0.2时，已知时，已知x0=0，y0=1，有，有 y(0.2) y1=0.21(401)0.8当当 k=1, x2=0.4时，已知时，已知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 School of Automation Engineering自自

19、 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1.2 梯形公式梯形公式为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii( 3.4 ) 代入代入(3.4)式式,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy( 3.5 )

20、由于数值积分的梯形公式比矩形公式的精度高，因此梯由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（形公式（3.5）比欧拉公式）比欧拉公式( 3.2 )的精度高一个数值方法。的精度高一个数值方法。 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 ),(),(2111iiiiiiyxfyxfhyy( 3.5 ) (3.5)式的右端含有未知的式的右端含有未知的yi+1,它是一个关于它是一个关于yi+1的函数的函数方程方程,这类数值方法称为隐式方法。相反地这类数值方法称为隐式方法。

21、相反地,欧拉法是关于欧拉法是关于yi+1的一个直接的计算公式，的一个直接的计算公式， 这类数值方法称为显式方法。这类数值方法称为显式方法。 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1.3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间上的两端在区间上 积分得积分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 3.6 ) 改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 )(,)(,1111iiiixxxy

22、xfxxdxxyxfii代入上式代入上式,并用并用yi近似代替式中近似代替式中y(xi)即可得到两步欧拉公式即可得到两步欧拉公式 ),(211iiiiyxhfyy ( 3.7 ) School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方法无论是欧拉方法,还是梯形方还是梯形方法，它们都是单步法法，它们都是单步法,其特点是在计算其特点是在计算yi+1时只用到前一步时只用到前一步的信息的信息yi;可是公式可是公式(3.7)中除了中除了yi外

23、外,还用到更前一步的信息还用到更前一步的信息yi-1,即调用了前两步的信息即调用了前两步的信息,故称其为两步欧拉公式。故称其为两步欧拉公式。School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.1.4 欧拉法的局部截断误差欧拉法的局部截断误差 衡量求解公式好坏的一个主要标准是求解公式的精度衡量求解公式好坏的一个主要标准是求解公式的精度, 因因此引入局部截断误差和阶数的概念。此引入局部截断误差和阶数的概念。 定义定义3.1 在在yi准确的前提下准确的前提下, 即即 时时, 用数值方法

24、计算用数值方法计算yi+1的误差的误差 , 称为该数值方法计算时称为该数值方法计算时yi+1的局部的局部截断误差。截断误差。 对于欧拉公式，假定对于欧拉公式，假定 ，则有，则有)(iixyy 11)(iiiyxyR)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy而将真解而将真解y(x)在在xi处按二阶泰勒展开：处按二阶泰勒展开： ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(!2)(211yhyxyii 因此有因此有 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常

25、微分方程的差分方法常微分方程的差分方法 定义定义3.2 数值方法的局部截断误差为数值方法的局部截断误差为 ,则称这种数值方则称这种数值方法的阶数是法的阶数是P。步长。步长(h N 结束。结束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxxSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 （2）改进欧拉法的流程图）改进欧拉法的流程图 开 始 输 入x0, y0,h , N 1 n x0 + h x1 y0+ h

26、f( x0,y0 ) yp y0+ h f( x1,yp) yc ( yp+ yc) /2 y1 输 出x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结 束 n y School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例3.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 1)0(2yyxyy区间为区间为 0,1 ,取步长取步长h=0.1。解解: 改进欧拉法的具体形式改进欧拉法的具体形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyy

27、xyyyyxyyy本题的精确解为本题的精确解为 ,计算请对比计算请对比P98列表所示。列表所示。 xxy21)(School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例3.3 对初值问题对初值问题 1)0(0yyy证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为 nnhhynhx 证明证明: 解初值问题的梯形公式为解初值问题的梯形公式为 ),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理成显式整理成显式 nnyhhy反复迭代反复迭代,得到

28、得到 yhhyhhyhhyhhynnnnn.10ynnhhy22 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2 龙格龙格-库塔（库塔（Runge-Kutta）法）法3.2.1 龙格龙格-库塔库塔(Runge-Kutta)法的基本思想法的基本思想 Euler公式可改写成公式可改写成 ),(111iiiiyxfKhKyy则则yi+1的表达式的表达式y(xi+1)与的与的Taylor展开式的前两项完全相同展开式的前两项完全相同,即局部截断误差为即局部截断误差为 。 改进的改进

29、的Euler公式又可改写成公式又可改写成 )(2hO),(),()(21121211hKyxfKyxfKKKhyyiiiiiiSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点:都是用都是用f(x,y)在某在某些点上值的线性组合得出些点上值的线性组合得出y(xi+1)的近似值的近似值yi+1,而且增加计算的而且增加计算的次数次数f(x,y)的次数的次数,可提高截断误差的阶。如欧拉公式可提高截断误差的阶。如欧拉公式:每步

30、计每步计算一次算一次f(x,y)的值的值,为一阶方法。改进欧拉公式需计算两次为一阶方法。改进欧拉公式需计算两次f(x,y)的值，它是二阶方法。它的局部截断误差为的值，它是二阶方法。它的局部截断误差为 。)(3hOSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 于是可考虑用函数于是可考虑用函数f(x,y)在若干点上的函数值的线性组合在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在来构造近似公式，构造时要求近似公式在(xi,yi)处的处的Taylor展开展开式与解式与

31、解y(x)在在xi处的处的Taylor展开式的前面几项重合，从而使近展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导似公式达到所需要的阶数。既避免求偏导,又提高了计算方法又提高了计算方法精度的阶数。或者说精度的阶数。或者说,在在 这一步内多预报几个点的斜率这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高精度的值，然后将其加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格计算格式，这就是龙格库塔（库塔（Runge-Kutta）法的基本思想。）法的基本思想。 1,iixxSchool of Automation Engineering自自 动

32、动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.2 二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi和和 ,以该两点处的斜率值以该两点处的斜率值k1和和k2的加权平均的加权平均(或称为线性组合或称为线性组合)来求取平均斜率来求取平均斜率k*的近似值的近似值K，即，即 1,iixxphxxipi2211kkK式中式中:k1为为xi点处的切线斜率值，点处的切线斜率值， k2为为 点处的切线斜率值点处的切线斜率值,比照改进的欧拉法比照改进的欧拉法,将将 视为视为 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12

33、phkyphxfkii对常微分方程初值问题对常微分方程初值问题(3.1)式的解式的解 y=y(x),根据微分中值定理，根据微分中值定理，存在点存在点 ，使得，使得 ),(1iixx)()()(11iiiixxyxyxySchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 )(,()(yfyK式中式中 K可看作是可看作是y=y(x)在区间在区间 上的平均斜率。所以可上的平均斜率。所以可得计算公式为：得计算公式为： 1,iixxhKxyxyii)()(1)()(2211kkhxyi（3.1

34、4） 将将y(xi)在在x=xi处进行二阶处进行二阶Taylor展开：展开： )()(! 2)()()(321hOxyhxyhxyxyiiii （3.15） 也即也即 hKxyxyii)()(1（3.13）School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 将将 在在x=xi处进行一阶处进行一阶Taylor展开：展开： ),()(12phkyphxfphxykiii)(),(),(),(),(22hOyxfyxfyxfphyxfkiiyiiiixii)()()(2hOxyphxyii

35、 将以上结果代入（将以上结果代入（3.14）得：）得： )()()(22111kkhxyxyii)()()()()(221hOxyphxyxyhxyiiii )()()()()(32221hOxyphxyhxyiii （3.16） 对式对式(3.15)和和(3.16)进行比较系数后可知进行比较系数后可知,只要只要 211221p（3.17） 成立成立,格式格式(3.14)的局部截断误差就等于的局部截断误差就等于)(3hO有有2阶阶精度精度School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分

36、方法 式式(3.17)中具有三个未知量中具有三个未知量,但只有两个方程但只有两个方程,因而有无穷多解。因而有无穷多解。若取若取 ,则则p=1，这是无穷多解中的一个解，将以上所，这是无穷多解中的一个解，将以上所解的值代入式解的值代入式(3.14)并改写可得并改写可得 2121),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不难发现，上面的格式就是改进的欧拉格式。凡满足不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（条件式（3.17）有一簇形如上式的计算格式，这些格式统称）有一簇形如上式的计算格式，这些格式统称为二阶龙格为二阶龙格库塔格式。因此改进的欧拉格式是众多

37、的二库塔格式。因此改进的欧拉格式是众多的二阶龙格阶龙格库塔法中的一种特殊格式。库塔法中的一种特殊格式。 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 若取若取 ,则则 ，此时二阶龙格，此时二阶龙格-库塔库塔法的计算公式为法的计算公式为 0121, 12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2 , 1 , 0ni 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。 21ix1,iix

38、xSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.3 三阶龙格三阶龙格-库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 pix) 1(qpqhxxiqi并用三个点并用三个点 , , 的斜率的斜率k1,k2,k3加权平均加权平均得出平均斜率得出平均斜率k*的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式: ixpixqix),(),()1 (1213211phkyphxfkyxfkkkkhyyiiiiii（3.18） 为

39、了预报点为了预报点 的斜率值的斜率值k3,在区间在区间 内有两个斜率值内有两个斜率值k1和和k2可以用可以用,可将可将k1,k2加权平均得出加权平均得出 上的平均斜率上的平均斜率,从而从而得到得到 的预报值的预报值 qixqiixx,qiixx,)(qixyqiy21)1 (kkqhyyiqiSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 于是可得于是可得 ),(3qiqiyxfk运用运用Taylor展开方法选择参数展开方法选择参数 ,可以使格式可以使格式(3.18)的局部截断误

40、差为的局部截断误差为 ,即具有三阶精度，这类格式统称为三即具有三阶精度，这类格式统称为三阶龙格阶龙格库塔方法。下列是其中的一种，称为库塔（库塔方法。下列是其中的一种，称为库塔（Kutta）公式。公式。 ,qp)(4hO)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（3.19） School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.4 四阶龙格四阶龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，

41、只需在如果需要再提高精度，用类似上述的处理方法，只需在区间区间 上用四个点处的斜率加权平均作为平均斜率上用四个点处的斜率加权平均作为平均斜率k*的近的近似值，构成一系列四阶龙格似值，构成一系列四阶龙格库塔公式。具有四阶精度，即局库塔公式。具有四阶精度，即局部截断误差是部截断误差是 。 由于推导复杂，这里从略，只介绍最常用的一种四阶经典龙由于推导复杂，这里从略，只介绍最常用的一种四阶经典龙格格库塔公式。库塔公式。 qiixx,)(5hO)22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（3.20）

42、 School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 3.2.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1) 计算步骤计算步骤 输入输入 ,h,N ; 使用龙格使用龙格库塔公式（库塔公式（3.20）计算出）计算出y1; 输出输出 ，并使，并使 转到转到 直至直至n N 结束。结束。 10,xx11, yx0101,yyxxSchool of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差

43、分方法 (2) 四阶龙格四阶龙格库塔算法流程图库塔算法流程图 开 始 输 入 x0, y0,h , N 1 n x0 + h x1 f(x0,y0 ) k1, f(x0+ h /2 ,y0 + h k1/2 ) k2 f(x0+ h /2 ,y0 + h k2/2 ) k3, f(x1,y0+ h k3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 输 出 x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结 束 n y School of Automation Engineering自自 动动 化化 工工 程程 学学 院院第第 三三 章章 常微分方程的差分方法常微分方程的差分方法 例例7.4 取步长取步长h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 1)0(2yxyy10 x解解: 由四阶龙格由四阶龙格-库塔公式可得库塔公式可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh4