龙格库塔法推导

8
引入记号
y ( xi +1 ) = y ( xi ) + K
K = hy′(ξi ) = hf [ξi , y (ξi )]
y i +1 = y i + K
K可以认为是y = y ( x)在区间[ xi , xi +1 ]上的平均斜率
y
只要使用适当的方法求 出y ( x)在区 间[ xi , xi +1 ]上平均斜率的近似值 K
y i +1 K1 K2
= = =
1 y i + (K 1 + K 2 ) 2 hf ( x i , y i ) hf ( x i + h , y i + K 1 )
y = y(x)
K1
K
xi
xi+1
x
因此，凡满足条件式有一簇形如上式的计算格 式，这些格式统称为二阶龙格—库塔格式。因此改 进的欧拉格式是众多的二阶龙格—库塔法中的一种 特殊格式。
[
]
将 K 2 = hy′( xi + a2h) = hf ( xi + a2h, yi + b21K1) 在x=xi处进行Taylor展开：
2011-12-8 12
′ ′ K 2 = h f ( xi , y i ) + a2 hf x + b21K1 f y + O( h 2 )
[
]
′ ′ = h f ( xi , yi ) + a2 hf x + b21hff y + O ( h 3 )
[
]
K1 =hf(xi, yi)
∴
yi +1 = yi + (c1K1 + c2 K 2 )
= y ( xi ) + c1hf ( xi , yi )
′ ′ + c2 h f ( xi , yi ) + a2 hf x + b21hff y + O ( h 3 )
[
]
= y ( xi ) + (c1 + c2 ) hf ( xi , yi )
′ ′ + a2c2 h 2 f x + b21c2 h 2 ff y + O ( h 3 )
2011-12-8 13
令 y ( x i + 1 ) = y i + 1 对应项的系数相等，得到
1 1 c1 + c 2 = 1 , a 2 c 2 = , b 21 c 2 = 2 2
这里有 4 个未知 个方程。 数，3 个方程。
2011-12-8 4
同理，改进Euler公式可改写成
1 1 yi +1 = yi + 2 K 1 + 2 K 2 K 1 = hf ( x i , y i ) K = hf ( x + h , y + K ) i i 1 2
局部截断误差为O(h3) 局部截断误差
上述两组公式在形式上共同点:都是用f(x,y)在某 些点上值的线性组合得出y(xi+1)的近似值yi+1, 且增 加计算的次数f(x,y)的次数,可提高截断误差的阶。如 欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进 欧拉法需计算两次f(x,y)的值，为二阶方法。
称为P阶龙格－库塔方法。 其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作 Tailor展开，通过相同项的系数确定参数。
2011-12-8 7
Runge-Kutta方法的推导思想 对于常微分方程的初值问题
y′ = f ( x , y ) y ( a ) = y0 a≤ x≤b
例如： 例如：取前两项可得到
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 2 )
2011-12-8 1
若取前三项，可得到截断误差为 若取前三项，可得到截断误差为O(h3)的公式 的公式
2011-12-8 15
若取 a 2 = b21 库塔公式。 阶龙格 - 库塔公式。
1 = , c1 = 0, c2 = 1 ，就是另一种形式的二 2
y i +1 K1 K2
= = =
yi + K 2 hf ( x i , y i ) 1 1 hf ( x i + h , y i + K 1 ) 2 2
′ ′′ ′′ ′′ ′ ′ yi′′′ = ( f x′ + ff y )′x = f xx + 2 f xy f + f yy f 2 + f x′ f y + ( f y ) 2 f
2011-12-8 2
显然p=1时, y i+1=y i+hf(xi,y i) 它即为我们熟悉的Euler方法。 当p≥2时,要利用泰勒方法就需要计算f(x,y)的高 阶微商。这个计算量是很大的，尤其当f(x,y)较复 杂时，其高阶导数会很复杂。因此，利用泰勒公 式构造高阶公式是不实用的。但是泰勒级数展开 法的基本思想是许多数值方法的基础。 R-K方法不是直接使用Taylor级数,而是利用它的思想
的解y=y(x)，在区间[xi, xi+1]上使用微分中值定理，有
y ( x i +1 ) − y ( x i ) = y ′(ξ i )( x i +1 − x i )
其中 ξ i ∈ ( x i , x i + 1 )
即
2011-12-8
y ( xi +1 ) = y ( xi ) + hy ′(ξ i )
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 )
y = y(x)
K2
K3
K1
K
xi
2011-12-8
xi +a2
xi +a3
18
x
同理推导二阶公式，将y(xi+1)和yi+1在x=xi处进 行Taylor展开，使局部截断误差达到O(h4)，使对 应项的系数相等，得到系数方程组：
c1 + c 2 + c3 = 1 a c + a c = 1 , 3 3 2 2 2 1 2 2 a 2 c 2 + a 3 c3 = , 3 1 2 2 b21c 2 + (b31 + b32 ) c3 = 3 a b c = 1 , 2 32 2 6
y = y(x)
K
就可得到相应的Runge-Kutta方法
xi xi +1
2011-12-8
x
9
如果以y ( x)在xi 处的斜率作为y ( x)在[ xi , xi +1 ]上的平均斜率K
即
K = hy′( xi ) = hf [ xi , y ( xi )]
= hf ( xi , yi )
如下图
2011-12-8
6
一般龙格－库塔方法的形式为
yi +1 = yi + c1K1 + c2K2 + L + c p K p K1 = hf ( xi , yi ) K2 = hf ( xi + a2h, yi + b21K1) L L L L L L K p = hf ( xi + a ph, yi + bp1K1 + L + bp, p−1K p−1)
§ 9.4 龙格－库塔方法 龙格－
得到高精度方法的一个直接想法是利用Taylor展开 y' =f(x,y) (a≤x≤b) 假设式 中的 f(x,y) 充分光滑,将y(xi+1)在x i点作Taylor展开，若 取右端不同的有限项作为y(xi+1)的近似值，就可得到 计算y(xi+1)的各种不同截断误差的数值公式。
2011-12-8 3
Runge-Kutta 方法是一种高精度的单步法,简称 方法是一种高精度的单步法,简称R-K法 法 龙格-库塔( 9.4.1 龙格-库塔(R-K)法的基本思想 ) Euler公式可改写成
yi+1 = yi + K K = hf ( x i , y i )
则yi+1 的表达式与y(xi+1)的Taylor展开式的前两项 完全相同,即局部截断误差 局部截断误差为O(h2)。 局部截断误差
2011-12-8
17
9.4.3 三阶龙格 库塔法 三阶龙格—库塔法 为进一步提高精度，在区间[xi, xi+1]上除两点xi和 xi+a2= xi +a2h,以外，再增加一点xi+a3= xi +a3h ，用这 三点处的斜率值K1、K2和K3的加权平均得出平均斜率 K*的近似值K，这时计算格式具有形式:
类似地，若取前P+1项作为 项作为y(x )的近似值 的近似值， 类似地，若取前P+1项作为y(xi+1)的近似值，便得到
2
[
]
yi +1
h2 h P ( P) = yi + hyi′ + yi′′ + L + yi 2! P!
P阶泰勒方法
其中 yi′ = f , yi′′ = f ( xi , yi )′ x = f x′ + ff y′
此计算公式称为变形的二阶龙格—库塔法。式中
1 xi + h 2
为区间 [xi , xi +1 ] 的中点。也称中点公式。 也称中点公式。 也称中点公式
Q:为获得更高的精度，应该如何进一步推广？ 为获得更高的精度，应该如何进一步推广？ 为获得更高的精度
2011-12-8 16
二级R-K方法是显式单步式,每前进一步需要计算两个 函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2, b21 ，可使每步计算两次函数值的二阶R-K方法达到二阶 精度。能否在计算函数值次数不变的情况下,通过选择不 同的参数值,使得二阶R-K方法的精度再提高呢? 答案是否定的！无论四个参数怎样选择,都不能使公 式的局部截断误差提高到三阶。 这说明每一步计算两个函数值的二阶R-K方法最高阶 为二阶。若要获得更高阶得数值方法,就必须增加计算函 若要获得更高阶得数值方法, 若要获得更高阶得数值方法 数值的次数。 数值的次数。
存在无穷多个解。所有满足上式的格式统称为2阶 无穷多个解。 无穷多个解 阶 龙格 - 库塔格式。 库塔格式。
2011-12-8 14
Leabharlann Baidu
1 注意到， 就是二阶龙格 注意到，a 2 = b21 = 1, c1 = c2 = 就是二阶龙格 - 库塔 2 y 公式，也就是改进的欧拉法 改进的欧拉法。 公式，也就是改进的欧拉法。 K2
y
则上式化为
y = y(x)
yi +1 = yi + hf ( xi , yi )
即Euler方法
K
K
xi
xi +1
x
Euler方法也称为一阶 一阶Runge-Kutta方法 一阶 方法
2011-12-8 10
二阶龙格— 9.4.2 二阶龙格—库塔法 在[xi, xi+1]上取两点xi和xi+a2= xi +a2h,以该两点处 的斜率值K1和K2的加权平均(或称为线性组合)来求取 平均斜率k*的近似值K，即 K = c1 K 1 + c 2 K 2 式中:K1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy'(xi) K K2为xi +a2h点处的切线斜率值,比照改进的欧 拉法,将xi+a2视为xi+1，即可得
K 2 = hf ( xi + a2h, yi + b21K1 )
确定系数 c1、c2、a2、b21 ，可得到有2阶精度的算法格式
2011-12-8 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 ( h 3 ) 2! h2 ′ + ff y + O ( h 3 ) ′ = y ( xi ) + hf ( xi , yi ) + fx 2!
h2 y ( xi +1 ) = y ( xi ) + hy ′( xi ) + y ′′( xi ) + O (h 3 ) 2 = yi + hf ( xi , yi )
h ′ + f x′ ( xi , yi ) + f ′( xi , yi ) f y ( xi , yi ) + O (h 3 ) 2
2011-12-8 5
于是可考虑用函数f(x,y)在若干点上的函数值的 线性组合来构造近似公式，构造时要求近似公式在 (xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式 的前面几项重合，从而使近似公式达到所需要的阶数。 既避免求高阶导数,又提高了计算方法精度的阶数。 或者说,在[xi,xi+1]这一步内多计算几个点的斜率值， 然后将其进行加权平均作为平均斜率，则可构造出更 高精度的计算格式，这就是龙格—库塔（Runge龙格— 龙格 库塔（ Kutta）法的基本思想 ）法的基本思想。