计算方法 81 - 83 Euler格式和Runge-Kutta格式

1、则初值问题则初值问题( (* *) )的连续可微解的连续可微解y(x)存在存在, ,且唯一且唯一. . 第八章第八章 常微分方程初值问题的数值解法常微分方程初值问题的数值解法讨论常微分方程（讨论常微分方程（ordinary differential equation）的定解问题，这类问题）的定解问题，这类问题的最简单型是如下一阶方程的初值问题：的最简单型是如下一阶方程的初值问题：00( , ), (), .yf x yaxby xyy R ,: ybxaD12,y y, ),(),(2121yyLyxfyxf 满足满足Lipschitz (Lipschitz (李普希兹李普希兹) )条件，即存

2、在正数条件，即存在正数L L，使得对所有的，使得对所有的和任何和任何 均有均有定理定理 设函数设函数f(x,y)在区域在区域 上连续上连续,且在区域且在区域D内关于内关于y由于数值解是找精确解由于数值解是找精确解 的近似值，因此，的近似值，因此，总假设总假设方程的解方程的解 存在且唯一，并具有充分的光滑性！存在且唯一，并具有充分的光滑性！( )y x( )y x , xa b(*) 求常微分方程数值解的必要性求常微分方程数值解的必要性1 1、方程本身很复杂，不能给出解析解，或难以求出解析解；、方程本身很复杂，不能给出解析解，或难以求出解析解；2 2、即使可以获得解析解，计算量太大或者计算过程太

3、复杂而不实用；、即使可以获得解析解，计算量太大或者计算过程太复杂而不实用；例例如：高阶常系数线性常微分方程如：高阶常系数线性常微分方程. .3 3、在实际应用中，只需求得解在某些特殊点上的近似值。、在实际应用中，只需求得解在某些特殊点上的近似值。数值解数值解是指：在是指：在解的存在区间解的存在区间上取一系列离散节点上取一系列离散节点012nxxxx0 1 2(, , ,)iy i 逐个给出逐个给出精确解精确解 的近似值的近似值()iy x 定义定义：相邻两个节点的间距相邻两个节点的间距 称为称为步长步长.1,nnnhxx 如何求出？如何求出？考虑考虑等距等距节点：节点：0 0 1, ,nxxn

4、h n 从初始条件从初始条件 出发，按照节点的排序，依次逐个计算出发，按照节点的排序，依次逐个计算 的的 值，称为值，称为步进法步进法. 一般有两种类型：一般有两种类型：单步法单步法、多步法多步法.00( )y xy iy注：注：bahn 为了考察数值解是否具有使用价值，必须解决的基本问题：为了考察数值解是否具有使用价值，必须解决的基本问题：解析解不能用初等函数及其积分表示！解析解不能用初等函数及其积分表示！当步长当步长h h取得充分小时，取得充分小时，所得的近似解所得的近似解y yn n能否以足够的能否以足够的精度精度逼近逼近初值初值问题的精确解问题的精确解y(y(x xn n) ). .这

5、就是这就是收敛问题收敛问题。即当。即当h 0h 0时时, , y yn n y y( (x xn n) ?) ?在数值求解的过程中，会产生若干类型的误差，具体分类如下：在数值求解的过程中，会产生若干类型的误差，具体分类如下： （1）局部截断误差局部截断误差 （2 2）局部舍入误差）局部舍入误差 （3 3）整体截断误差）整体截断误差 （4 4）整体舍入误差）整体舍入误差 （5 5）总误差）总误差（= =整体截断误差整体截断误差+ +整体舍入误差）整体舍入误差）因此，必须因此，必须估计估计精确解与近似解之间的精确解与近似解之间的误差误差。这就是这就是误差估计问题误差估计问题。由于初始值由于初始值y

6、 y0 0和右端项和右端项f(f(x,yx,y) )常常是通过测量得到的，所以必须考虑常常是通过测量得到的，所以必须考虑 它们的微小扰动，引起数值解的变化问题。即最初产生的误差在以它们的微小扰动，引起数值解的变化问题。即最初产生的误差在以 后各步的计算中是否会无限制扩大的问题，这就是后各步的计算中是否会无限制扩大的问题，这就是稳定性问题稳定性问题。在计算过程中无舍入误差，只有当问题的数值解对初始值具有某种连续在计算过程中无舍入误差，只有当问题的数值解对初始值具有某种连续依赖性时，方法才实用！依赖性时，方法才实用！可证明：可证明：若不考虑初始值误差，整体截若不考虑初始值误差，整体截断误差的阶由局

7、部截断误差的阶决定！断误差的阶由局部截断误差的阶决定！本章着重讨论一阶本章着重讨论一阶ODE初值问题的数值解初值问题的数值解.对于高阶方程（组）对于高阶方程（组）的数值解，其基本思想是完全一样的的数值解，其基本思想是完全一样的.计算一步所产生的误差。是算法中所固有的，与舍入误差无关计算一步所产生的误差。是算法中所固有的，与舍入误差无关初值问题的初值问题的解析解解析解表示过点表示过点 的一条（的一条（光滑光滑）曲线曲线. .00(,)x y 解析解与其数值解的解析解与其数值解的几何意义几何意义：oxynx 0 x0 00 0( (x x, ,y y ) ) 1x 初值问题的初值问题的数值解数值解

8、表示一组表示一组离散点列离散点列 （或一组数据点）（或一组数据点）(,)iix y可用可用拟合拟合方法求该组数据方法求该组数据 的的近似曲线近似曲线( , )iix y2x11(,)x y积分曲线积分曲线近似曲线近似曲线( , )yf x y (,)nnxy22(,)xy本章给出的几种方法本章给出的几种方法一、欧拉一、欧拉(Euler)(Euler)方法及其改进形式方法及其改进形式二、龙格二、龙格- -库塔库塔( (Runge-KuttaRunge-Kutta) )方法方法三、线性多步法三、线性多步法-Adams-Adams（埃德姆斯）方法（埃德姆斯）方法 得到得到数值解数值解有两个基本途径有

9、两个基本途径： 把近似解表示成有限个独立函数之和把近似解表示成有限个独立函数之和，例如：截断的幂（，例如：截断的幂（TaylorTaylor）级数或正交函数展开式中的前几项级数或正交函数展开式中的前几项. . 涉及到计算高阶导，尽管可用涉及到计算高阶导，尽管可用“数数值微分值微分”技术，但得到的公式太长、太复杂！通常比较适用于手算技术，但得到的公式太长、太复杂！通常比较适用于手算. . 离散化方法离散化方法（也称为（也称为差分方法差分方法），它），它提供了提供了用当前节点上的或前几用当前节点上的或前几个节点上的近似值来个节点上的近似值来计算计算下一个节点上的近似值的下一个节点上的近似值的规则规

10、则. .数值解所满足的数值解所满足的离散方程统称为离散方程统称为差分格式差分格式. . 它是本章中要研究的一种方法它是本章中要研究的一种方法. .微微分分方方程程区域剖分区域剖分递推计算或解线性递推计算或解线性代数方程组代数方程组微分方程离散微分方程离散初始和边界条件处理初始和边界条件处理解的存在性、唯一性解的存在性、唯一性解的收敛性和收敛速度解的收敛性和收敛速度解的稳定性解的稳定性得到数值解得到数值解离散系统的离散系统的性态研究性态研究现实问题现实问题数学模型数学模型离散格式离散格式模型模型误差误差舍入舍入误差误差观测观测误差误差截断截断误差误差数值解数值解数学模型数学模型离散化离散化计算计

11、算第第2 2节节 欧拉（欧拉（EulerEuler）方法）方法 最简单的一种方法，精度差，不推荐使用！最简单的一种方法，精度差，不推荐使用！ 欧拉格式的构造欧拉格式的构造解决问题的关键：解决问题的关键：如何处理方程中的导数项？如何处理方程中的导数项？00()(, () ,()nnny xf x y xy xy 在各节点在各节点 处处, ,有有求求 的近似值的近似值()ny xny方法：方法：将上述初值问题将上述初值问题化成化成节点离散方程，在节点上采用节点离散方程，在节点上采用离散化方法离散化方法（也也叫做差分方法叫做差分方法, ,通常用数值积分、微分、泰勒公式等）通常用数值积分、微分、泰勒公

12、式等），可逼近，可逼近节点离散方节点离散方程程，由此产生，由此产生可计算格式可计算格式，并用计算解，并用计算解 作为解析解作为解析解 的近似值的近似值. .iy( )iy x在在节点离散方程节点离散方程中直接用中直接用向前差商向前差商代替微商，得到代替微商，得到nx节点离散方程节点离散方程11()()(, ()nnnnnny xy xf x y xxx 1()()(, ()nnnny xy xhf x y x 101(,), , ,nnnnyyhf x yn 称为称为欧拉格式欧拉格式.当初值当初值 已知时，已知时，0y可递推求出可递推求出12,.,.nyyy( )iiyy x 令令bahn 该

13、切线与直线该切线与直线 的的交点交点的纵坐标的纵坐标现在从现在从 出发出发00(,)xyoxy),(00yx0 x),(00yx nx (,)nnxy 2x),(22yx 1x),(11yx 作解曲线的切线，作解曲线的切线，切线方程为：切线方程为：0000(,)()yyf xyxx 切线切线斜率斜率为为1y100010(,)()yyf xyxx 11( ,)x y1111(,)()yyf x yxx 欧拉方法的几何解释欧拉方法的几何解释00(,),f x y1xx 11( )y xy 再从再从 出发，以出发，以 为斜率作直线为斜率作直线11( ,)f x y该直线与该直线与直线直线 的交点的的

14、交点的纵坐标纵坐标2xx 2y211121(,)()yyf xyxx22()y xy 依次类推，依次类推， 的近似值的近似值 的计算公式：的计算公式：1()ny x 1ny 1(,)nnnnyyhf x y 即是即是欧拉公式欧拉公式.又称又称用一条折线近似代替积分曲线折线法折线法.欧拉方法欧拉方法的缺陷：误差比较大！的缺陷：误差比较大！( )y y x 令令令令( , )yf x y 方法可看做是方法可看做是TaylorTaylor级数前两项的近似！级数前两项的近似！y(x1)y(x2)y2y(x3)y3y1x1x2x3xn.x0y(x0)EulerEuler格式的整体截断误差与局部截断误差，

15、示意图格式的整体截断误差与局部截断误差，示意图整体截断误差局部截断误差在在节点离散方程节点离散方程中直接用中直接用向后差商向后差商代替微商，得到代替微商，得到11()()(, ()nnnnnny xy xf x y xxx 1()()(, ()nnnny xy xhf x y x 101(,), , ,nnnnyyhf x yn 称为称为隐式欧拉格式隐式欧拉格式. .在在节点离散方程节点离散方程中用中用中心差商中心差商代替微商，得到代替微商，得到1111()()(, ()nnnnnny xy xf x y xxx 112()()( , ( )nnnny xy xhf x y x 11201(,

16、), , ,nnnnyyhf x yn 称为称为欧拉两步格式欧拉两步格式.于是于是11101(,), , ,nnnnyyhf xyn 关于关于 的非线性方程，可用的非线性方程，可用迭代法求解！迭代法求解！1ny 欧拉格式的欧拉格式的局部截断误差局部截断误差假设假设 ，即第，即第n步的结果步的结果 是准确的，是准确的， 称为称为局部截断误差局部截断误差.()nnyy x ny11()nnyy x 定义定义2.1. 2.1. 如果一种方法的局部截断误差为如果一种方法的局部截断误差为 ，则称该数值方法的精度，则称该数值方法的精度 是是 阶阶的，或简称该方法是的，或简称该方法是 阶阶的的. .1()p

17、O h p11(,) ()(, ()() ()nnnnnnnnnnfyxy xyhf xyyhy xy xy xh 显式显式EulerEuler格式格式132112( )()(Tayl)()+or()nnnnny xy xhyxO hhxxy 根根据据公公式式，有有，13221112()( ) ( ) ()nnnnyy xhO hO hxxy 因因此此有有，计算一步所产生的误差计算一步所产生的误差一阶方法一阶方法( , )yf x y 局部截断误差的主项局部截断误差的主项p1111111(,) (,) ()()()nnnnnnnnnnf xyyhf xyyhyyxyhxxy 隐式隐式Euler

18、Euler格式格式Taylor根根据据公公式式，有有 1111232312111223()()()() ()+()()! =()nnnnnnyy xy xhyxy xh yxh yhOyh 因因此此有有1232211132()()+()!()nnnnnnyy xhyxhxxxxh yy ，2111112()+()()()nnnnnyxh yyxhyxxx ，111(),()(,) nnnnnyy xyxf xy 令令，则则一阶方法一阶方法 Euler两步格式两步格式11111222(,) (, () )()()nnnnnnnnnnf xy xyhf xyyhy xhxyyyx Taylor根根

19、据据公公式式，有有3343111133()()()() ! (=)nnnnh yxh yxO hOyy xh 因因此此有有13241213()()+()()!()()nnnnnhy xhyxh yxO hxy xy 13421312()+()()(!)(nnnnnh yy xhyxxy xh yxO h 11(),( )nnnnyy xyy x 令令，则则二阶方法二阶方法在区间在区间 上对微分方程上对微分方程 积分得积分得 欧拉格式的欧拉格式的积分学解释积分学解释( )( , ( )y xf x y x 1,kkxx 11( ) d( , ( ) dkkkkxxxxy xxf x y xx 即

20、即11()()( , ( ) dkkxkkxy xy xf x y xx 为了给出迭代格式，只要对积分项提供一种算法！而选择不同的算法，就为了给出迭代格式，只要对积分项提供一种算法！而选择不同的算法，就会得到不同的迭代格式会得到不同的迭代格式. . 利用左矩形公式，有利用左矩形公式，有21()()(, ()()kkkky xy xhf xy xO h 略去高阶项，便得略去高阶项，便得（显式）欧拉公式（显式）欧拉公式. 利用右矩形公式，有利用右矩形公式，有2111()()(, ()()kkkky xy xhf xy xO h 略去高阶项，便得略去高阶项，便得（隐式）欧拉公式（隐式）欧拉公式.略去

21、高阶项，然后用略去高阶项，然后用 代替代替 ，得，得利用梯形公式，有利用梯形公式，有11()()( , )kkxkkxy xy xf x y dx 1132(, ()(, ()()kkkkhfO hxy xf xy x()iy xiy 1112(,)(,)kkkkkkhyyf xyf xy 单步隐单步隐格式格式称为称为梯形格式梯形格式. .二阶方法二阶方法关于关于 的非线性方程，可用的非线性方程，可用简单迭代法求解！简单迭代法求解！1ky 例如：例如： 可构造可构造不动点迭代格式不动点迭代格式1111 20 1()( )(,),(,) ,nnnnnnkkhyyf xyf xky 其中迭代函数为

22、其中迭代函数为 1112(,)()(,)nnnnnnhyf xfyxyy 为了保证迭代法收敛，我们分析为了保证迭代法收敛，我们分析111111111112( )()( )()( )|(,)(,)| |()()|nnnnnnnnkkkkkkyyhf xf xyyyy 1112( )()|nnkkL yyh 于是，对任意正整数于是，对任意正整数p，有，有111111111121()( )()()()()()| |() ()(| )nnnnnnnnk pk pkk pkk pk pkyyyyyyyy 1111111121()()()()()( )|nnnnnnk pk pkk ppkkyyyyyy

23、1111111112112111221222( )()( )()()( )()( )|nnnnnnnnkkkkkppkkppkhhLyLyhhhLLyyyyyLy 111112()( )|nnkkyhyL 11101122( )( )|nnkLyyhhL 012hL 只只要要这表明：这表明：序列序列 是是CauchyCauchy列！列！ 10( )nkkky 将将显式欧拉公式显式欧拉公式和和梯形公式梯形公式联合使用联合使用11112 (,)(,)(,)nnnnnnnnnnyhf xyhyyyyf xyf x 即先用欧拉格式即先用欧拉格式预估预估一个粗糙的近似值，然后用梯形公式对其进行一个粗糙的

24、近似值，然后用梯形公式对其进行校正校正.这种这种预估预估-校正格式校正格式称为称为改进的改进的Euler格式格式.它可表为如下等价形式：它可表为如下等价形式：12112()nnkykyh 1(,)nnf xyk 0 1, ,n 12(,)nnf xh ykkh 单步显单步显格式格式是是 的斜率值的斜率值.nx是是 的斜率值，它是利用的斜率值，它是利用已知信息已知信息 通过通过Euler格式格式得到的得到的.1nx ny EulerEuler格式的格式的微分学解释微分学解释依据微分中值定理，依据微分中值定理，有有11 ()( ),( ,().( )nnnny xy xhxxhyfy 1*( ,

25、( ),( ),nnKfyy xxx 称称之之为为在在区区记记平平间间上上的的均均斜斜率率. .是是 两点的斜率值两点的斜率值 的算术平均作为的算术平均作为平均平均斜率斜率 的近似值的近似值.1,nnxx 12,kk*K局部截断误差局部截断误差为为O( (h3) )二阶方法二阶方法一般地，我们考虑：一般地，我们考虑：如果设法在如果设法在 内多预报几个点的斜率值，然内多预报几个点的斜率值，然后将它们加权平均作为后将它们加权平均作为平均斜率平均斜率 的近似值，则有可能构造出更高精度的近似值，则有可能构造出更高精度的计算格式的计算格式 这就是这就是Runge-KuttaRunge-Kutta方法的基

26、本思想方法的基本思想. .1,nnxx *K第第3 3节节 Runge-KuttaRunge-Kutta方法方法 - - 一种提高截断误差阶的方法一种提高截断误差阶的方法 二阶二阶Runge-KuttaRunge-Kutta格式格式121122111221 01 *,()n pnnn pnnnnxxphpxxKKxKKKyyhxKK 在在区区间间内内选选取取一一点点并并用用两两个个点点的的斜斜率率值值的的加加权权平平均均平平均均斜斜率率的的 近近似似值值，使使所所构构造造出出的的计计算算格格式式问问题题的的精精确确解解问问，题题且且 与与与与作作为为逼逼近近3 3具具有有二二阶阶精精度度. .

27、12(),(,)：n pn pnn pn py xyyhKKf xyp 改改进进的的E Eu ul le er r格格式式仿仿照照，用用E Eu ul le er r格格式式预预测测的的值值并并用用它它来来估估计计斜斜率率，于于是是得得到到如如下下计计算算格格式式：11211221 ()(,),(,).nnnnn pnyyhKKKf xyKf xyphK 1)()n pnn pny xy xKxx 问题：问题：适当选取待定适当选取待定参数，使得格式具有参数，使得格式具有二阶精度二阶精度. .工具：工具：用用TaylorTaylor公式，公式，分析格式的局部截断分析格式的局部截断误差误差nx1n

28、x n px ( , )yf x y 为此，考虑格式的局部截断误差为此，考虑格式的局部截断误差. . 1211212 =(,)(,)(,)(,)(,)(,)()(, ()()(,)()nnn pnnnnnnnnxnnynnnnnnnnKf xyKf xyhKf xph yphf xyph fxyfxyOyy xKKxy xpKf xyyxhf xy x 注注意意到到：设设，分分别别将将与与在在处处作作T Ta ay yl lo or r展展开开，有有 112122312122122 = 12=(,)()(,)()()()()()()()()()()()()nnxnnnynnnnnnnnnnnn

29、nfxyy xfxyKKy xy xphyxO hy xy xyxy xy xyxyyhy xhhhphO hy xhh 于于是是和和3)()O h 比较系数可知，若要格式具有二阶精度，只需比较系数可知，若要格式具有二阶精度，只需122112,.p 易易见见，满满足足条条件件的的参参数数不不止止一一组组，而而是是一一族族. .二二阶阶把把满满R Ru u足足这这一一条条n ng ge e- -件件的的一一族族格格式式统统称称K Ku ut tt t为为a a格格式式. .( , )yf x y 特别，若取特别，若取 有有 格式就是格式就是改进的改进的EulerEuler格式格式. .12211

30、2,.p 1,p 1212. 若取若取 有有 格式称为格式称为变形的变形的EulerEuler格式（也称为中点格式）格式（也称为中点格式）. .其形式为其形式为12,p 1201,. 1212112 2(,),(,).nnnnnnyyhKKf xyhKf xyK 三阶三阶Runge-KuttaRunge-Kutta格式格式12311223311122313 1 *,(,)，nn pn qnnn pn qnnnnxxxxqh pqxxxKKKKKKKyyxxhKKK 在在区区间间内内除除了了和和外外再再选选取取一一点点并并用用三三个个点点的的斜斜率率值值的的加加权权平平均均平平均均斜斜率率 ,

31、,问问题题4 4的的近近似似值值，使使所所构构造造出出的的计计算算格格式式问问题题的的精精确确解解，且且具具有有三三, ,作作为为逼逼近近阶阶精精度度. .112112212, ,(,)(),(,)., ,()(),：n pnnn pnnn pn pnn pn pn qnnn qn qn qnxxhx xKf x yy xyyhKKf xyxxhx xK Ky xyyhpqKpqK ，点点此此时时，在在区区间间内内有有一一个个斜斜率率值值用用E Eu ul le er r格格式式预预测测的的值值并并用用它它来来估估计计斜斜率率，点点此此时时，在在区区间间内内有有两两个个斜斜率率值值，同同样样，

32、用用E Eu ul le er r格格式式预预测测注注意意的的值值, ,到到注注意意到到, ,并并用用它它来来3(,).n qn qKf xy 估估计计斜斜率率nx1nx n px n qx 123123123,K K KKKK 利利用用这这三三个个斜斜率率值值作作加加权权平平均均得得出出上上的的平平均均斜斜率率，得得到到如如下下计计算算格格式式：123121123121312 ()(,),(,),(,().nnnnn pnn qnyyhKKKKf xyKf xyhKKf xyhKKpq 问题：问题：适当选取适当选取7 7个个待定参数，使得格式待定参数，使得格式具有三阶精度具有三阶精度. .类似地，利用类似地，