[线性代数]-23 微分方程和$e^{At}$


《线性代数》学习笔记

本节介绍一阶线性常微分方程的矩阵解法,也就是将微分方程用矩阵抽象,通过“解耦”,计算出对应系数,最终得到解。这里会牵涉到𝑒𝐴𝑥计算问题, (A是矩阵),所以也会引出幂指数是矩阵时算式的计算问题。最后扩展介绍了高阶微分方程的降阶求解方法。


1. 微分方程$\frac{\mathrm{d}u}{\mathrm{d}t}=Au$

本讲主要讲解解一阶方程(first-order system)一阶倒数(first derivative)常系数(constant coefficient)线性方程,上一讲介绍了如何计算矩阵的幂,本讲将进一步涉及矩阵的指数形式。

我们通过解一个例子来详细介绍计算方法。
有方程组$\begin{cases}\frac{\mathrm{d}u_1}{\mathrm{d}t}&=-u_1+2u_2\\\frac{\mathrm{d}u_2}{\mathrm{d}t}&=u_1-2u_2\end{cases}$,则系数矩阵是$A=\begin{bmatrix}-1&2\\1&-2\end{bmatrix}$,设初始条件为在$0$时刻$u(0)=\begin{bmatrix}u_1\\u_2\end{bmatrix}=\begin{bmatrix}1\\0\end{bmatrix}$。

  • 这个初始条件的意义可以看做在开始时一切都在$u_1$中,但随着时间的推移,将有$\frac{\mathrm{d}u_2}{\mathrm{d}t}>0$,因为$u_1$项初始为正,$u_1$中的事物会流向$u_2$。随着时间的发展我们可以追踪流动的变化;

  • 根据上一讲所学的知识,我们知道第一步需要找到特征值与特征向量。$A=\begin{bmatrix}-1&2\\1&-2\end{bmatrix}$,很明显这是一个奇异矩阵,所以第一个特征值是$\lambda_1=0$,另一个特征向量可以从迹得到$tr(A)=-3$。当然我们也可以用一般方法计算$\left|A-\lambda I\right|=\begin{vmatrix}-1-\lambda&2\\1&-2-\lambda\end{vmatrix}=\lambda^2+3\lambda=0$。
(教授提前剧透,特征值$\lambda_2=-3$将会逐渐消失,因为答案中将会有一项为$e^{-3t}$,该项会随着时间的推移趋近于$0$。答案的另一部分将有一项为$e^{0t}$,该项是一个常数,其值为$1$,并不随时间而改变。通常含有$0$特征值的矩阵会随着时间的推移达到稳态);

  • 求特征向量,$\lambda_1=0$时,即求$A$的零空间,很明显$x_1=\begin{bmatrix}2\\1\end{bmatrix}$;$\lambda_2=-3$时,求$A+3I$的零空间,$\begin{bmatrix}2&2\\1&1\end{bmatrix}$的零空间为$x_2=\begin{bmatrix}1\\-1\end{bmatrix}$;

  • 则方程组的通解为:$u(t)=c_1e^{\lambda_1t}x_1+c_2e^{\lambda_2t}x_2$,通解的前后两部分都是该方程组的纯解,即方程组的通解就是两个与特征值、特征向量相关的纯解的线性组合。验证一下,比如取$u=e^{\lambda_1t}x_1$带入$\frac{\mathrm{d}u}{\mathrm{d}t}=Au$,对时间求导得到$\lambda_1e^{\lambda_1t}x_1=Ae^{\lambda_1t}x_1$,化简得$\lambda_1x_1=Ax_1$。
对比上一讲,解$u_{k+1}=Au_k$时得到$u_k=c_1\lambda^kx_1+c_2\lambda^kx_2$,而解$\frac{\mathrm{d}u}{\mathrm{d}t}=Au$我们得到$u(t)=c_1e^{\lambda_1t}x_1+c_2e^{\lambda_2t}x_2$;

  • 继续求$c_1,c_2$,$u(t)=c_1\cdot 1\cdot\begin{bmatrix}2\\1\end{bmatrix}+c_2\cdot e^{-3t}\cdot\begin{bmatrix}1\\-1\end{bmatrix}$,已知$t=0$时,$\begin{bmatrix}1\\0\end{bmatrix}=c_1\begin{bmatrix}2\\1\end{bmatrix}+c_2\begin{bmatrix}1\\-1\end{bmatrix}$($Sc=u(0)$),所以$c_1=\frac{1}{3}, c_2=\frac{1}{3}$;

  • 最终结果,$u(t)=\frac{1}{3}\begin{bmatrix}2\\1\end{bmatrix}+\frac{1}{3}e^{-3t}\begin{bmatrix}1\\-1\end{bmatrix}$。


