线性代数学习笔记7-3：特征值的应用——解微分方程、矩阵的指数函数

之前介绍了求解一阶差分方程，本文介绍求解一阶导数常系数微分方程

常系数微分方程的解是指数形式的$e^{\lambda t}$，基于这个事实，我们得以将问题转为线性代数的问题，求解其指数和系数
另外，在微分方程中将会看到特征值的另一应用：除了可以帮助我们求矩阵的幂，还可以求矩阵的指数

举例：一阶微分方程

对于一阶微分方程$\{\begin{array}{cc}\frac{d{u}_1}{dt}& =-u_1+2u_2\\ \frac{d{u}_2}{dt}& =u_1-2u_2\end{array}$，初始条件为$u_1=1,u_2=0$

类比之前的差分方程，微分方程可以表示为一个方程组（未知量组成一个整体的列向量$\mathbf{u}$，系数为$\mathbf{A}$）$$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$$其中$\mathbf{A}=\left[\begin{array}{cc}-1& 2\\ 1& -2\end{array}\right],\mathbf{u}(0)=\left[\begin{array}{l}1\\ 0\end{array}\right]$
我们的最终目标是找到$u(t)$，或者说追踪$u$随时间的变化及最终的值

[结论] 方程$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$的通解形式为$$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$$

• 其中，特征向量${\lambda }_i$和特征值${\mathbf{x}}_i$来自特征值分解$\mathbf{A}=\mathbf{X}\mathbf{\Lambda }{\mathbf{X}}^{-1}$
  系数$c_i$由初值条件给出（一旦确定后系数恒定不变），本质就是通过特征向量的线性组合表示初始状态的组合系数，即$\mathbf{X}\mathbf{c}=\mathbf{u}(0)$；
• 注意，之前讲的差分方程中，出现纯幂形式$c_1{\lambda }_1^k{\mathbf{x}}_1+c_2{\lambda }_2^k{\mathbf{x}}_2+...$；
• 而微分方程中出现纯指数形式$c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2+...$，两者作用类似
  例如将特解$c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1$代入方程$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$中验证：得到$c_1{\lambda }_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1=\mathbf{A}\mathbf{u}=c_1\mathbf{A}{e}^{{\lambda }_{1}t}{\mathbf{x}}_1$，即${\lambda }_1{\mathbf{x}}_1=\mathbf{A}{\mathbf{x}}_1$，可见，这里纯指数主要是用于保证“微分后形式仍保持不变”

证明：方程的通解$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$ / 探究为何微分方程中出现指数函数
微分方程$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$中的$\frac{d\mathbf{u}}{dt}$，就相当于差分方程中前后两项之差${\mathbf{u}}_{k+1}-{\mathbf{u}}_k$，只不过两项的距离无限小；
或者说，从$u(0)$到$u(t)$认为经历了$N\to \infty$个差分方程，每个方程对应了时长为${\Delta }_t/N$的小时隙
那么，将某个时刻的值$u(t){|}_{t=k{\Delta }_t}$记为$u_k$，则微分方程可以暂时视为$$\frac{d\mathbf{u}}{dt}=\frac{{\mathbf{u}}_{k+1}-{\mathbf{u}}_k}{\Delta t}=\frac{{\mathbf{u}}_{k+1}-{\mathbf{u}}_k}{t/N}=\mathbf{A}{\mathbf{u}}_k$$
改写得到$${\mathbf{u}}_{k+1}=\left(\frac{t}{N}\mathbf{A}+1\right){\mathbf{u}}_k$$这是上一节所学的内容，可以知道最终的值$\mathbf{u}(\mathrm{t})$，等价于差分方程的最终值${\mathbf{u}}_N$$$\mathbf{u}(t)={\mathbf{u}}_N=c_1{\left(\frac{t}{N}{\lambda }_1+1\right)}^N{\mathbf{x}}_1+c_2{\left(\frac{t}{N}{\lambda }_2+1\right)}^N{\mathbf{x}}_2$$，并且考虑到${\lim }_{k\to \infty }{\left(\frac{1}{k}+1\right)}^k\to e$，从差分变为微分，得到其通解形式：$$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$$

求解$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$，其中$\mathbf{A}=\left[\begin{array}{cc}-1& 2\\ 1& -2\end{array}\right],\mathbf{u}(0)=\left[\begin{array}{l}1\\ 0\end{array}\right]$：

• 首先求$\mathbf{A}$的特征值：${\lambda }_1=0,\mathbf{x}1=\left[\begin{array}{l}2\\ 1\end{array}\right]$和${\lambda }_1=-3,\mathbf{x}1=\left[\begin{array}{l}1\\ -1\end{array}\right]$
• 将初始状态表示为特征向量的线性组合，也就是求解$\mathbf{X}\mathbf{c}=\mathbf{u}(0)$：$$\mathbf{u}(0)=\left[\begin{array}{l}1\\ 0\end{array}\right]=c_1\left[\begin{array}{l}2\\ 1\end{array}\right]+c_2\left[\begin{array}{c}1\\ -1\end{array}\right]\text{,其中c }1=\mathrm{c}2=1/3$$
• 根据上述的通解$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$，代入系数$c$和特征值$\lambda$、特征向量$\mathbf{x}$解得$$\begin{array}{l}\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2\\ =c_1{e}^0\left[\begin{array}{l}2\\ 1\end{array}\right]+c_2{e}^{-3t}\left[\begin{array}{c}1\\ -1\end{array}\right]\\ =\frac{1}{3}\left[\begin{array}{l}2\\ 1\end{array}\right]+\frac{1}{3}{e}^{-3t}\left[\begin{array}{c}1\\ -1\end{array}\right]\end{array}$$

分析：随着$t$的增大，答案中第二项会消失，而第一项为稳态（当有0特征值时，就出现稳态SteadyState）’
如果初态$\mathbf{u}(0)=\left[\begin{array}{l}1\\ 0\end{array}\right]$，那么最终的稳态就是$$\begin{array}{l}\mathbf{u}(\mathrm{t}){|}_{t=\infty }=\frac{1}{3}\left[\begin{array}{l}2\\ 1\end{array}\right]\end{array}$$

---

关于“稳态”：

对于通解$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$
矩阵的特征值给出了$\mathbf{u}(\mathrm{t})$的发展趋势（稳定性）：

• 将$e^{{\lambda }_{1}t}$视为$A{e}^{j\varphi }$的形式，则实部 R e { λ } Re\{\lambda\} Re{λ}决定了稳定性（即决定幅值的增长速度，因为$|e^{a+jb}|=|e^a||e^{jb}|=|e^a|$），虚部 I m { λ } Im\{\lambda\} Im{λ}对应了单位圆上的相位旋转

• R e { λ } > 0 Re\{\lambda\}>0 Re{λ}>0，对应项发散； R e { λ } = 0 Re\{\lambda\}=0 Re{λ}=0，对应项幅值稳定不变； R e { λ } < 0 Re\{\lambda\}<0 Re{λ}<0，对应项消失（$t\to \infty 时\mathbf{u}(\mathrm{t})\to 0$）

对比复习：
①之前的矩阵的幂的情况
如果其所有特征值$|{\lambda }_i|<1$，则$k\to \infty 时{\mathbf{A}}^k\to 0$（因为${\mathbf{\Lambda }}^k\to 0$，故${\mathbf{A}}^k=\mathbf{S}{\mathbf{\Lambda }}^k{\mathbf{S}}^{-1}\to 0$）
②之前的差分方程的稳态情况（主要取决于幅值$|{\lambda }_i|$）
对于实数特征值，特征值$|{\lambda }_i|<1$的项最终会消失，特征值$|{\lambda }_i|=1$的项恒定，特征值$|{\lambda }_i|>1$的项最终不断增长
对于复数特征值，虚部引入了复平面上的“旋转”，故特征值的幅值仍然确定稳态，而相位则对应了每次做矩阵乘法时特征向量的旋转角度

• 综合考虑所有特征值对于解的稳态的影响：
  ①若所有 R e { λ } ≤ 0 Re\{\lambda\}\leq0 Re{λ}≤0（对应项消失/幅值稳定不变），则可以进入稳态
  ②一旦存在 R e { λ } > 0 Re\{\lambda\}>0 Re{λ}>0，则发散

另外，我们比较关注二阶系统的稳定性，一个推论是：

• 假如行列式$det(\mathbf{A})>0$而$trace(\mathbf{A})<0$（这等价于说二阶矩阵有两个实部为负的特征值），则微分方程的解可以进入稳态

扩展：从方程解耦（特征值分解）的角度求解微分方程

下文一切讨论的前提：$\mathbf{A}$的特征向量矩阵$\mathbf{S}$可逆/$\mathbf{A}$有n个无关的特征向量，因为此时才能保证可对角化（从而用于解耦的特征向量数量是足够的）
将矩阵对角化（特征值分解）为$\mathbf{A}=\mathbf{S}\mathbf{\Lambda }{\mathbf{S}}^{-1}$

下面从另一角度考虑一般的一阶微分方程的求解原理，探究为什么一阶微分方程的解是指数函数$e^{\lambda t}$的和的形式（之前已经证明，这里从另一角度理解）；

出发角度是：
将问题转换到另一坐标系（基向量为特征向量），从而解耦方程，得到解，再转换回原坐标系；
具体的解耦方法：
对于未知量相互耦合的方程组，用特征值和特征向量来对角化方程的系数矩阵，可以实现解耦（各未知量没有关系）

关于变量的耦合与解耦
对于一阶微分方程$\{\begin{array}{cc}\frac{d{u}_1}{dt}& =-u_1+2u_2\\ \frac{d{u}_2}{dt}& =u_1-2u_2\end{array}$，也就是$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$

• 原方程有两个相互耦合(coupled)的未知函数$u_1,u_2$；
• 找出$\mathbf{A}$的特征值和特征向量（即对角化），可以实现解耦，即将方程组每一行变为可以独立求解的单个方程，他们之间互不相干
  实际上后面会看到，解耦的两个方程的解就对应了 最终通解的两个分量$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$；

---

进一步的，下面研究如何将解表示为特征值矩阵$\mathbf{\Lambda }$和特征向量矩阵$\mathbf{S}$的形式

对于$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$方程，原方程未知数$u_1,u_2$耦合，即$\mathbf{A}$不是对角阵
希望解耦，就是说希望将$\mathbf{A}$对角化$\mathbf{A}=\mathbf{S}\mathbf{\Lambda }{\mathbf{S}}^{-1}$（出现对角的系数矩阵意味着各未知数互不干扰），方法是：

• 将$\mathbf{u}$写为特征向量的线性组合形式$\mathbf{u}=\mathbf{S}\mathbf{v}$
  其中$\mathbf{v}$为新的未知量（代替原来的未知量$\mathbf{u}$），$\mathbf{S}$矩阵是$\mathbf{A}$的特征向量矩阵

• $\mathbf{u}=\mathbf{S}\mathbf{v}$带入原方程$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$，得到$$\begin{array}{l}\mathbf{S}\frac{d\mathbf{v}}{dt}=\mathbf{A}\mathbf{S}\mathbf{v}\\ \Rightarrow \frac{d\mathbf{v}}{dt}={\mathbf{S}}^{-1}\mathbf{A}\mathbf{S}\mathbf{v}=\mathbf{\Lambda }\mathbf{v}\end{array}$$

可以理解为：

• 以特征向量为坐标系的基，求解问题（带入$\mathbf{u}=\mathbf{S}\mathbf{v}$，原未知量$\mathbf{u}$在新的基下的坐标为$\mathbf{v}$，即新的未知量为$\mathbf{v}$）
• 问题变为：求解关于新的未知量$\mathbf{v}$的对角化方程组$\frac{d\mathbf{v}}{dt}=\mathbf{\Lambda }\mathbf{v}$

新方程组的系数矩阵$\mathbf{\Lambda }$为对角阵，即方程组每一行都形如$\frac{d{v}_i}{dt}={\lambda }_i{v}_i$；
此时新方程组不存在耦合，或者说方程组各个未知量之间没有联系（系数矩阵为对角阵导致的）
新的 对角化方程组$\frac{d\mathbf{v}}{dt}=\mathbf{\Lambda }\mathbf{v}$的解为$$\mathbf{v}(t)=e^{\Lambda t}\mathbf{v}(0)$$

• 通过坐标变换/基变换，可以得到原方程$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$的解为$$\mathbf{u}(t)=\mathbf{S}\mathbf{v}(t)=\mathbf{S}{e}^{\Lambda t}\mathbf{v}(0)=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}\mathbf{u}(0)=e^{\mathbf{A}t}\mathbf{u}(0)$$也就是说方程$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$的解就是$\mathbf{u}(t)=e^{\mathbf{A}t}\mathbf{u}(0)$，其中$e^{\mathbf{A}t}=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}$，意义是坐标变换（下面将会证明）

