120第6章 常微分方程初值问题数值解法6.1 问题的描述和基本概念1、常微分方程初值问题● 一般形式0(,)()y f x y y a y '=⎧⎨=⎩式中(,)f x y 已知，0()y a y =称为初值条件. ● 初值问题的数值方法和数值解求函数()y y x =在若干离散点k x 上的近似值(0,1,)k y k = 的方法称为初值问题的数值方法，而称(0,1,)k y k = 为初值问题的数值解.1212. 建立数值解法的思想与方法用离散化方法将初值问题化为差分方程, 然后再求解.设节点为011n n a x x x x +=<<<<<距离1k k k h x x +=-称为步长.求数值解一般是从0y 开使逐次顺序求出12,,y y . 初值问题的解法有单步法和多步法两种： ● 单步法：计算1k y +时只用到k y 一个值；● 多步法：计算1k y +时要用1,,,k k k l y y y -- 多个值。

数值解法还有显格式和隐格式之分。

122微分方程离散化方法主要有数值微分法，数值积分法和Taylor 展开法1) 数值微分法由'()(,())k k k y x f x y x =，用数值微分的2点前差公式代替'()k y x ，得近似离散化方程记1k k h x x +=-，做k k ，“”，得差分方程1(,)k k k k y y f x y h+-= 即 1(,)k k k k y y hf x y +=+（Euler 公式）由初值条件0()y y a =及Euler 公式可求出数值解12,,,,n y y y .Euler 公式是显式单步法.1232）数值积分法在1[,]k k x x +上对'(,)y f x y =两边取定积分，得 111()()'(,())kk k k x x k k x x y x y x y dx f x y x dx +++-==⎰⎰右端积分用左矩形公式（数值积分公式）得1()()(,())k k k k y x y x hf x y x +-≈于是得到求初值问题的Euler 方法1(,)k k k k y y hf x y +=+124右端积分用右矩形公式（数值积分公式）得111()()(,())k k k k y x y x hf x y x +++-≈于是得到求初值问题的后退Euler 方法1+1+1(,)k k k k y y hf x y +=+后退Euler 方法是隐式的.右端积分用梯形公式（数值积分公式）得近似离散化方程：于是得到求初值问题的梯形方法该公式是隐式单步法.1251263）Taylor 展开法因为初值问题中函数(,)f x y 是已知函数，由(,)y f x y '=，可以计算''y ，'''y ，…， 于是有函数()y y x =在k x 处的Taylor 展式取上式右端前若干项，得近似离散化方程. 例如取前两项有1()()(,())k k k k y x y x hf x y x +≈+于是又得到Euler 公式：1(,)k k k k y y hf x y +=+.1273. 数值解法的误差、阶与绝对稳定性单步法数学描述为11(,,,)k k k k k y y h x y y h ϕ++=+显式：1(,,)k k k k y y h x y h ϕ+=+其中(,,)x y h ϕ称为增量函数.128显式单步法的一些概念定义1 称111()k k k e y x y +++=-为单步法在节点1k x +的整体截断误差，而称11()()(,(),)k k k k k T y x y x h x y x h ϕ++=--为在1k x +点的局部截断误差。

