1、2022-5-271第第9章章 常微分方程初值问题数值解法常微分方程初值问题数值解法1. Euler公式公式2. 改进的欧拉公式改进的欧拉公式3. 龙格龙格库塔法库塔法4. 亚当斯法亚当斯法5. 算法的稳定性及收敛性算法的稳定性及收敛性2022-5-2729.1 9.1 引言引言 包含自变量、未知函数及未知函数的导数或微包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中分的方程称为微分方程。在微分方程中, , 自变量的自变量的个数只有一个个数只有一个, , 称为常微分方程。自变量的个数为称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方两个或两个
2、以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数程的阶数。如果未知函数y y及其各阶导数及其各阶导数都是一次的都是一次的, ,则称它是线性的则称它是线性的, ,否则称为非线性的。否则称为非线性的。 )(,nyyy 2022-5-273 在高等数学中，对于常微分方程的求解，给出在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解
3、法等。线性方程的解法等。但能求解的常微分方程仍然是但能求解的常微分方程仍然是有限的，大多数的常微分方程是不能给出解析解。有限的，大多数的常微分方程是不能给出解析解。 譬如譬如 22yxy 这个一阶微分方程就不能用初等函数来表达这个一阶微分方程就不能用初等函数来表达它的解。它的解。 2022-5-274再如，方程再如，方程 1)0(yyy的解的解 , ,虽然有表可查虽然有表可查, ,但对于表但对于表上没有给出上没有给出 的值的值, ,仍需插值方法来仍需插值方法来计算计算xey xe2022-5-275 从实际问题当中归纳出来的微分方程，通常主从实际问题当中归纳出来的微分方程，通常主要依靠数值解法
4、来解决。本章主要讨论一阶常微分要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题方程初值问题 00)(),(yxyyxfy ( 9.1 ) 在区间在区间a x ba x b上的数值解法上的数值解法。 可以证明可以证明, ,如果函数在带形区域如果函数在带形区域 R: axb,R: axb,-y y 内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数 L(L(它与它与x,yx,y无关无关) )使使 2121),(),(yyLyxfyxf 对对R内任意内任意x x及两个及两个 都成立都成立, ,则方程则方程(
5、9.1 )的的解解 在在 a, b 上上存在且唯一存在且唯一。 21, yy)(xyy 2022-5-276数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题( (9.1)式的数值解法，就式的数值解法，就是要算出精确解是要算出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节上的一系列离散节点点 处的函数值处的函数值 的近似值的近似值相邻两个节点的间距相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数，为定数，称为称为定步长定步长，这时节点可表示为，这时节点可表示为 数值解法需要把连
6、续性的问题加以离散化，数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。从而求出离散节点的数值解。 bxxxxann 110)(,),(),(10nxyxyxy.,10nyyyiiixxh 1niihxxi, 2 , 1,0 2022-5-277 对常微分方程数值解法的基本出发点就是离散对常微分方程数值解法的基本出发点就是离散化。其数值解法有化。其数值解法有两个基本特点，两个基本特点，11它们都采用它们都采用“步进式步进式”，即求解过程顺着节点排列的次序一步，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知一步地向前推进，描述这类算法，要求给出用已知信
7、息信息 计算计算 的递推公式。的递推公式。22建建立这类递推公式的基本方法是在这些节点上用数值立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值积分、数值微分、泰勒展开等离散化方法，对初值问题问题中的中的导数导数 进行不同的离散化进行不同的离散化处理处理。 iyyy,201 iyy 00)(),(yxyyxfy2022-5-278对于初值问题对于初值问题的数值解法，首先要解决的问题就是的数值解法，首先要解决的问题就是如何如何对微分方对微分方程进行离散化，建立求数值解的递推公式。递推公程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算式通常
8、有两类，一类是计算yi+1时只用到时只用到xi+1, xi 和和yi，即前一步的值，因此有了初值以后就可以逐步往下即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为计算，此类方法称为单步法单步法；其代表是；其代表是龙格龙格库塔库塔法。另一类是计算法。另一类是计算y yi+1i+1时，除用到时，除用到x xi+1i+1,x,xi i和和y yi i以外，以外，还要用到，还要用到， 即前面即前面k步的值，步的值，此类方法称为此类方法称为多步法多步法；其代表是亚当斯法。；其代表是亚当斯法。 00)(),(yxyyxfy), 2 , 1(,kpyxpipi 2022-5-2799.2 简单
9、的数值方法与基本概念简单的数值方法与基本概念9.2.1 Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最）方法是解初值问题的最简单的数值方法。初值问题简单的数值方法。初值问题的解的解y=y(x)y=y(x)为通过点为通过点 的一条积分曲线。的一条积分曲线。积分曲线上每一点积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在这点的值。在这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf2022-5-2710 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 Euler法的求解
10、过程是法的求解过程是: :从初从初始点始点P0(即点即点(x(x0 0,y,y0 0)出发出发, ,作积分曲线作积分曲线y=y(x)y=y(x)在在P0点上点上切线切线 ( (其斜率为其斜率为 ),),与与x=xx=x1 1直线直线10PP),()(000yxfxy相交于相交于P1点点( (即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1) )的近似值的近似值, ,如上图所示。过点如上图所示。过点(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0) )为为斜率的切线斜率的切线方程为方程为 当当 时时, ,得得 )(,(0000 x
11、xyxfyy 1xx )(,(010001xxyxfyy 这样就获得了这样就获得了P P1 1点的坐标。点的坐标。 hyxfyy),(0001 2022-5-2711 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同样同样, 过过点点P1(x x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)的切线的切线交直线交直线x=xx=x2 2于于P2点点, ,切线切线 的斜率的斜率直线方程为直线方程为21PP),()(111yxfxy )(,(1111xxyxfyy )(,(121112xxyxfyy 当当 时时,
12、 ,得得 2xx hyxfyy),(1112 2022-5-2712当当 时时, ,得得 Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程, ,就可获得一系就可获得一系列的点列的点: :P P1 1, ,P P1 1, , ,P Pn n。对已求得点对已求得点以以 为斜率作直线为斜率作直线 ),(nnnyxP),()(nnnyxfxy )(,(nnnnxxyxfyy 1 nxx)(,(11nxnnnnxxyxfyynnyxy )(取取2022-5-2713 从
13、图形上看从图形上看, ,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折线的折线 。 Pi+ 1 Pn y= y(x) P1 Pi Pn Pi+ 1 P0 x0 x1 xi xi+ 1 xn Pi P1 这样这样, ,从从x x0 0 出发出发逐个逐个算出算出对应的数值解nxxx,21nyyy,21nPPPP3212022-5-2714通常取通常取 ( (常数常数),),则则hhxxiii 1)2 . 9(1, 1 , 0)(),(001 nixyyyxhfyyiiii 00)(),(yxyyxfy微分方程的初值问题微分方程的初值问题Euler法的计算格式为：法的计算格
14、式为： 用折线近似于曲线得到的计算公式称为欧拉公式用折线近似于曲线得到的计算公式称为欧拉公式2022-5-2715 Euler格式格式还可用还可用数值微分数值微分、数值积分数值积分和和泰勒展泰勒展开开等方法得到。等方法得到。，微分方程),(yxfy hxyxyxynnn)()()(1 的导数用差商代替节点nx),()()(1nnnnyxfhxyxy ),(1nnnnyxhfyy 2022-5-2716 Euler格式格式用用泰勒展开泰勒展开方法得到。方法得到。，微分方程),(yxfy 展开处在节点将Taylorxxynn)(1 )(2)()()()(21 yhxyhxyhxyxynnnn),(
15、1nnnnyxhfyy )(22 yhh的高次项略去)()()(1nnnxyhxyxy 2022-5-2717 Euler 格式格式用用数值微分数值微分方法得到。方法得到。),(yxfy 1, iixx 11),(iiiixxxxdxyxfdxy)3 . 9()(,),()()(111 iiiixxxxiidxxyxfdxyxfxyxy 选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项就会得到不同的计算公式。就会得到不同的计算公式。 1)(,iixxdxxyxf将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，2022-5-2718 用左矩形方法计算积分项用左矩形
16、方法计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii 代入代入(9.3)(9.3)式式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到即可得到向前欧拉（向前欧拉（EulerEuler）公式）公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度很低，所以由于数值积分的矩形方法精度很低，所以欧拉（欧拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。 2022-5-2719例例9.1 用欧拉法解初值问题用欧拉法解初值问题 )6 . 00( ,1)0(2 xyxyyy取步长取步长h=0.2 ,h=0.2 ,计算过程保留计算过
17、程保留4 4位小数位小数 解解: :20.2,( , )hf x yyxy )2 , 1 ,0()4(2 .01 kyxyykkkk当当 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 )(
18、2 . 0),(21kkkkkkkkyxyyyxhfyy 2022-5-2720.2 .02(1 .0,(1计算及结果如下解： .1.1) nnnnnnnnnnnyxyyxyyyxhfyy clear; y=1, x=0, %初始化for n=1:10 y=1.1*y-0.2*x/y, x=x+0.1,endy = 1 x = 0 y = 1.1000 x = 0.1000y = 1.1918 x = 0.2000y = 1.2774 x = 0.3000y = 1.3582 x = 0.4000y = 1.4351 x = 0.5000y = 1.5090 x = 0.6000y = 1.5
19、803 x = 0.7000y = 1.6498 x = 0.8000y = 1.7178 x = 0.9000y = 1.7848 x = 1.0000)10( ,1)0(2, 1 . 0 xyyxyyEulerh公公式式求求解解利利用用取取例例2022-5-2721对方程对方程 的两端在区间的两端在区间 上积分得，上积分得，),(yxfy 1, iixx)4 . 9()(,)()(11 iixxiidxxyxfxyxy 2)(,()(,()(,1111 iiiiiixxxyxfxyxfxxdxxyxfii代入代入(9.4)(9.4)式式, ,即可得到梯形公式即可得到梯形公式 )5 . 9(
20、),(),(2111 iiiiiiyxfyxfhyy 由于数值积分的梯形公式比矩形公式的精度高，由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（9.59.5）是比欧拉公式）是比欧拉公式( 9.2 )( 9.2 )精度高精度高的一个数值方法。的一个数值方法。 为了提高精度为了提高精度,改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 9.2.2 梯形公式梯形公式2022-5-2722),(),(2111iiiiiiyxfyxfhyy( 9.5 ) ( (9.2) 和和(9.5)式粗看差不多式粗看差不多， )(),(001xyyyxhfyyiiii( 9.2 ) Eu
21、ler公式梯形公式 但但( (9.5)式的右端含有未知的式的右端含有未知的y yi+1i+1, ,它是一个关它是一个关于于y yi+1i+1的函数方程的函数方程, ,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地相反地, ,欧拉法是关于欧拉法是关于y yi+1i+1的一个直接的计算公式，的一个直接的计算公式， 这类数值方法称为这类数值方法称为显式方法显式方法。 2022-5-27239.2.3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间上的两端在区间上 积分得积分得 ),(yxfy 11, iixx 11)(,)()(11iixxiidxxyxfxyxy ( 9.6 ) 改
22、用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 )(,)()(,1111iiiixxxyxfxxdxxyxfii 代入上式代入上式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到两步即可得到两步欧拉公式欧拉公式 )7 .9(),(211iiiiyxhfyy 11111),(),( iiiiiiiiyPyxPyxP2022-5-2724 前面介绍过的数值方法前面介绍过的数值方法, ,无论是欧拉无论是欧拉方法方法, ,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法, ,其其特点是在计算特点是在计算y yi+1i+1时只用到前一步的信息时只用到
23、前一步的信息y yi i; ;可是公式可是公式( (9.7)中除了中除了y yi i外外, ,还用到更前一步还用到更前一步的信息的信息y yi-1i-1, ,即调用了前两步的信息即调用了前两步的信息, ,故称其故称其为两步欧拉公式为两步欧拉公式 2022-5-2725 定义定义9.1 在在yi准确准确的前提下的前提下, 即即 时时, 用数值方用数值方法计算法计算yi+1的误差的误差 , 称为该数值方法计称为该数值方法计算时算时yi+1的局部截断误差。的局部截断误差。)(iixyy 11)( iiiyxyR9.2.4 欧拉法的局部截断误差欧拉法的局部截断误差 衡量求解公式好坏的一个主要标准是求解
24、公式的衡量求解公式好坏的一个主要标准是求解公式的精度精度, 因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。故称为局部截断误差。得到的误差步精确的前提下在假设前,1 i整体截断误差。iiyxyr )(12022-5-2726假定假定 由欧拉公式，则有由欧拉公式，则有)(iixyy )()(1iiixyhxyy ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(! 2)(211 yhyxyii因此有因此有 展开处二阶在将精确解Taylorxxyii)(1 近似解局部截断误差局部截断误差 1ke2022-5-2727定义定义9.2 如果数值方法的局部截断误
25、差为如果数值方法的局部截断误差为 ,则称这种数值方法的则称这种数值方法的阶数是阶数是P。)(1 phO 显然，步长显然，步长(h N 结束。结束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx2022-5-2736（2）改进欧拉法的流程图）改进欧拉法的流程图 开 始 输 入x0, y0,h , N 1 n x0 + h x1 y0+ h 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
26、 2022-5-2737例例9.2 9.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 区间为区间为 0,10,1 , ,取步长取步长h=0.1h=0.1 解解: : 改进欧拉法的具体形式改进欧拉法的具体形式 本题的精确解为本题的精确解为xxy21)( 0,1,9i 000,1xy ,1)0(2 yyxyy,)(21)2( 1 . 0)2( 1 . 011 cpipipiciiiipyyyyxyyyyxyyy2022-5-2738clearx=0,yn=1 %初始化for n=1:10yp=yn+0.1*(yn-2*x/yn); %预测x=x+0.1;yc=yn+0.1*(yp-2*x/yp
27、) ;yn=(yp+yc)/2 %校正end2022-5-2739例例9.3 对初值问题对初值问题 1)0(0yyy证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为 nnhhy 22并证明当步长并证明当步长h h0 0时时,y,yn n收敛于精确解收敛于精确解xe ),(),(2111 nnnnnnyxfyxfhyyyyxf ),(211 nnnnyyhyy整理成显式整理成显式 nnyhhy 221反复迭代反复迭代, ,得到得到 012312122.222222yhhyhhyhhyhhynnnnn 10 ynnhhy 22证明证明: : 解初值问题的梯形公式为解初值问题的梯形公式为20
28、22-5-2740由于由于 ，有，有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy 22xnhy elim0 证毕证毕 2022-5-27419.3 9.3 龙格龙格- -库塔（库塔（Runge-KuttaRunge-Kutta）法）法hxyxyii)()(1 只有对平均斜率提供一种算法便可得到一种计算格式只有对平均斜率提供一种算法便可得到一种计算格式)(),()(hxyhxfhxyiii )(,1 yhfyyii9.3.1 龙格龙格-库塔法的基本思想库塔法的基本思想上上的的平平均均斜斜率率区区间间,)(,1 iixxyf),(y
29、xfy 11),(iiiixxxxdxyxfdxy2022-5-2742 111),(hKyyyxfKiiii )(2),(),(2111121KKhyyhKyxfKyxfKiiiiiiEuler公式可改写成公式可改写成改进的改进的Euler公式又可改写成公式又可改写成用左端点的斜率作为平均斜率得到的计算格式用左端点的斜率作为平均斜率得到的计算格式用左、右端点斜率的平均值作为平均斜率得到的计用左、右端点斜率的平均值作为平均斜率得到的计算格式算格式2022-5-2743 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点: :都是用都是用f(x,y)f(x,y)在某些点上值在某些点
30、上值( (斜率斜率) )的线性组合得出的线性组合得出y(xy(xi+1i+1) )的近似值的近似值y yi+1i+1, ,而且增加计算而且增加计算f f( (x x, ,y y) )的次数的次数, ,可提高截可提高截断误差的阶。如欧拉公式断误差的阶。如欧拉公式: :每步计算一次每步计算一次f f( (x x, ,y y) )的的值值, ,为一阶方法。改进欧拉公式需计算两次为一阶方法。改进欧拉公式需计算两次f f( (x x, ,y y) )的值，它是二阶方法。它的局部截断误差为的值，它是二阶方法。它的局部截断误差为 。)(3hO2022-5-2744 于是可考虑在于是可考虑在 内多内多预报几个
31、点的预报几个点的斜率值，然后将其加权平均作为平均斜率斜率值，然后将其加权平均作为平均斜率，构，构造时要求近似公式在造时要求近似公式在( (x xi i, ,y yi i) )处的处的TaylorTaylor展开式展开式与解与解y y( (x x) )在在x xi i处的处的TaylorTaylor展开式的前面几项重展开式的前面几项重合，从而使近似公式达到所需要的阶数。则可合，从而使近似公式达到所需要的阶数。则可构造出更高精度的计算格式，这就是龙格构造出更高精度的计算格式，这就是龙格库库塔（塔（Runge-KuttaRunge-Kutta）法的基本思想。）法的基本思想。 1, iixx2022-
32、5-2745只有多计算几个点的函数值，平均高度就更精确只有多计算几个点的函数值，平均高度就更精确.龙格龙格-库塔法的基本思想库塔法的基本思想函数值的线性组合，可以表示为上的定积分区间 1),(,1iixxiidxyxfxx),(yxfy 1),(1iixxiidxyxfyy多预报几个点的斜率值，将其加权平均作为平均斜率多预报几个点的斜率值，将其加权平均作为平均斜率.2022-5-27469.3.2 二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi i和和 , ,以该两点处的以该两点处的斜率值斜率值k1 1和和k2 2的加权平均的加权平均( (或称为线性组合或称为线性组合) )来求取平来
33、求取平均斜率均斜率k* *的近似值的近似值K，即，即 1,iixxphxxipi2211kkK式中式中: :k1 1为为xi i点处的切线斜率值，点处的切线斜率值， k2 2为为 点处的切线斜率值点处的切线斜率值, ,比照改进的欧比照改进的欧拉法拉法, ,将将 视为视为 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12phkyphxfkii对常微分方程初值问题对常微分方程初值问题(9.1)(9.1)式的解式的解 y= =y( (x),),根据微根据微分中值定理，存在点分中值定理，存在点 ，使得，使得 ),(1 iixx2022-5-2747)(,()(yfyK式中
34、式中 K K可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率。上的平均斜率。所以可得计算公式为：所以可得计算公式为： 1,iixxhKxyyii )(1)()(2211kkhxyi（9.14） 将将 在在x=xx=xi i处进行二阶处进行二阶TaylorTaylor展开：展开： )()(! 2)()()(321hOxyhxyhxyxyiiii （9 9.15） )()()(11iiiixxyxyxy也即也即 hKxyxyii)()(1（9.13）)(1 ixy2022-5-2748),()(12phkyphxfphxykiii )(),(),(),(),(22hOyxfyxf
35、yxfphyxfkiiyiiiixii)()()(2hOxyphxyii 将以上结果代入将以上结果代入)14. 9()()(22111kkhxyyii )()()()()(2211hOxyphxyxyhxyyiiiii )16. 9()()()()()(32221hOxyphxyhxyiii 展开处在将二元函数Taylorxphkyphxfkiii),(12 2022-5-2749)16. 9()()()()()(322211hOxyphxyhxyyiiii 对式对式(9.15)(9.15)和和(9.16)(9.16)进行比较系数后可知进行比较系数后可知, ,只要只要 )17.9(211221
36、 p成立成立, ,格式格式(9.14)(9.14)的局部截断误差就等于的局部截断误差就等于)(3hO)15. 9()()(! 2)()()(321hOxyhxyhxyxyiiii 2022-5-2750式式(9.17)(9.17)中具有三个未知量中具有三个未知量, ,但只有两个方程但只有两个方程, ,因而因而有无穷多解。若取有无穷多解。若取 , ,则则P P =1=1，这是无穷多解，这是无穷多解中的一个解，将以上所解的值代入式中的一个解，将以上所解的值代入式(9.14)(9.14)并改写并改写可得可得 2121 ),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不难发
37、现，上面的格式就是改进的欧拉格式。不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（凡满足条件式（9.179.17）有一簇形如上式的计算格式，）有一簇形如上式的计算格式，这些格式统称为二阶龙格这些格式统称为二阶龙格库塔格式。因此改进的库塔格式。因此改进的欧拉格式是众多的二阶龙格欧拉格式是众多的二阶龙格库塔法中的一种特殊库塔法中的一种特殊格式。格式。 2022-5-2751若取若取 , ,则则 ，此时二阶龙格，此时二阶龙格- -库塔库塔法的计算公式为法的计算公式为 0121, 12p 1,2 , 1 , 0 ni 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式
38、中 为区间为区间 的中点。的中点。 21 ix 1, iixx)2,(),(12121khyxfkyxfkiiii 21hkyyii 2022-5-27529.3.3 三阶龙格三阶龙格- -库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 pix )1( qpqhxxiqi并用三个点并用三个点 的斜率的斜率k1 1, ,k2 2, ,k3 3加权平均加权平均得出平均斜率得出平均斜率k* *的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式: : qipiixxx ,11 12233121(,)(,)iiiiiiyyhkkkkf x ykf xp
39、h yphk（9.18） 为了预报点为了预报点 的斜率值的斜率值k3 3, ,在区间在区间 内有两内有两个斜率值个斜率值k1 1和和k2 2可以用可以用, ,可将可将k1 1, ,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率, ,从而得到从而得到 的预报值的预报值 qix,qiixx )(qixyqiy 1122i qiyyqhkk,qiixx 2022-5-2753于是可得于是可得 ),(3qiqiyxfk 运用运用Taylor展开方法选择参数展开方法选择参数 , ,可以使格式可以使格式( (9.18)的局部截断误差为的局部截断误差为 , ,即具有即具有三阶精度三阶精度. .)
40、(4hO 613121112332223232121pqqpqp21321, qp满满足足：可可得得21321, qp。截断误差为截断误差为)(4hO2022-5-2754局部截断误差为局部截断误差为 , ,即具有三阶精度，这类格即具有三阶精度，这类格式统称为式统称为三阶龙格三阶龙格库塔方法库塔方法。)(4hO)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（9.19） 时，时，当当2, 1,61,64,61, 1,2121321 qp是其中的一种，称为库塔（是其中的一种，称为库塔（KuttaKutta）公式。）公式。 2
41、022-5-27559.3.4 四阶龙格四阶龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，如果需要再提高精度，用类似上述的处理方法，只需在区间只需在区间 上用四个点处的斜率加权平均作为上用四个点处的斜率加权平均作为平均斜率平均斜率k*的近似值，构成一系列四阶龙格的近似值，构成一系列四阶龙格库塔公库塔公式。具有四阶精度，即局部截断误差是式。具有四阶精度，即局部截断误差是 。 由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经典四阶经典RK公式公式。 ,qiixx )(5hO )22(6),()2,()2,(),(432113142213121
42、21kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（9.20） 2022-5-27569.3.5 9.3.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1) 计算步骤计算步骤 输入输入 , ,h,Nh,N 使用龙格使用龙格库塔公式（库塔公式（9.209.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n N N 结束。结束。 10,xx11, yx0101,yyxx2022-5-2757(2) (2) 四阶龙格四阶龙格库塔算法流程图库塔算法流程图 开 始 输 入 x0, y0,h , N 1 n x0 + h x1 f(
43、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 2022-5-2758(3)(3) 程序实现程序实现( (四阶龙格四阶龙格- -库塔法计库塔法计 算常微分方程初值问题算常微分方程初值问题) ) 例例9.4 9.4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 1)0(2y
44、xyy10 x解解: : 由四阶龙格由四阶龙格- -库塔公式可得库塔公式可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh41632. 0)0408. 1 , 2 . 0(),(3004fhkyxfkh2022-5-2759040811. 1)22(62 . 0432101kkkkyy可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为，从表中看到所求的数值解具有，从表中看到所求的数值解具有
45、4位有效数字。位有效数字。 2xey 龙格龙格库塔法的推导基于库塔法的推导基于Taylor展开方法，因而它展开方法，因而它要求所求的解具有较好的要求所求的解具有较好的光滑性。光滑性。如果解的光滑性如果解的光滑性差，那么，使用四阶龙格差，那么，使用四阶龙格库塔方法求得的数值解，库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的算法。时，应当针对问题的具体特点选择合适的算法。 2022-5-27609.3.6 变步长的龙格变步长的龙格-库塔法库塔法 在微分方程的数值解中，选择适当的步长是非常在微分方程的
46、数值解中，选择适当的步长是非常重要的。单从每一步看，重要的。单从每一步看，步长越小，截断误差就越步长越小，截断误差就越小；小；但随着步长的缩小，在一定的求解区间内所要但随着步长的缩小，在一定的求解区间内所要完成的完成的步数步数就就增加增加了了。这样会引起计算量的增大，这样会引起计算量的增大，并且会并且会引起舍入误差引起舍入误差的大量的大量积累与传播。积累与传播。因此微分因此微分方程数值解法也有选择步长的问题。方程数值解法也有选择步长的问题。 以经典的四阶龙格以经典的四阶龙格- -库塔法库塔法( (9.20)为例。从节点为例。从节点x xi i出发，先以出发，先以h为步长求出一个近似值，记为为步
47、长求出一个近似值，记为 ，由于局部截断误差为由于局部截断误差为 ，故有，故有 )(1hiy)(5hO5)(11)(chyxyhii当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。2022-5-2761然后将步长折半然后将步长折半, ,即以为即以为 步长步长, ,从节点从节点x xi i出发出发, ,跨跨两步到节点两步到节点x xi+1i+1, ,再求得一个近似值再求得一个近似值 , ,每跨一步的每跨一步的截断误差是截断误差是 , ,因此有因此有2h)2(1hiy52 hc5)2(1122)()(hcxyxyhii这样这样 161)()()(11)2
48、(11hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy由此可得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步的近似值，其误差可用步长折半前后两次计算结果的偏差长折半前后两次计算结果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy来判断所选步长是否适当来判断所选步长是否适当2022-5-2762当要求的数值精度为当要求的数值精度为时：时： （1 1）如果）如果，反复将步长折半进行计算，直，反复将步长折半进行计算，直至至为止为止, ,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如果）如果为止，并以上一次
49、步长的计算结果作为为止，并以上一次步长的计算结果作为 。 这种通过步长加倍或折半来处理步长的方法称为这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要变步长法。表面上看，为了选择步长，每一步都要反复判断反复判断，增加了计算工作量，但在方程的解，增加了计算工作量，但在方程的解y(x)y(x)变化剧烈的情况下，总的计算工作量得到减少，结变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。果还是合算的。1iy1iy2022-5-2763 9.4 算法的稳定性及收敛性算法的稳定性及收敛性9.4.1 稳定性稳定性 稳定性在微分方程的数值解法中是一个非常稳定性在微
50、分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误中，存在着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。稳定性问题。 2022-5-2764 当在某节点上当在某节点上x xi i的的y yi i值有大小为值有大小为的扰动
