第十一章 常微分方程数值解1

2008-10-23
20
三、龙格-库塔法（RK) 如何构造更高阶的单步法？
由 x(t n + 1 ) = x(t n ) + 设 h = t n +1 − t n ,
∫
t n +1 tn
f(x(t), t)dt
ri = t n + a i h, 0 ≤ a i ≤ 1, i = 0 : m 。 ∵
∴ errn +1 = Ο(h) !!
2008-10-23 14
评价 1、显式欧拉法是一阶收敛的算法； 2、显式欧拉法计算简单； 3、显式欧拉法对函数的光滑性要求不高； 4、显式欧拉法也形象地称为折线法。
2008-10-23
15
2、隐式欧拉法
利用算法构造模式1 ，积分用右矩公式近似有：
x(t n +1 ) = x(t n ) + h n f (x(t n +1 ), t n +1 ) + Τn +1 令 x(t n )的近似值为x n , f (x(t n ), t n )的近似值为f n = f (x n , t n )。 则得到：
2008-10-23
6
∗ 初始问题（ 11-1）的数值算法的构造原理：
1、初始问题（ 11-1）等价为积分方程： x(tn +1 ) = x(tn ) + ∫ 1, 2, n = 0， .
tn+1 tn
f ( x(t ), t )dt (11-3)
然后利用数值积分方法。
2008-10-23
7
∗ 初始问题（ 11-1）的数值算法的构造原理：
2008-10-23 24
Tn +1 = x(t n +1 ) − x(t n ) − c0 hk 0 − c1hk1 =x(t n +1 ) − x(t n ) − (c0 +c1 )x′(t n )h − (f x ⋅ x ′(t n ) ⋅ c1b11 + f t ⋅ c1a1 )h 2 + Ο(h 3 ) = (1 − c0 − c1 )x′(t n )h 1 ⎡ 1 ⎤ 2 3 ′ +（ c b f x (t ) c a f h (h )。 − ） ⋅ + （ − ） ⋅ + Ο 1 11 x n 1 1 t⎥ ⎢ 2 ⎣ 2 ⎦
2、采用台劳展开的方法，即 ′ (t n ) + x(tn +1 ) = x(tn ) + hn xn hn = tn +1 − tn 然后按精度要求截断。
(11-4)
2008-10-23
8
∗ 常微分方程初始问题的数值解法有 单步法和多步法。
定义：若计算xn +1只用到前面一步信息xn , 即 xn +1 = xn + hnφ ( xn , tn , hn ) 则称之为单步算法。否则为多步算法。
2008-10-23
17
二、改进欧拉法 利用算法构造模式1 公式近似有： ，其积分用梯形求积
x(t n +1 ) = x(t n ) + h n [ f (x(t n ), t n ) + f (x(t n +1 ), t n +1 ) ] + Τn +1 令 x(t n )的近似值为x n , f (x(t n ), t n )的近似值为f n = f (x n , t n ). 则得到：
2008-10-23
9
边值问题的算法构造方法有差分方法、 有限元法、打靶法。 本章只介绍初始问题的求解方法。
常微分方程初始问题在系统科学研 究中 占有最重要的地位，寻求其有效数值解法是解 决问题的关键。
2008-10-23
10
§11.2 初始问题的单步法
一、欧拉法(Euler)
利用算法构造模式 ，按模式2有
j i =1 j j−1 ) h k j−1
+ Ο (h 2 )
= x ( t n ) + h ∑ ( a i − a i −1 ) k i −1 + Ο ( h 2 ) = x ( t n ) + h ∑ b ji k i − 1
i =1
b ji 待 定 .
m ⎧ ⎪ x ( t n +1 ) = x ( t n ) + h ∑ c jk j + T n +1 j= 0 ⎪ ⎪ k 0 = f ( x ( t n ), t n ), ∴ ⎨ j ⎪k = f (x (t ) + h b ji k i − 1 , t n + a j h ) , ∑ n ⎪ j i =1 ⎪ ⎩ j = 1: m.
2008-10-23
22
利用台劳公式，将函数均在 t n处展开， 适当选取参数 c j , a j , b ji可构造出 m + 1阶的显 式单步法：
x
n +1
= x
n
+ h ∑ c jk
j= 0 j
m
j
Biblioteka Baidu
k 0 = f (x n , tn ) k
j
= f (x
n
+ h ∑ b ji k i − 1 , t n + a j h )
2008-10-23 12
显式Euler法的整体误差：errn +1 = x(t n +1 )− x n +1
分析： errn +1 = x(t n +1 ) − x(t n +1 ) + x(t n +1 ) − x n +1 ∵ x(t n +1 ) − x 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 ∴ errn +1 ≤ Tn +1 + (1 + h n L) errn 其中， L为 f的李布希兹常数。
x ′′( ξ ) 2 x (t n +1 ) = x ( t n ) + h n x ′ hn n (t n ) + 2 x ′′( ξ ) 2 = x ( t n ) + h n f ( x ( t n ), t n ) + hn 2 令 x(t n )的 近 似 值 记 为 x n , f (x (t n ), t n )的 近 似 值 记 为 f n = f ( x n , t n )。 则得到：
( n −1) ′， ， 。 x(t 0 ) = x0， x′(t 0 ) = x0 x (n-1) (t 0 ) = x0
令
y 1 = x , y 2 = x ′,
, y n = x ( n − 1 )。
可化为标准模式： dy = f ( y , t )， t > t 0 dt y (t0 ) = y 0
从上述三个算法的分析发现：
Ο(h p +1 ) ,则整体误差Ο( h p ) ！！ •局部截断误差若是
•收敛阶的提高需要更多函数值计算！
2008-10-23
19
定理：若单步法的局部截断误差Tn +1 = Ο(h p +1 ), 则此单步法的整体误差 errn +1 = Ο(h p )。 证明：略
以下，目的是构造更高阶的单步法！
∵ k 0 = f (x(t n ), t n ) = x′(t n ) k1 = f (x(t n ) + hb11k 0 , t n + a1h) = f (x(t n ), t n ) + f x (x(t n ), t n )ihb11x ′(t n ) + f t (x(t n ), t n )ia1h + Ο(h 2 ) = x′(t n ) + hb11 if x (x(t n ), t n )x ′(t n ) + ha1 if t (x(t n ), t n ) + Ο(h 2 )
2008-10-23 11
1、 显 式 欧 拉 法 ⎧ x n +1 = x n + h n f n ⎨ ⎩ n = 0,1, 2, 其 中 ， h n = t n + 1 − t n。
记 x(t n +1 ) = x(t n ) + h n f (x(t n ), t n ) 局部截断误差： x ′′(ξ) 2 Tn +1 = x(t n +1 ) − x(t n +1 ) = hn 2 2 Tn +1 = Ο(h n )
第十一章 常微分方程数值解
¾ §11.1 常微分方程定解问题
¾ §11.2 初始问题的单步法 ¾ §11.3 初始问题的线性多步法 ¾ §11.4 微分方程组和 刚性问题的解法
2008-10-23
1
§11.1 常微分方程定解问题
一般地： F (t , x, x′, , x ( n ) ) = 0, t > t0
⎧ x n +1 = x n + h n [f (x n , t n ) + f (x n +1 , t n +1 )] ⎨ ⎩n = 0,1, 2,
2008-10-23 18
可以证明：
局部截断误差： Τ n +1 = Ο (h 3 n ), 整体误差： errn +1 = Ο (h 2 )。 改进 欧拉法是 隐式 的 二阶收敛 的方法。
∫
t n +1 tn
f(x(t), t)dt ≈ h ∑ c j k j
j= 0
m
k j = f(x(r j ), r j ), j = 0 : m . ∴ x(t n + 1 ) = x(t n ) + h ∑ c j k j + T n + 1
j= 0
2008-10-23 21
m
取 r0 = t n， 则 k 0 = f ( x ( t n ) , t n ) , ∵ x ( r j ) = x ( r j− 1 ) + ( a j − a
2008-10-23 2
1、 初 始 问 题 ⎧ dx = f(x,t),t > t 0 ⎪ (11 - 1) ⎨ dt ⎪ ⎩ x(t 0 ) = x 0 n x(t)∈ R ,t ≥ t 0 是 其 解 析 解 。
2008-10-23
3
2、边值问题 ⎧ dx ⎪ = f ( x , t )， a < t < b (11-2) ⎨ dt ⎪ ⎩ x ( a ) = xa , x (b ) = xb n x (t ) ∈ R , a ≤ t ≤ b是其解析解。 还有第二、第三类边值条件（略）
2008-10-23
4
d 2x 2 dx ′ 例如 : 2 + (1 − x ) + x = f (t ), x (0) = x0 , x′(0) = x0 dt dt dx 引入 y1 = x, y2 = ， y = ( y1 , y2 ) Τ , 则 dt dy = f ( y , t ), y (0) = y0。 dt ′ )Τ , 其中,y0 = ( x0 , x0 f = ( f1 , f 2 ) = ( y2 , − (1 − y12 ) y2 − y1 + f (t )) Τ .
T
2008-10-23
5
非 线 性 问 题 很 难 求 得 解 析 解 x ( t )， 寻 求其数值解是现实并可行的途径。 数 值 解 法 是 ： 求 在 给 定 时 刻 t 0 , t1 , 的 近 似 值 x 0 , x1 , , tn
, x n。 然 后 在 利 用 数 值 逼
近 的 方 法 得 到 x ( t )的 近 似 x ( t )。
2008-10-23 13
设 h = max(h k , k = 1: n), M 2 = max( x′′(t) , t 0 ≤ t ≤ t n +1 ), M2 2 则 Tn +1 ≤ h ,及 2 M2 2 errn +1 ≤ h + (1 + hL) errn 2 M2 2 n h (∑ (1 + hL) k ) + (1 + hL) n +1 err0 ≤ 2 k =0 M 2 2 (1 + hL) n +1 − 1 M 2 (n +1)hL h ⋅ (e ≤ ≤ − 1)h 2 hL 2L
⎧ x n +1 = x n + h n f (x n +1 , t n +1 ) ⎨ ⎩n = 0,1, 2,
2008-10-23 16
可以证明：
隐式欧拉法的局部截断误差： Τ n +1 = Ο (h 2 n ), 隐式欧拉法的整体误差： errn +1 = Ο (h n ).
评价 1、隐式欧拉法是一阶收敛的算法； 2、隐式欧拉法每步需解非线性方程； 3、隐式欧拉法对函数的光滑性要求不高。
i =1
j= 1:m
称此结构为m+1级龙格-库塔法。 可以证明：它最多是m+1阶的。
2008-10-23 23
例一：推导二级二阶的龙格-库塔方法。
二级龙格 − 库塔法的形式为： ⎧ x(t n +1 ) = x(t n ) + c 0 hk 0 + c1 hk 1 + Tn +1 ⎪ ⎨ k 0 = f (x (t n ), t n ), ⎪ k = f (x (t ) + hb k , t + a h ) n 11 0 n 1 ⎩ 1