()k y x 表示解()y x 在k x 的值，是准确值，没有误差； k y 表示由数值解公式得出()k y x 的近似值，是数值解，有截断误差.129 局部截断误差1k T +的理解假设在计算()k y x 时没有误差（()k k y y x =）下,计算出的1k y +（1()(,(),)k k k k y y x h x y x h ϕ+=+）与1()k y x +的误差()111k k k T y x y +++=-(计算一步的误差）.定义2 如果数值解法的局部截断误差为11()P k T O h ++=则称该方法具有p 阶精度或该方法是p 阶方法.方法的阶越高，方法越好.130局部截断误差的主项如果某方法是p 阶方法，11()P k T O h ++=按h 可展为 1121()(,())()P P P k k k T O h g x y x h O h ++++==+ 则称1(,())P k k g x y x h +为局部截断误差的主项.在同阶方法中，局部截断误差的主项越小，方法越好.对Euler 方法1(,)k k k k y y hf x y +=+，有 1()()(,())k k k k k T y x h y x hf x y x +=+-- 将()k y x h +在k x 点展开，有2()()'()''()2!k k k k hy x h y x hy x y x +=+++2()(,())''()2!k k k k h y x hf x y x y x =+++ 故有231()().2k k y x T h O h +''=+Euler 方法是一阶方法.131例1 试求梯形方法的阶和局部截断误差主项. 解 该单步公式的局部截断误差是111()()((,())(,()))2k k k k k k k hT y x h y x f x y x f x y x +++=+--+ 1()()(()())2k k k k h y x h y x y x y x +''=+--+ 23()()()(()()23!2k k k k k h h hy x h y x y x y x y x ''''''''=++-+ 24()())()2k k h y x h y x O h '''''+++34()().12k h y x O h '''=-+ 故局部截断误差主项是3()12k h y x '''-，方法是二阶的.132定义3 设某种数值方法在k y 上大小为δ的扰动，于以后各()n y n k >上产生的偏差均不超过δ，则称该数值方法是稳定的。

通常用试验方程'y y λ= （λ为复数）来讨论求解数值方法绝对稳定性.Euler 方法稳定性将Euler 公式用于试验方程'y y λ=，得到1(1)k k k k y y h y h y λλ+=+=+设计算k y 时有误差,k ε则有11(1)()k k k k y h y ελε+++=++得1(1)k k h ελε+=+要想1k k εε+≤，只须11h λ+≤，因此Euler 方法在11h λ+≤时是绝对稳定的，其绝对稳定域为复平面h λ上以(-1,0)为中心的单位圆盘.绝对稳定区间为20.h λ-≤≤1336.2 Runge-Kutta 方法()11111,,(2,3,,)m k k i i i k k r rk r k rj j j y y h c K K f x y K f x a h y h b K r m +=-=⎧⎪=+⎪⎪=⎨⎪⎛⎫⎪=++= ⎪⎪⎝⎭⎩∑∑称为m 级R-K 方法.增量函数是()()1,,,,mi i i x y h c K x y h ϕ==∑134构造过程以2m =来说明Runge-Kutta 方法的构造方法和过程，对一般的Runge-Kutta 方法可类似处理.2m =的Runge-Kutta 公式为()11122k k y y h c K c K +=++式中 ()1,k k K f x y =，()22211,k k K f x a h y hb K =++.由(),y f x y '=，可得()()()(),(),(),()x y y x f x y x f x y x f x y x ''''=+ 在k x 处做Taylor 展开，有()()()()()()(2322!(,())(,())2!k k k k k k k x k k y x y x h y x y x h h O h h y x hf x y x f x y x '''+=+++'=++)()3(,())(,())y k k k k f x y x f x y x O h '++对()(),,k k x y x h ϕ在(,())k k x y x 做二元Taylor 展开，有()()(12,,(,())(,())k k k k k k x y x h c f x y x c f x y x ϕ=+21(,())(,())y k k k k hb f x y x f x y x '+)22(,())()x k k a hf x y x O h '++12()(,())k k c c f x y x =+(221(,())(,())y k k k k c b f x y x f x y x '+)22(,())()x k k a f x y x h O h '++由135()()()()1,,k k k k k T y x h y x h x y x h ϕ+=+--， 有()1122211(,())()(,())2k k k x k k T c c f x y x h c a f x y x +⎛'=--+- ⎝ 232211()(,())(,())()2y kk k k c b f x y x f x y x h O h ⎫'+-+⎪⎭选1222221111,0,022c c c a c b +=-=-=有局部截断误差()3T O h =，这样可得到二阶Runge-Kutta 公式.取20c t =≠，则式（6.13）的解为11c t =-，22112a b t==取不同的t 可得出不同的二阶Runge-Kutta 公式.如取12t =时，得到改进的Euler 公式()()()1,,,2k k k k k k k k h y y f x y f x h y hf x y +⎡⎤=++++⎣⎦ 1t =时，得到中点公式1(,(,)).22k k k k k k h h y y hf x y f x y +=+++136经典Runge-Kutta 公式()()()112341213243226,,22,22,k k k k k k k k k k hy y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎛⎫=++⎨ ⎪⎝⎭⎪⎪⎛⎫=++⎪ ⎪⎝⎭⎪⎪=++⎩四阶方法.137例1 设初值问题为()100y y y '=-+⎧⎪⎨=⎪⎩()00.4x ≤≤ 分别用Euler 方法（0.025h =），改进Euler 方法（0.05h =）和经典Runge-Kutta 方法（0.1h =）计算。