结论：

• （前提：若$\mathbf{A}$有n个无关的特征向量方程）$\frac{d\mathbf{u}}{dt}=\mathbf{A}\mathbf{u}$的解解就是$\mathbf{u}(t)=e^{\mathbf{A}t}\mathbf{u}(0)$，其中$e^{\mathbf{A}t}=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}$，意义是坐标变换（这里的$\mathbf{\Lambda }$和$\mathbf{S}$来自对角化$\mathbf{A}=\mathbf{S}\mathbf{\Lambda }{\mathbf{S}}^{-1}$）
• 上面的通解$\mathbf{u}(\mathrm{t})=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$是这里的一般方法在特征值个数为2时的一个特例，可以带入关系$\mathbf{S}\mathbf{c}=\mathbf{u}(0)$和$\mathbf{S}=[{\mathbf{x}}_1,{\mathbf{x}}_2]$验证，$\mathbf{u}(t)=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}\mathbf{u}(0)=\mathbf{S}{e}^{\mathbf{\Lambda }t}\mathbf{c}=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2$
  这里我们将解表示成了特征值矩阵$\mathbf{\Lambda }$和特征向量矩阵$\mathbf{S}$的形式，这种矩阵形式更加通用

证明：矩阵的指数函数$e^{\mathbf{A}t}=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}$