稳定性:这个流动过程从$u(0)=\begin{bmatrix}1\\0\end{bmatrix}$开始,初始值$1$的一部分流入初始值$0$中,经过无限的时间最终达到稳态$u(\infty)=\begin{bmatrix}\frac{2}{3}\\\frac{1}{3}\end{bmatrix}$。所以,要使得$u(t)\to 0$,则需要负的特征值。但如果特征值为复数呢?如$\lambda=-3+6i$,我们来计算$\left|e^{(-3+6i)t}\right|$,其中的$\left|e^{6it}\right|$部分为$\left|\cos 6t+i\sin 6t\right|=1$,因为这部分的模为$\cos^2\alpha+\sin^2\alpha=1$,这个虚部就在单位圆上转悠。所以只有实数部分才是重要的。所以我们可以把前面的结论改为需要实部为负数的特征值。实部会决定最终结果趋近于$0$或$\infty$,虚部不过是一些小杂音。

收敛态:需要其中一个特征值实部为$0$,而其他特征值的实部皆小于$0$。

发散态:如果某个特征值实部大于$0$。上面的例子中,如果将$A$变为$-A$,特征值也会变号,结果发散。


再进一步,我们想知道如何从直接判断任意二阶矩阵的特征值是否均小于零。对于二阶矩阵$A=\begin{bmatrix}a&b\\c&d\end{bmatrix}$,矩阵的迹为$a+d=\lambda_1+\lambda_2$,如果矩阵稳定,则迹应为负数;

但是这个条件还不够,有反例迹小于$0$依然发散:$\begin{bmatrix}-2&0\\0&1\end{bmatrix}$,迹为$-1$但是仍然发散。还需要加上一个条件,因为$\det A=\lambda_1\cdot\lambda_2$,所以还需要行列式为正数。


总结:原方程组有两个相互耦合的未知函数,$u_1, u_2$相互耦合,而特征值和特征向量的作则就是解耦,也就是对角化(diagonalize)。回到原方程组$\frac{\mathrm{d}u}{\mathrm{d}t}=Au$,将$u$表示为特征向量的线性组合$u=Sv$,代入原方程有$S\frac{\mathrm{d}v}{\mathrm{d}t}=ASv$,两边同乘以$S^{-1}$得$\frac{\mathrm{d}v}{\mathrm{d}t}=S^{-1}ASv=\Lambda v$。以特征向量为基,将$u$表示为$Sv$,得到关于$v$的对角化方程组,新方程组不存在耦合,此时$\begin{cases}\frac{\mathrm{d}v_1}{\mathrm{d}t}&=\lambda_1v_1\\\frac{\mathrm{d}v_2}{\mathrm{d}t}&=\lambda_2v_2\\\vdots&\vdots\\\frac{\mathrm{d}v_n}{\mathrm{d}t}&=\lambda_nv_n\end{cases}$,这是一个各未知函数间没有联系的方程组,它们的解的一般形式为$v(t)=e^{\Lambda t}v(0)$,则原方程组的解的一般形式为$u(t)=e^{At}u(0)=Se^{\Lambda t}S^{-1}u(0)$。这里引入了指数部分为矩阵的形式。


2. 指数矩阵$e^{At}$

在上面的结论中,我们见到了$e^{At}$。这种指数部分带有矩阵的情况称为指数矩阵(exponential matrix)。

理解指数矩阵的关键在于,将指数形式展开称为幂基数形式,就像$e^x=1+\frac{x^2}{2}+\frac{x^3}{6}+\cdots$一样,将$e^{At}$展开成幂级数的形式为:
$$e^{At}=I+At+\frac{(At)^2}{2}+\frac{(At)^3}{6}+\cdots+\frac{(At)^n}{n!}+\cdots$$

