1、1第七章 常微分方程的数值解法1 引言1、一阶常微分方程初值问题（微分方程加初值条件）(1)0 )(,yxbafy特征： 为已知的二元函数，只有一个变元 ，一阶导数 ，要求满足微分,(f xdxy方程且过点 的曲线，一个未知函数 （微分方程的解函数） 。)0y )(y2、解的存在性若 连续且对 满足 Lipschitz 条件，即存在常数 ，使对 ，有),(yxfy 0LRy21,，则初值问题(1)的解存在唯一，本章均假设 对 满足2121,Lf )(xfLipschitz 条件。3、数值解法问题：求 精确解（解析表达式）极其困难，实际应用中只要求数值解。)(xy数值解：在 取 一系列等距离散节
2、点，步长,ba 1210 nxx，求 的近似值 。1nxhn )(ny0hx01h02 hn0 hnxn)1(01)(y)()(y )(x)(y 012 ny1n方法：建立 的递推公式， ，从而 按节ny0)(x 10 nyy点排列顺序（步进式） 。22 Euler 方法1、 显示、隐式 Euler 法与梯形法对（1）的方程两端由 到 定积分得：nx1 1)(,)(1nxn dxyfyx从而有： （*）)(,)(1yhfyxn,n（积分中值定理 ）11 ),()(,)1 nnfhfdn （1）若取 ，得 （左矩形公式）nx(,),(1 nx xyfhyfn由（*）得： )ny， ,则有： )(
3、nny)(11n ,10),(1nyhfn称式为解初值问题的 Euler 法.（显格式）10几何意义：初值问题(1)的解曲线 过点 ，从 出发，以 为斜率作一)(xy),(0yP0x),(0yxf段直线，与直线 交点于 ，显然有 ，再从 出发，以 为斜1x,(1xP1hf1P,1f率作直线推进到 上一点 ，其余类推，这样得到解曲线的一条近似曲线，它就是折线22（某点导数值的几何意义为该点切线的斜率） 。210P（2）若取 ，得 （右矩形公式）1nx )(,),(11 nxc xyfhdyfn类似可得： 01yn称式为隐式 Euler 法（向后 Euler 法，隐格式）（3）取 （梯形公式）)(
4、,)(,2),( 11 nnxc xyfxyfhdfn则得： ,011 ynnn称式为梯形法.（隐格式）例 1 用 Euler 法、隐式 Euler 法、梯形法及改进 Euler 法求解 ，取1)0(, yxy，计算到 ，并与精确解比较.0h50x解：由于 , ,1),(yf .h0x1y3Euler 法： )1(1nnxyhy 1.0.9.0( nnxyhxy时， .0.0.9.01隐式 Euler 法： )11nn解出 ).0.(.)(1 nn xyhxyhy当 时，090.0梯形法： ),(),(211 nnn yxfyxfynh解出 2)()(1 hxyynn )21.0.9.(1nn
5、xy当 时，00476.12.9.改进 Euler 法:),(,(),(211 nnnn yxhfxfyxfhy即 1nxy)2(2)()( xynn1.095.0.x当 时， n 05.195.01 精确解： ，(梯形法效果最好，改进 Euler 较好：具有相同的误差数量级，其它不xey)(好!) 0xEuler 法 n隐式 Euler 法 ny梯形法 ny改进 Euler 法ny精确解 )(nxy0.1 1.000 000 1.009 091 1.004 762 1.005 000 1.004 8370.2 1.010 000 1.026 446 1.018 594 1.019 025 1
6、.019 7310.3 1.029 000 1.051 315 1.040 633 1.041 218 1.040 8180.4 1.056 100 1.083 014 1.070 097 1.070 802 1.070 3200.5 1.090 490 1.120 922 1.106 278 1.107 076 1.106 5312、隐式法 的计算（Euler 法及梯形法）1ny方法 1：显示化（ 对 线性时，以 为未知量的一元线性方程，见上例） fy1ny方法 2：迭代法（ 对 非线性时，可看作一个关于 方程，利用迭代法求解 ）1ny1nyEuler 法： （*1）,0),(11 xhfn
7、n（以 为未知量的一元非线性方程）ny4, n （*2）nknknyxfhy)0(1)(1)1( ,0梯形法： （*3）,10),(),21 yxffnn， （*4）),(,)0(1 )()1( nnkkhfyxy ,计算步骤：（1）初值 ；),(00)(1xfy迭代 （当 ）1)(1()( yykk )(1)(kky（2） ；迭代,11)0(hf 322) 收敛性：只要步长 足够小，就可保证迭代收敛。（a）当 （*2）收敛 （b）当 （*4）收敛.LLh证明：（b） 迭代序列 收敛 当 有 ；由（*3） 、 （*4）式知：)(1kny,k01)(nky（满足李氏条件）1)(1)(11)( 2
8、),2 nknknnk yxfxfhy，则当 ， （ ）时，有)0(kyL Lh，从而收敛得证。01)(1)(nkny3、 改进 Euler 法（梯形法隐式迭代式仅限迭代一次）,为避免迭代，可先用 Euler 法计算出 的近似 1ny1ny显式预报，隐式校正),(),(2111nnn yxfyfhy或 （6）),(11n hxfy称以上两式为改进 Euler 法（校正就是迭代一次） 右端已不含 ，显式方法。1ny例 2：用梯形法的迭代格式求 的数值解， ，计算到)0(,2 yxy .0h3.0x解：梯形公式 nnnnnnn yyhy 1.5.5.(2 1111 迭代式 ，解为：nnnnkkkn
9、 yxyxhfy2.0.1),(.5.05.)0(1)(1)(1)( 4509.836.32y53、 Euler 法的局部截断误差定义 1 在 的前提下，称 为在 的局部截断误差。)(nxy1)(nyx1nx定义 2 若一种数值方法的局部截断误差： ，则称这种方法是精)(phO度为 阶的。 （含 的项称为局部截断误差主项.）p1ph注：按某种方法由 算出 这一步的误差。ny1n（1）Euler 法（*） )(2)()()( nnnnn xyhxyhxy , ffy（*）)()()()( 1 nnnnnn xyh由(*)（*）式有： 221Oxyhx根据定义 2，Euler 法中的 故此方法为一
10、阶方法.p（2）隐式 Euler 法（分析时隐式公式右边的 ）)(11nnxy,)()(,)(, n fxyxfyxy )()()(1111 nnnnn xyhhh .)( nn（*）)()(2 nxyxy由（*）及（*）式： )()(221 hOxyn故局部截断误差主项为 ，也是一阶方法.),(2pxyhn（3）梯形法 )(!3)(2)()()( 1 nnnnnn xyhxyxy,21xffh类似有 ，方法是二阶的.1)(nyx )()(233hOyhn（4）改进 Euler 法),() yxfy)(,(),() xyfyxf则 nnynn ,( 1 hfhff 6 ),(),(),(),(
11、 nnynxn yxhxfhfyfh),(,(),(211 nnnn fff= ) xyhxy(2( nnnh )(!3)(2) 1 nnxyhxyxyhxy,方法是二阶的.(3On73 Runge-Kutta 方法一、 显式 Runge-Kutta 法的一般形式 一般形式： ),(,)()1,122211 rrnrnr nrn khbkyhaxfkfxkckchy二、 2 级显式 R-K 方法（r=2）设想构造 RK 公式： （*）),(21221kbyhaxfkcnn试确定 ，使上式为二阶格式，即 ：211,bac )(31hOyxn假设 ，对(7.3.6)式在 处按 Taylor 公式展
12、开，由于nyx)( ),(n；,)(, nn yxffy )(,)(,)()( nnyn xyfxfdxn,(),( ynxfff)(!3)(2)( 1 hxhy nnnn ,(2), yxfyfyfxhfxn )(),()()( 21212212 kbhaOxfkbakyk nnn 将上述结果代入(*)得： )(,(),()(),()1(),()()( 3212211 hOyxffcyfcyxhfcxy nnynxnn 要使公式为二阶,必须：（三个方程四个求知数，无穷多解）2121bca8（1）取 ，有 （改进欧拉法）1,221bac),(2112211hkyxfkynn（2）取 ，有 （中
13、点法）21,021bac)2,(1212khyxfknn取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.三、 四阶 R-K 方法及步长的自动选择类似可得到，经典的四阶 R-K 方法是：、)2(64311 kkhyn ),(1nyxf,2xfkn )23kh)(34它的局部截断误差 ，故 p=4,这是最常用的四阶 R-K 方法，数学库中都有用此方法)5hO求解初值问题的软件.这种方法的优点是精度较高，缺点是每步要算 4 个右端函数值，计算量较大.例 7.3 用经典四阶 R-K 方法解例 7.1 的初值问题 ，取 ，1)0(, yxy 1.0h，计算到 ，并与改进 Euler 法、梯形法
14、在 处比较其误差大小.解 用四阶 R-K 方法公式 (7.3.12)，此处 ，于是当n=0 时于是 ，按公式(7.3.12)可算出9此方法误差： 改进 Euler 法误差： 梯形法误差： 可见四阶 R-K 方法的精度比二阶方法高得多.用四阶 R-K 方法求解初值问题(1)精度较高，但要从理论上给出误差 的估计式nyx)(则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长 h?通常是通过不同步长在计算机上的计算结果近似估计.设 在 处的值 ,当 时，的近似为 ，于是由四阶 R-K 方法有5)(1(hcyxnn若以 为步长，计算两步到 ，则有2h15)2(1)(cyxnhn于是得： 6
15、)()2(1hn即 )5(1)2()( hnyyx或(7.3.13)(1)2()2(1( hnhn它给出了误差的近似估计.如果 ( 为给定精度)，则认为以 为步长的计算)(1)2(5hny 2h结果 满足精度要求，若 ，则还可放大步长 .因此(7.3.13) 提供了自动选)2(1hny )(1)2(hn择步长的方法.104 线性多步法一、 线性多步法的一般公式前面给出了求解初值问题(1)的单步法，其特点是计算 时只用到 的值（要计算多点1nyny处的函数值，计算量较大） ，此时 的值均已算出.如果在计算 时除用 的值nyy,110 1ny外，还用到 的值，这就是多步法.若记 , 为步长， 21
16、,rnny khxk0，则线性多步法可表示为（一般形式）： ),)()(kxfykk )( 1101110 rnnnrnn fffhaa 其中 为常数，若 ，称(7.5.1)为线性 步法.计算时用到前面已算出的i, 022rnr个值. 当 时，为显式方法，当 则称为隐式多步法.隐式方法r11,nyy 1与梯形方法一样，计算时要用迭代法求 多步法(7.5.1)的局部截断误差定义也与单步法类1ny似.线性多步法在 处的局部截断误差定义为：假定 时1nx 1,)(rnixyi 若 ，称其为是 p 阶的。)()(1pnhOy二、 Adams 显式与隐式方法（即 ）)( 11011 rnnnn fffy
17、 1,2,01 rnii，的 步法称为 Adams 方法，当 时为 Adams 显式方法，当 时，称为 Adamsr 1隐式方法.对初值问题(1)的方程两端从 到 积分得nx1 )(,)(1nn dxyfy选取不同的节点，对上式中被积函数进行不同的多项式插值，从而得到不同的的线性多步法。1、 Adams 四步显式公式（外插公式）选 ， ， ，）（ )(,nnxyfx ）（ )(,11nnxyf ）（ )(,22nnxyf四节点作三次多项式插值：）（ (,33 )(33RPx从而 1)(331nxn dRP取 与 为插值节点做 的插值多项式：23,x)(,xyf)(,)() 321 nnnn fP ,)(111n xyxxx)(,)( 23222 nnnn f
