[常微分方程的數值解法系列五] 龍格-庫塔(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})$。
但這裏在求解過程當中$x_{i+p}$是未知量，因此須要先求解出，能夠仿照改進歐拉公式的思想，得
$$x_{i+p}=x_i+p\Delta tf(t_i,x_i)=x_i+p\Delta t{k}_1$$
因此整理得
$$\{\begin{array}{l}{x}_{i+1}=x_i+\Delta t({\lambda }_1{k}_1+{\lambda }_2{k}_2)\\ k_1=f(t_i,x_i)\\ k_2=f(t_i+p\Delta t,x_i+p\Delta t{k}_1)\end{array}\text{\hspace{1em}(3)}$$
其中公式$(3)$有三個待確認參數${\lambda }_1,{\lambda }_2,p$，利用泰勒展開選擇合適的參數使上訴公式具備二階精度，即局部截斷偏差爲$\mathbf{O}(\Delta t^3)$，知足這類條件的公式$(3)$就是二階龍格-庫塔公式。

求解參數

下面利用泰勒展開求解待確認參數${\lambda }_1,{\lambda }_2,p$，使公式$(3)$的局部截斷偏差爲$\mathbf{O}(\Delta t^3)$。
分別對公式$3$的$k_1,k_2$做泰勒展開爲
$$\begin{array}{rl}{k}_1& =f(t_i,x_i)={\mathbf{x}}^{\prime }(t_i)\\ k_2& =f(t_i+p\Delta t,x_i+p\Delta t{k}_1)=f(t_i+p\Delta t,x_i+p\Delta t\cdot f(t_i,x_i))\\ & =f(t_i,x_i)+p\Delta t{f}_t^{\prime }(t_i,x_i)+(p\Delta t\cdot f(t_i,x_i))\cdot f_x^{\prime }(t_i,x_i)+\mathbf{O}(\Delta t^2)\\ & ={\mathbf{x}}^{\prime }(t_i)+p\Delta t{\mathbf{x}}^{\prime \prime }(t_i)+\mathbf{O}(\Delta t^2)\end{array}$$
因此可得
$$\begin{array}{rl}{x}_{i+1}& =x_i+\Delta t({\lambda }_1{k}_1+{\lambda }_2{k}_2)\\ & =x_i+\Delta t({\lambda }_1{\mathbf{x}}^{\prime }(t_i)+{\lambda }_2({\mathbf{x}}^{\prime }(t_i)+p\Delta t{\mathbf{x}}^{\prime \prime }(t_i)+\mathbf{O}(\Delta t^2)))\\ & =x_i+\Delta t({\lambda }_1+{\lambda }_2){\mathbf{x}}^{\prime }(t_i)+{\lambda }_{2}p\Delta t^2{\mathbf{x}}^{\prime \prime }(t_i)+\mathbf{O}(\Delta t^3)\end{array},\text{\hspace{1em}(4)}$$
又$\mathbf{x}(t_{i+1})$的二階泰勒展開爲
$$\mathbf{x}(t_{i+1})=\mathbf{x}(t_i+\Delta t)=\mathbf{x}(t_i)+\Delta t{\mathbf{x}}^{\prime }(t_i)+\frac{\Delta t^2}{2!}{\mathbf{x}}^{\prime \prime }(t_i)+\mathbf{O}(\Delta t^3),\text{\hspace{1em}(5)}$$
比較公式$(4)$和公式$(5)$，若是要使公式$(3)$的局部截斷偏差知足$x_{i+1}-\mathbf{x}(t_{i+1})=\mathbf{O}(\Delta t^3)$，即要求公式具備二階精度，須要$\mathbf{O}(\Delta t^3)$前面的項相等，即須要知足
$$\{\begin{array}{l}{\lambda }_1+{\lambda }_2=1\\ {\lambda }_{2}p=\frac{1}{2}\end{array},\text{\hspace{1em}(6)}$$
知足上訴條件的公式$(3)$統稱爲二階龍格-庫塔公式。

特殊二階

若是當$p=1,{\lambda }_1=\frac{1}{2},{\lambda }_2=\frac{1}{2}$，二階龍格-庫塔公式就成了改進歐拉公式，具體詳解見[常微分方程的數值解法系列三] 改進歐拉法(預估校訂法)。簡單來講改進歐拉公式就是以$\mathbf{x}$函數在$t_i$點和$t_{i+1}$點處的斜率$k_1$和$k_2$的算數平均值做爲$\mathbf{x}$在$[t_i,t_{i+1}]$上的平均斜率$\bar{k}$。
若是當$p=\frac{1}{2},{\lambda }_1=0,{\lambda }_2=1$，二階龍格-庫塔公式就成了變形歐拉公式，具體詳解見[常微分方程的數值解法系列四] 中值法。簡單來講變形歐拉公式或者說是中值法就是以$\mathbf{x}$函數在$t_i$點和$t_{i+1}$中點處的斜率$k$做爲$\mathbf{x}$在$[t_i,t_{i+1}]$上的平均斜率$\bar{k}$。

三階

在$[t_i,t_{i+1}]$上取三點$t_i,t_{i+p},t_{i+q}(0<p\le q\le 1)$的斜率值分別爲$k_1,k_2,k_3$，把它們的線性組合${\lambda }_1{k}_1+{\lambda }_2{k}_2+{\lambda }_3{k}_3$做爲$\bar{k}$的近似值，因此根據公式$(1)$可得
$$x_{i+1}=x_i+\Delta t({\lambda }_1{k}_1+{\lambda }_2{k}_2+{\lambda }_3{k}_3)$$
其中
$$\{\begin{array}{l}{k}_1=f(t_i,x_i)\\ k_2=f(t_i+p\Delta t,x_i+p\Delta t{k}_1)\end{array}$$
爲了求$x_{i+q}$點的斜率$k_3$，最直接方法就是利用改進歐拉法進行預報得
$$k_3=f(t_{i+q},{\hat{x}}_{i+q})=f(t_i+q\Delta t,x_i+q\Delta t{k}_1)$$
可是這樣作效率比較差，由於在區間$[t_i,t_{i+q}]$區間內已經有$k_1,k_2$兩個斜率可使用了，因此能夠把$k_1,k_2$的線性組合當作$[x_i,x_{i+q}]$上平均斜率的近似值，這確定比上面求${\hat{x}}_{i+q}$精度要好。由此，可得
$${\hat{x}}_{i+q}=x_i+q\Delta t({\mu }_1{k}_1+{\mu }_2{k}_2)$$
因此
$$k_3=f(t_{i+q},{\hat{x}}_{i+q})=f(t_i+q\Delta t,x_i+q\Delta t({\mu }_1{k}_1+{\mu }_2{k}_2))$$
因此整理得
$$\{\begin{array}{l}{x}_{i+1}=x_i+\Delta t({\lambda }_1{k}_1+{\lambda }_2{k}_2+{\lambda }_3{k}_3)\\ k_1=f(t_i,x_i)\\ k_2=f(t_i+p\Delta t,x_i+p\Delta t{k}_1)\\ k_3=f(t_i+q\Delta t,x_i+q\Delta t({\mu }_1{k}_1+{\mu }_2{k}_2))\end{array},\text{\hspace{1em}(7)}$$
其中公式$(7)$有七個待確認參數${\lambda }_1,{\lambda }_2,{\lambda }_3,{\mu }_1,{\mu }_2,p,q$，相似前面的推導利用泰勒展開選擇合適的參數使上訴公式具備三階精度，即局部截斷偏差爲$\mathbf{O}(\Delta t^4)$，知足這類條件的公式$(7)$就是三階龍格-庫塔公式。
推導過程和二階相似，再此省略，只給出比較經常使用的一致形式
$$\{\begin{array}{l}{x}_{i+1}=x_i+\frac{1}{6}\Delta t(k_1+4k_2+k_3)\\ k_1=f(t_i,x_i)\\ k_2=f(t_i+\frac{1}{2}\Delta t,x_i+\frac{1}{2}\Delta t\cdot k_1)\\ k_3=f(t_i+\Delta t,x_i-\Delta t\cdot k_1+2\Delta t\cdot k_2))\end{array},\text{\hspace{1em}(7)}$$

高階

爲了進一步提升精度，在$[t_i,t_{i+1}]$上能夠選取更多點的，利用改進歐拉法中預報法，預報這些點的斜率值(以前在三階也討論過，爲了預報更準確，利用前面全部點的斜率來預報後面點的估計值進一步求得斜率，而不只僅只利用第一個點)，而後對這些斜率值加權平均做爲做爲$\mathbf{x}$在$[t_i,t_{i+1}]$上的平均斜率$\bar{k}$。而後在利用泰勒展開，對比相應的係數，從而肯定在知足特定$n$階精度下全部未知參數須要知足的條件。
前面也提到過，不考慮計算量的狀況下，理論上是能夠構造任意高階的龍格-庫塔公式。但在實際中發現，高於四階的龍格-庫塔公式，不但計算量很是大，並且精度進一步提高的有限。因此在實際使用中，四階的龍格-庫塔公式是精度和計算量都比較理想的公式。
推導過程和二階相似，再此省略，只給出最經典的四階的龍格-庫塔公式：
$$\{\begin{array}{l}{x}_{i+1}=x_i+\Delta t(\frac{1}{6}{k}_1+\frac{1}{3}{k}_2+\frac{1}{3}{k}_3+\frac{1}{6}{k}_4)=\frac{1}{6}\Delta t(k_1+2k_2+2k_3+k_4)\\ k_1=f(t_i,x_i)\\ k_2=f(t_i+\frac{1}{2}\Delta t,x_i+\frac{1}{2}\Delta t\cdot k_1)\\ k_3=f(t_i+\frac{1}{2}\Delta t,x_i+\frac{1}{2}\Delta t\cdot k_2)\\ k_4=f(t_i+\Delta t,x_i+\Delta t\cdot k_3)\end{array},\text{\hspace{1em}(8)}$$

步長選擇

以前討論的全部的龍格-庫塔方法都是以$\Delta t$定步長來展開的，但從$x_i\Rightarrow x_{i+1}$單步遞推過程來講，步長$\Delta t$越小，局部截斷偏差越小(方法肯定狀況下)，可是隨着步長的縮小，不但會引發計算量的增長，並且也有可能引發舍入偏差的嚴重積累；但步長$\Delta t$太大又不能達到預期的精度要求，因此選擇合適的步長$\Delta t$，在實際計算中也是比較重要的。其實有時候在實際使用中步長並不須要算法肯定，而是須要根據數據幀率來肯定的，好比imu數據。
下面給出求解步長的步驟：

1. 以步長$\Delta t$開始，利用龍格-庫塔公式計算$x_i\Rightarrow x_{i+1}$獲得一個近似值$x_{i+1}^{\Delta t}$；
2. 而後步長減半爲$\Delta t/2$，利用龍格-庫塔公式分兩步計算$x_i\Rightarrow x_{i+\frac{1}{2}}\Rightarrow x_{i+1}$獲得一個近似值$x_{i+1}^{\Delta t/2}$；
3. 計算$|x_{i+1}^{\Delta t/2}-x_{i+1}^{\Delta t}<ϵ|$是否成立，若是成立直接步長選擇$\Delta t$，不然繼續步長減半$\Delta t/4$，重複上訴步驟直到知足精度要求。

例子

待續~~

本文同步分享在 博客「無比機智的永哥」（CSDN）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。