线性动态系统建模

状态外插方程的推导

我们的目标是推导出如下形式的状态外插方程：

为达到这个目的，我们需要对动态系统进行建模。也就是说，思考出动态系统的 状态空间表达。下面两个方程就是LTI系统的状态空间表达：

为了找到 状态转移矩阵 $\boldsymbol{F}$ 和 输入转移矩阵 $\boldsymbol{G}$，我们需要求解状态空间微分方程。

下图总结了状态外插方程的推导过程。

状态空间表达

你也许会好奇为什么状态空间表达必须是下述的形式：

$$\boldsymbol{\dot{x}(t) = Ax(t) + Bu(t)}$$ $$\boldsymbol{y(t) = Cx(t) + Du(t)}$$

许多计算机软件包求解微分方程时需要这样的描述。

最好的讲述状态空间表达的方式是通过示例。

示例 - 匀速运动的物体

由于没有外部力作用在物体上，系统没有输入：

$$\boldsymbol{u(t)} = 0$$

状态空间向量 $\boldsymbol{x(t)}$ 是物体的位移量 $p(t)$ 和速度 $v(t)$.

$$\boldsymbol{x(t)} = \left[ \begin{matrix} p\\ v\\ \end{matrix} \right]$$

$$\boldsymbol{\dot{x}(t) = Ax(t) + Bu(t)}$$

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \boldsymbol{A} \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \boldsymbol{B} \times 0$$

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right]$$

$$\boldsymbol{A} = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right]$$

我们获取了微分形式的第一个方程。

动态系统输出 $\boldsymbol{y(t)}$ 是物体的位移 $\boldsymbol{p(t)}$，也是系统的状态变量。

$$\boldsymbol{y(t) = Cx(t) + Du(t)}$$

$$p = \boldsymbol{C} \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \boldsymbol{D} \times 0$$

$$p = \left[ \begin{matrix} 1 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right]$$

$$\boldsymbol{C} = \left[ \begin{matrix} 1 & 0\\ \end{matrix} \right]$$

到这里就完成了。

建模高阶动态系统

许多动态系统模型由高阶的微分方程描述。微分方程的阶次就是微分方程中具有最高阶导数的变量的阶次。

为了处理高阶方程，我们应该通过引入新变量并将其代入最高阶项的方式来将方程逐渐降为一阶微分方程。

微分方程降阶算法

一个 $n$ 阶线性微分方程可以由 $n$ 个一阶线性微分方程组来表示。

动态系统的高阶 支配方程 Governing Equation 有如下的形式：

$$a_{n}\frac{d^{n}y}{dt^{n}}+ a_{n-1}\frac{d^{n-1}y}{dt^{n-1}}+ \ldots + a_{2}\frac{d^{2}y}{dt^{2}}+a_{1}\frac{d^{1}y}{dt^{1}}+ a_{0}y = u$$

这个 支配方程 完全描述了系统的动态状态。

降阶：

1. 隔离出最高阶项：
$$\frac{d^{n}y}{dt^{n}}= -\frac{a_{0}}{a_{n}}y- \frac{a_{1}}{a_{n}}\frac{d^{1}y}{dt^{1}}- \frac{a_{2}}{a_{n}}\frac{d^{2}y}{dt^{2}}- \ldots - \frac{a_{n-1}}{a_{n}}\frac{d^{n-1}y}{dt^{n-1}}+\frac{1}{a_{n}}u$$

$y$ 及其前 $n - 1$ 阶导是这个系统的状态。

2. 定义新变量 $x_{1}, x_{2}, \ldots, x_{n}$:

令 $x_{i} = \frac{d^{i-1} y}{dt^{i-1}},$，有：

$$x_{1} \left( t \right) =y$$ $$x_{2} \left( t \right) = \frac{dy}{dt}$$ $$x_{3} \left( t \right) = \frac{d^{2}y}{dt^{2}}$$ $$\vdots$$ $$x_{n} \left( t \right) = \frac{d^{n-1}y}{dt^{n-1}}$$

此时 $x_{i}(t)$ 就是状态变量。