我们希望研究$e^{\mathbf{A}t}$，这是个矩阵，但是其元素没有显式计算表示，而是需要通过下面的幂级数公式来计算！！！

目标：证明$e^{\mathbf{A}t}$与$\mathbf{A}$的特征值矩阵$\mathbf{\Lambda }$和特征向量矩阵$\mathbf{S}$的关系为$e^{\mathbf{A}t}=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}$

利用指数函数的幂级数公式$e^x={\sum }_{n=0}^{\infty }\frac{x^n}{n!}=1+x+\frac{x^2}{2}+\frac{x^3}{6}+\cdots$，可以将指数部分$e^x$变为幂次项$x^n$之和的形式，类比得到矩阵的指数函数$e^{\mathbf{A}t}$，可以写为$$e^{\mathbf{A}t}=I+\mathbf{A}t+\frac{(\mathbf{A}t{)}^2}{2}+\frac{(\mathbf{A}t{)}^3}{6}+\cdots$$利用${\mathbf{A}}^k=\mathbf{S}{\mathbf{\Lambda }}^k{\mathbf{S}}^{-1}$，可以得到$$\begin{array}{rl}{e}^{\mathbf{A}t}& =I+\mathbf{A}t+\frac{(\mathbf{A}t{)}^2}{2}+\frac{(\mathbf{A}t{)}^3}{6}+\cdots \\ & =\mathbf{S}{\mathbf{S}}^{-1}+\mathbf{S}\mathbf{\Lambda }{\mathbf{S}}^{-1}t+\frac{\mathbf{S}{\mathbf{\Lambda }}^2{\mathbf{S}}^{-1}}{2}{t}^2+\frac{\mathbf{S}{\mathbf{\Lambda }}^3{\mathbf{S}}^{-1}}{6}{t}^3+\cdots \\ & =\mathbf{S}\left(I+\mathbf{\Lambda }t+\frac{{\Lambda }^2}{2}{t}^2+\frac{{\Lambda }^3}{6}{t}^3+\cdots \right){\mathbf{S}}^{-1}\\ & =\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}\end{array}$$
注意，这里出现了${\mathbf{S}}^{-1}$，那么下面的一切成立的前提是，矩阵A必须具有n个线性无关的特征向量，从而矩阵才能对角化

