粗谈常系数线性常微分方程的求解

芋圆公式

最近发现大一同学学习《力学》课程时，对老师直接写出常系数高阶线性常微分方程的解感到困惑。不如搬运《常微分方程教程》一书的内容展示给大家。本帖内容会在数学专业春季学期开设的专业选修课《微分方程》中涉及。

称 $$ y^{(n)} + a_{1} y^{(n - 1)} + \dotsb + a_{n - 1} y' + a_{n} y = f(x) \tag{6.65} $$ 为 $n$ 阶线性常系数微分方程，相应的齐次线性方程为 $$ y^{(n)} + a_{1} y^{(n - 1)} + \dotsb + a_{n - 1} y' + a_{n} y = 0. \tag{6.66} $$ 其中，

$y^{(n)}$ 是指 $y$ 的 $n$ 阶导函数。当 $n = 1, 2$ 时，也可以换用我们更熟悉的记号 $y', y''$。
线性是指等号左边是 $y, y', \dotsc, y^{(n)}$ 的线性组合。如果还是不理解线性组合，则可以粗略地这么看：左边不会出现 $y$ 及其导数的乘积。
常系数是指 $a_{1}, \dotsc, a_{n}$ 是实常数。
$f(x)$ 是某个区间上的实值连续函数。

以后学到线性空间和仿射空间，就会发现：

方程 $(6.65)$ 的解构成仿射空间，而方程 $(6.66)$ 的解构成 $n$ 维线性空间。
方程 $(6.65)$ 的通解可表达为方程 $(6.65)$ 的某个解加上方程 $(6.66)$ 的通解。

结论的证明应该不难，有空补上。

齐次方程的通解

最让大一同学困惑的应当是方程 $(6.66)$ 的通解的模样。可能有的同学已经知道方程 $(6.66)$ 的特解由其特征方程的根决定。将 $(6.66)$ 中 $y^{(k)}$ 换成 $\lambda^{k}$ 就得到 $(6.66)$ 的特征方程（这里 $k = 0, 1, \dotsc, n$）。这些操作看起来毫无理由，因为涉及到的线性代数需要一定的篇幅来解释。

令 $y_{1} = y$ 、 $y_{2} = y'$ 、……、 $y_{n} = y^{(n - 1)}$ ，则方程 $(6.66)$ 被等价成常系数齐次线性方程组
$$
\frac{\mathrm{d} \boldsymbol{y}}{\mathrm{d} x} = A \boldsymbol{y},
\tag{6.67}
$$
其中 $$ \boldsymbol{y} = \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n - 1} \\ y{n} \end{pmatrix}, \qquad A = \begin{pmatrix} 0 & 1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 1 \\ -a_{n} & -a_{n - 1} & -a_{n - 2} & \cdots & -a_{2} & -a_{1} \end{pmatrix}. \tag{6.68} $$

考虑方程组 $(6.67)$ 在 $n = 1$ 时的情形：
$$
\frac{\mathrm{d} y}{\mathrm{d} x} = a y,
\tag{6.26}
$$
容易算得通解 $y = C \mathrm{e}^{a x}$（分离变量并单独讨论 $y \equiv 0$ 的情形即可），其中 $C \in \mathbb{R}$。对于任意自然数 $n$ 的方程组 $(6.67)$，这个结论便推广为 $\boldsymbol{y} = \mathrm{e}^{x A} \boldsymbol{C}$，其中涉及的矩阵指数定义为
$$
\mathrm{e}^{x A} := \sum_{k = 0}^{\infty} \frac{(x A)^{k}}{k!},
$$
$\boldsymbol{C}$ 为任意常数列向量。这里省略了大量的严格性讨论，包括