3. 这些变量对 $t$ 求导：
$$\frac{dx_{1}}{dt}=x_{2} \left( t \right)$$ $$\frac{dx_{2}}{dt}=x_{3} \left( t \right)$$ $$\vdots$$ $$\frac{dx_{n-1}}{dt}=x_{n} \left( t \right)$$ $$\frac{dx_{n}}{dt}=\frac{d^{n}y}{dt^{n}}$$
4. 把一开始隔离出的 $\frac{d^{n}y}{dt^{n}}$ 项代入最后一个方程：
$$\frac{dx_{1}}{dt}=x_{2} \left( t \right)$$ $$\frac{dx_{2}}{dt}=x_{3} \left( t \right)$$ $$\vdots$$ $$\frac{dx_{n-1}}{dt}=x_{n} \left( t \right)$$ $$\frac{dx_{n}}{dt}=-\frac{a_{0}}{a_{n}}x_{1}- \frac{a_{1}}{a_{n}}x_{2}- \frac{a_{2}}{a_{n}}x_{3}- \ldots - \frac{a_{n-1}}{a_{n}}x_{n}+\frac{1}{a_{n}}u$$
5. 用向量矩阵形式改写上述方程组：
$$\left[ \begin{matrix} \frac{dx_{1}}{dt}\\ \frac{dx_{2}}{dt}\\ \vdots \\ \frac{dx_{n-1}}{dt}\\ \frac{dx_{n}}{dt}\\ \end{matrix} \right] = \left[ \begin{matrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \vdots \\ \dot{x}_{n-1}\\ \dot{x}_{n}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 1\\ -\frac{a_{0}}{a_{n}} & -\frac{a_{1}}{a_{n}} & -\frac{a_{2}}{a_{n}} & \cdots & -\frac{a_{n-2}}{a_{n}} & -\frac{a_{n-1}}{a_{n}}\\ \end{matrix} \right] \left[ \begin{matrix} x_{1}\\ x_{2}\\ \vdots \\ x_{n-1}\\ x_{n}\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ 0\\ \vdots \\ 0\\ \frac{1}{a_{n}}\\ \end{matrix} \right] u$$

搞定了。这个方程就是符合下述形式的状态空间方程：

$$\boldsymbol{\dot{x}(t) = Ax(t) + Bu(t)}$$

其中：

$$\boldsymbol{A} = \left[ \begin{matrix} 0 & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 1\\ -\frac{a_{0}}{a_{n}} & -\frac{a_{1}}{a_{n}} & -\frac{a_{2}}{a_{n}} & \cdots & -\frac{a_{n-2}}{a_{n}} & -\frac{a_{n-1}}{a_{n}}\\ \end{matrix} \right]$$

$$\boldsymbol{B} = \left[ \begin{matrix} 0\\ 0\\ \vdots \\ 0\\ \frac{1}{a_{n}}\\ \end{matrix} \right]$$

再来看两个例子。

示例 - 匀加速运动物体

本例中，物体上有外力作用。

匀加速运动物体的支配方程就是牛顿第二定律：

$$m\ddot{p} = F$$

式中：

$p$ 是物体位移

$m$ 是物体质量

$F$ 是外部作用在物体上的力

我们对上述方程应用前面提到的降阶算法。

1. 隔离出最高阶项：
$$\ddot{p} = \frac{1}{m}F$$
2. 定义新变量 $x_{1}$ 和 $x_{2}$:
$$x_{1} = p$$ $$x_{2} = \dot{p}$$

$x_{1}$ 和 $x_{2}$ 就是状态变量。

3. 状态变量各自对时间求导：
$$\dot{x}_{1} = x_{2}$$ $$\dot{x}_{2} = \ddot{p}$$
4. 把最高阶项 $\ddot{p}$ 代入第二个方程：
$$\dot{x}_{1} = x_{2}$$ $$\dot{x}_{2} = \frac{F}{m}$$
5. 用矩阵形式描述上述方程组：
$$\left[ \begin{matrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} x_{1}\\ x_{2}\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \frac{1}{m}\\ \end{matrix} \right] F$$

记住 $x_{2} = \dot{p}$. 记 $x_{2}$ 为 $v$ 会更有意义，因为 $x_{2}$ 本身就是物体的速度。

上述方程可以重写为：

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \frac{1}{m}\\ \end{matrix} \right] F$$

或：

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ 1\\ \end{matrix} \right] a$$

其中 $a$ 是物体被施加作用力 $F$ 后产生的加速度。

于是我们得到了这个形式的方程：

$$\boldsymbol{\dot{x}(t) = Ax(t) + Bu(t)}$$

状态向量 $\boldsymbol{x(t)}$ 是物体的位移 $p(t)$ 和速度 $v(t)$.

动态系统的输出 $\boldsymbol{y(t)}$ 是物体的位移 $p(t)$，也是系统状态变量的分量。

$$\boldsymbol{y(t) = Cx(t) + Du(t)}$$

