龙格库塔法推导

类似地，若取前P+1项作为y(xi+1)的近似值，便得到
[
]
yi +1
h2 h P ( P) = yi + hyi′ + yi′′ + L + yi P! 2!
P阶泰勒方法
其中 yi′ = f , yi′′ = f ( xi , yi )′ x = f x′ + ff y′
′ )′x = f xx ′′ + 2 f xy ′′ f + f yy ′′ f 2 + f x′ f y ′ +( fy ′ )2 f yi′′′ = ( f x′ + ff y
的解y=y(x)，在区间[xi, xi+1]上使用微分中值定理，有
y ( xi +1 ) − y ( xi ) = y ′(ξ i )( xi +1 − xi )
其中ξ i ∈ ( xi , xi +1 )
即
2014/9/2
y ( xi +1 ) = y ( xi ) + hy ′(ξ i )
8
引入记号
′ + b21c2 h 2 ff y ′ + O ( h3 ) + a2c2 h 2 f x
2014/9/2 13
令 y ( xi +1 ) = yi +1 对应项的系数相等，得到
1 1 c1 + c2 = 1 , a2c2 = , b21c2 = 2 2
这里有 4 个未知 数，3 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶 龙格 - 库塔格式。
h2 y ( xi +1 ) = y ( xi ) + hy ′( xi ) + y ′′( xi ) + O(h 3 ) 2 = yi + hf ( xi , yi )
h2 ′ ( xi , y i ) + O ( h 3 ) f x′ ( xi , yi ) + f ′( xi , yi ) f y + 2
2014/9/2 14
1 注意到，a2 = b21 = 1, c1 = c2 = 就是二阶龙格 - 库塔 2 y 公式，也就是改进的欧拉法。 K2
yi +1 K1 K2
= = =
1 yi + (K1 + K 2 ) 2 hf ( xi , yi ) hf ( xi + h, yi + K1 )
y = y( x )
y ( xi +1 ) = y ( xi ) + K
K = hy′(ξi ) = hf [ξi , y (ξi )]
yi +1 = yi + K
K可以认为是y = y ( x)在区间[ xi , xi +1 ]上的平均斜率
y
只要使用适当的方法求出y ( x)在区 间[ xi , xi +1 ]上平均斜率的近似值K
K 2 = hf ( xi + a2 h, yi + b21K1 )
确定系数 c1、c2、a2、b21 ，可得到有2阶精度的算法格式
2014/9/2 11
因此
y ( xi +1 ) = y ( xi ) + K
= y ( xi ) + (c1K1 + c2 K 2 )
将y(xi+1)在x=xi处进行Taylor展开： h2 y ( xi +1 ) = y ( xi ) + hy′( xi ) + y′′( xi ) + O(h3 ) 2! h2 ′ + ff y ′ + O ( h3 ) fx = y ( xi ) + hf ( xi , yi ) + 2!
2014/9/2 4
同理，改进Euler公式可改写成
1 1 yi +1 = yi + 2 K1 + 2 K 2 K1 = hf ( xi , yi ) K = hf ( x + h, y + K ) 1 i i 2
局部截断误差为O(h3)
上述两组公式在形式上共同点:都是用f(x,y)在某 些点上值的线性组合得出y(xi+1)的近似值yi+1, 且增 加计算的次数f(x,y)的次数,可提高截断误差的阶。如 欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进 欧拉法需计算两次f(x,y)的值，为二阶方法。
2014/9/2 16
二级R-K方法是显式单步式,每前进一步需要计算两个 函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2, b21 ，可使每步计算两次函数值的二阶 R-K 方法达到二阶 精度。能否在计算函数值次数不变的情况下,通过选择不 同的参数值,使得二阶R-K方法的精度再提高呢? 答案是否定的！无论四个参数怎样选择 ,都不能使公 式的局部截断误差提高到三阶。 这说明每一步计算两个函数值的二阶R-K方法最高阶 为二阶。若要获得更高阶得数值方法,就必须增加计算函 数值的次数。
y = y( x )
K
就可得到相应的Runge-Kutta方法
2014/9/2
xi xi +1
x
9
如果以y ( x)在xi 处的斜率作为y ( x)在[ xi , xi +1 ]上的平均斜率K
即
K = hy′( xi ) = hf [ xi , y ( xi )]
= hf ( xi , yi )
如下图
2014/9/2 3
Runge-Kutta 方法是一种高精度的单步法,简称R-K法 9.4.1 龙格-库塔(R-K)法的基本思想 Euler公式可改写成
yi +1 = yi + K K = hf ( xi , yi )
则 yi+1 的表达式与 y(xi+1) 的 Taylor 展开式的前两项 完全相同,即局部截断误差为O(h2)。
K1
K
xi
xi +1
x
因此，凡满足条件式有一簇形如上式的计算格 式，这些格式统称为二阶龙格—库塔格式。因此改 进的欧拉格式是众多的二阶龙格—库塔法中的一种 特殊格式。
2014/9/2 15
1 若取 a2 = b21 = , c1 = 0, c2 = 1 ，就是另一种形式的二 2
阶龙格 - 库塔公式。
2014/9/2
17
9.4.3 三阶龙格—库塔法 为进一步提高精度，在区间[xi, xi+1]上除两点xi和 xi+a2= xi +a2h,以外，再增加一点xi+a3= xi +a3h ，用这 三点处的斜率值K1、K2和K3的加权平均得出平均斜率 K*的近似值K，这时计算格式具有形式:
y
yi +1 = yi + c1K1 + c2 K 2 + c3 K 3 K1 = hf ( xi , yi ) K 2 = hf ( xi + a2 h, yi + b21K1 ) K 3 = hf ( xi + a3h, yi + b31K1 + b32 K 2 )
称为P阶龙格－库塔方法。 其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作 Tailor展开，通过相同项的系数确定参数。
2014/9/2 7
Runge-Kutta方法的推导思想 对于常微分方程的初值问题
y′ = f ( x , y ) y(a ) = y0 a≤ x≤b
y
则上式化为
yi +1 = yi + hf ( xi , yi )
y = y( x )
K
K
即Euler方法
xi xi +1
x
Euler方法也称为一阶Runge-Kutta方法
2014/9/2 10
9.4.2
二阶龙格—库塔法
在[xi, xi+1]上取两点xi和xi+a2= xi +a2h,以该两点处 的斜率值K1和K2的加权平均(或称为线性组合)来求取 平均斜率k*的近似值K，即 K = c1K1 + c2 K 2 式中:K1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy'(xi) K2为xi +a2h点处的切线斜率值,比照改进的欧 拉法,将xi+a2视为xi+1，即可得
2014/9/2
6
一般龙格－库塔方法的形式为
yi +1 = yi + c1K1 + c2 K 2 + L + c p K p K1 = hf ( xi , yi ) K 2 = hf ( xi + a2 h, yi + b21K1 ) L L L L L L K p = hf ( xi + a p h, yi + b p1K1 + L + b p, p −1K p −1 )
例如：取前两项可得到
y ( xi +1 ) = y ( xi ) + hy ′( xi ) + O(h 2 )
= y ( xi ) + hf ( xi , y ( xi )) + O(h 2 ) = yi + hf ( xi , yi ) + O(h )
2014/9/2 1
2
若取前三项，可得到截断误差为O(h3)的公式
[
]
将 K 2 = hy′( xi + a2h) = hf ( xi + a2h, yi + b21K1 ) 在x=xi处进行Taylor展开：