 |
RMVL
1.4.0-dev
Robotic Manipulation and Vision Library
|
载入中...
搜索中...
未找到
|
RMVL
1.4.0-dev
Robotic Manipulation and Vision Library
|
涉及常用的 Euler 公式 与 Runge-Kutta 算法
上一篇教程：非线性方程（组）数值解与迭代法 ↑
下一篇教程：基于 TOPSIS 模型的熵权法 ↓
\[ \def\transparent#1{\color{transparent}{#1}} \]
一个质量-弹簧-阻尼系统的运动微分方程可以表示为
\[m\frac{\mathrm d^2x}{\mathrm dt^2}+c\frac{\mathrm dx}{\mathrm dt}+kx=p(t)\tag{1-1a}\]
即
\[m\ddot x+c\dot x+kx=p(t)\tag{1-1b}\]
可写成一阶方程组的形式，令 \(x_1=x,\quad x_2=\dot x\)，有
\[\def\mat#1#2{\begin{bmatrix}#1\\#2\end{bmatrix}} \left\{\begin{align}\dot x_1&=x_2\\ \dot x_2&=-\frac kmx_1-\frac cmx_2+\frac1mp(t) \end{align}\right.\tag{1-2}\]
记 \(\dot{\pmb x}=\begin{bmatrix}\dot x_1\\\dot x_2\end{bmatrix},\quad A=\begin{bmatrix}0&1\\-\frac km&-\frac cm\end{bmatrix},\quad \pmb b(t)=\begin{bmatrix}0\\\frac1mp(t)\end{bmatrix},\quad \pmb x^{(0)}=\begin{bmatrix}x_1(t_0)\\x_2(t_0)\end{bmatrix}\)，有
\[\dot{\pmb x}=A\pmb x+\pmb b(t),\quad\pmb x(t_0)=\pmb x^{(0)}\tag{1-3}\]
\(\text{(1-3)}\)在控制系统中能够经常遇见，这种表示一个在时刻 \(t_0\)带有初始条件的 2 阶线性系统，对于一般的（非线性）方程组，我们可以表示为
\[\def\rkf#1{\dot x_{#1}=f_{#1}(t,x_1,x_2,\cdots,x_k),\quad x_{#1}(t_0)=x_{#1}^{(0)}}\rkf1\\\rkf2\\\vdots\\\rkf k\]
即
\[\dot{\pmb x}=\pmb F(t,\pmb x),\quad\pmb x(t_0)=\pmb x^{(0)}\tag{1-4a}\]
仅有 1 条的常微分方程则表示为
\[\dot x=f(t,x),\quad x(t_0)=x^{(0)}\tag{1-4b}\]
这就是常微分方程的初值问题，即给定一个一阶常微分方程和对应的初始条件，求解出指定点或者指定范围的函数值。对于 \(\text{(1-4b)}\)形式的初值问题，我们通常截取等步长，如取
\[\begin{align} &t_0<t_1<\cdots<t_n<t_{n+1}<\cdots\\&h=t_{n+1}-t_n\quad或\quad t_{n+1}=t_n+h\quad(h=0,1,\cdots) \end{align}\]
先考虑一阶常微分方程 \(\text{(1-4b)}\)的形式，在取定的 \(t_n\)处 Taylor 展开，有
\[\begin{align} x_{n+1}\approx x(t_{n+1})&=x(t_n)+x'(t_n)h+x''(t_n)\frac{h^2}2+\cdots\\&=x_n+h\dot x_n+\frac{h^2}2\ddot x_n+\cdots\\ 忽略h^2以上的高阶项，有\quad&\approx x_n+hf(t_n,x_n) \end{align}\tag{2-1}\]
可以得到
\[x_{n+1}=x_n+hf(t_n,x_n)\tag{2-2}\]
这便是求解常微分方程的基本方法，公式 \(\text{(2-2)}\)被称为 Explicit Euler 显式欧拉公式。
① Implicit Euler 隐式欧拉公式
\[x_{n+1}=x_n+hf(t_{n+1},x_{n+1})\tag{2-3a}\]
② 梯形公式
\[x_{n+1}=x_n+\frac h2\left[f(t_n,x_n)+f(t_{n+1},x_{n+1})\right]\tag{2-3b}\]
③ 改进 Euler 公式
\[\left\{\begin{align}&\tilde x_{n+1}=x_n+hf(t_n,x_n)\\ &x_{n+1}=x_n+\frac h2\left[f(t_n,x_n)+f(t_{n+1},\tilde x_{n+1})\right] \end{align}\right.\tag{2-3c}\]
同样的，对于一阶方程组，公式 \(\text{(2-2)}\)可以改写成
\[\pmb x_{n+1}=\pmb x_n+h\pmb F(t_n,\pmb x_n)\tag{2-4}\]
示例
用 Euler 方法和步长 \(h=0.1\)求解常微分方程的初值问题
\[\left\{\begin{align}y'&=x^3+y^3+1\quad(0\leq x\leq0.8)\\y(0)&=0\end{align}\right.\]
计算结果保留小数点后 6 位
解答
这里， \(f(x,y)=x^3+y^3+1,\ x_n=nh+x_0=0.1n+0=0.1n\quad(n=0,1,\cdots,8),\ y_0=0\)。由 Euler 公式计算可得
\[\def\hf#1#2{y_{#1}+h(x_{#1}^3+y_{#1}^3+1)=#2} y(0.1)\approx y_1=\hf0{0.100000}\\ y(0.2)\approx y_2=\hf1{0.200200}\\ y(0.3)\approx y_3=\hf2{0.301802}\\ y(0.4)\approx y_4=\hf3{0.407249}\\ y(0.5)\approx y_5=\hf4{0.520403}\\ y(0.6)\approx y_6=\hf5{0.646995}\\ y(0.7)\approx y_7=\hf6{0.795680}\\ y(0.8)\approx y_8=\hf7{0.980355}\]
上文的 4 种 Euler 单步法公式，使用哪种精度更高？这里要引入局部阶段误差的概念，因为每一次求解 \(x_i\)都会引入误差，并且误差会进行累积，局部截断误差不考虑迭代求解 \(x_n\)及之前的累积误差，仅考虑从 \(x_n\)到 \(x_{n+1}\)产生的误差，即认为 \(x_n=x(t_n)\)。对于显式 Euler 单步法，可以计算出其局部截断误差 \(T_{n+1}\)。
\[\begin{align} T_{n+1}&=x(t_{n+1})-x_{n+1}\\ &=x(t_{n+1})-x_n-hf(t_n,x_n)\\ &=x(t_{n+1})-x(t_n)-hx'(t_n)\\ x(t_{n+1})在t_n处\text{ Taylor }展开，有\quad&=x(t_n)+hx'(t_n)+\frac{h^2}2x''(t_n)+o(h^3)-x(t_n)-hx'(t_n)\\ &=\frac12x''(t_n)h^2+o(h^3) \end{align}\tag{2-5}\]
可以看出局部截断误差的主项为 \(0.5x''(t_n)h^2\)，我们称显式 Euler 单步法具有 1 阶精度。
再来考虑梯形公式 \(\text{(2-3b)}\)的局部截断误差 \(T_{n+1}\)。
\[\begin{align} T_{n+1}&=x(t_{n+1})-x(t_n)-\frac12hf(t_n,x(t_n))-\frac12hf(t_{n+1},x(t_{n+1}))\\ &=x(t_{n+1})-x(t_n)-\frac12hx'(t_n)-\frac12hx'(t_{n+1})\\ &=\left[x(t_n)+hx'(t_n)+\frac{h^2}2x''(t_n)+\frac{h^3}6x'''(t_n)+o(h^4)\right]-x(t_n)-\\ &\qquad\frac12hx'(t_n)-\frac12h\left[x'(t_n)+hx''(t_n)+\frac{h^2}2x'''(t_n)+o(h^3)\right]\\ &=-\frac1{12}x'''(t_n)h^3+o(h^4) \end{align}\tag{2-6}\]
局部截断误差的主项为 \(-\frac1{12}x'''(t_n)h^3\)，我们称梯形公式具有 2 阶精度。
对于 \(\text{(2-3c)}\)的改进 Euler 公式，可以改写成
\[\left\{\begin{align} x_{n+1}&=x_n+h\left(\frac12k_1+\frac12k_2\right)\\ k_1&=f(t_n,x_n)\\k_2&=f(t_n+h,x_n+hk_1) \end{align}\right.\tag{3-1}\]
模仿这一写法，我们继续构造新的 2 阶公式
\[\begin{align}\frac{\mathrm df}{\mathrm dx}&=\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}·\frac{\mathrm dy}{\mathrm dx}\\ f'&=f_x+f_yy'\end{align}\tag{i}\]
以及多元函数 \(f(x,y)\)的 Taylor 展开，令 \(\pmb x=(x-x_0,\ y-y_0)^T\)，则多元函数 Taylor 展开如下\[\begin{align}f(x,y)&=f(x_0,y_0)+\begin{bmatrix}f_x(x_0,y_0)&f_y(x_0,y_0)\end{bmatrix}\pmb x+\\ &\transparent=\frac1{2!}\pmb x^T\begin{bmatrix}f_{xx}(x_0,y_0)&f_{xy}(x_0,y_0)\\f_{yx}(x_0,y_0)&f_{yy}(x_0,y_0)\end{bmatrix}\pmb x+o^n\end{align}\tag{ii}\]
若仅想了解最终结果，请跳过此小节\[\left\{\begin{align} x_{n+1}&=x_n+h(\lambda_1k_1+\lambda_2k_2)\\ k_1&=f(t_n,x_n)\\k_2&=f(t_n+ph,x_n+phk_1) \end{align}\right.\tag{3-2}\]
其中 \(\lambda_1,\lambda_2,p\)为待定参数，先分别将 \(x(t_{n+1})\)和 \(x_{n+1}\)在 \(t_n\)处 Taylor 展开。
对 \(x(t_{n+1})\)，有
\[\begin{align} x(t_{n+1})&=x(t_n)+hx'(t_n)+\frac{h^2}2x''(t_n)+o(h^3)\\ &=x(t_n)+hf(t_n,x(t_n))+\frac{h^2}2\frac{\mathrm d}{\mathrm dt}f(t_n,x(t_n))+o(h^3)\\ &=x(t_n)+hf(t_n,x(t_n))+\\ &\transparent=\frac{h^2}2\left[\frac{\partial f(t_n,x(t_n))}{\partial t}+\frac{\partial f(t_n,x(t_n))}{\partial x}·\frac{\mathrm dx}{\mathrm dt}\right]+o(h^3)\\ 令f(t_n,x(t_n))=(f)_{(n)}\quad&=x(t_n)+h(f)_{(n)}+\frac{h^2}2\left[\frac{\partial (f)_{(n)}}{\partial t}+\frac{\partial (f)_{(n)}}{\partial x}\frac{\mathrm dx}{\mathrm dt}\right]+o(h^3)\\ 由f(t_n,x(t_n))=\frac{\mathrm dx}{\mathrm dt}\quad&=x(t_n)+h(f)_{(n)}+\frac{h^2}2\left[\frac{\partial (f)_{(n)}}{\partial t}+\frac{\partial (f)_{(n)}}{\partial x}(f)_{(n)}\right]+o(h^3)\\ &=x(t_n)+h(f)_{(n)}+\frac{h^2}2(f_t+f_xf)_{(n)}+o(h^3)\tag{3-3} \end{align}\]
对 \(x_{n+1}\)，有
\[\begin{align}x_{n+1}&=x_n+h\left[\lambda_1f(t_n,x_n)+\lambda_2f(t_n+ph,x_n+phf(t_n,x_n))\right]\\ &=x_n+h\lambda_1f(t_n,x_n)+h\lambda_2\left[f(t_n,x_n)+ph\frac{\partial f(t_n,x_n)}{\partial t}+\right.\\ &\transparent=\left.phf(t_n,x_n)\frac{\partial f(t_n,x_n)}{\partial x}+o(h^2)\right]\\ 令f(t_n,x_n)=(f)_n\quad&=x_n+h\lambda_1(f)_n+\\ &\transparent=h\lambda_2\left[(f)_n+ph\frac{\partial(f)_n}{\partial t}+ph(f)_n\frac{\partial(f)_n}{\partial x}+o(h^2)\right]\\ &=x_n+h\lambda_1(f)_n+h\lambda_2\left[(f)_n+ph(f_t)_n+ph(f)_n(f_x)_n+o(h^2)\right]\\ &=x_n+h(\lambda_1+\lambda_2)(f)_n+ph^2\lambda_2(f_t+f_xf)_n+o(h^3) \end{align}\tag{3-4}\]
为求局部截断误差，则满足 \(x(t_n)=x_n,\ (f)_{(n)}=(f)_n,\ (f_t+f_xf)_{(n)}=(f_t+f_xf)_n\)，因此 \(T_{n+1}=\text{(3-3)}-\text{(3-4)}\)，即
\[\begin{align} T_{n+1}&=x(t_{n+1})-x_{n+1}\\ &=(1-\lambda_1-\lambda_2)(f)_n+\left(\frac12-p\lambda_2\right)h^2(f_t+f_xf)_n+o(h^3) \end{align}\tag{3-5}\]
我们希望这个公式具有 2 阶精度，即局部截断误差的主项为 0，因此有
\[\left\{\begin{align} &\lambda_1+\lambda_2=1\\ &p\lambda_2=\frac12 \end{align}\right.\tag{3-6}\]
\(\text{(3-2)}\)和 \(\text{(3-6)}\)联立构成 2 阶精度的单步显式公式族，由于使用到了 2 个斜率值 \(k_1\)和 \(k_2\)，因此称为 2 级 2 阶 Runge-Kutta 公式族，最常用的是当 \(\lambda_1=0\)即 \(\lambda_2=1,\quad p=\frac12\)，这时得到所谓 2 级 2 阶中点公式。
\[\begin{align}x_{n+1}&=x_n+hk_2\\k_1&=f(t_n,x_n)\\ k_2&=f\left(t_n+\frac h2,x_n+\frac h2k_1\right)\end{align}\tag{3-7}\]
为了便于记录上文的 Runge-Kutta 中点公式，我们将公式 \(\text{(3-2)}\)改写为
\[\left\{\begin{align}x_{n+1}&=x_n+h(\lambda_1k_1+\lambda_2k_2)\\k_1&=f(t_n+p_1h,x_n+h(a_{11} k_1+a_{12}k_2))\\k_2&=f(t_n+p_2h,x_n+h(a_{21}k_1+a_{22}k_2))\end{align}\right.\tag{3-8}\]
令 \(\pmb p=\mat{p_1}{p_2},\quad\pmb\lambda=(\lambda_1,\lambda_2),\quad R=\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}\)，则
\[\begin{array}{c|c}\pmb p&R\\\hline&\pmb\lambda\end{array}\Rightarrow\begin{array}{c|cc}p_1& a_{11}&a_{12}\\p_2&a_{21}&a_{22}\\\hline&\lambda_1&\lambda_2\end{array}\tag{3-9}\]
被称为 Butcher 表，例如，上文的中点公式可以表示为
\[\begin{array}{c|cc}0&0&0\\1&1&0\\\hline&\frac12&\frac12\end{array}\tag{3-10}\]
一般的，对于一阶方程 \(x'=f(t,x),\ x(t_0)=x^{(0)}\)，有以下 \(n\)阶公式
\[\left\{\begin{align} x_{n+1}&=x_n+h(\lambda_1k_1+\lambda_2k_2+\cdots+\lambda_nk_n)\\ k_1&=f(t_n+p_1h,x_n+h(a_{11}k_1+a_{12}k_2+\cdots+a_{1n}k_n))\\ k_2&=f(t_n+p_2h,x_n+h(a_{21}k_1+a_{22}k_2+\cdots+a_{2n}k_n))\\&\vdots\\ k_n&=f(t_n+p_nh,x_n+h(a_{n1}k_1+a_{n2}k_2+\cdots+a_{nn}k_n))\\ \end{align}\right.\tag{3-11}\]
同样可以用 Butcher 表来表示：
\[\begin{array}{c|cccc} p_1&a_{11}&a_{12}&\cdots&a_{1n}\\ p_2&a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ p_n&a_{n1}&a_{n2}&\cdots&a_{nn}\\ \hline&\lambda_1&\lambda_2&\cdots&\lambda_n\end{array}\tag{3-12}\]
这里直接给出对应的 Butcher 表。
一种 3 级 3 阶 Runge-Kutta 方法的具体形式对应的 Butcher 表
\[\begin{array}{c|ccc}0&0&0&0\\\frac12&\frac12&0&0\\1&-1&2&0\\ \hline&\frac16&\frac23&\frac16\end{array}\tag{3-13a}\]
经典 4 级 4 阶 Runge-Kutta 方法对应的 Butcher 表
\[\begin{array}{c|cccc}0&0&0&0&0\\\frac12&\frac12&0&0&0\\\frac12&0&\frac12&0& 0\\1&0&0&1&0\\\hline&\frac16&\frac13&\frac13&\frac16\end{array}\tag{3-13b}\]
对于一阶方程组 \(\pmb x'=\pmb F(t,\pmb x),\ \pmb x(t_0)=\pmb x^{(0)}\)，公式 \(\text{(3-11)}\)可以改写为
\[\left\{\begin{align} \pmb x_{n+1}&=\pmb x_n+h(\lambda_1\pmb k_1+\lambda_2\pmb k_2+\cdots+\lambda_n\pmb k_n)\\ \pmb k_1&=\pmb F(t_n+p_1h,\pmb x_n+h(a_{11}\pmb k_1+a_{12}\pmb k_2+\cdots+a_{1n}\pmb k_n))\\ \pmb k_2&=\pmb F(t_n+p_2h,\pmb x_n+h(a_{21}\pmb k_1+a_{22}\pmb k_2+\cdots+a_{2n}\pmb k_n))\\&\vdots\\ \pmb k_n&=\pmb F(t_n+p_nh,\pmb x_n+h(a_{n1}\pmb k_1+a_{n2}\pmb k_2+\cdots+a_{nn}\pmb k_n))\\ \end{align}\right.\tag{3-14}\]
公式 \(\text{(3-14)}\)与 \(\text{(3-11)}\)基本一致，因此同样可以使用 Butcher 表来描述常微分方程组的 Runge-Kutta 公式。
RMVL 的相关类请参考 rm::RungeKutta
示例
使用 2 阶中点公式求解 \([0,1]\)上常微分方程的初值问题，并在 \(t=1\)处与实际值进行比较
\[\left\{\begin{align}\dot x_1&=2x_2+t\\\dot x_2&=-x_1-3x_2\end{align}\right.\\ 且满足x_1(0)=1,\quad x_2(0)=-1\]
解答
① 数值解
使用中点公式：
\[\begin{align}\pmb x_{n+1}&=\pmb x_n+h\pmb k_2\\\pmb k_1&=\pmb F(t_n,\pmb x_n)\\ \pmb k_2&=\pmb F(t_n+\frac h2,\pmb x_n+\frac h2\pmb k_1)\end{align}\]
即
\[\begin{align}x_{1_{n+1}}&=x_{1_n}+hk_2\\x_{2_{n+1}}&=x_{2_n}+hl_2\\ k_1&=f_1(t_n,x_{1_n},x_{2_n})\\l_1&=f_2(t_n,x_{1_n},x_{2_n})\\ k_2&=f_1(t_n+\frac h2,x_{1_n}+\frac h2k_1,x_{2_n}+\frac h2l_1)\\ l_2&=f_2(t_n+\frac h2,x_{1_n}+\frac h2k_1,x_{2_n}+\frac h2l_1)\end{align}\]
迭代运行结果如下：
因此，在 \(t=1\)时， \(x\approx(0.587286,-0.219401)^T\)
② 精确解
对待求的初值问题写成矩阵形式
\[\dot X=\begin{bmatrix}0&2\\-1&-3\end{bmatrix}X+\mat10t,\quad X(0)=\mat1{-1}\]
由矩阵分析的知识，可以知道，对于形如
\[\dot X=AX+B(t)\tag{a}\]
的矩阵微分方程，其解为
\[X=e^{At}X(0)+\int_0^te^{A(t-\tau)}B(\tau)\mathrm d\tau\tag{b}\]
其中 \(e^{At}\)为矩阵函数，可使用
的方式进行求解，这里使用最小多项式的 Hamilton-Cayley 定理进行求解。
\(A\) 的特征多项式为 \(\varphi(\lambda)=\lambda^2+3\lambda+2=(\lambda+1)(\lambda+2)\)，无重根，因此 \(A\) 的最小多项式为 \(\mu(\lambda)=(\lambda+1)(\lambda+2)\)，由 Hamilton-Cayley 定理可知 \(A^2+3A+2I=0\)，对 \(e^{At}\)进行 Taylor 展开，并使用带余除法，一定可以得到
\[e^{At}=a_0I+a_1A\tag{c}\]
将 \(\lambda=-1,-2\)代入 \(e^{At}\)，可得
\[\left\{\begin{align}a_0+a_1(-1)&=e^{-t}\\a_0+a_1(-2)&=e^{-2t}\end{align}\right.\tag{d}\]
解得 \(\left\{\begin{align}a_0&=2e^{-t}-e^{-2t}\\a_1&=e^{-t}-e^{-2t}\end{align}\right.\)，因此
\[\begin{align}e^{At}&=(2e^{-t}-e^{-2t})I+(e^{-t}-e^{-2t})A\\ &=\begin{bmatrix}2&2\\1&1\end{bmatrix}e^{-t}+\begin{bmatrix}-1&-2\\1&2\end{bmatrix}e^{-2t} \end{align}\tag{e}\]
可以算出
\[e^{At}X(0)=\mat00e^{-t}+\mat1{-1}e^{-2t}=\mat1{-1}e^{-2t}\tag{f}\]
同理有
\[\int_0^te^{A(t-\tau)}B(\tau)\mathrm d\tau=\mat2{-1}\left(e^{-t}+t-1\right)+ \mat{-1}1\left(\frac14e^{-2t}+\frac12t-\frac14\right)\tag{g}\]
整理得
\[\begin{align}X&=e^{At}X(0)+\int_0^te^{A(t-\tau)}B(\tau)\mathrm d\tau\\ &=\mat{0.75}{-0.75}e^{-2t}+\mat2{-1}e^{-t}+\mat{1.5}{-0.5}t+\mat{-1.75}{0.75}\end{align}\tag{h}\]
把 \(t=1\) 代入上式，可得 \(X=\mat{0.75e^{-2}+2e^{-1}-0.25}{-0.75e^{-2}-e^{-1}+0.25}\approx\mat{0.587260}{-0.219381}\)，数值解 \(X=\mat{0.587286}{-0.219401}\)与之相符。
代码
1.12.0