第五章 常微分方程数值解

tn1 tn
f (x(t),t)dt ， n 0,1 ,2, .
(5-3)
然后数值积分。 或采用台劳展开的 ，方 即 hn 法 tn1 tn (tn) x(tn1) x(tn) hnxn (5- 4)
然后按精度要求截断。 定义：若 x 只 计用 算到 前 面 x 一 ,即 步 算 信 法 息 是 n 1 n x x h ( x ,t ,h ) n 1 n n n n n
本章只介绍初始问题的解法，原因是初始问题系统科学中 研究的最重要的问题。 §5.2 初始问题的单步法
一、欧拉法 1、显式欧拉法 x n 1 x n h n f n n 0 ,1, 2 , 其中， h n t n 1 t n , f n f ( x n , t n ) 局部截断误差： T x ( t ) ~ x (t
二、龙格-库塔法（RK)
我们的目的是构造更高阶的单步法！
p 1 结论：若单 截 步 断 法 T 误 的 o 差 ( h 局 ), 则 部 此单 n 1
p 步法的整 err 体 o ( 误 h ). 即 差 此单 p 阶 步方 法法 为 n 1
分 析 ： 由 x (t n 1 ) x ( t n )
j i 1 j
j 1

x ( t n ) h ( a i a i 1 ) k i 1 o ( h 2 ) x ( t n ) h b ji k i 1
i 1 j
b ji 待 定
k j f ( x ( t n ) h b ji k i 1 , t n a j h )
这样做的根据是：1）对实验模型数值不稳定的方法，不可用； 2）一般的初始问题在其解的存在区域内，可局部线性化转化 为试验方程。
几个常用算法的数值稳定性
1、显式欧拉法的绝对稳 定区域 显式 Euler 用于试验模型的计算公 式为 x n 1 x n h x n (1 h ) x n n 1 (1 h ) n 绝对稳定区域 : 1 h 1 当 为实数时，绝对稳定区 间： h （ 2， 0） 步长 h 不仅要满足算法收敛（ o ( h ))限制，而且 还要满足算法数值稳定 性（0 h 2 )限制.
i1
j 1: m
称此结构为m+1级龙格-库塔法。它最多是m+1阶的。
一些常用的显式RK方法 1、二级二阶RK法 改进的欧拉法（p277),中点公式（p278),休恩方法（p278)等。 2、三级三阶RK法 休恩三阶方法（P279),库塔三阶方法（P279)等
3、四级四阶RK法
经典四级四阶RK法（P280),基尔四级四阶RK法(P280)
当 为实数时，绝对稳定区 步长 h 不仅考虑满足算法收敛 还要考虑数值稳定性（
n 1 n 1
~ x ( t n 1 ) x ( t n ) h n f ( x ( t n ), t n ) T n 1 x ( ) 2 hn 2 !!
n 1 )
总体误差： err n 1 x ( t n 1 ) x n 1 x (t n 1 ) ~ x (t n 1 ) ~ x (t n 1 ) x n 1 ~ x (t ) x
4 、经典四级四阶
RK 法的绝对稳定区域 RK 法用于试验模型的计算 公式为
经典四级四阶
(h )2 (h )3 (h )4 x n 1 （ 1 h ) xn 2! 3! 4! (h )2 (h )3 (h )4 n 1 （ 1 h ) n 2! 3! 4! 绝对稳定区域 (h )2 (h )3 (h )4 :（ 1 h ) 1 2! 3! 4! 间： h （ 2 . 78 ， 0） （ o ( h 2 )) 限制，而且
一 般F (t , x, x,, x( n) ) 0 均 可 化为 dx f ( x, t ) dt 形 式 ， 称 此 为标 准 模。 式
例如 d 2 x ( 1 x 2 ) dx x f ( t ) dt dt 2 引入 x 1 x , x 2 dx ，则 dt dx 1 x2 dt dx 2 ( 1 x 12 ) x 2 x 1 f ( t ) dt
第五章 常微分方程数值解 §5.1 问题的提出
常微分方程定解问题 1 、初始问题 dx f ( x , t ), t t 0 dt (5 - 1) x (t0 ) x 0 x ( t ) R n , t t 0 是其解析解。 2 、边值问题 dx f ( x , t ) ， a t b (5 - 2) dt x (a ) x a , x (b ) xb x ( t ) R n , a t b 是其解析解。
的绝对稳定区域包含了
A - 稳定的算法，只要
3、改进的欧拉法的绝对 稳定区域 改进的 Euler 用于试验模型的计算公 式为 h 2 h x n 1 x n （ x n x n 1） x n 1 xn 2 2 h n 1 2 h n 2 h 2 h 绝对稳定区域 : 1, 是 A 稳定的。 2 h 当 为实数时，绝对稳定区 间： h （ ， 0） 步长 h 只需考虑满足算法收敛 （ o ( h 2 ))限制即可！
对于非线性问题很难求 得解析解 x ( t) 。数值解法 求在给定时刻 t ,t , ,t 的近似值 x ,x , ,x . 然后在利 0 1 n 0 1 n ~ 值逼近的方法得到 x ( t) 得近似 x ( t).
* 初始问题算法的构 微造 分是 方程 5（ 1 ） 积 分 方 程 ： x(tn1) x(tn)