$$p = \boldsymbol{C} \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \boldsymbol{D} a$$

$$p = \left[ \begin{matrix} 1 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \end{matrix} \right] a$$

示例 - 质量-弹簧-阻尼系统

质量-弹簧-阻尼系统包含三个基本元件：

• 质量 - 惯性元件
• 弹簧 - 弹性元件
• 阻尼 - 摩擦元件

Each of the elements has one of two possible energy behaviors: 每个元素具有两种可能的能量转移行为之一：

• 将输入的能量存储在自身
• 将输入的能量转化为“摩擦”热量

质量将输入的能量以动能形式存储。

当弹簧被拉伸超过原始长度时，输入能量以势能形式存储。

阻尼以热能形式将输入能量耗散。

这三个元件：质量、弹簧和阻尼，能够以一般形式描述任何动态系统的响应。

本系统的受力分析图如下所示。

弹簧弹力正比于质量的位移量 $p$.

粘滞阻力正比于质量的速度 $\dot{p}$.

牛顿第二定律表明：

$$\sum _{}^{}F = ma = m\frac{d^{2}p}{dt^{2}} = m\ddot{p}$$

We proceed by summing the forces and applying Newton's second law: 我们将上述所有作用力加和在一起，并运用牛顿第二定律：

$$\sum_{}^{}F = F(t) - c\dot{p} - kp = m\ddot{p}$$

其中：

$p$ 是物体位移

$m$ 是物体质量

$F$ 是外部作用在物体上的力

$k$ 是弹簧弹力系数

$c$ 是阻尼系数

这个方程就是完全描述上述系统的动态状态的支配方程。

接下来我们对这个支配方程应用降阶算法。

1. 隔离出最高阶项：
$$\ddot{p} = -\frac{k}{m}p - \frac{c}{m}\dot{p} + \frac{1}{m}F$$
2. 定义新变量 $x_{1}$ 和 $x_{2}$:
$$x_{1} = p$$ $$x_{2} = \dot{p}$$

$x_{1}$ 和 $x_{2}$ 就是状态变量。

3. 状态变量各自对时间求导：
$$\dot{x}_{1} = x_{2}$$ $$\dot{x}_{2} = \ddot{p}$$
4. 把一开始隔离出的 $\ddot{p}$ 代入第二个方程：
$$\dot{x}_{1} = x_{2}$$ $$\dot{x}_{2} = -\frac{k}{m}p - \frac{c}{m}\dot{p} + \frac{1}{m}F$$
5. 用向量矩阵形式改写上述方程组：
$$\left[ \begin{matrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ -\frac{k}{m} & -\frac{c}{m}\\ \end{matrix} \right] \left[ \begin{matrix} x_{1}\\ x_{2}\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \frac{1}{m}\\ \end{matrix} \right] F$$

记住 $x_{2} = \dot{p}$. 记 $x_{2}$ 为 $v$ 会更有意义，因为 $x_{2}$ 本身就是物体的速度。

上述方程可以重写为：

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ -\frac{k}{m} & -\frac{c}{m}\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \frac{1}{m}\\ \end{matrix} \right] F$$

于是我们得到了这个形式的方程：

$$\boldsymbol{\dot{x}(t) = Ax(t) + Bu(t)}$$

状态向量 $\boldsymbol{x(t)}$ 是物体的位移 $p(t)$ 和速度 $v(t)$.

动态系统的输出 $\boldsymbol{y(t)}$ 是物体的位移 $p(t)$，也是系统状态变量的分量。

$$\boldsymbol{y(t) = Cx(t) + Du(t)}$$

$$p = \boldsymbol{C} \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \boldsymbol{D} F$$

$$p = \left[ \begin{matrix} 1 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \end{matrix} \right] F$$

更多示例

你可以在ShareTechnote找到更多有意义的状态空间建模示例。

求解微分方程

没有输入的动态系统

没有输入的LTI系统可以用一阶微分方程组描述：

$$\boldsymbol{\dot{x} = Ax}$$

其中 $\boldsymbol{A}$ 是 动态矩阵。

我们的目标是找到 状态转移矩阵 $\boldsymbol{F}$.

需要求解微分方程来获取 $\boldsymbol{F}$.

在一维情况下，微分方程具有如下形式：

$$\frac{dx}{dt} = kx$$ $$\frac{dx}{x} = kdt$$

两侧同时积分得到：

$$\int _{x_{0}}^{x_{1}}\frac{1}{x}dx = \int _{0}^{ \Delta t}kdt$$

