第五章 常微分方程数值解
5.2 Runge-Kutta法§
5.2 Runge-Kutta法§
),( 111 ??? += kkkk yxhfyy
)],(),([2 111 kkkkkk yxfyxfhyy ++= ???
考虑改进 Euler法
如果将其改成 )(2 211 KK
hyy
kk ++= ?
),( 111 ??= kk yxfK
),( 112 hKyxfK kk += ?
)( 00 xyy =?
?
?
???
?
----------(1)
改进 Euler法是由梯形公式和 Euler公式复合而成
梯形公式具有 2阶精度
形如 (1)式的求解公式称为 二阶 Runge-Kutta法
同样可以证明 ,改进 Euler法也具有 2阶精度
对于 Simpson求解公式 ：
)],())2(,2(4),([6 11111 kkkkkkkk yxfhxyhxfyxfhyy +++++= ?????
这是隐式多步法 选取适当的显化方法 ,可得类似 (1)的高阶 Runge-Kutta方法
以下使用中值定理进行推导
一 、 Runge-Kutta方法的导出
??
?
=
≤≤=′
0)(
),(
yay
bxayxfy
对于常微分方程的边值问题
的解 ),(xyy = 有上使用微分中值定理在区间 ,],[ 1 nn xx ?
))(()()( 111 ??? ?′=? nnnnn xxyxyxy x ),( 11 nnn xx ?? ∈x
)()()( 111 ??? ′+=+ nnn yhxyhxy x即 ----------(3)
为了同学们课后复习的方便 ,以下的内容将 k写成 n.
上的平均斜率在区间可以认为是 ],[)( 1 nn xxxyyK ?=
hKxyhxy nn +=+ ?? )()( 11 ----------(3)引入记号
)( 1?′= nyK x )](,[ 11 ??= nn yf xx
Kxx
xy
nn 上平均斜率的近似值间
在区出只要使用适当的方法求
],[
)(
1?
就可得到相应的 Runge-Kutta方法
1?nx nx x
y
)(xyy =
hKyy nn += ?1 ----------(4)
即 (4)式
K
二 、 低阶 Runge-Kutta方法
1?nx nx x
y
)(xyy =
如下图
Kxxxyxxy nnn 上的平均斜率在处的斜率作为在如果以 ],[)()(.1 11 ??
即 )( 1?′= nxyK )](,[ 11 ??= nn xyxf
则 (4)式化为
),( 111 ??? += nnnn yxhfyy
),( 11 ??= nn yxf
即 Euler方法
Euler方法也称为 一阶 Runge-Kutta方法
)()( 2hOhen =由于
----(5)
K K
1?nx nx x
y
)(xyy =
上的平均斜率在
的算术平均值作为和处的斜率和在如果以
],[)(
)(.2
1
211
nn
nn
xxxy
KKxxxy
?
?
)( 11 ?′= nxyK ),( 11 ??= nn yxf
)(2 nxyK ′= )](,[ nn xyxf=
),( nn yxf≈
),( 11 hKyxf nn += ?(由 (5)式 )
令 2 21 KKK +=
则 (4)式化为
K1K
2K
)(2 211 KKhyy nn ++= ?
),( 111 ??= nn yxfK
),( 112 hKyxfK nn += ?
)( 00 xyy =?
?
?
???
?
-----------(6)
和 (1)式 一致 ,即改进 Euler公式 ,也称为 二 阶 Runge-Kutta法
)()( 3hOhen =
三 、 高阶 Runge-Kutta方法
2],[ 1212111
hxxxxx
nnnnn +== ??+?? 上增加一点如果
上的平均斜率在作为
的加权平均值和、处的斜率和、在且以
],[)(
)(
1
321
2
11
nn
nnn
xxxy
KKKxxxxy
?
??
)( 11 ?′= nxyK ),( 11 ??= nn yxf
)(
2
12 ?′= nxyK )](,[
2
1
2
1 ??= nn xyxf
1?nx nx21
hx
n +? x
y
)(xyy =
)(3 nxyK ′= )](,[ nn xyxf= 未知
K1K
2K 3K
)(
2
1?nxy
)2,2( 111 Khyhxf nn ++= ??令
11 2 K
hy
n +≈ ?
)(
2
12 ?′= nxyK
)2( 121 KKhyn ?+≈ ?)( nxy
))2(,( 1211 KKhyhxf nn ?++= ??)(3 nxyK ′=
)(
2
111 ?? nn xyKx 预测处的斜率如果以
)(2 2111 nnn xyKKhxx 预测、处的斜率、同样以 +??
令
参照 Simpson求解公式
)],())2(,2(4),([6 11111 nnnnnnnn yxfhxyhxfyxfhyy +++++= ?????
公式前进 Euler
取 321 616461 KKKK ++=
则 (4)式 化为
)4(6 3211 KKKhyy nn +++= ?
),( 111 ??= nn yxfK
)2,2( 1112 KhyhxfK nn ++= ??
))2(,( 12113 KKhyhxfK nn ?++= ??
-----------(7)
)( 00 xyy =?
?
?
???
?
(7)式称为 三阶 Runge-Kutta方法
展开式处的在、、分别作 TailorxKKK n 1321 ?
),( 111 ??= nn yxfK )( 1?′= nxy
)2,2( 1112 KhyhxfK nn ++= ??
),(2),(2),( 1111111 ?????? ′+′+= nnynnxnn yxfhKyxfhyxf
)()],(4),(4),(4[!21 311
2
1
2
11
1
2
11
2
hOyxfKhyxfKhyxfh nnynnxynnx +′′+′′+′′+ ??????
)(2)( 11 ?? ′′+′= nn xyhxy )()(8 31
2
hOxyh n +′′′+ ?
))2(,( 12113 KKhyhxfK nn ?++= ??
),()2(),(),( 11121111 ?????? ′?+′+= nnynnxnn yxfKKhyxfhyxf
)()],()2(
),()2(),([!21
3
11
2
12
2
1112
2
11
2
hOyxfKKh
yxfKKhyxfh
nny
nnxynnx
+′′?
+′′?+′′+
??
????
),()22()()( 111211 ???? ′?+′′+′= nnynn yxfKKhxyhxy
)()],()44(
),()22()([2
3
11212
11121
2
hOyxfKKK
yxfKKxyh
nny
nnxyn
+′′?
+′′?+′′′+
??
???
)()(2)()( 31
2
11 hOxy
hxyhxy
nnn +′′′+′′+′= ???
展开式处的在再作 Tailorxxyy n 1)( ?=
)()(!3)(2)()()( 41
3
1
2
111 hOxy
hxyhxyhxyhxy
nnnnn +′′′+′′+′+=+ ?????
)4(6 3211 KKKhyy nn +++= ?
因此
)()(!3)(2)( 41
3
1
2
11 hOxy
hxyhxyhy
nnnn +′′′+′′+′+= ???? ----(8)
----(9)比较 (7)(8)两式 ,可知
)()( 4hOhen = 因而三阶 R-K方法 (7)具有 3阶精度
)22(6 43211 KKKKhyy nn ++++= ?
),( 111 ??= nn yxfK
)2,2( 1112 KhyhxfK nn ++= ??
),( 3114 hKyhxfK nn ++= ??
)( 00 xyy =?
?
?
???
?
类似于 (7)式 ,还可构造 四阶 (经典 )Runge=Kutta方法
)2,2( 2113 KhyhxfK nn ++= ??
)()( 5hOhen =
-----------(10)
因而方法 (10)有 4阶精度
例 1. 使用高阶 R-K方法计算初值问题
?
?
?
=
≤≤=′
1)0(
5.002
y
xyy
.1.0=h取
解 : (1) 使用三阶 R-K方法 (7)
时1=n 201 yK = 1=
2
102 )2
1.0( KyK += 1025.1=
2
1203 ))2(1.0( KKyK ?+= 2555.1=
)4(61.0 32101 KKKyy +++= 1111.1=
RK.m
其余结果如下 :
(2) 如果使用四阶 R-K方法 (10)
时1=n 201 yK = 1=
2
102 )2
1.0( KyK += 1025.1=
n xn k1 k2 k3 yn
1.0000 0.1000 1.0000 1.1025 1.2555 1.1111
2.0000 0.2000 1.2345 1.3755 1.5945 1.2499
3.0000 0.3000 1.5624 1.7637 2.0922 1.4284
4.0000 0.4000 2.0404 2.3423 2.8658 1.6664
5.0000 0.5000 2.7768 3.2587 4.1634 1.9993
2
203 )2
1.0( KyK += 1133.1=
2
304 )1.0( KyK += 2351.1=
)22(61.0 432101 KKKKyy ++++= 1111.1=
其余结果如下 :
n xn k1 k2 k3 k4 yn
1.0000 0.1000 1.0000 1.1025 1.1133 1.2351 1.1111
2.0000 0.2000 1.2346 1.3756 1.3921 1.5633 1.2500
3.0000 0.3000 1.5625 1.7639 1.7908 2.0423 1.4286
4.0000 0.4000 2.0408 2.3428 2.3892 2.7805 1.6667
5.0000 0.5000 2.7777 3.2600 3.3476 4.0057 2.0000
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
1
1.2
1.4
1.6
1.8
2
2.2
精 确 解
Euler法
预 测 -校 正 系 统
改 进 Euler法
3阶 RK法
4阶 RK法
0.128 0.13 0.132 0.134 0.136 0.138 0.14 0.142 0.144 0.146 0.148
1.13
1.14
1.15
1.16
1.17
1.18
1.19
精 确 解
Euler法
预 测 -校 正 系 统
改 进 Euler法
3阶 RK法
4阶 RK法
0.2965 0.297 0.2975 0.298 0.2985 0.299 0.2995 0.3 0.3005
1.418
1.42
1.422
1.424
1.426
1.428
精 确 解
改 进 Euler法
3阶 RK法
4阶 RK法
