L23-微分方程和exp(At)

参考

Unit II: Least Squares, Determinants and Eigenvalues - Differential Equations and exp(At) | pdf
课本章节：Read Section 6.3 in the 4^th or 5^th edition.
练习题：L23-微分方程和exp(At)-习题集

微分方程组

两个变量 $u_{1}$ 和 $u_{2}$ 随着时间 $t$ 变化（即它们都是关于 $t$ 的函数），它们的变化率如下

\begin{aligned} \left\{\begin{matrix} \cfrac{\mathrm{d} u_{1}}{\mathrm{d} t}=-u_{1}+2u_{2} \\ \cfrac{\mathrm{d} u_{2}}{\mathrm{d} t}=u_{1}-2u_{2} \end{matrix}\right. \end{aligned}

已知在初始状态 $t=0$ 时 $u_{1}=1$ 、 $u_{2}=0$

提示

由于 $u_{1}$ 的变化率是 $-u_{1}+2u_{2}$ ， $u_{2}$ 的变化率是 $u_{1}-2u_{2}$ ，所以两者的变化是相互影响的

并结合初始值，可以推测出 $u_{1}$ 是随着时间 $t$ 的增加而减小的，由于它的变化率是负值；而 $u_{2}$ 的变化则是相反的。这个变化规律可以通过求出 $u_{1}$ 和 $u_{2}$ 具体的函数表达式来验证。

根据以上条件求解 $u_{1}$ 和 $u_{2}$ （函数的具体表达式）

说明

以上两个变化率等式都是一阶线性微分方程

微分方程是指具有函数关系以及含有一个或以上的导数的方程
一阶是指等式中的导数部分（对于这个例子就是指等式的左侧）只包含一阶导数 $\cfrac{\mathrm{d} y}{\mathrm{d} x}$ 而没有更高阶的导数
线性是指等式中的其他部分（对于这个例子就是指等式的右侧）是线性的 ❓，没有高阶的变量

使用线性代数求解微分方程组

首先令 $u=\begin{bmatrix} u_{1} \\ u_{2} \end{bmatrix}$ 将以上微分方程组转换为矩阵的形式（将系数矩阵记作 $A$ ）

\cfrac{\mathrm{d} u}{\mathrm{d} t}= \begin{bmatrix} -1 & 2 \\ 1 & -2 \end{bmatrix}u= Au

则在 $t=0$ 时 $u(0)=\begin{bmatrix} u_{1} \\ u_{2} \end{bmatrix}=\begin{bmatrix} 1 \\ 0 \end{bmatrix}$

提示

经过转换后的等式 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 是一阶常系数线性方程，这里的常系数是指矩阵 $A$ ，而且没有常数项

对于这种形式的微分方程，它的解 $u(t)$ 为指数形式

转换后的微分方程的通解形式是 $u(t)=c_{1}e^{\lambda_{1}t}x_{1}+c_{2}e^{\lambda_{2}t}x_{2}$

其中 $c_{1}$ 和 $c_{2}$ 是系数， $\lambda_{1}$ 和 $\lambda_{2}$ 是矩阵 $A$ 的特征值， $x_{1}$ 和 $x_{2}$ 是矩阵 $A$ 的特征向量

提示

在课堂上并没有直接推导证明这个通解，而是通过代入特解进行验证

以上的通解 $u(t)$ 可以看作是由两个向量 $e^{\lambda_{1}t}x_{1}$ 和 $e^{\lambda_{2}t}x_{2}$ 通过线性组合所构成的（线性组合的系数为 $c_{1}$ 和 $c_{2}$ ），可以代入其中一个特解进行验证

当 $u=e^{\lambda_{1}t}x_{1}$ 时

对于 $\cfrac{\mathrm{d} u}{\mathrm{d} t}$ 可得

\cfrac{\mathrm{d} u}{\mathrm{d} t}=\cfrac{\mathrm{d}}{\mathrm{d} t}(e^{\lambda_{1}t}x_{1})=\lambda_{1}e^{\lambda_{1}t}x_{1}

而对于 $Au$ 可得

Au=A(e^{\lambda_{1}t}x_{1})

对比两个等式

\begin{aligned} A{\color{Red}e^{\lambda_{1}t} }x_{1} &\leftrightarrow \lambda_{1}{\color{Red}e^{\lambda_{1}t} }x_{1} \\ &\Downarrow \\ Ax_{1} &\leftrightarrow \lambda_{1}x_{1} \end{aligned}

再根据特征向量的定义可知 $Ax_{1} = \lambda_{1}x_{1}$ 即以上左右两部分是相等的

所以 $u=e^{\lambda_{1}t}x_{1}$ 是 $u(t)$ 的特解，使用同样的方法可以验证另一个特解

所以 $u(t)=c_{1}e^{\lambda_{1}t}x_{1}+c_{2}e^{\lambda_{2}t}x_{2}$ 是形式如 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 的微分方程的通解

在微分方程的通解中 $u(t)=c_{1}e^{\lambda_{1}t}x_{1}+c_{2}e^{\lambda_{2}t}x_{2}$ 特征值是指数形式的

可以与上一节课的结论进行横向对比，在差分方程的通解中 $u_{k}=c_{1}\lambda_{1}^{k}x_{1}+c_{2}\lambda_{2}^{k}x_{2}+\dots+c_{n}\lambda_{n}^{k}x_{n}$ 特征值是幂形式的

根据通解的形式，则需要求出矩阵 $A$ 的特征值和特征向量

A= \begin{bmatrix} -1 & 2 \\ 1 & -2 \end{bmatrix}

观察矩阵 $A$ 可知它是一个奇异矩阵/不可逆矩阵，则 $0$ 必然是该矩阵的特征值之一 $\lambda_{1}=0$ ，再根据特征值的特性可知矩阵的迹 $trace(A)=\lambda_{1}+\lambda_{2}=(-1)+(-2)=-3$ ，可得 $\lambda_{2}=-3$

所以矩阵 $A$ 的特征值是 $\lambda_{1}=0$ 和 $\lambda_{2}=-3$

提示

也可以通过求解特征方程 $det(A-\lambda I)=0$ 来得到矩阵的所有特征值

\begin{aligned} det(A-\lambda I)&= \begin{bmatrix} -1-\lambda & 2 \\ 1 & -2-\lambda \end{bmatrix} \\ &=(-1-\lambda)(-2-\lambda) - 2 \times 1 \\ &=\lambda^{2}+3\lambda \\ &=\lambda(\lambda+3) \\ &=0 \end{aligned}

解得 $\lambda_{1}=0$ 或 $\lambda_{2}=-3$

将以上求得的特征值分别代入方程 $(A-\lambda I)x=0$ 求出相应的特征向量

当 $\lambda_{1}=0$ 时 $\begin{aligned} (A-\lambda I)x&=Ax \\ &= \begin{bmatrix} -1 & 2 \\ 1 & -2 \end{bmatrix}x \\ &=0 \end{aligned}$
其中一个特解 $x_{1}=\begin{bmatrix} 2 \\ 1 \end{bmatrix}$
当 $\lambda_{1}=-3$ 时 $\begin{aligned} (A-\lambda I)x&=Ax \\ &= \begin{bmatrix} 2 & 2 \\ 1 & 1 \end{bmatrix}x \\ &=0 \end{aligned}$
其中一个特解 $x_{2}=\begin{bmatrix} 1 \\ -1 \end{bmatrix}$

将以上所得的矩阵 $A$ 的特征值和特征向量代入通解 $u(t)$ 的表达式可得

\begin{aligned} u(t)&=c_{1}e^{0t}\begin{bmatrix} 2 \\ 1 \end{bmatrix}+c_{2}e^{-3t}\begin{bmatrix} 1 \\ -1 \end{bmatrix} \\ &=c_{1}\begin{bmatrix} 2 \\ 1 \end{bmatrix}+c_{2}e^{-3t}\begin{bmatrix} 1 \\ -1 \end{bmatrix} \end{aligned}

再根据初始状态求出系数

\begin{aligned} u(0)&=c_{1}\begin{bmatrix} 2 \\ 1 \end{bmatrix}++c_{2}e^{-3 \times 0}\begin{bmatrix} 1 \\ -1 \end{bmatrix} \\ &=c_{1}\begin{bmatrix} 2 \\ 1 \end{bmatrix}+c_{2}\begin{bmatrix} 1 \\ -1 \end{bmatrix} \\ &=\begin{bmatrix} 1 \\ 0 \end{bmatrix} \end{aligned}

提示

其实可以将以上等式写成矩阵形式

\begin{aligned} u(0)&=c_{1}\begin{bmatrix} 2 \\ 1 \end{bmatrix}+c_{2}\begin{bmatrix} 1 \\ -1 \end{bmatrix} \\ &=\begin{bmatrix} 2 & 1 \\ 1& -1 \end{bmatrix} \begin{bmatrix} c_{1} \\ c_{2} \end{bmatrix} \\ &=\begin{bmatrix} 1 \\ 0 \end{bmatrix} \end{aligned}

那么通解在初始状态用矩阵表示就是 $u(0)=Sc$ ，其中 $S$ 是特征值构成的矩阵， $c$ 是系数构成的列向量

将以上等式转换为方程组形式可得

\begin{aligned} \left\{\begin{matrix} 2c_{1}+c_{2}=1 \\ c_{1}-c_{2}=0 \end{matrix}\right. \end{aligned}

解得

\begin{aligned} \left\{\begin{matrix} c_{1}=\cfrac{1}{3} \\ c_{2}=\cfrac{1}{3} \end{matrix}\right. \end{aligned}

所以微分方程 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 的通解为 $u(t)=\cfrac{1}{3}\begin{bmatrix} 2 \\ 1 \end{bmatrix}+\cfrac{1}{3}e^{-3t}\begin{bmatrix} 1 \\ -1 \end{bmatrix}$

提示

将微分方程 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 的通解 $u(t)$ 进行「拆分」就可以得到原微分方程组 $u_{1}$ 和 $u_{2}$ 的通解

变化趋势与稳定性

观察微分方程的通解 $u(t)=\cfrac{1}{3}\begin{bmatrix} 2 \\ 1 \end{bmatrix}+\cfrac{1}{3}e^{-3t}\begin{bmatrix} 1 \\ -1 \end{bmatrix}$ 可知它随时间的变化情况：

开始状态为 $u(0)=\begin{bmatrix} 1 \\ 0 \end{bmatrix}$ 即 $u_{1}=1$ 、 $u_{2}=0$
然后随着时间的增大 $u_{1}$ 的值减小， $u_{2}$ 的值增大
当 $t \to \infty$ 时， $u(\infty)$ 趋于稳态 $u(\infty) \to \cfrac{1}{3}\begin{bmatrix} 2 \\ 1 \end{bmatrix}$ ，即 $u_{1} \to \cfrac{2}{3}$ 、 $u_{2} \to \cfrac{1}{3}$

由于在时间增大的全过程中 $u_{1}$ 和 $u_{2}$ 之和都为 $1$ ，从开始的 $u_{1}=1$ 、 $u_{2}=0$ 变成最终的 $u_{1} \to \cfrac{2}{3}$ 、 $u_{2} \to \cfrac{1}{3}$ 可以将此变化看作是 $u_{1}$ 的部分值「流向」了 $u_{2}$ ，并最后趋于稳态 steady

稳定性

对于微分方程 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 的通解 $u(t)=c_{1}e^{\lambda_{1}t}x_{1}+c_{2}e^{\lambda_{2}t}x_{2}$ 当 $t \to \infty$ 时，其稳定性与系数矩阵 $A$ 的特征值相关：

Stability 稳定：如果所有特征值的实数部分都小于零 $Re(\lambda)<0$ ，则当 $t \to \infty$ 函数 $u(t) \to 0$ 都会趋于零
Steady State 稳态：如果其中一个特征值的实数部分等于零 $Re(\lambda_{i})=0$ ，而其余的特征值的实数部分都小于零 $Re(\lambda)<0$ ，则当 $t \to \infty$ 函数会趋于一个数 $u(t) \to c_{i}x_{i}$ （和 Stability 稳定一样，也是会趋于一个固定值，但这个值和哪个特征值为 $0$ 相关）
Blow up 无法收敛：如果所有特征值的实数部分都大于零 $Re(\lambda)>0$ ，则当 $t \to \infty$ 函数 $u(t)$ 无法收敛

对于微分方程 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 如果其系数矩阵 $A$ 是一个 $2 \times 2$ 矩阵，由于它的特征值只有两个，所以可以直接通过矩阵的元素来判断通解 $u(t)$ 的稳定性

对于 $A_{2 \times 2}$ 矩阵

A= \begin{bmatrix} a & b \\ c & d \end{bmatrix}

根据特征值的特性可知

\begin{aligned} \left\{\begin{matrix} trace(A)=\lambda_{1}+\lambda_{2}=a+d \\ detA=\lambda_{1}\lambda_{2} \end{matrix}\right. \end{aligned}

如果矩阵的迹小于零 $trace(A)<0$ 并且矩阵的行列式大于零 $detA>0$ ，则可以推导出矩阵 $A$ 的两个特征值 $\lambda_{1}$ 和 $\lambda_{2}$ 都是负数，则当 $t \to \infty$ 微分方程的通解会处于Stability 稳定 $u(t) \to 0$

注意前面在讨论特征值时，都只是考虑特征值的实数部分 $Re(\lambda)$ （特征值可能是复数，具有虚数部分），由于将虚数映射到直角坐标系后，（在几何形状上）它是绕着一个单位圆变化的 ❓ 所以该部分的值对于函数的收敛性并没有很大的影响，可以忽略不计

解耦

使用矩阵 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 表达微分方程组时，矩阵 $A$ 的作用相当于将变量 $u_{1}$ 和 $u_{2}$ 进行耦合，这可能会让求解过程变得更复杂

可以利用矩阵的对角化进行解耦，即将以上等式中的列向量 $u$ 表达为矩阵 $A$ 的特征向量的某种线性组合，通过换元可以得到变量之间解耦的微分方程，便于求解

如果使用矩阵 $A$ 的特征向量 $x_{i}$ 为基，则列向量 $u$ 可以表达为

u=Sv

其中矩阵 $S$ 是由特征向量构成的列向量，列向量 $v$ 则是线性组合的系数构成的

说明

由于矩阵 $A$ 是根据微分方程组构造的，而矩阵 $S$ 是由矩阵 $A$ 的特征向量所构成的，所以矩阵 $S$ 是可确定的，这里看作为常量

而列向量 $v$ 是变量，它跟随原向量 $u$ 的不同而变化

经过换元后，微分方程组从原来针对 $u$ 求导，变成了对 $v$ 求导，其形式如下

\cfrac{\mathrm{d} (Sv)}{\mathrm{d} t}=S\cfrac{\mathrm{d} v}{\mathrm{d} t}=ASv

可以在等式两边同时乘上逆矩阵 $S^{-1}$ 可得

\cfrac{\mathrm{d} v}{\mathrm{d} t}=S^{-1}ASv

再根据矩阵对角化的公式 $S^{-1}AS=\Lambda$ 进一步化简

\cfrac{\mathrm{d} v}{\mathrm{d} t}=S^{-1}ASv=\Lambda v

对角矩阵的解耦

由于对角矩阵只有对角线上存在元素，其他位置的元素都是 $0$

\begin{bmatrix} \lambda_{1} & 0 & \dots & 0 \\ 0 & \lambda_{2} & \dots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 & \dots & \lambda_{n} \end{bmatrix}

如果这种矩阵与列向量相乘

\begin{bmatrix} \lambda_{1} & 0 & \dots & 0 \\ 0 & \lambda_{2} & \dots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 & \dots & \lambda_{n} \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{bmatrix}= \begin{bmatrix} \lambda_{1}x_{1} \\ \lambda_{2}x_{2} \\ \vdots \\ \lambda_{n}x_{n} \end{bmatrix}

得到的结果列向量的每个元素中的变量，都只包含原来列向量的相应行的元素（变量），即对角矩阵并没有让不同行进行耦合，这样可以让结果列向量便于后续的计算

由于换元后微分方程的系数矩阵 $\Lambda$ 为对角矩阵，所以如果将矩阵形式转换为方程组的形式，则每个微分方程等式中都只会包含一个变量，即 $\cfrac{\mathrm{d} v_{i}}{\mathrm{d} t}=\lambda_{i}v_{i}$ （而不像原微分方程组中，每个方程等式都包含 $u_{1}$ 和 $u_{2}$ ），这样就实现了解耦

注意

进行以上变换操作需要矩阵 $A$ 满足对角化的前提，即矩阵 $A$ 具有 $n$ 个线性不相关的特征向量

指数矩阵

对于换元后所得的微分方程 $\cfrac{\mathrm{d} v}{\mathrm{d} t}=\Lambda v$ 参考前面的示例，其通解为

v(t)=e^{\Lambda t}v(0)

其中 $e^{\Lambda t}$ 称为矩阵指数

那么对于换元前的微分方程 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 其（使用矩阵指数）对应的通解为

u(t)=e^{At}u(0)

提示

在课堂上并没有证明以上的通解表达式，它的形式也是由指数构成的，可以通过参考前面的微分方程组的示例，类推而得到该表达式

矩阵指数 Matrix exponential $e^{At}$ （其中 $A$ 为矩阵），其定义可以由指数的幂级数 power series 展开形式给出

e^{At}=I+At+\cfrac{(At)^{2}}{2}+\cfrac{(At)^{3}}{3}+\dots

幂级数展开式

两个经典的幂级数展开形式（泰勒级数 ❓）

\begin{aligned} {\color{Orange}e^{x} }&=\sum_{n=0}^{\infty} \cfrac{x^{n}}{n!} \\ &=1+x+\cfrac{x^{2}}{2}+\cfrac{x^{3}}{6}+\dots \end{aligned}

\begin{aligned} {\color{Orange}\cfrac{1}{1-x} }&=\sum_{n=0}^{\infty}x^{n} \\ &=1+x+x^{2}+x^{3}+\dots \end{aligned}

对于矩阵也同样适用

根据第一个经典幂级数展开式，对矩阵指数 $e^{At}$ 进行定义

{\color{Orange}e^{At} }=I+At+\cfrac{(At)^{2}}{2}+\cfrac{(At)^{3}}{3}+\dots

对于以上展开式，随着 $n$ 增大（展开式靠后面的项）对应项的值会趋于 $0$ ，所以该展开式更容易收敛

另一个经典的展开式则是 $(I-At)^{-1}$

{\color{Orange}(I-At)^{-1} }=I+At+(At)^{2}+(At)^{3}+\dots

对于以上展开式，如果从变量 $t$ 来考虑，只有当它很小时后面的高阶项（如 $t^{2}$ 、 $t^{3}$ 等，会趋于 $0$ ）才可以忽略不计，则可以利用该展开式求得逆矩阵的近似值 $(I-At)^{-1} \approx I+At$

如果从矩阵 $A$ 来考虑，展开式的收敛需要满足特定的条件

根据上一节课可知矩阵的幂可以分解为 $(A)^{k}=S\Lambda^{k}S$

则当矩阵 $At$ 的各个特征值都满足 $\left | \lambda(At) \right |<1$ 以上展开式才会收敛

对于矩阵指数 $e^{At}$ （其中 $A$ 为矩阵），可以基于幂级数展开式和矩阵对角化进行转换 $e^{At}=Se^{\Lambda t}S^{-1}$

其中矩阵 $S$ 是特征向量矩阵（由矩阵 $A$ 的特征向量作为列向量构成的），矩阵 $\Lambda$ 是特征值矩阵（它是一个对角矩阵，对角线上的元素是矩阵 $A$ 的特征值，其他元素都是 $0$ ）

注意

进行以上变换操作需要矩阵 $A$ 满足对角化的前提，即矩阵 $A$ 具有 $n$ 个线性不相关的特征向量

证明

基于幂级数展开式和矩阵对角化对矩阵指数 $e^{At}$ （其中 $A$ 为矩阵）进行转换

e^{At}=I+At+\cfrac{(At)^{2}}{2}+\cfrac{(At)^{3}}{6}+\dots

根据矩阵的对角化分解公式 $A={\color{Red}S \Lambda S^{-1} }$ 以及高阶矩阵（矩阵的幂）的分解公式 $A^{k}={\color{Red}S \Lambda^{k} S^{-1}}$ 可以将以上式子进行变换

\begin{aligned} e^{At}&=I+At+\cfrac{(At)^{2}}{2}+\cfrac{(At)^{3}}{6}+\dots \\ &=SS^{-1}+{\color{Red}S \Lambda S^{-1}}t+\cfrac{{\color{Red}S \Lambda^{2} S^{-1}}t^{2}}{2}+\cfrac{{\color{Red}S \Lambda^{3} S^{-1} }t^{3}}{6}+\dots \end{aligned}

将以上展开式的每项中的公因子 $S$ 和 $S^{-1}$ 提出

e^{At}=S(I+\Lambda t+\cfrac{\Lambda^{2}t^{2}}{2}+\cfrac{\Lambda^{3}t^{3}}{6}+\dots)S^{-1}

由于 $e^{\Lambda t}$ 的幂级数展开式为 $I+\Lambda t+\cfrac{\Lambda^{2}t^{2}}{2}+\cfrac{\Lambda^{3}t^{3}}{6}+\dots$ 所以以上展开式可以进一步简化为

\begin{aligned} e^{At}&=S{\color{Red}(I+\Lambda t+\cfrac{\Lambda^{2}t^{2}}{2}+\cfrac{\Lambda^{3}t^{3}}{6}+\dots) }S^{-1} \\ &=S{\color{Red}e^{\Lambda t} }S^{-1} \end{aligned}

所以对于微分方程 $\cfrac{\mathrm{d} u}{\mathrm{d} t}=Au$ 的矩阵指数形式的通解可以进一步化简

u(t)=e^{At}u(0)=Se^{\Lambda t}S^{-1}u(0)

由于 $e^{\Lambda t}$ 中的 $\Lambda$ 是对角矩阵（除了对角线上的元素，其他位置上的元素都是 $0$ ），所以该矩阵指数除了幂级数展开的形式，也具有矩阵形式

e^{\Lambda t}= \begin{bmatrix} e^{\lambda_{1}t} & 0 & \dots & 0 \\ 0 & e^{\lambda_{2}t} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & e^{\lambda_{n}t} \end{bmatrix}

与幂级数展开形式（有无数项相加）相比，它的矩阵形式仅含有限的 $n \times n$ 个元素，更便于计算

应用

对于二阶微分方程，可以通过构造微分方程组并写成矩阵的形式，将其转换为一阶微分方程的问题，然后再根据前面小节的步骤进行求解

提示

二阶是指等式中的导数部分具有的最高导数是二阶导 $y''$ 而没有更高阶的导数

已知二阶微分方程 $y''+by'+ky=0$

令 $u=\begin{bmatrix} y' \\ y \end{bmatrix}$

则以下微分方程组

\begin{aligned} \left\{\begin{matrix} y''+by'+ky=0 \\ y'=y' \end{matrix}\right. \end{aligned}

可以写成矩阵形式

u'= \begin{bmatrix} y'' \\ y' \end{bmatrix}= \begin{bmatrix} -b & -k \\ 1 & 0 \end{bmatrix} \begin{bmatrix} y' \\ y \end{bmatrix}

将系数矩阵记为 $A$ 则以上等式为 $u'=Au$ 它就是一阶微分方程

提示

可以将此方法进行推广，应用于含有 $k$ 阶导的微分方程，通过类似的方法构造得到一个 $k \times k$ 的系数矩阵，其中原方程的系数在第一行，然后从第二行的第一个元素开始，沿着在对角线的元素为 $1$ ，其他位置的元素为 $0$

A= \begin{bmatrix} \ast & \ast & \dots & \ast \\ 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & \dots \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 1 & 0 \end{bmatrix}