龙格库塔方法的由来和推导
进的欧拉格式是众多的二阶龙格—库塔法中的一种
特殊格式。
2018/10/18
15
若取 a2 阶龙格
-= b库21塔= 公12 ,式c1。=
0,
c2 = 1 ,就是另一种形式的二
yi+1 = yi + K2
K1 = hf (xi , yi )
K2
=
hf
( xi
+
1 2
h,
yi
+
1 2
K1)
此计算公式称为变形的二阶龙格—库塔法。式中
2
=
hf (xi
+
h 2
,
yi
+
1 2
K1)
K3 = hf (xi + h, yi - K1 + 2K 2 )
2018/10/18
20
9.4.4 四阶(经典)龙格—库塔法
如果需要再提高精度,用类似上述的处理 方法,只需在区间[xi,xi+1]上用四个点处的斜 率加权平均作为平均斜率K*的近似值,构成一 系列四阶龙格—库塔公式。具有四阶精度,即 局部截断误差是O(h5)。
yi+1 = yi + K
K可以认为是y = y(x)在区间[xi , xi+1]上的平均斜率
y
只要使用适当的方法求出y(x)在区 间[xi , xi+1]上平均斜率的近似值K
就可得到相应的Runge-Kutta方法
y = y(x)
K
2018/10/18
xi
xi +1
x
9
如果以y(x)在xi处的斜率作为y(x)在[xi , xi+1]上的平均斜率K
2018/10/18
6
一般龙格-库塔方法的形式为
yi+1 = yi + c1K1 + c2K2 + + cpK p
K1
=
hf
(
xi
,
yi
)
K
2
=
hf
( xi
+
a2h,
yi
+
b21K1)
• •• •• •• •• ••
K p = hf (xi + a ph, yi + bp1K1 + + bp, p-1K p-1)
2018/10/18
5
于是可考虑用函数f(x,y)在若干点上的函数值的 线性组合来构造近似公式,构造时要求近似公式在 (xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式 的前面几项重合,从而使近似公式达到所需要的阶数。 既避免求高阶导数,又提高了计算方法精度的阶数。 或者说,在[xi,xi+1]这一步内多计算几个点的斜率值, 然后将其进行加权平均作为平均斜率,则可构造出更 高精度的计算格式,这就是龙格—库塔(RungeKutta)法的基本思想 。
应项的系数相等,得到系数方程组:
c1 + c2 + c3 = 1
a2c2
+
a3c3
=
1 2
,
a22c2
+
a32c3
=
1 3
,
b221c2 + (b31 +
a 2b32 c2
=
1 6
,
b32 )2 c3
=
1 3
b21c2
+
(b31
+
b32 )c3
=
1 2
a 2 b21c 2
= y(xi ) + c1hf (xi , yi )
[ ] + c2h f (xi , yi ) + a2hfx¢ + b21hff y¢ + O(h3)
= y(xi ) + (c1 + c2 )hf (xi , yi ) + a2c2h2 fx¢ + b21c2h2 ff y¢ + O(h3)
2018/10/18
就是二阶ห้องสมุดไป่ตู้格
-库塔
公式,也就是改进的欧拉法。
y
K2
yi +1
=
yi
+
1 2
(K1
+
K2 )
K1 = hf ( xi , yi )
y = y(x)
K1
K
K2 = hf ( xi + h, yi + K1)
xi
xi+1
x
因此,凡满足条件式有一簇形如上式的计算格
式,这些格式统称为二阶龙格—库塔格式。因此改
y(xi+1) - y(xi ) = y¢(xi )(xi+1 - xi )
其中xi Î ( xi , xi+1)
即
y(xi+1) = y(xi ) + hy¢(xi )
2018/10/18
8
引入记号
y(xi+1) = y(xi ) + K
K = hy¢(xi ) = hf [xi , y(xi )]
2018/10/18
2
显然p=1时, y i+1=y i+hf(xi,y i)
它即为我们熟悉的Euler方法。
当p≥2时,要利用泰勒方法就需要计算f(x,y)的高 阶微商。这个计算量是很大的,尤其当f(x,y)较复 杂时,其高阶导数会很复杂。因此,利用泰勒公 式构造高阶公式是不实用的。但是泰勒级数展开 法的基本思想是许多数值方法的基础。
推导过程与前面类似,由于过程复杂,这 里从略,只介绍最常用的一种四阶经典龙格— 库塔公式。
2018/10/18
21
设 yi+1=yi+c1K1+c2K2+c3K3+c4K4
这里K1、K2、K3、K4为四个不同点上的函数值, 分别设其为
K1=hf (xi, yi) K2=hf (xi+a2h, yi+b21K1) K3=hf (xi+a3h, yi+b31K1+b32K2) K4=hf (xi+a4h, yi+b41K1+b42K2+b43K3) 其中c1、c2、c3、c4、a2、a3、a4、b21、b31、 b32、b41、b42、b43均为待定系数。
在[xi, xi+1]上取两点xi和xi+a2= xi +a2h,以该两点处 的斜率值K1和K2的加权平均(或称为线性组合)来求取 平均斜率k*的近似值K,即
K = c1K1 + c2K2
式中:K1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy'(xi) K2为xi +a2h点处的切线斜率值,比照改进的欧
K3
K2 = hf (xi + a2h, yi + b21K1)
K1
K
K3 = hf (xi + a3h, yi + b31K1 + b32K2)
2018/10/18
xi
xi +a 2
xi+a3
x
18
同理推导二阶公式,将y(xi+1)和yi+1在x=xi处进 行Taylor展开,使局部截断误差达到O(h4),使对
yi +1
=
yi
+ hyi¢
+
h2 2!
yi¢¢ +
+
hP P!
yi( P)
P阶泰勒方法
其中 yi¢ = f , yi¢¢ = f (xi , yi )¢x = f x¢ + ff y¢
yi¢¢¢= ( fx¢ + ff y¢ )¢x = f x¢¢x + 2 f x¢¢y f + f y¢¢y f 2 + f x¢ f y¢ + ( f y¢ )2 f
2018/10/18
4
同理,改进Euler公式可改写成
yi+1 = yi K1 = hf (
+ xi
1 2
K1
, yi )
+
1 2
K2
K
2
=
hf
( xi
+
h, yi
+
K1)
局部截断误差为O(h3)
上述两组公式在形式上共同点:都是用f(x,y)在某 些点上值的线性组合得出y(xi+1)的近似值yi+1, 且增 加计算的次数f(x,y)的次数,可提高截断误差的阶。如 欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进 欧拉法需计算两次f(x,y)的值,为二阶方法。
2018/10/18
12
[ ] K2 = h f (xi, yi) + a2hfx¢ + b21K1 f y¢ + O(h2)
[ ] = h f (xi, yi ) + a2hfx¢ + b21hff y¢ + O(h3)
K1 =hf(xi, yi)
\
yi+1 = yi + (c1K1 + c2K2 )
13
令 y(xi+1) = yi+1 对应项的系数相等,得到
c1 + c2 = 1 ,
a2c2
=
1 2
,
b21c2
=
1 2
这里有4 个未知 数,3 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶 龙格 -库塔格式。
2018/10/18