第二节 龙格-库塔方法

h( k1 + k2 ) yn +1 = yn + 2 hk1 1 1 3 3 k1 = f [ xn + ( + )h, yn + )hk2 ] +( + 2 6 4 4 6 hk2 1 3 1 3 k2 = f [ xn + ( − )h, yn + ( − )hk1 + ] 2 6 4 6 4
1 库塔三阶方法 库塔三阶方法
h h k1 = f ( xn , yn ); k2 = f ( xn + , yn + k1 ) 2 2
h yn +1 = yn + ( k1 + 4k2 + k3 ) 6
k3 = f ( xn + h, yn − hk1 + 2hk2 )
四级方法： 四级方法：N = 4 常见的2种四阶方法： 常见的 种四阶方法： 种四阶方法 经典龙格 库塔方法 龙格1经典龙格-库塔方法
1 中点法（修正的Euler法） 取 c1 = 0, c2 = 1, a2 = 1 中点法（修正的 法 2 h h
y n +1 = y n + hf ( x n +
常见的3种二级方法： 常见的 种二级方法： 种二级方法
二阶龙格库塔方法 2 二阶龙格库塔方法
h yn +1 = yn + [ f ( xn , yn ) + f ( xn + h, yn + hf ( xn , yn ))] 2
1 取 c1 = c2 = , a2 = 1 2
2
, yn +
2
f ( x n , y n ))
类似于N 的推导方法， 的推导方法 三级方法： 三级方法：N = 3 类似于 = 2的推导方法，可得到 1 c1 + c 2 + c 3 = 1; c 2 a 2 + c 3 a 3 = ; 2 O ( h4 ) 1 1 2 2 c 2 a 2 + c 3 a 3 = ; c 3 a 2 b32 = 6 3 常见的2种三阶方法： 常见的 种三阶方法： 种三阶方法
yi + p = yi + phK1
并用它来估计斜率 K 2 ：
K 2 = f ( x i + p , yi + p )
于是得到如下形式的算法： 于是得到如下形式的算法：
yi +1 = yi + h(λ1 K1 + λ2 K 2 ), K1 = f ( xi , yi ), K = f ( x , y + phK ). i +1 i 1 2
K * = λ1 K1 + λ2 K 2
即取
yi +1 = yi + h(λ1 K1 + λ2 K 2 )
是待定常数。 其中 λ1 和 λ2 是待定常数。若取 K1 = f ( xi , yi ) ，则 问题在于如何确定 xi + p 处的斜率 K 2 和常数 λ1 和 λ2 。
仿照改进的欧拉方法， 的值， 仿照改进的欧拉方法，用欧拉方法预测 y( xi + p ) 的值，
en +1 = en + h[φ ( xn , y( xn ), h) − φ ( xn , yn , h)] + Tn +1
p阶的：∃h0 , 0 < h ≤ h0 满足| Tn+1 |≤ Ch p+1 因为单步法 单步法是 因为单步法是 阶的： p +1 | en +1 |≤| en | + hL | en | + Ch ≤α | en | +β
( h / 2) n +1 ( h) n +1
1 ≈ 16
记 ∆ =| y
(h / 2) n +1
−y
(h) n +1
|
16( y( xn +1 ) − y
(h / 2) n +1
) ≈ y( xn+1 ) − y
( h) n +1
y ( x n +1 ) − y
如果 ∆ 此时取 y
( h / 2) n +1
ki = f ( xn + ai h, yn + h∑ bij k j ); i = 1, ..., N
j =1
N
(1)一 二阶的隐式中点方法 的隐式中点方法： (1)一级二阶的隐式中点方法：
yn +1 = yn + hk1
hk1 h k1 = f ( xn + , yn + ) 2 2
(2)二 四阶的隐式 的隐式R-K方法： 方法： (2)二级四阶的隐式 方法
下面以2阶龙格-库塔方法为例来阐述这种思想 下面以 阶龙格-库塔方法为例来阐述这种思想 阶龙格 考察区间 [ xi , xi +1 ] 上的一点 xi + p = xi + ph, 0 < p ≤ 1 ， 用 xi 和 xi + p 的斜率 K1 和 K 2 的加权平均作为平均
K * 的近似值： 斜率 的近似值：
p(N)
≤ N −2
若要得到N阶以上方法，则使用N级隐式 级隐式R-K方法 若要得到N阶以上方法，则使用N级隐式R-K方法 阶以上方法 N级隐式 级隐式R-K方法的一般形式： 方法的一般形式： 级隐式 方法的一般形式
y n +1 = y n + h ∑ c i k i
i =1
N
N级隐式 级隐式R-K法 级隐式 法 可以达到2N阶 可以达到 阶
2x dy 龙格例2：用经典的龙格-库塔方法 ： 经典的龙格 库塔方法 = y− y x∈(0,1) ∈ 求解下列初值问题 h = 0.1 。 dx y (0) = 1 解：经典的四阶龙格-库塔公式： 经典的四阶龙格-库塔公式： 龙格 h yn+1 = yn + ( k1 + 2k2 + 2k3 + k4 ) 6 h 2( xn + ) 2x 2xn h 2 k1 = yn − ; k2 = yn + k1 − h yn 2 yn + k1 2 h 2( xn + ) h 2 ; k = y + hk − 2( xn + h) k3 = yn + k2 − 4 n 3 h 2 yn + hk3 yn + k2 2
可将改进的欧拉格式改写成
1 y i + 1 = 2 ( K 1 + K 2 ), K 1 = y i + h f ( x i , y i ), K = y + h f ( x , K ). i i +1 1 2
该公式可以看作是用 xi 和 xi +1 两个点处的斜率 K1 和 K 2 的算术平均值作为平均斜率。 的算术平均值作为平均斜率。 由改进型欧拉公式我们可以猜想， 由改进型欧拉公式我们可以猜想，如果在[ xi , xi +1 ] 内多预测几个点的斜率，再对他们进行加权平均， 内多预测几个点的斜率，再对他们进行加权平均， 可能得到精度更好的平均斜率！ 可能得到精度更好的平均斜率！
其中 α
en = O ( h ) = 1 + hL, β = Ch 2 | en |≤ α | en −1 | + β ≤ α | en − 2 | +αβ + β
2xபைடு நூலகம்+1
xn yn xn yn
0.1 0.2 0.3 0.4 0.5 1.0954 1.1832 1.2649 1.3416 1.4142 0.6 0.7 0.8 0.9 1.0 1.4832 1.5492 1.6125 1.6733 1.7321
高阶和隐式 和隐式Runge-Kutta方法 二、高阶和隐式 方法 对于显式N级 方法， 级方法； 注： 对于显式 级R-K方法，最多只能得到 级方法； 方法 最多只能得到N级方法 已经证明N级 方法的阶 具有下列关系： 已经证明 级R-K方法的阶 p( N )具有下列关系： 方法的 N 1,2,3,4 N 5,6,7 N-1 8,9 N-2 10,11,…
局部截断误差 O ( h )
5
h yn+1 = yn + ( k1 + 2k2 + 2k3 + k4 ) 6 k1 = f ( xn , yn )
h h k2 = f ( xn + , yn + k1 ) 2 2 h h k3 = f ( x n + , y n + k 2 ) 2 2
k4 = f ( xn + h, yn + hk3 )
y ( x n +1 ) − y
( h) n +1
≈ Ch
5
h Step2:取 为步长，再从 xn 出发，经过两步计算得到 出发，经过两步 两步计算得到 取 为步长， 2
y( xn +1 ) − y
(h ) 2 n +1
h 5 ≈ 2C ( ) 2
y( xn+1 ) − y y( xn+1 ) − y
则将步长折半进行计算， 折半进行计算 > ε ，则将步长折半进行计算，直到 ∆ < ε 为止
(h / 2) 为最终结果； n +1 为最终结果；
1 (h / 2) ( h) ≈ [ yn +1 − yn +1 ] 15
则将步长加倍进行计算， 加倍进行计算 如果 ∆ < ε ，则将步长加倍进行计算，直到 ∆ > ε 为止 此时将步长折半一次计算，得到的为最终结果。 折半一次计算 此时将步长折半一次计算，得到的为最终结果。