6.1 初值问题的Euler方法
设一阶常微分方程初值问题 dy f ( x, y ) dx y ( x0 ) y0 记 xn a nh (n 0,1,2,...), h 0为步长，一般总 假定为常数。该式的数值解是指通过某种方法 去获得解y ( x)在点xn 上的近似值yn , 即 y ( xn ) y n (n 0,1,2,...)
误差概述
一般说来，一个算法，局部截断 误差阶ｐ 越大， 则精度相对的越高。 h2 对于显式Euler法：Tn 1 y( xn ) O(h 3 ) 2 h2 对于隐式Euler法：Tn 1 y( xn ) O( h 3 ) 2 h3 ( xn ) O (h 4 ) 对于梯形公式：Tn 1 y 2
这就是隐式的`Euler 公式或向后Euler方法，它与显式 的不同在于，它每算一步要解函数方程(2)才能得到 yn 1。
初值问题的Euler方法
如果取以上两式的算术平均值的结果，则得 h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] （n 0,1,2,...） 2 称为梯形公式。 计算yn时常用以下迭代式：
~ ~
误差概述
定义6.1. 如果给果给定的算法的截断 2 误差为 Tn 1 O(h p 1 )
p 1
则称该算法具有p阶精度。如果 Tn 1 g ( xn ,y ( xn )) h O(h
p2
)
则非零项g ( xn , y ( xn )) h p 1称为为局部截断误差主。
1 yn 1 yn (k1 k 2 ) 2 2 xn 改进的Euler法： k1 0.1( yn ) yn 2( xn 0.1) k 2 0.1( yn k1 ) yn k1 计算结果如下表所示：
(n 0,1,2,...)
初值问题的Euler方法
对于( 2)计算yn 1 ,由于迭代工作量较大, 一般只 迭代一次, 构成一类预估 校正算法, 即
( p) yn 1 yn hf ( xn , yn ) h (c) ( p) yn ` yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] 2 ( 并取yn 1 ync )1。
6.1.2 误差概述
显式单步法一般形式为 yn 1 yn h ( xn , yn , h) 而隐式单步法一般形式为 yn 1 yn h ( xn , yn , xn 1 , yn 1 , h) 函数与f ( x, y )有关，称为增量函数。
误差概述
定义6.1.1 从初值y ( x0 ) y0出发，由单步法显式
或隐式逐步计算，得xn 1的值yn 1 , 则en 1 y ( xn 1 ) yn 1 称为在点xn 1上的整体截断误差。如果第n步在点xn的 值计算没有误差，即yn y ( xn ),由单步法计算出 y n 1 , 则Tn 1 y ( xn 1 ) y n 1 , 称为点xn 1上的局部截断误差。
于是 yn 1 y ( xn ) (1 2 )hy( xn ) [ 2 f x( xn , yn ) 2 f ( xn , yn ) f y ( xn , yn )]h 2 O(h3 )
Runge-Kutta方法
与Taylor展式相比较得 : 1 2 1 1 2 2 2 1 2 由于四个参数，三个方程，因此有一个 自由参数，即解答不唯一。
这就是显式的`Euler公式，它可以从y0出发，逐次
初值问题的Euler方法
如果用xn 1代替x0 , 于是该式可离散为： y ( xn h ) y ( xn ) f ( xn 1 , y ( xn 1 )) h 以yn 表示y ( xn )的近似值，则有 yn 1 yn hf ( xn 1 , yn 1 ) （n 0,1,2,...） ( 2)
初值问题的Euler方法
例6.1.1 设初值问题
2x dy y y dx y (0) 1 试分别用Euler法和改进Euler法求解, 并与 精确解y 1 2 x 进行比较。
初值问题的Euler方法
解：取h 0.1, 计算x [0,1]上结果，此时 2 xn Euler法：yn 1 yn 0.1( yn ) yn (n 0,1,2,...)
( yn0 )1 yn hf ( xn , yn ) （3） h ( k 1) (k ) yn 1 yn 2 [ f ( xn , yn ) f ( xn 1 , yn 1 )] （k 0,1,2,...） ( ( 当 | ynk 1) ynk ) | 1 ( 时, 取yn 1 ynk 1) 1
证明：由式（2）和（3）有 | y n 1 y
(k 1) n 1
h ( ) | | f ( xn 1 , yn 1 ) f ( xn 1 , ynk 1 ) | 2 hL ( ) | yn 1 ynk 1 | 2 ......
hL k 1 ( ( ) | yn 1 yn0 )1 | 2 hL k 1 ( 由假设知 : lim ( ) 0 , 故有 lim ynk 1) yn 1。 1 k 2 k
初值问题的Euler方法
定理6.1.1 设函数f ( x, y )对变量y满足Lipschitz
hL 条件，L为Lipschitz常数。如果步长h满足0 1, 2 2 ( ) 即h 时，则由（3）产生的序列{ ynk 1}( k 0,1,2...) L 收敛。
初值问题的Euler方法
1 1 (1) 取1 , 可得2 , 1, 此时算式为 2 2 1 yn 1 yn 2 (k1 k 2 ) k1 hf ( xn , yn ) k hf ( x h, y k ) n n 1 2 这是改进的Euler方法。
~ ~
数值稳定性分析


定义6.1.3 若某数值算法的绝对稳定性区 域包含hλ平面上的左半平面Re(hλ)＜0， 则称该方法是A稳定的。 隐式Euler法是A稳定的。
6.2 Runge-Kutta方法
受改进的Euler方法启发，更一般算式可设为 yn 1 yn 1k1 2 k 2 (n 0,1,2,...) k1 hf ( xn , yn ) k hf ( x h, y k ) n n 1 2 适当选择参数1，2，，，使局部截断误差 Tn 1 y ( xn 1 ) yn 1 O(h ), 这里仍假定yn y ( xn )。
1 3 2 (3) 取1 , 可得2 , , 又有算式 4 4 3 1 （ yn 1 yn 4 k1 3k 2 ) k1 hf ( xn , yn ) 2 2 k 2 hf ( xn h, yn k1 ) 3 3 这也是二阶R - K方法。
初值问题的Euler方法
上式还常写成 1 yn 1 yn 2 f (k1 k 2 ) k1 hf ( xn , yn ) k hf ( x h, y k ) (n 0,1,2,...) n n 1 2 该式称为改进Euler方法, 亦可写成 h yn 1 yn [ f ( xn yn ) f ( xn 1 , yn hf ( xn , yn ))] 2
Runge-Kutta方法
1 (2) 取1 0, 可得2 1, , 此时算式为 2 yn 1 yn k 2 k1 hf ( xn , yn ) 1 1 k 2 hf ( xn h, yn k1 ) 2 2 这是二阶R - K方法.
初值问题的Euler方法
为实现这一目标，Euler方法首先将微分算子离 散化，并用xn 代替x0 , 于是该式可离散为： y ( xn h ) y ( xn ) f ( xn , y ( xn )) h 以yn 表示y ( xn )的近似值，则有 yn 1 yn hf ( xn , yn ) （n 0,1,2,...） 算出y1 , y2 , y3 ...。 (1)
3
Runge-Kutta方法
由二元函数Taylor展开式 : k 2 hf ( xn , yn ) h 2 f x( xn , yn ) hk1 f y ( xn , yn ) O(h 3 ) hy( xn ) h 2 (f x( xn , yn ) f ( xn , yn ) f y ( xn , yn )) O(h 3 )
x
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.000 1.000000 1.191818 1.277438 1.358213 1.435133 1.508966 1.580338 1.649783 1.717779 1.784770
改进的Euler法y 精确解
1.000000 1.095909 1.184097 1.266201 1.343360 1.416402 1.485956 1.552514 1.616475 1.678166 1.737867 1.000000 1.095445 1.183216 1.264911 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051
第6章 常微分方程数值解法
绪论
在工程和科学计算中，所建立的各 种常微分方程的初值或边值问题，除很 少几类的特殊方程能给出解析解，绝大 多数的方程是很难甚至不可能给出解析 解的，其主要原因在于积分工具的局限 性。因此，人们转向用数值方法去解常 微分方程，并获得相当大的成功，讨论 和研究常微分方程的数值解法是有重要 意义的。