引入矩阵的范数，证明 $n$ 阶实方阵构成的线性空间装配范数后是完备的赋范空间（通常是展开讨论微积分的前提）。
矩阵的指数的定义是合适的。
性质：当矩阵 $A$ 和 $B$ 可交换时，$\mathrm{e}^{A + B} = \mathrm{e}^{A} \mathrm{e}^{B}$；矩阵指数总是可逆的，且 $\bigl( \mathrm{e}^{A} \bigr)^{-1} = \mathrm{e}^{-A}$；若 $P$ 是可逆矩阵，则 $\mathrm{e}^{P A P^{-1}} = P \mathrm{e}^{A} P^{-1}$。

容易证明，矩阵 $\mathrm{e}^{x A}$ 本身也满足方程 $(6.67)$。实际上这是因为矩阵 $\mathrm{e}^{x A}$ 的每一列都满足方程 $(6.67)$。若 $P$ 是可逆的 $n \times n$ 常数矩阵，则易证 $\mathrm{e}^{x A} P$ 也满足方程 $(6.67)$。

由矩阵指数的定义，可知计算矩阵指数等价于计算矩阵的任意幂次 $A^{k}$，但是这往往不好算。若 $A$ 相似于某个对角矩阵 $D = \operatorname{diag}(\lambda_{1}, \dotsc, \lambda_{n})$ ，也就是说存在可逆矩阵 $P$ 使得 $P^{-1} A P = D$，则容易证明 $A^{k} = P D^{k} P^{-1} = P \operatorname{diag}\bigl(\lambda_{1}^{k}, \dotsc, \lambda_{n}^{k} \bigr) P^{-1}$ 。但是一般的矩阵在实数范围内无法对角化，最多相似于它的 Jordan 标准型 $J$（某个和 $A$ 同样尺寸的矩阵），但是这对于本帖的问题已经足够了。

求 Jordan 标准型首先要求矩阵 $A$ 的特征值，因为这些特征值会依次排列在 $A$ 的 Jordan 标准型的对角线上。特征值满足特征方程 $\det(\lambda I_{n} - A) = 0$ ，其中 $$ \det(\lambda I_{n} - A) = \begin{vmatrix} \lambda & -1 & 0 & \cdots & 0 & 0 \\ 0 & \lambda & -1 & \ddots & \ddots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \lambda & -1 & 0 \\ 0 & 0 & 0 & \cdots & \lambda & -1 \\ a_{n} & a_{n - 1} & a_{n - 2} & \cdots & a_{2} & \lambda + a_{1} \end{vmatrix}. \tag{6.69} $$ 将行列式按最后一行展开就得到特征方程 $$ \lambda^{n} + a_{1} \lambda^{n - 1} + \dotsb + a_{n - 1} \lambda + a_{n} = 0. \tag{6.70} $$ 设特征方程在 $\mathbb{C}$ 中互不相同的解为 $\lambda_{1}, \dotsc, \lambda_{s}$ ，相应的重数分别为 $n_{1}, \dotsc, n_{s}$ ，即特征方程等号左边有因式分解：$$ (\lambda - \lambda_{1})^{n_{1}} \dotsm (\lambda - \lambda_{s})^{n_{s}} = 0. $$ 再经过一点“超纲”的内容的讨论，就得到矩阵 $A$ 的 Jordan 标准型 $$ J = \begin{pmatrix} J_{1} \\ ~ & J_{2} \\ ~ & ~ & \ddots \\ ~ & ~ & ~ & J_{s} \end{pmatrix}, $$ 其中 $J_{k}$ 是 $n_{k} \times n_{k}$ 的 Jordan 块：$$ J_{k} = \begin{pmatrix} \lambda_{k} & 1 \\ ~ & \lambda_{k} & \ddots \\ ~ & ~ & \ddots & 1 \\ ~ & ~ & ~ & \lambda_{k} \end{pmatrix}. $$ 接下来需要计算 Jordan 标准型的指数，这里省略详细过程（但是计算过程还挺有趣的）。$$ \mathrm{e}^{x J} = \begin{pmatrix} \mathrm{e}^{x J_{1}} \\ ~ & \mathrm{e}^{x J_{2}} \\ ~ & ~ & \ddots \\ ~ & ~ & ~ & \mathrm{e}^{x J_{s}} \end{pmatrix}, $$ 其中 $$ \mathrm{e}^{x J_{i}} = \mathrm{e}^{\lambda_{i} x} = \begin{pmatrix} 1 & x & \cdots & \dfrac{x^{n_{i} - 2}}{(n_{i} - 2)!} & \dfrac{x^{n_{i} - 1}}{(n_{i} - 1)!} \\ ~ & 1 & \ddots & \ddots & \dfrac{x^{n_{i} - 2}}{(n_{i} - 2)!} \\ \phantom{\dfrac{1}{2}} & ~ & \ddots & \ddots & \vdots \\ \phantom{\dfrac{1}{2}} & ~ & ~ & 1 & x \\ \phantom{\dfrac{1}{(2 - 2)}} & \phantom{\dfrac{1}{(2 - 2)}} & \phantom{\dfrac{1}{(2 - 2)}} & ~ & 1 \end{pmatrix}. \tag{6.32} $$ BTW，Jordan 块有如下分解式：$$ J_{i} = \begin{pmatrix} \lambda_{i} \\ ~ & \lambda_{i} \\ ~ & ~ & \ddots \\ ~ & ~ & ~ & \lambda_{i}, \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ ~ & 0 & \ddots \\ ~ & ~ & \ddots & 1 \\ ~ & ~ & ~ & 0 \end{pmatrix}, $$ 其中第二个矩阵是 $n_{i}$ 次幂零矩阵。这是计算 Jordan 块的指数的关键。

