龙格库塔算法

前言: 本文从简单欧拉法讲起，详细探究了龙格库塔算法的原理。

文章目录

龙格库塔算法
- 0. 问题提出
- 1. 欧拉法
- - 1.1 基本原理
  - 1.2 误差分析
- 2. 龙格库塔算法
- - 2.1 基本原理
  - 2.2 精度问题

0. 问题提出

求解下面的一阶微分方程
{y′(x)=dydx=f(x,y)y(x0)=y0\left\{\begin{aligned} &y'(x)=\frac{dy}{dx} = f(x,y) \\ &y(x_0) = y_0 \end{aligned} \right. ⎩⎨⎧y′(x)=dxdy=f(x,y)y(x0)=y0
即已知上述条件，则 y(x) 如何求出？

1. 欧拉法

1.1 基本原理

上述方程组给出了某点的坐标，以及该点的斜率值。可以使用欧拉法对 y(x) 进行迭代的近似求解。求解分析如下：

已知 x=x0x = x_0x=x0 时 y=y0y=y_0y=y0 ，且斜率为 y′(x0)=f(x0,y0)y'(x_0)=f(x_0,y_0)y′(x0)=f(x0,y0) ，那么就可以用过该点且斜率为y′(x0)y'(x_0)y′(x0)的直线近似求解x0x_0x0附近一点x1x_1x1的函数值。

构建直线：
y(x)=y0+f(x0,y0)(x−x0)y(x) = y_0 + f(x_0,y_0)(x-x_0) y(x)=y0+f(x0,y0)(x−x0)
那么，在x=x1x=x_1x=x1时，求得
y(x1)=y0+f(x0,y0)(x1−x0)y(x_1) = y_0+f(x_0,y_0)(x_1-x_0) y(x1)=y0+f(x0,y0)(x1−x0)
即在 x=x1x=x_1x=x1 时，y=y(x1)y=y(x_1)y=y(x1)，且斜率为 y′(x1)=f(x1,y1)y'(x_1)=f(x_1,y_1)y′(x1)=f(x1,y1)。重复以上步骤，可以求得 y(x2),y(x3),⋯,y(xn)y(x_2),y(x_3),\cdots,y(x_n)y(x2),y(x3),⋯,y(xn)

如果等间隔取 x 值，即 x1−x0=x2−x1=x3−x2=⋯=xn−xn−1=hx_1-x_0 = x_2-x_1=x_3-x_2=\cdots=x_{n}-x_{n-1}=hx1−x0=x2−x1=x3−x2=⋯=xn−xn−1=h，也即xn=x0+n⋅hx_n=x_0+n\cdot hxn=x0+n⋅h可以得到
y(x1)=y0+f(x0,y0)⋅hy(x2)=y1+f(x1,y1)⋅h⋮y(xn)=yn−1+f(xn−1,yn−1)⋅hy(x_1)=y_0+f(x_0,y_0)\cdot h \\ y(x_2)=y_1+f(x_1,y_1)\cdot h \\ \vdots\\ y(x_n)=y_{n-1}+f(x_{n-1},y_{n-1})\cdot h y(x1)=y0+f(x0,y0)⋅hy(x2)=y1+f(x1,y1)⋅h⋮y(xn)=yn−1+f(xn−1,yn−1)⋅h
总结为：
{y(xn)=yn−1+f(xn−1,yn−1)⋅hxn=x0+n⋅h\left\{\begin{aligned} &y(x_n)=y_{n-1}+f(x_{n-1},y_{n-1})\cdot h \\ &x_n = x_0 + n\cdot h \end{aligned} \right. {y(xn)=yn−1+f(xn−1,yn−1)⋅hxn=x0+n⋅h
该公式就是欧拉公式。

1.2 误差分析

值得注意的是，如上求解 y(x1)y(x_1)y(x1) 是有误差的，而在求解 y(x2)y(x_2)y(x2) 的时候，虽然斜率是准确的，但函数值是基于 y(x1)y(x_1)y(x1) 的，即在第一次误差的基础上计算的，因此会累积误差。

可以使用泰勒展开分析每次迭代产生的误差。对 y(x)y(x)y(x) 进行泰勒展开
y(x)=y(xk)+y′(xk)(x−xk)+y′′(xk)2!(x−xk)2+⋯+yn+1(ϵ)(n+1)!(x−xk)n+1y(x) = y(x_k) +y'(x_k)(x-x_k)+\frac{y''(x_k)}{2!}(x-x_k)^2+\cdots+\frac{y^{n+1}(\epsilon)}{(n+1)!}(x-x_k)^{n+1} \\ y(x)=y(xk)+y′(xk)(x−xk)+2!y′′(xk)(x−xk)2+⋯+(n+1)!yn+1(ϵ)(x−xk)n+1
则对于 y(xk+1)y(x_{k+1})y(xk+1) 有
y(xk+1)=y(xk)+y′(xk)(xk+1−xk)+O((xk+1−xk)2)=y(xk)+y′(xk)⋅h+O(h2)\begin{aligned} y(x_{k+1}) &= y(x_k) + y'(x_k)(x_{k+1}-x_k)+O((x_{k+1}-x_k)^2) \\ &=y(x_k)+y'(x_k)\cdot h+O(h^2) \\ \end{aligned} y(xk+1)=y(xk)+y′(xk)(xk+1−xk)+O((xk+1−xk)2)=y(xk)+y′(xk)⋅h+O(h2)
等式右侧y(xk)+y′(xk)⋅hy(x_k)+y'(x_k)\cdot hy(xk)+y′(xk)⋅h 即为用欧拉法的 yk+1y_{k+1}yk+1的估计值，所以可以看出，每步误差即为 O(h2)O(h^2)O(h2) ，也即使用欧拉法具有一阶精度。该误差也被称为局部截断误差。当然，步长 h 越小，误差越小。

2. 龙格库塔算法

2.1 基本原理

该算法有些类似于欧拉法的思想，算是对其的改进，都是使用过(xk,yk)(x_k,y_k)(xk,yk)的直线来近似下一点的函数值yk+1y_{k+1}yk+1，区别在于斜率的选择。

在欧拉法中，直接使用该点的切线y′(xk)y'(x_k)y′(xk) 当做直线的斜率，这样会造成一个明显的结果，即对于 y(xk+1)y(x_{k+1})y(xk+1) 的求解一定会有误差，除非y(x)为直线时，结果才是准确的。

基于上面的考虑，最终需要解决的便是，如何得到一条通过(xk,yk)(x_k,y_k)(xk,yk)和下一点(xk+1,yk+1)(x_{k+1},y_{k+1})(xk+1,yk+1)的直线呢？或者说在知道(xk,yk)(x_k,y_k)(xk,yk)的情况下，如何确定直线斜率，确保恰好经过待估测点(xk+1,yk+1)(x_{k+1},y_{k+1})(xk+1,yk+1) 呢？

答案是利用拉格朗日中值定理：

在(xk,xk+1)(x_k,x_{k+1})(xk,xk+1) 区间内存在一点 ϵ(xk<ϵ<xk+1)\epsilon (x_{k}<\epsilon<x_{k+1})ϵ(xk<ϵ<xk+1) ，使得等式 f(xk+1)−f(xk)=f′(ϵ)(xk+1−xk)f(x_{k+1})-f(x_k)=f'(\epsilon)(x_{k+1}-x_{k})f(xk+1)−f(xk)=f′(ϵ)(xk+1−xk) 成立。

所以准确的直线方程即为：
y(x)=y(xk)+y′(ϵ)(x−xk)=y(xk)+y′(xk+θh)(x−xk)y(x) = y(x_k) + y'(\epsilon)(x-x_k) = y(x_k)+y'(x_k+\theta h)(x-x_k) y(x)=y(xk)+y′(ϵ)(x−xk)=y(xk)+y′(xk+θh)(x−xk)
如果等间距取值，即xk+1−xk=hx_{k+1}-x_k=hxk+1−xk=h ，则有
y(xk+1)=y(xk)+y′(xk+θh)⋅hy(x_{k+1}) = y(x_k)+y'(x_k+\theta h)\cdot h y(xk+1)=y(xk)+y′(xk+θh)⋅h
且θ∈(0,1)\theta\in(0,1)θ∈(0,1) ，这种方法可以确保估测的y(xk+1)y(x_{k+1})y(xk+1)是极其准确，而另 θ=0\theta=0θ=0 的时候，其实就是欧拉法。

2.2 精度问题

原理已经明确了，下面就是如何确定 y′(ϵ)y'(\epsilon)y′(ϵ) 的问题了，这就涉及到了龙格库塔算法的精度问题。

改进欧拉法

既然直接使用端点处y′(xk)y'(x_k)y′(xk) 的斜率不准确，那使用两个端点斜率取平均，是否准确度就提升了一点呢，即取
y′(ϵ)=y′(xk)+y′(xk+1)2y'(\epsilon) = \frac{y'(x_k)+y'(x_{k+1})}{2} y′(ϵ)=2y′(xk)+y′(xk+1)
直观感受，应该会如此。为了方便书写，将y′(xk)=κ1，y′(xk+1）=κ2y'(x_k)=\kappa_1，y'(x_{k+1}）=\kappa_2y′(xk)=κ1，y′(xk+1）=κ2 ，κ\kappaκ 即表示斜率。则上式
y′(ϵ)=12(κ1+κ2)y'(\epsilon)= \frac{1}{2}(\kappa_1+\kappa_2) y′(ϵ)=21(κ1+κ2)
根据已知条件方程 y′(x)=f(x,y)y'(x)=f(x,y)y′(x)=f(x,y)，则
κ1=f(xk,yk)κ2=f(xk+1,yk+1)\kappa_1 = f(x_k,y_k) \\ \kappa_2 = f(x_{k+1},y_{k+1}) κ1=f(xk,yk)κ2=f(xk+1,yk+1)
f(x,y)是已知的，而值得注意的是，f(xk+1,yk+1)f(x_{k+1},y_{k+1})f(xk+1,yk+1) 中 yk+1y_{k+1}yk+1 是未知的，待求的，此处的yk+1y_{k+1}yk+1可以用欧拉法来估计，即
yk+1=yk+y′(xk)⋅h=yk+κ1hy_{k+1}=y_k+y'(x_k)\cdot h=y_k+\kappa_1 h yk+1=yk+y′(xk)⋅h=yk+κ1h
当然这是不准确的，但精度影响肯定小于上述的欧拉法（如何减少影响，下述），此时
κ1=f(xk,yk)κ2=f(xk+1,yk+κ1⋅h)\kappa_1 = f(x_k,y_k) \\ \kappa_2 = f(x_{k+1},y_k+\kappa_1\cdot h) κ1=f(xk,yk)κ2=f(xk+1,yk+κ1⋅h)
总结一下，在使用端点斜率的均值当做直线斜率思路下，递推求解方程如下
{yk+1=yk+(12κ1+12κ2)⋅hκ1=f(xk,yk)κ2=f(xk+1,yk+κ1⋅h)\left\{\begin{aligned} &y_{k+1} = y_k+(\frac{1}{2}\kappa_1+\frac{1}{2}\kappa_2)\cdot h \\ &\kappa_1 = f(x_k,y_k) \\ &\kappa_2 = f(x_{k+1},y_k+\kappa_1\cdot h) \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧yk+1=yk+(21κ1+21κ2)⋅hκ1=f(xk,yk)κ2=f(xk+1,yk+κ1⋅h)
该种方法称为改进欧拉法，精度可以达到 O(h3)O(h^3)O(h3) 。

事实上，对于两端点斜率取平均完全是出于习惯，更一般描述是对两斜率进行加权，取平均只是习惯性的把两者权重取相等而已，但实际应用中，不一定权重相等更精确，因此更一般描述为
{yk+1=yk+(λ1κ1+λ2κ2)⋅hκ1=f(xk,yk)κ2=f(xk+1,yk+κ1⋅h)\left\{\begin{aligned} &y_{k+1} = y_k+(\lambda_1\kappa_1+\lambda_2\kappa_2)\cdot h \\ &\kappa_1 = f(x_k,y_k) \\ &\kappa_2 = f(x_{k+1},y_k+\kappa_1\cdot h) \end{aligned} \right. ⎩⎪⎨⎪⎧yk+1=yk+(λ1κ1+λ2κ2)⋅hκ1=f(xk,yk)κ2=f(xk+1,yk+κ1⋅h)
通用方法

依据上述思路，如何更精确一点呢？方法是用更多点加权来估算f′(ϵ)f'(\epsilon)f′(ϵ) 的值。

假设在(xk,xk+1)(x_k,x_{k+1})(xk,xk+1) 区间取n个点，对应为 xk+α1h,xk+α2h,⋯,xk+αnhx_k+\alpha_1 h,\ x_k+\alpha_2 h,\ \cdots,\ x_k+\alpha_n hxk+α1h, xk+α2h, ⋯, xk+αnh，且α∈[0,1]\alpha\in[0,1]α∈[0,1]，在每个点的斜率记为κ1,κ2,⋯,κn\kappa_1,\kappa_2,\cdots,\kappa_nκ1,κ2,⋯,κn，对应加权为λ1,λ2,⋯,λn\lambda_1,\lambda_2,\cdots,\lambda_nλ1,λ2,⋯,λn那么有
f′(ϵ)=λ1κ1+λ2κ2+⋯+λnκnf'(\epsilon) = \lambda_1\kappa_1+\lambda_2\kappa_2+\cdots+\lambda_n\kappa_n f′(ϵ)=λ1κ1+λ2κ2+⋯+λnκn
则
yk+1=yk+(λ1κ1+λ2κ2+⋯+λnκn)⋅hy_{k+1} = y_k+(\lambda_1\kappa_1+\lambda_2\kappa_2+\cdots+\lambda_n\kappa_n)\cdot h yk+1=yk+(λ1κ1+λ2κ2+⋯+λnκn)⋅h
λ\lambdaλ 为自取值，对于κ\kappaκ 有
κ1=f′(xk+α1h)=f(xk+α1h,y(xk+α1h))κ2=f′(xk+α2h)=f(xk+α2h,y(xk+α2h))⋮kn=f′(xk+αnh)=f(xk+αnh,y(xk+αnh))\kappa_1 = f'(x_k+\alpha_1 h)=f(x_k+\alpha_1h,y(x_k+\alpha_1h)) \\ \kappa_2 = f'(x_k+\alpha_2h)=f(x_k+\alpha_2h,y(x_k+\alpha_2h)) \\ \vdots \\ k_n=f'(x_k+\alpha_nh) = f(x_k+\alpha_nh,y(x_k+\alpha_nh)) κ1=f′(xk+α1h)=f(xk+α1h,y(xk+α1h))κ2=f′(xk+α2h)=f(xk+α2h,y(xk+α2h))⋮kn=f′(xk+αnh)=f(xk+αnh,y(xk+αnh))
通常α1\alpha_1α1取0，即对应xkx_kxk点；此时，κ1=f(xk,yk)\kappa_1=f(x_k,y_k)κ1=f(xk,yk)

同样，对于上面其余式子 f 自变量中的 y(xk+αnh)y(x_k+\alpha_nh)y(xk+αnh) 使用欧拉法来求，即
κ1=f(xk,yk)y(xk+α2h)=y(xk)+α2h⋅f′(xk)=yk+α2h⋅κ1y(xk+α3h)=y(xk+α2h)+(α3−α2)h⋅κ2=yk+α2h⋅κ1+(α3−α2)h⋅κ2⋮y(xk+αnh)=yk+α2h⋅h+(α3−α2)h⋅κ2+⋯+(αn−αn−1)h⋅κn−1\kappa_1=f(x_k,y_k) \\ y(x_k+\alpha_2h) = y(x_k)+\alpha_2h\cdot f'(x_k)=y_k+\alpha_2h\cdot\kappa_1 \\ y(x_k+\alpha_3h) = y(x_k+\alpha_2h)+(\alpha_3-\alpha_2)h\cdot \kappa_2=y_k+\alpha_2h\cdot\kappa_1 + (\alpha_3-\alpha_2)h\cdot\kappa_2 \\ \vdots \\ y(x_k+\alpha_n h)=y_k+\alpha_2h\cdot h+(\alpha_3-\alpha_2)h\cdot \kappa_2+\cdots+(\alpha_n-\alpha_{n-1})h\cdot\kappa_{n-1} κ1=f(xk,yk)y(xk+α2h)=y(xk)+α2h⋅f′(xk)=yk+α2h⋅κ1y(xk+α3h)=y(xk+α2h)+(α3−α2)h⋅κ2=yk+α2h⋅κ1+(α3−α2)h⋅κ2⋮y(xk+αnh)=yk+α2h⋅h+(α3−α2)h⋅κ2+⋯+(αn−αn−1)h⋅κn−1
值得注意的是，对于欧拉法来说，两点间距越小，误差越小。对应于该处，取样点越多，间距越小，此处欧拉法误差越小。

总结以上，即
{yk+1=yk+(λ1κ1+λ2κ2+⋯+λnκn)⋅hκ1=f(xk,yk)κ2=f(xk+α1,yk+α1hκ1)κ3=f(xk+α2,yk+α1hκ1+(α2−α1)hκ2)⋮κn=f(xk+αn,yk+α1hκ1+(α2−α1)hκ2+⋯+(αn−αn−1)hκn−1)\left\{\begin{aligned} &y_{k+1} = y_k+(\lambda_1\kappa_1+\lambda_2\kappa_2+\cdots+\lambda_n\kappa_n)\cdot h \\ &\kappa_1 = f(x_k,y_k) \\ &\kappa_2 = f(x_k+\alpha_1,y_k+\alpha_1 h\kappa_1) \\ &\kappa_3 = f(x_k+\alpha_2,y_k+\alpha_1h\kappa_1+(\alpha_2-\alpha_1)h\kappa_2)\\ \vdots \\ &\kappa_n = f(x_k+\alpha_n,y_k+\alpha_1h\kappa_1+(\alpha_2-\alpha_1)h\kappa_2+\cdots+(\alpha_n-\alpha_{n-1})h\kappa_{n-1}) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧⋮yk+1=yk+(λ1κ1+λ2κ2+⋯+λnκn)⋅hκ1=f(xk,yk)κ2=f(xk+α1,yk+α1hκ1)κ3=f(xk+α2,yk+α1hκ1+(α2−α1)hκ2)κn=f(xk+αn,yk+α1hκ1+(α2−α1)hκ2+⋯+(αn−αn−1)hκn−1)
其中：

λ1,λ2,⋯,λn\lambda_1,\lambda_2,\cdots,\lambda_nλ1,λ2,⋯,λn 为对应采样点斜率权值，满足 λ1+λ2+⋯+λn=1\lambda_1+\lambda_2+\cdots+\lambda_n=1λ1+λ2+⋯+λn=1；

α1,α2,⋯,αn\alpha_1,\alpha_2,\cdots,\alpha_nα1,α2,⋯,αn 为采样点位置比率，满足0<α1<α2<⋯<αn0<\alpha_1<\alpha_2<\cdots<\alpha_n0<α1<α2<⋯<αn；

实际应用

采样两个点来估算斜率，其中一个为端点（xk,yk）（x_k,y_k）（xk,yk），另一个点取α1\alpha_1α1 。

可以得到二阶龙格库塔公式
{yk+1=yk+(λ1κ1+λ2κ2)hκ1=f(xk,yk)κ2=f(xk+α1h,yk+α1hκ1)\left\{\begin{aligned} &y_{k+1} = y_k+(\lambda_1\kappa_1+\lambda_2\kappa_2)h \\ &\kappa_1 = f(x_k,y_k) \\ &\kappa_2 = f(x_k+\alpha_1h,y_k+\alpha_1h\kappa_1) \end{aligned} \right. ⎩⎪⎨⎪⎧yk+1=yk+(λ1κ1+λ2κ2)hκ1=f(xk,yk)κ2=f(xk+α1h,yk+α1hκ1)
依次类推，采样适当数量的点，设置适当的α和λ\alpha和\lambdaα和λ 就可以得到三阶、四阶等更高精度的龙格库塔公式。具体推导不赘述。

仅需要说明的是，对于更高阶常用龙格库塔公式，对于采样点的选取以及采样点函数值的计算并非一定按照上述方式计算。如采样点选取可能会选取重复点，但采用不同函数值计算方案（以上都是基于前一点，有些就基于更前面的点），总之十分灵活，但原理和思想都是一样的。