第九章 常微分方程数值解

1、石家庄铁道学院 数理系Ch9常微分方程数值解常微分方程数值解 在自然界与工程界中,很多问题的数学表述都可以归结为解常微分方程和偏微分方程问题.但是,只有极少数特殊的方程有解析解。对于绝大部分的微分方程求其解析解是很困难的。因此研究微分方程的数值解是很重要的。 本章讨论常微分方程的数值解法欧拉法欧拉法,后退的欧拉法后退的欧拉法,梯形公式梯形公式;显式的显式的,隐式的方法隐式的方法,单步法多步法单步法多步法,局部截断误差局部截断误差;常用的四阶常用的四阶R-K方法方法另外,还有收敛性分析,稳定性分析,刚性方程,高性能算法等.石家庄铁道学院 数理系对于一个常微分方程：对于一个常微分方程：, , ),

2、(baxyxfdxdyy通常会有无穷个解。例如：通常会有无穷个解。例如：Rccxyxdxdy,)sin( )cos(因此，我们要加入一个限定条件。通常会在端点出给出，如下面的初值问题：因此，我们要加入一个限定条件。通常会在端点出给出，如下面的初值问题：0)(, , ),(yaybaxyxfdxdy为了使解存在唯一，一般，要加限制条件在为了使解存在唯一，一般，要加限制条件在 f 上，要求上，要求 f 对对y 满足满足Lipschitz条件：条件：2121),(),(yyLyxfyxf如何求解如何求解石家庄铁道学院 数理系解析解法解析解法：（常微分方程理论）：（常微分方程理论）只能求解极少一类常微

3、分方程；实际中给定的问题不一只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。定是解析表达式，而是函数表，无法用解析解法。计算解函数计算解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyyii 节点间距节点间距 为步长，通常采用为步长，通常采用等距节点等距节点，即取即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii数值解法数值解法: 求解所有的常微分方程求解所有的常微分方程步进式步进式：根据已知的或已求出的节点上的函数值计算：根据已知的或已求出的节点上的函数值

4、计算当前节点上的函数值，一步一步向前推进。因此只需当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。函数值的递推公式即可。石家庄铁道学院 数理系例：我们对区间做等距分割：例：我们对区间做等距分割： , ()/ixh ihbam设解函数在节点的近似为设解函数在节点的近似为iy由数值微分公式，我们有由数值微分公式，我们有iixxxxyxfdxdy),(，则：，则：),(1iiiiyxfhyy向前差商公式向前差商公式),( 1iiiiyxfhyy可以看到，给出初值，就可以用上式求出所有的可以

5、看到，给出初值，就可以用上式求出所有的iy这种方法这种方法 ，称为，称为数值离散方法数值离散方法。在一系列离散点列上，求未知。在一系列离散点列上，求未知函数函数y在这些点上的值的近似。在这些点上的值的近似。石家庄铁道学院 数理系为了考察数值方法提供的数值解，是否为了考察数值方法提供的数值解，是否有实用价值，需要知道如下几个结论：有实用价值，需要知道如下几个结论： 步长充分小时，所得到的数值解步长充分小时，所得到的数值解能否逼近问题得真解；即收敛性问题能否逼近问题得真解；即收敛性问题 误差估计误差估计 产生得舍入误差，在以后得各步产生得舍入误差，在以后得各步计算中，是否会无限制扩大；即稳计算中，

6、是否会无限制扩大；即稳定性问题定性问题石家庄铁道学院 数理系单步法：在计算yi+1 时只利用y i多步法：在计算yi+1 时不仅利用y i , 还要利用 yi1, yi2,k步法：在计算yi+1 时要用到yi,yi1,yik+1显式计算公式可写成:yk+1=yk+hf(xk，yk;h）隐式格式:yk+1=yk+hf（xk，yk,yk+1;h）它每步求解yk+1需要解一个隐式方程石家庄铁道学院 数理系1 .Euler公式公式做等距分割，利用数值微分代替导数项，做等距分割，利用数值微分代替导数项，imabxiI :向前差商公式向前差商公式)( 2)( )()(1nnnnyhxyhxyxy)( 2)

7、(,()()(1nnnnnyhxyxfhxyxy)( 2)(,()()(21nnnnnyhxyxhfxyxy于是得到于是得到Euler公式公式),(1nnnnyxhfyy称为局部截断误差。称为局部截断误差。显然，这个误差在显然，这个误差在逐步计算过程中会逐步计算过程中会传播，积累。因此传播，积累。因此还要估计这种积累还要估计这种积累nx1nx2nx石家庄铁道学院 数理系定义定义 在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断误差局部截断误差 /* local tr

8、uncation error */。定义定义 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。),(1nnnnyxhfyyEuler公式公式一阶一阶显式的显式的石家庄铁道学院 数理系2. 后退的后退的Euler公式公式)( 2)( )()(11nnnnyhxyhxyxy)( 2)(,()()(111nnnnnyhxyxfhxyxy)( 2)(,()()(2111nnnnnyhxyxhfxyxy),(111nnnnyxhfyy是隐式的，要迭代求解是隐式的，要迭代求解后退的后退的Euler公式公式一阶隐式的一阶隐式的预测预测-校正系统校

9、正系统: 用显式算法预测用显式算法预测, 再用隐式算法校正再用隐式算法校正),(),(111nnnnnnnnyxhfyyyxhfyynx1nxPQR石家庄铁道学院 数理系),(1nnnnyxhfyy),(111nnnnyxhfyy两式相加平均得到梯形公式梯形公式),(),(2111nnnnnnyxfyxfhyy预测预测-校正系统校正系统: 用显式算法预测用显式算法预测, 再用隐式算法校正再用隐式算法校正),(),(2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy校正预测2阶隐式的阶隐式的改进的改进的欧拉格式欧拉格式3.梯形公式梯形公式2阶隐式的阶隐式的单步法单步法石家庄铁道学

10、院 数理系1+1(,) (,)12pnnncnnpnpcyyh f xyyyh f xyyyy或几何解释几何解释xnxn+1ABPn+1=(A+B)/2尤拉法尤拉法后退尤拉法后退尤拉法梯形法梯形法石家庄铁道学院 数理系 基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到),(yxfy ,1npnxx1)(,()()(1ipnxxpnndxxyxfxyxy只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通，则可通过过 近似近似y(xn+1) 。而。而选用不同近似式选用不同近似式 Ik，可得到不，可得到不同的计算公式同的计算公式。1)(,(npnxxkdxxyxf

11、IkpnnIyy1Another point of view石家庄铁道学院 数理系所以，有格式为：),(),(2111nnnnnnyxfyxfhyy11 ()()( , ( )nnxnnxy xy xf x y x dx)( 12)(,()(,(2)()(2111fhxyxfxyxfhxyxynnnnnn局部截断误差梯形公式梯形公式隐式的隐式的对微分方程, , ),(baxyxfdxdyy做积分，则：1),(nnxxyxfdxdy例如例如石家庄铁道学院 数理系公式公式局部截断误差局部截断误差精精度度显显隐隐稳稳定定性性步数步数尤拉显尤拉显式公式式公式1 1阶阶显显差差单步单步尤拉隐尤拉隐式公式

12、式公式1 1阶阶隐隐好好单步单步梯形公梯形公式式2 2阶阶隐隐差差单步单步中点法中点法2 2阶阶显显好好二步二步3(3)3nhyx3(3)12nhyx2(2)2nhyx2(2)2nhyxsummary石家庄铁道学院 数理系4. RungeKutta法法由由Taylor展开展开)()!1()(!)( )()()1(1)(1nkknkknnnykhxykhxhyxyxy记为1nhT),()( yxfxy),(),()( yyxfyxfxyyx)( xy所以，可以构造格式所以，可以构造格式),(),(),(! 2),(21nnnnynnxnnnnyxfyxfyxfhyxhfyy这种格式使用到了各阶偏

13、导数，使用不便。这种格式使用到了各阶偏导数，使用不便。石家庄铁道学院 数理系从另一个角度看，从另一个角度看，11),(,()()(nnnnnhTfxyxhhFxyxy取取(x,y)及其附近的点做线性组合，表示及其附近的点做线性组合，表示F，问题就好办了。，问题就好办了。当然，要求此时的展开精度相同。这种方法称为当然，要求此时的展开精度相同。这种方法称为RungeKutta法。法。石家庄铁道学院 数理系在在(x,y)处展开，处展开，)(),(),(),(),(),(),(221221hOyxfyxhfbyxhfayxfcyxfcfyxhFyx比较比较11),(),(),(! 2),()()(nn

14、nnnynnxnnnnhTyxfyxfyxfhyxfhxyxy以以2阶为例，设阶为例，设),(,(),(),(21221yxhfbyhaxfcyxfcfyxhF石家庄铁道学院 数理系有：有：2/12/112122221bcaccc112/121221bacc改进的改进的Euler公式公式),(),( 2/121211hKyhxfKyxfKKKhyynnnnnn石家庄铁道学院 数理系常用的四阶常用的四阶RungeKutta方法方法),()2,2()2,2(),(22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn四阶四阶RungeK

15、utta方法的每一步需要四次计算函数值方法的每一步需要四次计算函数值f,可以证明其截断误差为可以证明其截断误差为O(h5)石家庄铁道学院 数理系2 Runge-Kutta Method注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的

16、导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。石家庄铁道学院 数理系方程组和高阶方程的数值解法方程组和高阶方程的数值解法bxaayayyyxfdxdyyyxfdxdymmmmmm , )()(),(),(111111写成向量的形式：写成向量的形式：)(),(aYYxFdxdY石家庄铁道学院 数理系各种方法都可以直接运用过来。各种方法都可以直接运用过来。bxazazyayzyxgdxdzzyxfdxdy , )()(),(),(0

17、0Euler公式公式),(),(11nnnnnnnnnnzyxhgzzzyxhfyy以两个方程的方程组为例以两个方程的方程组为例石家庄铁道学院 数理系Runge-Kutta公式公式112341(22)6nnnnyyhKKKKzz)2,2(12KhYhxFK1(1)(2)112(1)(2)11(,)(,)(,)222(,)222nnnnnnnnnnnnf xyzKg xyzhhhf xyKzKKhhhg xyKzK 石家庄铁道学院 数理系(1)(1)223(1)(1)22(1)(1)334(1)(1)33(,)222(,)222(,)(,)nnnnnnnnnnnnhhhf xyKzKKhhhg

18、xyKzKf xh yhKzhKKg xh yhKzhK石家庄铁道学院 数理系0.05 (1)0.002200.09 (1)0.1515(0)0.193(0)0.083duuuuvdtdvvvuvdtuv1、( , , )0.05 (1/20) 0.002( , , )0.09 (1/15) 0.15f u v tuuuvFg u v tvvuv2、确定方法，然后求解、确定方法，然后求解（0.20276 0.0881157）（0.213007 0.0934037）（0.223763 0.0988499）（0.235052 0.104437）（0.246902 0.110146）4阶阶Runge-Kutta法，法，h=1例例石家