接下来考察 $P^{-1} A P = J$ 中的 $P = (p_{ij})_{n \times n}$。由于 $$ A P = P J = P \begin{pmatrix} J_{1} \\ ~ & J_{2} \\ ~ & ~ & \ddots \\ ~ & ~ & ~ & J_{s} \end{pmatrix}, \tag{6.72} $$ 则可断言，对 $$ m_{1} = 1, \quad m_{2} = n_{1} + 1, \quad \dotsc, \quad m_{s} = n_{1} + \dotsb + n_{s - 1} + 1, $$ 有 $p_{1, m_{j}} \neq 0$ ，$j = 1, 2, \dotsc, s$。否则，假设存在 $k$ 使得 $p_{1, m_{k}} = 0$ ，结合式 $(6.32)$ 对照式 $(6.72)$ 等号两边矩阵的前 $n - 1$ 行第 $m_{k}$ 列元素，有 $$ 0 = p_{1, m_{k}} = p_{2, m_{k}} = \dotsb = p_{n, m_{k}}, $$ 说明 $P$ 中第 $m_{k}$ 列是零向量，这与 $P$ 可逆矛盾。

由前述性质 $\mathrm{e}^{P^{-1} A P} = P^{-1} \mathrm{e}^{A} P$ 可知 $$ \mathrm{e}^{x A} P = P \begin{pmatrix} \mathrm{e}^{x J_{1}} \\ ~ & \mathrm{e}^{x J_{2}} \\ ~ & ~ & \ddots \\ ~ & ~ & ~ & \mathrm{e}^{x J_{s}} \end{pmatrix}. \tag{6.35} $$ 结合式 $(6.32)$ 可知 $\mathrm{e}^{x A} P$ 的第一行元素为 $$ \begin{aligned} & p_{1, m_{1}} \mathrm{e}^{\lambda_{1} x}, \quad * \mathrm{e}^{\lambda_{1} x} + p_{1, m_{1}} x \mathrm{e}^{\lambda_{1} x}, \quad \dotsc, \quad * \mathrm{e}^{\lambda_{1} x} + \dotsb + \frac{p_{1, m_{1}}}{(n_{1} - 1)!} x^{n_{1} - 1} \mathrm{e}^{\lambda_{1} x}, \\ & \cdots\cdots \\ & p_{1, m_{s}} \mathrm{e}^{\lambda_{s} x}, \quad * \mathrm{e}^{\lambda_{s} x} + p_{1, m_{s}} x \mathrm{e}^{\lambda_{s} x}, \quad \dotsc, \quad * \mathrm{e}^{\lambda_{s} x} + \dotsb + \frac{p_{1, m_{s}}}{(n_{s} - 1)!} x^{n_{s} - 1} \mathrm{e}^{\lambda_{s} x}, \end{aligned} $$ 其中 $*$ 表示常值系数。对 $\mathrm{e}^{x A} P$ 施加一系列初等列变换（等价地右乘可逆矩阵 $Q$，而且 $\mathrm{e}^{x A} P Q$ 也满足方程 $(6.67)$），可以将第一行元素整理为 $$ \begin{cases} \mathrm{e}^{\lambda_{1} x}, x \mathrm{e}^{\lambda_{1} x}, \dotsc, x^{n_{1} -1} \mathrm{e}^{\lambda_{1} x}, \\ \cdots\cdots \\ \mathrm{e}^{\lambda_{s} x}, x \mathrm{e}^{\lambda_{1} x}, \dotsc, x^{n_{s} - 1} \mathrm{e}^{\lambda_{s} x}. \end{cases} \tag{6.71} $$ 给 $\mathrm{e}^{x A} P Q$ 右乘任意实数列向量 $\boldsymbol{C}$ 可得 $\boldsymbol{y}$ 的第一个分量的通解。这样我们就证明了

