在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程。但之前求解微分方程的解析方法主要是应用于一些简单和特殊的微分方程求解中,对于一般形式的微分方程,一般很难用解析方法求出精确解,只能用数值方法求解。该系列主要介绍一些常用的常微分方程的数值解法,主要包括:
- [常微分方程的数值解法系列一] 常微分方程
- [常微分方程的数值解法系列二] 欧拉法
- [常微分方程的数值解法系列三] 改进欧拉法(预估校正法)
- [常微分方程的数值解法系列四] 中值法
- [常微分方程的数值解法系列五] 龙格-库塔(RK4)法
这个系列后面文章会用到前面文章的理论和技术,所以建议按照顺序查看。
简介
而在具体求解微分方程中,一般来说给定的条件是
{
x
′
(
t
)
=
f
(
t
,
x
(
t
)
)
,
a
≤
t
≤
b
x
(
t
0
)
=
x
0
\begin{cases} \mathbf{x}^{\prime}(t)=f(t, \mathbf{x}(t)), \quad a \leq t \leq b \\ \mathbf{x}({t_0}) = x_0 \end{cases}
{x′(t)=f(t,x(t)),a≤t≤bx(t0)=x0
即给定微分方程以及原方程在初始点的值,求原方程在某个
t
t
t下的
x
(
t
)
x(t)
x(t)原方程的值,这类问题就是(一阶)常微分方程初值求解问题。
(这里不对常微分方程或者偏微分方程等概念,以及求解微分方程的其他条件如边界条件情况做详细介绍,需要了解的话可以自己google,不影响本系列介绍。)
解法介绍
一般来说,(一阶)常微分方程数值解法基本思想是:
在区间
[
a
,
b
]
[a, b]
[a,b]中插入一系列间隔相同为
Δ
t
\Delta t
Δt的离散点
a
≤
t
0
<
t
1
<
⋯
<
t
i
<
⋯
<
t
n
−
1
<
t
n
≤
b
a \leq t_0 < t_1< \cdots < t_i < \cdots < t_{n-1} < t_n \leq b
a≤t0<t1<⋯<ti<⋯<tn−1<tn≤b
常微分方程的数值解法的目的就是在给定合适初始值前提下,建立求解
x
(
t
)
x(t)
x(t)的近似值
x
t
x_t
xt的递推方式,这样求得
x
(
t
)
x(t)
x(t)在各个离散点的近似值。
给定一个具体的时间间隔
Δ
t
\Delta_t
Δt,可以把该方程转换成差分方程
x
(
t
+
Δ
t
)
=
x
(
t
)
+
∫
τ
∈
[
t
,
t
+
Δ
t
]
f
(
τ
,
x
(
τ
)
)
d
τ
,
a
≤
t
≤
b
,
(1)
\mathbf{x}(t+\Delta t)=\mathbf{x}(t)+\int_{\tau \in [t, t+\Delta t]} f(\tau, \mathbf{x}(\tau)) d \tau, \quad a \leq t \leq b, \tag{1}
x(t+Δt)=x(t)+∫τ∈[t,t+Δt]f(τ,x(τ))dτ,a≤t≤b,(1)
这样利用之前说的离散点可以把
t
t
t离散化,那么
t
i
=
a
+
i
Δ
t
t_{i}=a+i \Delta t
ti=a+iΔt以及
x
i
≈
x
(
t
i
)
\mathbf{x}_{i} \approx \mathbf{x}(t_{i})
xi≈x(ti),所以上面公式
(
2
)
(2)
(2)也可以写成
x
(
t
i
+
1
)
≈
x
i
+
1
=
x
i
+
∫
τ
∈
[
a
+
i
Δ
t
,
a
+
(
i
+
1
)
Δ
t
]
f
(
τ
,
x
(
τ
)
)
d
τ
,
(2)
\mathbf{x}(t_{i+1}) \approx x_{i+1}=\mathbf{x}_{i}+\int_{\tau \in [a+i \Delta t, a+(i+1) \Delta t]} f(\tau, \mathbf{x}(\tau)) d \tau, \tag{2}
x(ti+1)≈xi+1=xi+∫τ∈[a+iΔt,a+(i+1)Δt]f(τ,x(τ))dτ,(2)
截断误差
上节介绍常微分方程的数值解法就是利用初值和离散点获得近似值
x
t
x_t
xt,但是和真实值
x
(
t
)
x(t)
x(t)之前还是存在误差,即
e
i
=
x
i
−
x
(
t
i
)
,
(3)
\boldsymbol{e}_i = x_i - \mathbf{x}(t_i), \tag{3}
ei=xi−x(ti),(3)
e
i
\boldsymbol{e}_i
ei则表示该数值解法在离散点
t
i
t_i
ti处的局部截断误差。局部截断误差在一定程度上反应了该数值解法的精度。
一般来说常用泰勒展开式来计算讨论局部截断误差,后面具体方法会给出具体的对应的截断误差的讨论。这里简单给出定义:
如果数值解法的局部截断误差为
e
i
=
O
(
Δ
t
P
+
1
)
\boldsymbol{e}_i = \boldsymbol{O}(\Delta t^{P+1})
ei=O(ΔtP+1),则称该解法具有P阶精度或P阶解法。
这里 O ( Δ t P + 1 ) \boldsymbol{O}(\Delta t^{P+1}) O(ΔtP+1)为泰勒展开式的余项