[常微分方程的数值解法系列五] 龙格-库塔(RK4)法

在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态，其中主要就是求解微分方程。但之前求解微分方程的解析方法主要是应用于一些简单和特殊的微分方程求解中，对于一般形式的微分方程，一般很难用解析方法求出精确解，只能用数值方法求解。该系列主要介绍一些常用的常微分方程的数值解法，主要包括：

• [常微分方程的数值解法系列一] 常微分方程
• [常微分方程的数值解法系列二] 欧拉法
• [常微分方程的数值解法系列三] 改进欧拉法(预估校正法)
• [常微分方程的数值解法系列四] 中值法
• [常微分方程的数值解法系列五] 龙格-库塔(RK4)法

这个系列后面文章会用到前面文章的理论和技术，所以建议按照顺序查看。

简介

在之前常微分方程的数值解法系列中，已经介绍了欧拉法，改进欧拉法以及中值法等多种常微分方程的数值解法。但是之前讲解的方法的局部截断误差相对来说比较大，精度有限，在某些情况下需要更高精度的方法。
本来主要介绍的龙格-库塔法，可以提供更高精度的常微分方程的数值解法。而且龙格-库塔法是常微分方程的数值解法的基本理论，后面讲解的过程中会发现，之前讲解的多种方法都可以归类到龙格-库塔法的体系中。而且不考虑计算量的情况下，理论上是可以构造任意高阶的龙格-库塔公式，这里阶数可以认为和 [常微分方程的数值解法系列一] 常微分方程介绍的截断误差的阶数概念是一致的，阶数越高代表截断误差越小，精度越高。但比较常用的是四阶龙格-库塔公式，后面会给出具体详细的讲解。

基本思想

由微分中值定理可知，存在点$\xi$使得
$$\frac{\mathbf{x}(t_{i+1})-\mathbf{x}(t_i)}{\Delta t}={\mathbf{x}}^{\prime }(\xi ),\xi \in (t_i,t_{i+1})$$
即
$$\mathbf{x}(t_{i+1})=\mathbf{x}(t_i)+\Delta t{\mathbf{x}}^{\prime }(\xi )=\mathbf{x}(t_i)+\Delta t\cdot f(\xi ,\mathbf{x}(\xi ))=\mathbf{x}(t_i)+\Delta t\bar{k},\text{\hspace{1em}(1)}$$
其中$\bar{k}=f(\xi ,\mathbf{x}(\xi ))$称为$\mathbf{x}(t)$在$[t_i,t_{i+1}]$上的平均斜率。所以基于这点考虑，只要对平均斜率$\bar{k}$提供一种近似算法，那么就可以利用$x_i$去递推求解$x_{i+1}$，这就是龙格-库塔方法的基本思想。
当然不同的近似算法求解精度不同，像之前讲解的向前欧拉法、改进欧拉法和中值法利用龙格-库塔方法思想解释就是分别利用$t_i$点斜率、$t_i$点和$t_{i+1}$点斜率平均值和$t_i$点和$t_{i+1}$点间中点斜率来近似$\bar{k}$。

具体方法

一阶

以$\mathbf{x}(t)$在$t_i$点斜率作为$\mathbf{x}(t)$在$[t_i,t_{i+1}]$上的平均斜率$\bar{k}$，并令$x_i\approx \mathbf{x}(t_i)$，可得
$$\bar{k}={\mathbf{x}}^{\prime }(t_i)=f(t_i,\mathbf{x}(t_i))\approx f(t_i,x_i)$$
所以
$$x_{i+1}=x_i+\Delta t\cdot f(t_i,x_i),\text{\hspace{1em}(2)}$$
这就是[常微分方程的数值解法系列二] 欧拉法文中讲解的向前欧拉法，也叫作一阶龙格-库塔方法。公式$(2)$就是一阶龙格-库塔公式。

$x_i$是递推过程中$\mathbf{x}(t_i)$的近似值，后面都用$x_n$表示$\mathbf{x}(t_n)$。

二阶

在$[t_i,t_{i+1}]$上取两点$t_i,t_{i+p}(0<p\le 1)$的斜率值分别为$k_1,k_2$，把它们的线性组合${\lambda }_1{k}_1+{\lambda }_2{k}_2$作为$\bar{k}$的近似值，其中${\lambda }_1,{\lambda }_2,p$都是待确定常数，所以根据公式$(1)$可得
$$x_{i+1}=x_i+\Delta t({\lambda }_1{k}_1+{\lambda }_2{k}_2)$$
其中$k_1=f(t_i,x_i)$，$k_2$是$(t_i,t_{i+1}]$内任意一点$t_{i+p}=t_i+p\Delta t(0<p\le 1)$的斜率 k 2 = f ( t i + p , x i + p ) k_2 = f(t_{i+p}, x_{i+p})