求解积分：

$$ln \left( x_{1} \right) - ln \left( x_{0} \right) =k \Delta t$$ $$ln \left( x_{1} \right) = ln \left( x_{0} \right) +k \Delta t$$ $$x_{1}= e^{ln \left( x_{0} \right) +k \Delta t}$$ $$x_{1}= e^{ln \left( x_{0} \right) }e^{k \Delta t}$$ $$x_{1}= x_{0}e^{k \Delta t}$$

类似地，对多维情况：

$$\boldsymbol{\dot{x} = Ax}$$

其解为：

$$x_{n+1}= x_{n}e^{\boldsymbol{A} \Delta t}$$

于是我们就找到了状态转移矩阵 $\boldsymbol{F}$:

$$\boldsymbol{F} = e^{\boldsymbol{A} \Delta t}$$

$e^{\boldsymbol{A}t}$ 是 矩阵指数.

矩阵指数可以通过泰勒展开计算：

$$e^{\boldsymbol{X}}= \sum _{k=0}^{\infty}\frac{1}{k!}\boldsymbol{X}^{k}$$

因此：

$$\boldsymbol{F} = e^{\boldsymbol{A} \Delta t} = \boldsymbol{I} + \boldsymbol{A} \Delta t+ \frac{ \left( \boldsymbol{A} \Delta t \right) ^{2}}{2!}+ \frac{ \left( \boldsymbol{A} \Delta t \right) ^{3}}{3!}+ \frac{ \left( \boldsymbol{A} \Delta t \right) ^{4}}{4!}+ \ldots$$

其中 $\boldsymbol{I}$ 是单位矩阵。

示例 续 - 匀速运动的物体

现在我们可以得到匀速运动物体的状态转移矩阵 $\boldsymbol{F}$ 了。

下面的微分方程组可以描述匀速运动模型：

$$\begin{cases} \frac{dp}{dt}=v\\ \frac{dv}{dt}=0\\ \end{cases}$$

用矩阵形式表示：

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right]$$ $$\boldsymbol{A} = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right]$$

我们来计算 $\boldsymbol{A}^{2}$

$$\boldsymbol{A}^{2} = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 0\\ 0 & 0\\ \end{matrix} \right]$$

由于 $\boldsymbol{A}^{2} = 0$，$\boldsymbol{A}$ 的更高阶次幂也全部是0.

于是我们可以得到匀速运动物体的状态转移矩阵 $\boldsymbol{F}$:

$$F= e^{\boldsymbol{A} \Delta t} = \boldsymbol{I} + \boldsymbol{A} \Delta t = \left[ \begin{matrix} 1 & 0\\ 0 & 1\\ \end{matrix} \right] + \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \Delta t = \left[ \begin{matrix} 1 & \Delta t\\ 0 & 1\\ \end{matrix} \right]$$ $$\boldsymbol{x}_{n+1} = \boldsymbol{F}\boldsymbol{x}_{n} = \left[ \begin{matrix} 1 & \Delta t\\ 0 & 1\\ \end{matrix} \right] \boldsymbol{x}_{n}$$ $$\left[ \begin{matrix} x_{n+1}\\ \dot{x}_{n+1}\\ \end{matrix} \right] = \left[ \begin{matrix} 1 & \Delta t\\ 0 & 1\\ \end{matrix} \right] \left[ \begin{matrix} x_{n}\\ \dot{x}_{n}\\ \end{matrix} \right]$$

有输入的动态系统

假设零阶保持采样，即把输入看作一个分段常量，那么一个下述形式的状态空间微分方程：

$$\boldsymbol{\dot{x}(t) = Ax(t) + Bu(t)}$$

具有通解为：

$$\boldsymbol{x ( t+ \Delta t )} = \color{red}{\underbrace{e^{\boldsymbol{A} \Delta t}}_\textrm{F}} \boldsymbol{x(t)} + \color{blue}{\underbrace{\int _{0}^{\Delta t} e^{\boldsymbol{A}t}dt\boldsymbol{B}}_\textrm{G}} \boldsymbol{u(t)}$$

去掉输入的话 $(\boldsymbol{u(t)} = 0)$，我们就得到了前面例子里推导出的解。

这里的证明就不展开了 - 你可以在 Tom M. Apostol 的 “微积分” 教科书里的定理8.3处找到证明，也可以在任何其他一本微积分书上找到。

我们再来求解匀加速运动物体和质量-弹簧-阻尼系统示例。