定理 6.6*： 式 $(6.71)$ 中的 $n$ 个函数是方程 $(6.66)$ 的一组线性独立的特解，它们的线性组合是方程 $(6.66)$ 通解。

另外，特征方程 $(6.70)$ 虽然是实方程，但是可以成对地解出共轭复根 $\lambda, \overline{\lambda}$。共轭复根对应的特解是复函数 $$ \begin{aligned} z_{+}(x) &= \mathrm{e}^{\lambda x} = \mathrm{e}^{x (\operatorname{Re} \lambda + \mathrm{i} \operatorname{Im} \lambda)} = \mathrm{e}^{x \operatorname{Re} \lambda} \bigl( \cos (x \operatorname{Im} \lambda) + \mathrm{i} \sin (x \operatorname{Im} \lambda) \bigr), \\ z_{-}(x) &= \mathrm{e}^{\overline{\lambda} x} = \mathrm{e}^{x (\operatorname{Re} \lambda - \mathrm{i} \operatorname{Im} \lambda)} = \mathrm{e}^{x \operatorname{Re} \lambda} \bigl( \cos (x \operatorname{Im} \lambda) - \mathrm{i} \sin (x \operatorname{Im} \lambda) \bigr). \end{aligned} $$ 它们的线性组合仍然是原方程的特解，比如最终想要得到的线性独立的实函数 $$ \begin{aligned} u(x) &= \frac{z_{+}(x) + z_{-}(x)}{2} = \mathrm{e}^{x \operatorname{Re} \lambda} \cos (x \operatorname{Im} \lambda), \\ v(x) &= \frac{z_{+}(x) - z_{-}(x)}{2 \mathrm{i}} = \mathrm{e}^{x \operatorname{Re} \lambda} \sin (x \operatorname{Im} \lambda). \end{aligned} $$

非齐次方程的通解

非齐次情形 $(6.65)$ 普适的特解由下面定理给出，定理不证。