有两个极具美感的泰勒级数:$e^x=\sum \frac{x^n}{n!}$与$\frac{1}{1-x}=\sum x^n$,如果把第二个泰勒级数写成指数矩阵形式,有$(I-At)^{-1}=I+At+(At)^2+(At)^3+\cdots$,这个式子在$t$非常小的时候,后面的高次项近似等于零,所以可以用来近似$I-At$的逆矩阵,通常近似为$I+At$,当然也可以再加几项。第一个级数对我们而言比第二个级数好,因为第一个级数总会收敛于某个值,所以$e^x$总会有意义,而第二个级数需要$A$特征值的绝对值小于$1$(因为涉及矩阵的幂运算)。我们看到这些泰勒级数的公式对矩阵同样适用。


回到正题,我们需要证明$Se^{\Lambda t}S^{-1}=e^{At}$,继续使用泰勒级数:
$$ e^{At}=I+At+\frac{(At)^2}{2}+\frac{(At)^3}{6}+\cdots+\frac{(At)^n}{n!}+\cdots\\ e^{At}$$

$$=SS^{-1}+S\Lambda S^{-1}t+\frac{S\Lambda^2S^{-1}}{2}t^2+\frac{S\Lambda^3S^{-1}}{6}t^3+\cdots+\frac{S\Lambda^nS^{-1}}{n!}t^n+\cdots\\ e^{At}$$

$$=S\left(I+\Lambda t+\frac{\Lambda^2t^2}{2}+\frac{\Lambda^3t^3}{3}+\cdots+\frac{\Lambda^nt^n}{n}+\cdots\right)S^{-1}\\ e^{At}=Se^{\Lambda t}S^{-1} $$

需要注意的是,$e^{At}$的泰勒级数展开是恒成立的,但我们推出的版本却需要矩阵可对角化这个前提条件。

最后,我们来看看什么是$e^{\Lambda t}$,我们将$e^{At}$变为对角矩阵就是因为对角矩阵简单、没有耦合

$e^{\Lambda t}=\begin{bmatrix}e^{\lambda_1t}&0&\cdots&0\\0&e^{\lambda_2t}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&e^{\lambda_nt}\end{bmatrix}$

有了$u(t)=Se^{\Lambda t}S^{-1}u(0)$,再来看矩阵的稳定性可知,所有特征值的实部均为负数时矩阵收敛,此时对角线上的指数收敛为$0$。

如果我们画出复平面,则要使微分方程存在稳定解,则特征值存在于复平面的左侧(即实部为负);要使矩阵的幂收敛于$0$,则特征值存在于单位圆内部(即模小于$1$),这是幂稳定区域。(上一讲的差分方程需要计算矩阵的幂); 

同差分方程一样,我们来看二阶情况如何计算,有$y''+by'+k=0$。我们也模仿差分方程的情形,构造方程组$\begin{cases}y''&=-by'-ky\\y'&=y'\end{cases}$,写成矩阵形式有$\begin{bmatrix}y''\\y'\end{bmatrix}=\begin{bmatrix}-b&-k\\1&0\end{bmatrix}\begin{bmatrix}y'\\y\end{bmatrix}$,令$u'=\begin{bmatrix}y''\\y'\end{bmatrix}, \ u=\begin{bmatrix}y'\\y\end{bmatrix}$。

继续推广,对于$5$阶微分方程$y'''''+by''''+cy'''+dy''+ey'+f=0$,则可以写作

$\begin{bmatrix}y'''''\\y''''\\y'''\\y''\\y'\end{bmatrix}=\begin{bmatrix}-b&-c&-d&-e&-f\\1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\end{bmatrix}\begin{bmatrix}y''''\\y'''\\y''\\y'\\y\end{bmatrix}$

这样就把一个五阶微分方程化为$5\times 5$一阶方程组了,然后就是求特征值、特征向量了步骤了。




在实际情况下使用矩阵对角化,特征值等方法求解微分方程,给出了一种使用矩阵求解微分方程的通用规律,即高阶降阶,一阶用特征值和特征向量将原系数矩阵 A解耦,最后得到结果。并介绍了在我们解耦A时使用矩阵对角化将其与特征向量联系起来运算的方法。另外介绍了判断收敛性的方法,即看特征值实部绝对值与 1的大小关系。

2024-08-06

声明:ZY|版权所有,违者必究|如未注明,均为原创|本网站采用BY-NC-SA协议进行授权

转载:转载请注明原文链接 - [线性代数]-23 微分方程和$e^{At}$