示例 续 - 匀加速运动物体

回忆匀加速运动物体的状态空间方程：

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ 1\\ \end{matrix} \right] a$$

其中：

$$\boldsymbol{A} = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right]$$

$$\boldsymbol{B} = \left[ \begin{matrix} 0\\ 1\\ \end{matrix} \right]$$

我们来求解方程：

$$\boldsymbol{x ( t+ \Delta t )} = \color{red}{\underbrace{e^{\boldsymbol{A} \Delta t}}_\textrm{F}} \boldsymbol{x(t)} + \color{blue}{\underbrace{\int _{0}^{\Delta t} e^{\boldsymbol{A}t}dt\boldsymbol{B}}_\textrm{G}} \boldsymbol{u(t)}$$

解 F

$$\boldsymbol{F} = e^{\boldsymbol{A} \Delta t}$$

我们已经得到了在 $A = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right]$ 时的解：

$$\boldsymbol{F} = \left[ \begin{matrix} 1 & \Delta t\\ 0 & 1\\ \end{matrix} \right]$$

解 G

$$\boldsymbol{G} = \int _{0}^{ \Delta t}e^{\boldsymbol{A}t}dt \boldsymbol{B}$$

$\int _{0}^{ \Delta t}e^{\boldsymbol{A}t}dt$ 的通解公式是幂级数：

$$\int _{0}^{ \Delta t}e^{\boldsymbol{A}t}dt= \Delta t \left( I+ \frac{\boldsymbol{A} \Delta t}{2!}+ \frac{ \left( \boldsymbol{A} \Delta t \right) ^{2}}{3!}+ \frac{ \left( \boldsymbol{A} \Delta t \right) ^{3}}{4!}+ \ldots \right)$$

$$\boldsymbol{A}^{2} = \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 0\\ 0 & 0\\ \end{matrix} \right]$$

$\boldsymbol{A}$ 的更高阶次幂都是0.

$$\int _{0}^{ \Delta t}e^{\boldsymbol{A}t}dt= \Delta t \left( \left[ \begin{matrix} 1 & 0\\ 0 & 1\\ \end{matrix} \right] + \left[ \begin{matrix} 0 & 1\\ 0 & 0\\ \end{matrix} \right] \frac{ \Delta t}{2} \right) = \left[ \begin{matrix} \Delta t & \frac{1}{2}\Delta t^{2}\\ 0 & \Delta t\\ \end{matrix} \right]$$

$$\boldsymbol{G} = \int _{0}^{ \Delta t}e^{\boldsymbol{A}t}dt\boldsymbol{B} = \left[ \begin{matrix} \Delta t & \frac{1}{2}\Delta t^{2}\\ 0 & \Delta t\\ \end{matrix} \right] \left[ \begin{matrix} 0\\ 1\\ \end{matrix} \right] = \left[ \begin{matrix} \frac{1}{2}\Delta t^{2}\\ \Delta t\\ \end{matrix} \right]$$

现在我们得到状态外插方程：

$$\left[ \begin{matrix} p_{n+1}\\ \dot{p}_{n+1}\\ \end{matrix} \right] = \left[ \begin{matrix} 1 & \Delta t\\ 0 & 1\\ \end{matrix} \right] \left[ \begin{matrix} p_{n}\\ \dot{p}_{n}\\ \end{matrix} \right] + \left[ \begin{matrix} \frac{1}{2}\Delta t^{2}\\ \Delta t\\ \end{matrix} \right] a$$

示例 续 - 质量-弹簧-阻尼系统

回忆质量-弹簧-阻尼系统的状态空间方程：

$$\left[ \begin{matrix} \dot{p}\\ \dot{v}\\ \end{matrix} \right] = \left[ \begin{matrix} 0 & 1\\ -\frac{k}{m} & \frac{c}{m}\\ \end{matrix} \right] \left[ \begin{matrix} p\\ v\\ \end{matrix} \right] + \left[ \begin{matrix} 0\\ \frac{1}{m}\\ \end{matrix} \right] F$$

其中：

$$\boldsymbol{A} = \left[ \begin{matrix} 0 & 1\\ -\frac{k}{m} & \frac{c}{m}\\ \end{matrix} \right]$$

$$\boldsymbol{B} = \left[ \begin{matrix} 0\\ \frac{1}{m}\\ \end{matrix} \right]$$

本例中，矩阵指数的计算并不简单，因为此时矩阵A的高阶次幂不再是0了。

这个微分方程的求解细节超出了本教程的范围。