定理 6.3*： 设 $\varphi_{1}, \dotsc, \varphi_{n}$ 张成方程 $(6.66)$ 的解空间，即方程 $(6.66)$ 的通解为 $$ y = C_{1} \, \varphi_{1}(x) + \dotsb + C_{n} \, \varphi_{n}(x), $$ 其中 $C_{1}, \dotsc, C_{n}$ 为任意常数。那么方程 $(6.65)$ 的通解为 $$ y = C_{1} \, \varphi_{1}(x) + \dotsb + C_{n} \, \varphi_{n}(x) + \varphi^{*}(x), \tag{6.57} $$ 其中 $$ \varphi^{*}(x) = \sum_{k = 1}^{n} \varphi_{k}(x) \int_{x_{0}}^{x} \frac{W_{k}(s)}{W(s)} \, f(s) \,\mathrm{d}{s}. \tag{6.58} $$ $W(x)$ 是 $\varphi_{1}, \dotsc, \varphi_{n}$ 的 Wroński 行列式：$$ W(x) = \begin{vmatrix} \varphi_{1}(x) & \varphi_{2}(x) & \cdots & \varphi_{n}(x) \\ \varphi_{1}'(x) & \varphi_{2}'(x) & \cdots & \varphi_{n}'(x) \\ \vdots & \vdots & \ddots & \vdots \\ \varphi_{1}^{(n - 1)}(x) & \varphi_{2}^{(n - 1)}(x) & \cdots & \varphi_{n}^{(n - 1)}(x) \\ \end{vmatrix}, \tag{6.53} $$ $W_{k}(x)$ 是 $W(x)$ 的第 $n$ 行第 $k$ 列元素的代数余子式（即，用 $(0, 0, \dotsc, 0, 1)^{\text{T}}$ 替换 $W(x)$ 的第 $k$ 列后得到的行列式）。

当方程 $(6.65)$ 中的 $f(x)$ 取特殊的形式时，可以按照经验先假设特解 $\varphi^{*}(x)$ 的形式，然后用待定系数法确定特解。

若 $f(x) = P_{m}(x) \, \mathrm{e}^{\mu x}$，其中 $P_{m}(x)$ 是 $x$ 的 $m$ 次多项式，$\mu$ 不是特征方程的根。方程的特解应为 $\varphi^{*}(x) = Q_{m}(x) \, \mathrm{e}^{\mu x}$ ，其中多项式 $Q_{m}(x)$ 的系数待定。
在上一条中，修改条件：设 $\mu$ 是特征方程的 $k$ 重根，则特解应为 $\varphi^{*}(x) = x^{k} \, Q_{m}(x) \, \mathrm{e}^{\mu x}$。
若 $f(x) = \bigl( A_{m}(x) \cos(\beta x) + B_{l}(x) \sin(\beta x) \bigr) \mathrm{e}^{\alpha x}$ ，其中 $A_{m}(x)$ 和 $B_{l}(x)$ 分别为 $x$ 的 $m$ 次和 $l$ 次多项式。那么特解应为 $\varphi^{*}(x) = x^{k} \bigl( C_{n}(x) \cos(\beta x) + D_{n}(x) \sin(\beta x) \bigr) \mathrm{e}^{\alpha x}$ ，其中非负整数 $k$ 是特征根 $\alpha \pm \mathrm{i} \beta$ 的重数、 $n = \max\{m, l\}$ 、 $C_{n}(x)$ 和 $D_{n}(x)$ 为 $n$ 次多项式。

回到物理：理论力学

现在可以求解理论力学中的自由振动、阻尼振动、强迫振动的运动方程了（例子可以不提）。

自由振动

重新定义系数后运动方程为
$$
\ddot{x} + \omega^{2} x = 0,
$$
特征方程为 $\lambda^{2} + \omega^{2} = 0$，有一对共轭复根 $\lambda_{\pm} = \pm \mathrm{i} \omega$，因此通解为
$$
x(t) = A \cos \omega t + B \sin \omega t.
$$

阻尼振动

重新定义系数后运动方程为 $$ \ddot{x} + 2 \gamma \dot{x} + \omega_{0}^{2} x = 0, $$ 其中 $\gamma < \omega_{0}$ 。特征方程为 $\lambda^{2} + 2 \gamma \lambda+ \omega_{0}^{2} = 0$ ，有一对共轭复根 $\lambda_{\pm} = -\gamma \pm \mathrm{i} \underbrace{\sqrt{\omega_{0}^{2} - \gamma^{2}}}_{{} =: \omega}$ ，因此通解为
$$
x(t) = \mathrm{e}^{-\gamma t} (A \cos \omega t + B \sin \omega t).
$$

