变步长的龙格库塔法

变步长的龙格—库塔方法
以经典四阶龙格—库塔公式为例。从节点xn出发，以h为
(h) 步长求一近似值yn 1 , (h) 5 y ( xn 1 ) yn 1 ch ,
将步长折半，即取 值y
(h ／ 2) n 1
h 为步长从xn 跨两步到xn 1，求一近似 2
5
h , 每跨一步的截断误差是c ,因此有 2
5 (h ／ 2) n 1
h y ( xn 1 ) y 2c , 2 (h ／ 2) y ( xn 1 ) yn 1 1 由上两式 . (h) y ( xn 1 ) yn 1 16
(h ／ 2) y ( xn 1 ) yn 1
1 (h ／ 2) ( h) [ yn yn 1 1 ]. 15
公 式
Euler公式 改进Euler公式 yn 1 yn hK1 K1 f ( xn , yn ) 1 1 y y h ( K K2 ) n 1 n 1 2 2 K f (x , y ) 1 n n K 2 f ( xn h, yn hK1 )
一般地，RK方法设近似公式为
p yn 1 yn h ci K i i 1 K1 f ( xn , yn ) i 1 K i f ( xn ai h, yn h bij K j ) (i 2, 3 , p) j 1 确定原则是使近似公式在( xn , yn )处的Taylor展开式与
y ( x )在xn处的Taylor展开式的前面项尽可能多地重合。
7.4.6 变步长的龙格-库塔法
在微分方程的数值解中，选择适当的步长是非常 重要的。单从每一步看，步长越小，截断误差就越 小；但随着步长的缩小，在一定的求解区间内所要 完成的步数就增加了。这样会引起计算量的增大， 并且会引起舍入误差的大量积累与传播。因此微分 方程数值解法也有选择步长的问题。 以经典的四阶龙格-库塔法(7.20)为例。从节点 (h) xi出发，先以h为步长求出一个近似值，记为 yi 1 ， 由于局部截断误差为 O(h 5 ) ，故有
y( xi 1 ) y
( h) i 1
ch
5
当h值不大时，式中的系数c可近似地看作为常数。
然后将步长折半,即以为 步长,从节点xi出发,跨 h ( ) 两步到节点xi+1,再求得一个近似值 yi 21 ,每跨一步的 h ,因此有 截断误差是 c
5
h 2
2
y ( xi 1 ) y ( x
h ) 2 i 1 (
这样
由此可得
2 y ( x i 1 ) y i 1 1 (h) 16 y ( x i 1 ) y i 1
h ) 2c 2 h ( )
5
y( xi 1 ) y
h ( ) 2 i 1
h ( ) 2 i 1
1 (y 15
h ( ) 2 i 1
h) yi( 1)
这表明以 y 作为 y( xi 1 ) 的近似值，其误差可用步 h 长折半前后两次计算结果的偏差 ( ) 来判断所选步长是否适当
(h) 2 yi y 1 i 1
当要求的数值精度为ε时：
（1）如果Δ>ε，反复将步长折半进行计算，直
y i 1 至Δ<ε为止,并取其最后一次步长的计算结果作为
（ 2 ）如果Δ<ε，反复将步长加倍，直到Δ>ε
为止，并以上一次步长的计算结果作为 y i 1 。
这种通过步长加倍或折半来处理步长的方法称为
变步长法。表面上看，为了选择步长，每一步都要 反复判断Δ，增加了计算工作量，但在方程的解y(x)
变化剧烈的情况下，总的计算工作量得到减少，结
果还是合算的。
高阶Runge—Kutta Method
y
i 1
§2 Runge-Kutta Method
y h [ K K ... K ]
i 1 1 2 2 m m i i
K f (x , y )
1
K f ( x h, y
2 i 2 i
21
hK )
1
K f ( x h, y
3 i 3 i
31
hK
1
32
hK )
2
... ... K
m
f ( x h, y
i m
m1
hK
1
m2
hK ...
2
m m 1
hK
m 1
)
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确 定这些系数的步骤与前面相似。
§2 Runge-Kutta Method
最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： h
y i 1 yi 6 ( K 1 2 K 2 2 K 3 K 4 ) f ( xi , yi ) h f ( xi h , y K1 ) i 2 2 f ( xi h , yi h K2 ) 2 2 f ( x i h, yi hK 3 ) K1 K2 K3 K4
Gill公式：4阶经典龙格-库塔公式的一种改进
yi 1 yi h 6 K1 2 2 K 2 2 2 K 3 K 4 K1 K2 K3 K4 f ( xi , yi ) h f ( xi h , y i 2 2 K1 ) f ( xi h 2 , yi f ( xi h , yi
2 1 2 2 2





hK1 1
hK 2 1


2 2 2 2
hK )
3
hK )
2
§2 Runge-Kutta Method
注：
龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的
值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系：
每步须算Ki 的个数
2
3
4
5
6
7
O(h6 )
n8
O(hn2 )
4 可达到的最高精度 O(h2 ) O(h3 ) O(h4 ) O(h ) O(h5 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
变步长的Runge—Kutta Method
§2 Runge-Kutta Method
Rn1 ch p1 y( p1) xn
Q: 由局部截断误差可以看出，步长 h 越小，局部截断 误差越小；但步长减小，在一定求解范围（区间）内 要完成的步数就增加了，步数增加会引起计算量增大， 导致舍入误差积累。因此要选取适当的步长。 选择步长时要考虑两个问题： 1.如何衡量和检验计算结果的精度？ 2.如何根据所获得的精度处理步长？
HW: p.201 #6-8