最终，我们可以将$e^{\mathbf{A}t}$分解为$$e^{\mathbf{A}t}=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1},其中e^{\Lambda t}=\left[\begin{array}{cccc}{e}^{{\lambda }_{1}t}& 0& \cdots & 0\\ 0& e^{{\lambda }_{2}t}& & 0\\ ⋮& & \ddots & ⋮\\ 0& \cdots & 0& e^{{\lambda }_{n}t}\end{array}\right]$$

可见，若$\mathbf{A}$有n个无关的特征向量方程，则$e^{\mathbf{A}t}$与$\mathbf{A}$的特征值矩阵$\mathbf{\Lambda }$和特征向量矩阵$\mathbf{S}$的关系为$e^{\mathbf{A}t}=\mathbf{S}{e}^{\mathbf{\Lambda }t}{\mathbf{S}}^{-1}$（给出了$e^{\mathbf{A}t}$的显式计算式）

二阶微分方程

像之前的Fibonacci数列的例子，可以从二阶差分方程构造一阶的差分方程
同样的，给出二阶微分方程（同时出现了前后三项）$y^{\prime \prime }+b{y}^{\prime }+ky=0$，我们也可以增加一个方程，得到一个方程组（可表示为矩阵向量乘法），从而将整体列向量变为新的变量，得到一阶微分方程

构造方程组$$\{\begin{array}{c}{y}^{\prime \prime }=-b{y}^{\prime }-ky\\ y^{\prime }=y^{\prime }\end{array}$$并令新的未知量$\mathbf{u}=\left[\begin{array}{l}{y}^{\prime }\\ y\end{array}\right]$，则得到新方程组$$\left[\begin{array}{l}{y}^{\prime \prime }\\ y^{\prime }\end{array}\right]=\left[\begin{array}{cc}-b& -k\\ 1& 0\end{array}\right]\left[\begin{array}{l}{y}^{\prime }\\ y\end{array}\right]即{\mathbf{u}}^{\prime }=\left[\begin{array}{cc}-b& -k\\ 1& 0\end{array}\right]\mathbf{u}$$

可以推广到k阶微分方程，我们仍然添加(k-1)个额外方程，则${\mathbf{u}}^{\prime }=\mathbf{A}\mathbf{u}$的系数矩阵$\mathbf{A}$是一个k x k矩阵，这样k阶微分方程可以转化为一阶微分方程