非周期阻尼

运动方程仍然为 $$ \ddot{x} + 2 \gamma \dot{x} + \omega_{0}^{2} x = 0, $$ 其中 $\gamma \geqslant \omega_{0}$ 。容易算得

$\gamma > \omega_{0}$ 时， $x(t) = A \exp\left( -\left( \gamma - \sqrt{\gamma^{2} - \omega_{0}^{2}} \right) t \right) + B \exp\left( -\left( \gamma + \sqrt{\gamma^{2} - \omega_{0}^{2}} \right) t \right)$ ；
$\gamma = \omega_{0}$ 时，$x(t) = (A + B t) \mathrm{e}^{-\gamma t}$。

（无摩擦）强迫振动

重新定义系数后运动方程为
$$
\ddot{x} + \omega^{2} x = \frac{f}{m} \cos(\gamma t + \beta),
$$
其相应的齐次方程的通解为 $A \cos\omega t + B \sin\omega t$。

非共振

此时 $\gamma \neq \omega$，假设特解 $x^{*}(t) = b \cos(\gamma t + \beta)$（这里的假设和前述的理论部分不一样，但是容易证明这样假设是对的），代入方程算得 $b = \dfrac{f}{m \bigl( \omega^{2} - \gamma^{2} \bigr)}$。于是
$$
x(t) = A \cos\omega t + B \sin\omega t + \frac{f}{m \bigl( \omega^{2} - \gamma^{2} \bigr)} \cos(\gamma t + \beta).
$$

共振

此时 $\gamma = \omega$，而 $\pm\mathrm{i} \omega$ 是特征方程的 $1$ 重根，因此假设特解 $x^{*}(t) = t (a \cos\omega t + b \sin\omega t)$，代入方程算得系数 $a = \dfrac{f}{2 m \omega} \sin\beta$ 和 $b = \dfrac{f}{2 m \omega} \cos\beta$。于是整理可得
$$
x(t) = A \cos\omega t + B \sin\omega t + \frac{f}{2 m \omega} t \sin(\omega t + \beta).
$$

有摩擦强迫振动

只讨论 $\gamma < \omega_{0}$ 的情形。运动方程为 $$ \ddot{x} + 2 \eta \dot{x} + \omega_{0}^{2} x = \frac{f}{m} \cos\gamma t, $$ 其相应的齐次方程的通解为 $\mathrm{e}^{-\gamma t} (A \cos\omega t + B \sin\omega t)$，其中 $\omega = \sqrt{\omega_{0}^{2} - \gamma^{2}}$ 。假设特解 $x^{*}(t) = a \cos\gamma t + b \sin\gamma t$，代入方程算得系数 $a = \dfrac{f}{m} \dfrac{\omega_{0}^{2} - \gamma^{2}}{\bigl( \omega_{0}^{2} - \gamma^{2} \bigr) + 4 \eta^{2} \gamma^{2}}$ 和 $b = \dfrac{f}{m} \dfrac{2 \eta \gamma}{\bigl( \omega_{0}^{2} - \gamma^{2} \bigr) + 4 \eta^{2} \gamma^{2}}$ 。整理可得 $$ x(t) = \mathrm{e}^{-\gamma t} (A \cos\omega t + B \sin\omega t) +\frac{f}{m} \frac{1}{\sqrt{\bigl( \omega_{0}^{2} - \gamma^{2} \bigr) + 4 \eta^{2} \gamma^{2}}} \cos(\gamma t - \delta), $$ 其中 $\tan\delta = \dfrac{2 \eta \gamma}{\omega_{0}^{2} - \gamma^{2}}$。

cai_b

太棒了！！！

芋圆公式

数学物理方法

题目来自明平兵老师班的作业。

Image description

叠加态

《利用$\text{Mathematica}$计算》，作为作业确实有点复杂了

芋圆公式

积分因子法

下面的例子来自舒幼生《力学》。

Image description