2、 隐 式 欧 拉 法 的 绝 对 稳
定区域 式为 1 xn 1 h
隐 式 Euler 用 于 试 验 模 型 的 计 算 公 x n 1 x n h x n 1 x n 1
1 n 1 n 1 h 绝对稳定区域 : 1 h 1 当 为实数时，绝对稳定区 步长 h 只需考虑满足算法收敛 定义：一个算法如果它 面的整个左半平面，则 显然，凡是 是绝对稳定的。 称此算法是 间： h （ ， 0） （ o ( h )) 限 制 即 可 ！ h 复平 A 稳定的。 Re( ) 0, h 0 都
则称之为单步 则 算 为 法 多 。 步 否 算法。
* 边值问题的算法构造方 有限元法等。差分法的
法很多，有差分方法、 例子：
x (1 x 2 ) x x f ( t ), 0 t 1 x ( 0 ) 0 , x (1 ) 1 x ( t i ) 1 ( x i 1 x i ), x ( t i ) 12 ( x i 1 2 x i x i 1 ) h h 12 ( x i 1 2 x i x i 1 ) (1 x i2 ) 1 ( x i 1 x i ) x i f ( t i ) h h i 1:n 1 x 0 0 , 1 ( x n x n 1 ) 1 h 由此解出 x 0 , x1 , , x n .
设 h max( hk , k 1 : n ), M 2 max( x ( t ) , t 0 t t n 1 ), M2 2 则 T n 1 h ,及 2 M2 2 err n 1 h (1 hL ) err n 2 M2 2 n h ( (1 hL ) k ) (1 hL ) n 1 err 0 2 k 0 M 2 2 (1 hL ) n 1 1 M 2 ( n 1) hL err n 1 h (e 1) h 2 hL 2L err n 1 o ( h ) !!
n 1 n 1
x ( t n ) x n h n ( f ( x ( t n ), t n ) f ( x n , t n )) ~ x (t n 1 ) x n 1 x (t n ) x n h n L x (t n ) x n ~ x ( t n 1 ) x n 1 (1 h n L ) x ( t n ) x n err n 1 T n 1 (1 h n L ) err n 其 中 ， L 为 f的 李 布 希 兹 常 数
三、单步法的数值稳定性 若计算 x n 产生误差 n ( n )，即实际得到近似值
ˆn , x
ˆ n , 则用单步法 n xn x x n 1 x n h ( x n , t n , h ) 计算 x n 1时，将产生误差 ˆ n 1 n 1 x n 1 x ˆ n h [ ( x n , t n , h ) ( x ˆ n , t n , h )] xn x ˆ n , t n , h ) 太依赖于 f ( x , t ), 估计其大 由于 ( x n , t n , h ) ( x 小很困难。为摆脱这种 依赖性，在讨论方法的 数值 稳定性时，都针对同一 试验模型： x x , 为复常数
以上说明显式欧拉法是一阶收敛的算法
2、隐式欧拉法
xn 1 xn hn f ( xn 1 , t n 1 ) n 0,1,2, 一阶收敛的方法，即 errn 1 o ( h ) 3、改进的欧拉法（隐式） xn 1 xn hn [ f ( xn , t n ) f ( xn 1 , t n 1 )] n 0,1,2, 二阶收敛的方法，即 errn 1 o ( h 2 )。
i 1
j 1:m
(**)
由 式* （ ）和式 * * （ ）可知：适当选 cj ,取 aj ,b 参 可 数 构造 m 出 1 阶 ji 的显式单步法：