龙格库塔算法原理详解_龙格库塔法的基本原理

龙格库塔算法

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

0. 问题提出

求解下面的一阶微分方程
{ y ′ ( x ) = d y d x = f ( x , y ) y ( x 0 ) = y 0 \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=x_0$ 时 $y=y_0$ ，且斜率为 $y^{\prime }(x_0)=f(x_0,y_0)$ ，那么就可以用过该点且斜率为$y^{\prime }(x_0)$的直线近似求解$x_0$附近一点$x_1$的函数值。

构建直线：
$$y(x)=y_0+f(x_0,y_0)(x-x_0)$$
那么，在$x=x_1$时，求得
$$y(x_1)=y_0+f(x_0,y_0)(x_1-x_0)$$
即在 $x=x_1$ 时，$y=y(x_1)$，且斜率为 $y^{\prime }(x_1)=f(x_1,y_1)$。重复以上步骤，可以求得 $y(x_2),y(x_3),\cdots ,y(x_n)$

如果等间隔取 x 值，即 $x_1-x_0=x_2-x_1=x_3-x_2=\cdots =x_n-x_{n-1}=h$，也即$x_n=x_0+n\cdot h$可以得到
$$\begin{gathered} y(x_1)=y_0+f(x_0,y_0)\cdot h \\ y(x_2)=y_1+f(x_1,y_1)\cdot h \\ ⋮ \\ y(x_n)=y_{n-1}+f(x_{n-1},y_{n-1})\cdot h \end{gathered}$$
总结为：
{ y ( x n ) = y n − 1 + f ( x n − 1 , y n − 1 ) ⋅ h x n = x 0 + 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(x_1)$ 是有误差的，而在求解 $y(x_2)$ 的时候，虽然斜率是准确的，但函数值是基于 $y(x_1)$ 的，即在第一次误差的基础上计算的，因此会累积误差。

可以使用泰勒展开分析每次迭代产生的误差。对 $y(x)$ 进行泰勒展开
$$y(x)=y(x_k)+y^{\prime }(x_k)(x-x_k)+\frac{y^{\prime \prime }(x_k)}{2!}(x-x_k{)}^2+\cdots +\frac{y^{n+1}(ϵ)}{(n+1)!}(x-x_k{)}^{n+1}$$
则对于 $y(x_{k+1})$ 有
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 ) ⋅ h + O ( h 2 )

$$\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(x_k)+y^{\prime }(x_k)\cdot h$ 即为用欧拉法的 $y_{k+1}$的估计值，所以可以看出，每步误差即为 $O(h^2)$ ，也即使用欧拉法具有一阶精度。该误差也被称为局部截断误差。当然，步长 h 越小，误差越小。

2. 龙格库塔算法

2.1 基本原理

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

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

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

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

在$(x_k,x_{k+1})$ 区间内存在一点 $ϵ(x_k<ϵ<x_{k+1})$ ，使得等式 $f(x_{k+1})-f(x_k)=f^{\prime }(ϵ)(x_{k+1}-x_k)$ 成立。

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

2.2 精度问题

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

改进欧拉法

既然直接使用端点处$y^{\prime }(x_k)$ 的斜率不准确，那使用两个端点斜率取平均，是否准确度就提升了一点呢，即取
$$y^{\prime }(ϵ)=\frac{y^{\prime }(x_k)+y^{\prime }(x_{k+1})}{2}$$
直观感受，应该会如此。为了方便书写，将$y^{\prime }(x_k)={\kappa }_1，y^{\prime }(x_{k+1}）={\kappa }_2$ ，$\kappa$ 即表示斜率。则上式
$$y^{\prime }(ϵ)=\frac{1}{2}({\kappa }_1+{\kappa }_2)$$
根据已知条件方程 $y^{\prime }(x)=f(x,y)$，则
$$\begin{gathered} {\kappa }_1=f(x_k,y_k) \\ {\kappa }_2=f(x_{k+1},y_{k+1}) \end{gathered}$$
f(x,y)是已知的，而值得注意的是，$f(x_{k+1},y_{k+1})$ 中 $y_{k+1}$ 是未知的，待求的，此处的$y_{k+1}$可以用欧拉法来估计，即
$$y_{k+1}=y_k+y^{\prime }(x_k)\cdot h=y_k+{\kappa }_{1}h$$
当然这是不准确的，但精度影响肯定小于上述的欧拉法（如何减少影响，下述），此时
$$\begin{gathered} {\kappa }_1=f(x_k,y_k) \\ {\kappa }_2=f(x_{k+1},y_k+{\kappa }_1\cdot h) \end{gathered}$$
总结一下，在使用端点斜率的均值当做直线斜率思路下，递推求解方程如下
{ y k + 1 = y k + ( 1 2 κ 1 + 1 2 κ 2 ) ⋅ h κ 1 = f ( x k , y k ) κ 2 = f ( x k + 1 , y k + κ 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(h^3)$ 。

事实上，对于两端点斜率取平均完全是出于习惯，更一般描述是对两斜率进行加权，取平均只是习惯性的把两者权重取相等而已，但实际应用中，不一定权重相等更精确，因此更一般描述为
{ y k + 1 = y k + ( λ 1 κ 1 + λ 2 κ 2 ) ⋅ h κ 1 = f ( x k , y k ) κ 2 = f ( x k + 1 , y k + κ 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^{\prime }(ϵ)$ 的值。

假设在$(x_k,x_{k+1})$ 区间取n个点，对应为 $x_k+{\alpha }_{1}h,x_k+{\alpha }_{2}h,\cdots ,x_k+{\alpha }_{n}h$，且$\alpha \in [0,1]$，在每个点的斜率记为${\kappa }_1,{\kappa }_2,\cdots ,{\kappa }_n$，对应加权为${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$那么有
$$f^{\prime }(ϵ)={\lambda }_1{\kappa }_1+{\lambda }_2{\kappa }_2+\cdots +{\lambda }_n{\kappa }_n$$
则
$$y_{k+1}=y_k+({\lambda }_1{\kappa }_1+{\lambda }_2{\kappa }_2+\cdots +{\lambda }_n{\kappa }_n)\cdot h$$
$\lambda$ 为自取值，对于$\kappa$ 有
$$\begin{gathered} {\kappa }_1=f^{\prime }(x_k+{\alpha }_{1}h)=f(x_k+{\alpha }_{1}h,y(x_k+{\alpha }_{1}h)) \\ {\kappa }_2=f^{\prime }(x_k+{\alpha }_{2}h)=f(x_k+{\alpha }_{2}h,y(x_k+{\alpha }_{2}h)) \\ ⋮ \\ {k}_n=f^{\prime }(x_k+{\alpha }_{n}h)=f(x_k+{\alpha }_{n}h,y(x_k+{\alpha }_{n}h)) \end{gathered}$$
通常${\alpha }_1$取0，即对应$x_k$点；此时，${\kappa }_1=f(x_k,y_k)$

同样，对于上面其余式子 f 自变量中的 $y(x_k+{\alpha }_{n}h)$ 使用欧拉法来求，即
$$\begin{gathered} {\kappa }_1=f(x_k,y_k) \\ y(x_k+{\alpha }_{2}h)=y(x_k)+{\alpha }_{2}h\cdot f^{\prime }(x_k)=y_k+{\alpha }_{2}h\cdot {\kappa }_1 \\ y(x_k+{\alpha }_{3}h)=y(x_k+{\alpha }_{2}h)+({\alpha }_3-{\alpha }_2)h\cdot {\kappa }_2=y_k+{\alpha }_{2}h\cdot {\kappa }_1+({\alpha }_3-{\alpha }_2)h\cdot {\kappa }_2 \\ ⋮ \\ y(x_k+{\alpha }_{n}h)=y_k+{\alpha }_{2}h\cdot h+({\alpha }_3-{\alpha }_2)h\cdot {\kappa }_2+\cdots +({\alpha }_n-{\alpha }_{n-1})h\cdot {\kappa }_{n-1} \end{gathered}$$
值得注意的是，对于欧拉法来说，两点间距越小，误差越小。对应于该处，取样点越多，间距越小，此处欧拉法误差越小。

总结以上，即
{ y k + 1 = y k + ( λ 1 κ 1 + λ 2 κ 2 + ⋯ + λ n κ n ) ⋅ h κ 1 = f ( x k , y k ) κ 2 = f ( x k + α 1 , y k + α 1 h κ 1 ) κ 3 = f ( x k + α 2 , y k + α 1 h κ 1 + ( α 2 − α 1 ) h κ 2 ) ⋮ κ n = f ( x k + α n , y k + α 1 h κ 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​+α1​hκ1​)κ3​=f(xk​+α2​,yk​+α1​hκ1​+(α2​−α1​)hκ2​)κn​=f(xk​+αn​,yk​+α1​hκ1​+(α2​−α1​)hκ2​+⋯+(αn​−αn−1​)hκn−1​)​
其中：

${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$ 为对应采样点斜率权值，满足 ${\lambda }_1+{\lambda }_2+\cdots +{\lambda }_n=1$；

${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n$ 为采样点位置比率，满足$0<{\alpha }_1<{\alpha }_2<\cdots <{\alpha }_n$；

实际应用

采样两个点来估算斜率，其中一个为端点$（x_k,y_k）$，另一个点取${\alpha }_1$ 。

可以得到二阶龙格库塔公式
{ y k + 1 = y k + ( λ 1 κ 1 + λ 2 κ 2 ) h κ 1 = f ( x k , y k ) κ 2 = f ( x k + α 1 h , y k + α 1 h κ 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​+α1​h,yk​+α1​hκ1​)​
依次类推，采样适当数量的点，设置适当的$\alpha 和\lambda$ 就可以得到三阶、四阶等更高精度的龙格库塔公式。具体推导不赘述。

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