- 1苹果Vision pro 头显已发售，库克多次表示AR是终极方向？
- 2Python-VBA函数之旅-any函数
- 3成功解决UnicodeDecodeError: ‘utf-8‘ codec can‘t decode byte 0xc2 in position 0: invalid continuation byt_unicodedecodeerror: 'utf-8' codec can't decode byt
- 4安卓AS：简简单单就实现了一个qq登录界面_as登录界面
- 5致力微服务 SpringCloud 插曲 篇五 数据缓存-整合redis_springcloud结合缓存整合
- 6vivado 设置 VIO 核以执行测量、查看 VIO 核状态_set vio short activity persistence
- 7huggingface实战bert-base-chinese模型(训练+预测)
- 8git提交错误：non-fast-forward_git non-fast-forward
- 9免费社工库_@didisgk_bot
- 10记录一个Arduino调用MPU6050的姿态解算算法代码_mpu6050姿态解算代码 arduino
Runge-Kutta（龙格-库塔）方法 | 基本思想 + 二阶格式 + 四阶格式
赞
踩
Euler方法有各种格式，但其精度最高不超过2阶，一般难以满足实际计算的精度要求。因此，有必要构造精度更高的数值计算公式求解微分方程。Runge-Kutta方法就是一种高精度的经典的解常微分方程的单步方法。
1. Runge-Kutta 方法的基本思想
对于函数
y
=
y
(
x
)
y=y(x)
y=y(x)，设
y
n
=
y
(
x
n
)
y_n=y(x_n)
yn=y(xn)，将
y
(
x
n
+
1
)
y(x_{n+1})
y(xn+1)在点
x
n
x_n
xn处展开为Taylor级数：
y
(
x
n
+
1
)
=
y
(
x
n
)
+
h
y
′
(
ϵ
)
x
n
<
ϵ
<
x
n
+
1
y(x_{n+1})=y(x_n)+hy'(\epsilon) \quad x_n<\epsilon<x_{n+1}
y(xn+1)=y(xn)+hy′(ϵ)xn<ϵ<xn+1
利用
y
′
(
x
n
)
=
f
(
x
n
,
y
n
)
y'(x_n)=f(x_n,y_n)
y′(xn)=f(xn,yn)可得：
y
(
x
n
+
1
)
=
y
(
x
n
)
+
h
f
(
ϵ
,
y
(
ϵ
)
)
y(x_{n+1})=y(x_n)+hf(\epsilon,y(\epsilon))
y(xn+1)=y(xn)+hf(ϵ,y(ϵ))
其中，
f
(
ϵ
,
y
(
ϵ
)
)
f(\epsilon,y(\epsilon))
f(ϵ,y(ϵ))为区间
(
x
n
,
x
n
+
1
)
(x_n,x_{n+1})
(xn,xn+1)上的平均斜率，记为
k
∗
k^*
k∗。因此，只要能设法提供一种算法求得
k
∗
k^*
k∗，就可相应地导出一种计算格式。
回顾一下前面所讲的Euler方法可知：
若取 x n x_n xn点处的斜率 k 1 = f ( x n , y n ) k_1=f(x_n,y_n) k1=f(xn,yn)作为 k ∗ k^* k∗，则可得到显示Euler格式，但只有一阶精度；
若取 x n x_n xn和 x n + 1 x_{n+1} xn+1两点处的斜率 k 1 = f ( x n , y n ) k_1=f(x_n,y_n) k1=f(xn,yn)和 k 2 = f ( x n + 1 , y n + h k 1 ) k_2=f(x_{n+1},y_{n}+hk_1) k2=f(xn+1,yn+hk1)的算术平均值作为平均斜率 k o k^o ko，即 k o = ( k 1 + k 2 ) / 2 k^o=(k_1+k_2)/2 ko=(k1+k2)/2，则得到改进的Euler格式，则有二阶精度。
由此推而广之，如果设法在区间 ( x n , x n + 1 ) (x_n,x_{n+1}) (xn,xn+1)内取若干个点的斜率，然后将它们加权平均作为平均斜率 k o k^o ko，则可以构造一种具有更高精度的计算格式。这就是Runge-Kutta方法的基本思想。
2. 二阶Runge-Kutta格式
在改进的Euler格式的基础上加以修改，即在区间
[
x
n
,
x
n
+
1
]
[x_n,x_{n+1}]
[xn,xn+1]内，除
x
n
x_n
xn外再取区间中任一点p，则有：
x
n
+
p
=
x
n
+
p
h
0
<
p
≤
1
x_{n+p}=x_n+ph \quad 0<p\leq 1
xn+p=xn+ph0<p≤1
取
x
n
x_n
xn和
x
n
+
p
x_{n+p}
xn+p两点处的斜率
k
1
k_1
k1和
k
2
k_2
k2加权平均作为平均斜率
k
o
k^o
ko，即：
k
o
=
(
1
−
λ
)
k
1
+
λ
k
2
λ
为
待
定
系
数
k^o=(1-\lambda)k_1+\lambda k_2 \quad \lambda 为待定系数
ko=(1−λ)k1+λk2λ为待定系数
这样构造出的计算格式为：
y
n
+
1
=
y
n
+
h
[
(
1
−
λ
)
k
1
+
λ
k
2
]
(1a)
y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2] \tag{1a}
yn+1=yn+h[(1−λ)k1+λk2](1a)
其中：
k
1
=
f
(
x
n
,
y
n
)
(
1
b
)
k
2
=
f
(
x
n
+
p
,
y
n
+
p
h
k
1
)
(
1
c
)
k_1=f(x_n,y_n) \quad (1b)\\ k_2=f(x_{n+p},y_n+phk_1) \quad (1c)
k1=f(xn,yn)(1b)k2=f(xn+p,yn+phk1)(1c)
（1）式称为二阶的Runge-Kutta格式，其中
p
,
λ
p,\lambda
p,λ为待定参数。在确定
p
,
λ
p,\lambda
p,λ的同时，使二阶Runge-Kutta格式具有较高的精度，具体推导如下：
根据假设
y
n
=
y
(
x
n
)
y_n=y(x_n)
yn=y(xn)，有：
k
1
=
f
(
x
n
,
y
n
)
=
y
′
(
x
n
)
k
2
=
f
(
x
n
+
p
,
y
n
+
p
h
k
1
)
=
f
(
x
n
+
p
h
,
y
n
+
p
h
k
1
)
k_1=f(x_n,y_n)=y'(x_n) \\ k_2=f(x_{n+p},y_n+phk_1)=f(x_n+ph,y_n+phk_1)
k1=f(xn,yn)=y′(xn)k2=f(xn+p,yn+phk1)=f(xn+ph,yn+phk1)
将
k
2
k_2
k2右端在点
(
x
n
,
y
n
)
(x_n,y_n)
(xn,yn)处展开为Taylor级数：
k
2
=
f
(
x
n
,
y
n
)
+
p
h
⋅
f
x
(
x
n
,
y
n
)
+
p
h
y
′
(
x
n
)
⋅
f
y
(
x
n
,
y
n
)
+
O
(
h
2
)
k_2=f(x_n,y_n)+ph·f_x(x_n,y_n)+phy'(x_n)·f_y(x_n,y_n)+O(h^2)
k2=f(xn,yn)+ph⋅fx(xn,yn)+phy′(xn)⋅fy(xn,yn)+O(h2)
而
p
h
⋅
f
x
(
x
n
,
y
n
)
+
p
h
⋅
y
′
(
x
n
)
⋅
f
y
(
x
n
,
y
n
)
=
p
h
⋅
f
′
(
x
n
,
y
n
)
=
p
h
⋅
y
′
′
(
x
n
)
ph·f_x(x_n,y_n)+ph· y'(x_n)·f_y(x_n,y_n)=ph·f'(x_n,y_n)=ph·y''(x_n)
ph⋅fx(xn,yn)+ph⋅y′(xn)⋅fy(xn,yn)=ph⋅f′(xn,yn)=ph⋅y′′(xn)
故
k
2
=
f
(
x
n
,
y
n
)
+
p
h
y
′
′
(
x
n
)
+
O
(
h
2
)
=
y
′
(
x
n
)
+
p
h
y
′
′
(
x
n
)
+
O
(
h
2
)
k_2=f(x_n,y_n)+phy''(x_n)+O(h^2)=y'(x_n)+phy''(x_n)+O(h^2)
k2=f(xn,yn)+phy′′(xn)+O(h2)=y′(xn)+phy′′(xn)+O(h2)
将
k
1
k_1
k1和
k
2
k_2
k2代入
y
n
+
1
=
y
n
+
h
[
(
1
−
λ
)
k
1
+
λ
k
2
]
y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2]
yn+1=yn+h[(1−λ)k1+λk2]中，有：
y
n
+
1
=
y
n
+
h
[
(
1
−
λ
)
k
1
+
λ
k
2
]
=
y
(
x
n
)
+
h
y
′
(
x
n
)
+
λ
p
h
2
y
′
′
(
x
n
)
+
λ
O
(
h
3
)
y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2] \\ =y(x_n)+hy'(x_n)+\lambda ph^2y''(x_n)+\lambda O(h^3)
yn+1=yn+h[(1−λ)k1+λk2]=y(xn)+hy′(xn)+λph2y′′(xn)+λO(h3)
再将
y
(
x
n
+
1
)
y(x_{n+1})
y(xn+1)在点
(
x
n
,
y
n
)
(x_n,y_n)
(xn,yn)处展开为Taylor级数：
y
(
x
n
+
1
)
=
y
(
x
n
+
h
)
=
y
(
x
n
)
+
h
y
′
(
x
n
)
+
1
2
h
2
y
′
′
(
x
n
)
+
O
(
h
3
)
y(x_{n+1})=y(x_n+h)=y(x_n)+hy'(x_n)+\frac{1}{2}h^2y''(x_n)+O(h^3)
y(xn+1)=y(xn+h)=y(xn)+hy′(xn)+21h2y′′(xn)+O(h3)
比较两式，得：
y
(
x
n
+
1
)
−
y
n
+
1
=
(
1
2
h
2
−
λ
p
h
2
)
y
′
′
(
x
n
)
−
λ
O
(
h
3
)
y(x_{n+1})-y_{n+1}=(\frac{1}{2}h^2-\lambda ph^2)y''(x_n)-\lambda O(h^3)
y(xn+1)−yn+1=(21h2−λph2)y′′(xn)−λO(h3)
因此，当
λ
p
=
1
2
\lambda p=\frac{1}{2}
λp=21
时，二阶Runge-Kutta格式具有为2阶精度。可见二阶Runge-Kutta格式有很多具体形式：
若取 λ = 1 2 , p = 1 \lambda=\frac{1}{2},p=1 λ=21,p=1，则可得改进的Euler格式；
若取 λ = 3 4 , p = 2 3 \lambda=\frac{3}{4},p=\frac{2}{3} λ=43,p=32，则可得Heun（休恩）格式；
若取
λ
=
1
,
p
=
1
2
\lambda=1,p=\frac{1}{2}
λ=1,p=21，则可得变形的Euler格式，即中点格式：
{
y
n
+
1
=
y
n
+
h
k
2
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
1
2
,
y
n
+
h
2
k
1
)
3. 四阶Runge-Kutta格式
在区间
[
x
n
,
x
n
+
1
]
[x_n,x_{n+1}]
[xn,xn+1]内，使用两个不同的点可以构造出二阶Runge-Kutta格式。依此规律，在区间
[
x
n
,
x
n
+
1
]
[x_n,x_{n+1}]
[xn,xn+1]内，取3个不同的点可以构造出三阶Runge-Kutta格式；取4个不同的点可以构造出四阶Runge-Kutta格式。对此可以加以证明。在实际中，应用最广泛的是四阶经典的Runge-Kutta格式：
{
y
n
+
1
=
y
n
+
h
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
1
2
,
y
n
+
h
2
k
1
)
k
3
=
f
(
x
n
+
1
2
,
y
n
+
h
2
k
2
)
k
4
=
f
(
x
n
+
1
,
y
n
+
h
k
3
)
四阶Runge-Kutta方法的优点如下：
（1）它是一种高精度的单步法，可达四阶精度 O ( h 5 ) O(h^5) O(h5)；
（2）数值稳定性较好；
（3）只需知道一阶导数，无需明确定义或计算其他高阶导数；
（4）只需给出 y n y_n yn就能计算出 y n + 1 y_{n+1} yn+1，所以能够自启动；
（5）编程容易。
四阶Runge-Kutta法除经典的格式外，还有其他的格式。例如：
（1）Kutta格式
{
y
n
+
1
=
y
n
+
h
8
(
k
1
+
3
k
2
+
3
k
3
+
k
4
)
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
1
3
,
y
n
+
h
3
k
1
)
k
3
=
f
(
x
n
+
2
3
,
y
n
−
h
3
k
2
+
h
k
2
)
k
4
=
f
(
x
n
+
1
,
y
n
+
h
k
1
−
h
k
2
+
h
k
3
)
（2）Gill格式：
{
y
n
+
1
=
y
n
+
h
6
(
k
1
+
(
2
−
2
)
k
2
+
(
2
+
2
)
k
3
+
k
4
)
k
1
=
f
(
x
n
,
y
n
)
k
2
=
f
(
x
n
+
1
2
,
y
n
+
1
2
h
k
1
)
k
3
=
f
(
x
n
+
1
2
,
y
n
+
2
−
1
2
h
k
1
+
2
−
2
2
h
k
2
)
k
4
=
f
(
x
n
+
1
,
y
n
−
2
2
h
k
2
+
2
+
2
2
h
k
3
)
