6.3 微分方程系统（1）

一、常系数微分方程

特征值与特征向量和 $A=X\Lambda X^{-1}$ 很适合求解 $A^k$，对于微分方程 $d\mathbf{u}/dt=A\mathbf{u}$ 的求解也很合适。本节大部分都是线性代数，但是需要知道微积分的一个知识：$e^{\lambda t}$ 的导数是 $\lambda e^{\lambda t}$。本节的重点是：把常系数微分方程转换成线性代数。
常微分方程 $\frac{du}{dt}=u$ 和 $\frac{du}{dt}=\lambda u$ 的解是指数形式：$$\frac{du}{dt}=u的解是u(t)=C{e}^t\frac{du}{dt}=\lambda u的解是u(t)=C{e}^{\lambda t}(\mathrm{6.3.1})$$在时间 $t=0$ 时这些解都包含 $e^0=1$，所以此时它们都简化成 $u(0)=C$，这个 “初始值（initial value）” 告诉我们 $C$ 是什么。在 $t=0$ 时从数字 $u(0)$ 开始的上式的解是 $u(t)=u(0)e^t$ 和 $u(t)=u(0)e^{\lambda t}$。
我们仅仅求解了 $1\times 1$ 的问题，线性代数会移到 $n\times n$。未知数是一个向量 $\mathbf{u}$（现在加粗了），它由给定的初始向量 $\mathbf{u}(0)$ 开始。这 $n$ 个方程包含一个方阵 $A$，我们希望从 $n$ 个 ${\lambda }^{\prime }s$ 得到 $\mathbf{u}(t)$ 的 $n$ 个指数 $e^{\lambda t}$：

$$n个方程的系统\frac{d\mathbf{u}}{dt}=A\mathbf{u}从t=0时的向量\mathbf{u}(0)=\left[\begin{array}{c}{u}_1(0)\\ \cdots \\ u_n(0)\end{array}\right]开始(\mathrm{6.3.2})$$

这些微分方程是线性的，如果 $\mathbf{u}(t)$ 和 $\mathbf{v}(t)$ 是解，则 $C\mathbf{u}(t)+D\mathbf{v}(t)$ 也是解。我们需要 $n$ 个像 $C$ 和 $D$ 一样的常数来匹配 $\mathbf{u}(0)$ 的分量，第一份工作是使用 $A\mathbf{x}=\lambda \mathbf{x}$ 求出 $n$ 个 “纯指数解（pure exponential solutions）” $\mathbf{u}=e^{\lambda t}\mathbf{x}$。
注意 $A$ 是常数矩阵，在其它的线性方程组中，$A$ 随着 $t$ 改变；非线性方程组中，$A$ 随着 $\mathbf{u}$ 改变。我们这里没有那么困难，$d\mathbf{u}/dt=A\mathbf{u}$ 是 “线性且是常系数”。只有这样的微分方程组我们才能直接转换成线性代数，下面是关键：$$当A\mathbf{x}=\lambda \mathbf{x}时，用指数形式的e^{\lambda t}\mathbf{x}求解线性常系数微分方程。$$

二、du/dt = Au 的解

我们的纯指数解是 $e^{\lambda t}$ 乘上一个固定的向量 $\mathbf{x}$，可以猜到这个 $\lambda$ 就是 $A$ 的特征值，$\mathbf{x}$ 是特征向量。将 $\mathbf{u}(t)=e^{\lambda t}\mathbf{x}$ 代入到方程 $d\mathbf{u}/dt=A\mathbf{u}$ 中可以验证这个猜测是正确的。因子 $e^{\lambda t}$ 会消去最终留下 $\lambda \mathbf{x}=A\mathbf{x}$：

$$\begin{array}{l}当A\mathbf{x}=\lambda \mathbf{x}时，\\ 令\mathbf{u}=e^{\lambda t}\mathbf{x}\end{array}\frac{d\mathbf{u}}{dt}=\lambda e^{\lambda t}\mathbf{x}和A\mathbf{u}=A{e}^{\lambda t}\mathbf{x}是一致的(\mathrm{6.3.3})$$

$\mathbf{u}=e^{\lambda t}\mathbf{x}$ 这个特解的所有分量都有相同的 $e^{\lambda t}$。当 $\lambda >0$ 时，解增长；当 $\lambda <0$ 时解衰减。如果 $\lambda$ 是一个复数，它的实部决定增长或衰减；虚部 $\omega$ 产生振荡，像正弦波一样。

【例1】求解 $\frac{d\mathbf{u}}{dt}=A\mathbf{u}=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\mathbf{u}$，从 $\mathbf{u}(0)=\left[\begin{array}{c}4\\ 2\end{array}\right]$ 开始。
这是一个关于 $\mathbf{u}$ 的向量方程，它含有两个分量分别是 $y$ 和 $z$ 的数量方程。这两个方程 “耦合（coupled together）” 在了一起，因为矩阵 $A$ 不是对角矩阵：$$\frac{d\mathbf{u}}{dt}=A\mathbf{u}\frac{d}{dt}\left[\begin{array}{c}y\\ z\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}y\\ z\end{array}\right]表示\frac{dy}{dt}=z和\frac{dz}{dt}=y$$特征向量的思想是将这些方程以某种方式组合起来，然后回到 $1\times 1$ 的问题。本题可以组合成 $y+z$ 和 $y-z$，将两个方程相加和相减：$$\frac{d}{dt}(y+z)=z+y和\frac{d}{dt}(y-z)=-(y-z)$$$y+z$ 的组合像 $e^t$ 一样增长，因为它的 $\lambda =1$；$y-z$ 的组合像 $e^{-t}$ 一样衰减，因为它的 $\lambda =-1$。重点是：我们没有必要用原始方程 $d\mathbf{u}/dt=A\mathbf{u}$ 来找这些特殊的组合，$A$ 和特征值和特征向量会给我们做好这件事。
矩阵 $A$ 的特征值是 $1$ 和 $-1$，特征向量 $\mathbf{x}$ 是 $(1,1)$ 和 $(1,-1)$。纯指数解 ${\mathbf{u}}_1$ 和 ${\mathbf{u}}_2$ 都是指数 $e^{\lambda t}$ 的形式，其中 ${\lambda }_1=1，{\lambda }_2=-1$：

$$\boxed{{\mathbf{u}}_1(t)=e^{{\lambda }_{1}t}{\mathbf{x}}_1=e^t\left[\begin{array}{c}1\\ 1\end{array}\right]}和\boxed{{\mathbf{u}}_2(t)=e^{{\lambda }_{2}t}{\mathbf{x}}_2=e^{-t}\left[\begin{array}{c}1\\ -1\end{array}\right]}(\mathrm{6.3.4})$$

注意：这些 ${\mathbf{u}}^{\prime }s$ 满足 $A{\mathbf{u}}_1={\mathbf{u}}_1$ 和 $A{\mathbf{u}}_2=-{\mathbf{u}}_2$，就如同 ${\mathbf{x}}_1$ 和 ${\mathbf{x}}_2$ 一样。因子 $e^t$ 和 $e^{-t}$ 随时间 $t$ 改变。这些因子得到 $d{\mathbf{u}}_1/dt={\mathbf{u}}_1=A{\mathbf{u}}_1$ 和 $d{\mathbf{u}}_2/dt=-{\mathbf{u}}_2=A{\mathbf{u}}_2$。对于 $d\mathbf{u}/dt=A\mathbf{u}$，我们得到了两个解。为了找到其它的所有解，用任意数字 $C$ 和 $D$ 乘上这些特解，然后加起来：$$\begin{gathered} 全解\mathbf{u}(t)=C{e}^t\left[\begin{array}{c}1\\ 1\end{array}\right]+D{e}^{-t}\left[\begin{array}{c}1\\ -1\end{array}\right]=\left[\begin{array}{c}C{e}^t+D{e}^{-t \\ }\\ C{e}^t-D{e}^{-t}\end{array}\right](\mathrm{6.3.5}) \end{gathered}$$有了这两个常数 $C$ 和 $D$，我们可以和任意的起始向量 $\mathbf{u}(0)=(u_1(0),u_2(0))$ 相匹配。令 $t=0$，则 $e^0=1$，例 $1$ 要求的初始值是 $\mathbf{u}(0)=(4,2)$：$$\mathbf{u}(0)决定C和DC\left[\begin{array}{c}1\\ 1\end{array}\right]+D\left[\begin{array}{c}1\\ -1\end{array}\right]得到C=3和D=1$$解式（6.2.5）中的 $C=3,D=1$，初始值问题就完全解决了。
解微分方程 $d\mathbf{u}/dt=A\mathbf{u}$ 和求解 ${\mathbf{u}}_{k+1}=A{\mathbf{u}}_k$ 一样有三步：

1. 将 $\mathbf{u}(0)$ 写成 $A$ 特征向量的组合 $c_1{\mathbf{x}}_1+c_2{\mathbf{x}}_2+\cdots +c_n{\mathbf{x}}_n$。
2. 每个特征向量 ${\mathbf{x}}_i$ 左边乘上它的增长因子 $e^{{\lambda }_{i}t}$.
3. 解就是这些纯指数解 $e^{\lambda t}\mathbf{x}$ 的组合，和步骤 $1$ 的组合相同：$$\frac{d\mathbf{u}}{dt}=A\mathbf{u}\mathbf{u}(t)=c_1{e}^{{\lambda }_{1}t}{\mathbf{x}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{x}}_2+\cdots +c_n{e}^{{\lambda }_{n}t}{\mathbf{x}}_n(\mathrm{6.3.6})$$

注：不包含下面的情况：如果有两个 ${\lambda }^{\prime }s$ 相等，且它们只有一个特征向量，那么还需要另一个解（是 $t{e}^{\lambda t}\mathbf{x}$）。步骤 $1$ 需要对角化 $A=X\Lambda X^{-1}$：$n$ 个特征向量构成一组基。

【例2】求解 $d\mathbf{u}/dt=A\mathbf{u}$，已知 $A$ 的特征值是 $\lambda =1,2,3$：$$\begin{array}{l}\mathbf{u}的方程的典型例题\\ 初始条件是\mathbf{u}(0)\end{array}\frac{d\mathbf{u}}{dt}=\left[\begin{array}{ccc}1& 1& 1\\ 0& 2& 1\\ 0& 0& 3\end{array}\right]\mathbf{u}初始值\mathbf{u}(0)=\left[\begin{array}{c}9\\ 7\\ 4\end{array}\right]$$特征向量是 ${\mathbf{x}}_1=(1,0,0)、{\mathbf{x}}_2=(1,1,0)$ 和 ${\mathbf{x}}_3=(1,1,1)$。

步骤1： 向量 $\mathbf{u}(0)=(9,7,4)$ 是 $2{\mathbf{x}}_1+3{\mathbf{x}}_2+4{\mathbf{x}}_3$，因此 $(c_1,c_2,c_3)=(2,3,4)$。
步骤2： 因子 $e^{\lambda t}$ 得到指数解 $e^t{\mathbf{x}}_1、e^{2t}{\mathbf{x}}_2$ 和 $e^{3t}{\mathbf{x}}_3$。
步骤3： 从 $\mathbf{u}(0)$ 开始的组合是 $\mathbf{u}(t)=2e^t{\mathbf{x}}_1+3e^{2t}{\mathbf{x}}_2+4e^{3t}{\mathbf{x}}_3$。
系数 $2,3,4$ 是解线性方程 $c_1{\mathbf{x}}_1+c_2{\mathbf{x}}_2+c_3{\mathbf{x}}_3=\mathbf{u}(0)$ 得出的：$$\left[\begin{array}{ccc}& & \\ {\mathbf{x}}_1& {\mathbf{x}}_2& {\mathbf{x}}_3\\ & & \end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{ccc}1& 1& 1\\ 0& 1& 1\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}2\\ 3\\ 4\end{array}\right]=\left[\begin{array}{c}9\\ 7\\ 4\end{array}\right]这个就是X\mathbf{c}=\mathbf{u}(0)(\mathrm{6.3.7})$$现在已经有了如何求解 $d\mathbf{u}/dt=A\mathbf{u}$ 的基本概念，后面会更深入些，求解的方程会包含二阶导数，因为这些在应用中经常出现。还有是如何决定 $\mathbf{u}(t)$ 是趋近于零还是爆炸或振荡。
后面会有矩阵指数（matrix exponential） $e^{At}$，短公式 $e^{At}\mathbf{u}(0)$ 求解方程 $d\mathbf{u}/dt=A\mathbf{u}$ 和 $A^k{\mathbf{u}}_0$ 求解方程 ${\mathbf{u}}_{k+1}=A{\mathbf{u}}_k$ 方式一样。例 3 会展示 “差分方程” 如何帮助求解微分方程。
这些步骤都使用到了 ${\lambda }^{\prime }s$ 和 ${\mathbf{x}}^{\prime }s$，这节是将求解常系数的问题转换成线性代数的问题，它阐明了最简单但是最重要的微分方程 —— 它们的解都基于增长因子 $e^{\lambda t}$。

三、二阶方程

力学中最重要的方程是 $m{y}^{\prime \prime }+b{y}^{\prime }+ky=0$。 第一项是质量 $m$ 乘上加速度 $a=y^{\prime \prime }$，$ma$ 这一项平衡了力 $F$（牛顿第二定律），这个方程包含了阻尼力 $-b{y}^{\prime }$ 和正比于位移的弹性力 $-ky$。这个方程是一个二阶方程，因为它包含二阶导数 $y^{\prime \prime }=d^{2}y/d{t}^2$，它仍然是一个常系数为 $m,b,k$ 的线性方程。
在微积分中，求解微分方程的方法是将 $y=e^{\lambda t}$ 代入进去，$y$ 的每一阶导数都会提出一个 $\lambda$，我们希望 $y=e^{\lambda t}$ 是这个方程的解：

$$m\frac{d^{2}y}{d{t}^2}+b\frac{dy}{dt}+ky=0变为(m{\lambda }^2+b\lambda +k)e^{\lambda t}=0(\mathrm{6.3.8})$$

一切都与 $m{\lambda }^2+b\lambda +k=0$ 有关，这个关于 $\lambda$ 的方程有两个根 ${\lambda }_1$ 和 ${\lambda }_2$，则关于 $y$ 的方程有两个纯指数解 $y_1=e^{{\lambda }_{1}t}$ 和 $y_2=e^{{\lambda }_{2}t}$，当 ${\lambda }_1\ne {\lambda }_2$ 时，它们的组合 $c_1{y}_1+c_2{y}_2$ 就是全解。
在线性代数中，我们希望有矩阵和特征值，因此我们将数量方程（带有 $y^{\prime \prime }$）转换成一个关于 $y$ 和 $y^{\prime \prime }$ 的向量方程：只含有一阶导数。假设质量 $m=1$，关于 $\mathbf{u}=(y,y^{\prime })$ 的两个方程可以得到 $d\mathbf{u}/dt=A\mathbf{u}$：$$\begin{array}{l}dy/dt=y^{\prime }\\ d{y}^{\prime }/dt=-ky-b{y}^{\prime }\end{array}转换成\frac{d}{dt}\left[\begin{array}{c}y\\ y^{\prime }\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -k& -b\end{array}\right]\left[\begin{array}{c}y\\ y^{\prime }\end{array}\right]=A\mathbf{u}(\mathrm{6.3.9})$$第一个方程 $dy/dt=y^{\prime }$ 是个平凡方程，但是是成立的，第二个方程将 $y^{\prime \prime }、y^{\prime }$ 和 $y$ 联系起来，实际上就是方程（6.3.8）。将它们结合起来，可以将 ${\mathbf{u}}^{\prime }$ 和 $\mathbf{u}$ 联系起来，所以我们可以通过 $A$ 的特征值来求解 ${\mathbf{u}}^{\prime }=A\mathbf{u}$：

$$\boxed{A-\lambda I=\left[\begin{array}{cc}-\lambda & 1\\ -k& -b-\lambda \end{array}\right]}的行列式\boxed{{\lambda }^2+b\lambda +k=0}$$

这个关于 ${\lambda }^{\prime }s$ 的方程和方程（6.3.8）一样！它仍然是 ${\lambda }^2+b\lambda +k=0$，因为 $m=1$。 根 ${\lambda }_1$ 和 ${\lambda }_2$ 现在是 $A$ 的特征值，所有也称为特征根。特征向量和解是：$${\mathbf{x}}_1=\left[\begin{array}{c}1\\ {\lambda }_1\end{array}\right]{\mathbf{x}}_1=\left[\begin{array}{c}1\\ {\lambda }_2\end{array}\right]\mathbf{u}(t)=c_1{e}^{{\lambda }_{1}t}\left[\begin{array}{c}1\\ {\lambda }_1\end{array}\right]+c_2{e}^{{\lambda }_{2}t}\left[\begin{array}{c}1\\ {\lambda }_2\end{array}\right]$$$\mathbf{u}(t)$ 的第一分量是 $y=c_1{e}^{{\lambda }_{1}t}+c_2{e}^{{\lambda }_{2}t}$，和前面的解是一样的，也不可能是其它的解。$\mathbf{u}(t)$ 的第二分量是速度 $dy/dt$。向量问题与数量问题是完全一致的，$2\times 2$ 的 矩阵$A$ 称为友矩阵（companion matrix）—— 一个带有 $y^{\prime \prime }$ 的二阶方程的伴随。

【例3】$绕一个圆运动的方程y^{\prime \prime }+y=0和y=\cos t$
这就是主方程，其中质量 $m=1$，刚度系数 $k=1$，且 $b=0$：没有阻尼 。将 $y=e^{\lambda t}$ 代入 $y^{\prime \prime }+y=0$ 得到 ${\lambda }^2+1=0$，根是 $\lambda =i$ 和 $\lambda =-i$，$e^{it}+e^{-it}$ 的一半得到解 $y=\cos t$。
作为一个一阶系统，初始值 $y(0)=1，y^{\prime }(0)=0$ 进入 $\mathbf{u}(0)=(1,0)$：$$使用y^{\prime \prime }=-y\frac{d\mathbf{u}}{dt}=\frac{d}{dt}\left[\begin{array}{c}y\\ y^{\prime }\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\left[\begin{array}{c}y\\ y^{\prime }\end{array}\right]=A\mathbf{u}(\mathrm{6.3.10})$$$A$ 的特征值和 $\lambda =i$ 和 $\lambda =-i$ 一样，这符合预期，$A$ 是一个反对称（anti-symmetric）矩阵，特征向量是 ${\mathbf{x}}_1=(1,i)$ 和 ${\mathbf{x}}_2=(1,-i)$，适配 $\mathbf{u}(0)=(1,0)$ 的组合是 $\frac{1}{2}({\mathbf{x}}_1+{\mathbf{x}}_2)$，步骤 $2$ 用 $e^{it}$ 和 $e^{-it}$ 分别乘上 ${\mathbf{x}}^{\prime }s$，步骤 $3$ 是组合纯振荡函数进入 $\mathbf{u}(t)$，如期望的得到 $y=\cos t$：$$\mathbf{u}(t)=\frac{1}{2}{e}^{it}\left[\begin{array}{c}1\\ i\end{array}\right]+\frac{1}{2}{e}^{-it}\left[\begin{array}{c}1\\ -i\end{array}\right]=\left[\begin{array}{c}\cos t\\ -\sin t\end{array}\right]这就是\left[\begin{array}{c}y(t)\\ y^{\prime }(t)\end{array}\right]$$向量 $\mathbf{u}=(\cos t,-\sin t)$ 环绕一个圆（Figure 6.3），由于 ${\cos }^{2}t+{\sin }^{2}t=1$，所以这个圆的半径为 $1$。

四、微分方程的注解

常系数线性方程很容易求解，这里只是微分方程的一部分，微分方程还有更多的内容：

1. 二阶方程 $m{u}^{\prime \prime }+b{u}^{\prime }+ku=0$ 在应用上非常重要，解 $u=e^{\lambda t}$ 的指数 $\lambda$ 是 $m{\lambda }^2+b\lambda +k=0$ 的解，阻尼系数 $b$ 很关键：
   欠阻尼（Underdamping）$b^2<4mk$临界阻尼（Critical damping）$b^2=4mk$过阻尼（Overdamping）$b^2>4mk$
   这个决定了 ${\lambda }_1$ 和 ${\lambda }_2$ 是实根、重根还是复数根。如果 ${\lambda }^{\prime }s$ 是复数，则解 $u(t)$ 会边振荡边衰减。
2. 我们的方程没有外力的项 $f(t)$，我们是求 “零空间解”。除了 ${\mathbf{u}}_n(t)$ 我们还需要一个平衡外力 $f(t)$ 的特解 $u_p(t)$：$$\begin{array}{l}时间s时输入f(s)\\ 增长因子e^{A(t-s)}\\ 累加在一起得到时间t时的输出\end{array}{\mathbf{u}}_{\text{particular}}={\int }_0^t{e}^{A(t-s)}f(s)ds$$

在实际应用中，非线性微分方程都是使用数值方法求解，有一个标准且准确的方法龙格-库塔法 “Rung-Kutta”，这个方法是以它的发现者命名的。通过分析可以发现 $du/dt=f(u)$ 的常数解，是 $u(t)=Y$ 且 $f(Y)=0$ 和 $du/dt=0$ 的解：没有移动。我们也能够了解接近 $u=Y$ 是稳定还是不稳定，远离 $Y$ 的，使用计算机来做。

五、差分方程

要将一个圆显示在屏幕上，我们用一个差分方程来代替 $y^{\prime \prime }=-y$，在使用 $Y(t+\Delta t)-2Y(t)+Y(t-\Delta t)$ 除以 $(\Delta t{)}^2$ 来近似 $y^{\prime \prime }$ 时有三种选择：

$$\begin{array}{ll}\text{F}& 从n-1向前\\ \text{C}& 中心在时间n\\ \text{B}& 从n+1向后\end{array}\frac{Y_{n+1}-2Y_n+Y_{n-1}}{(\Delta t{)}^2}=\begin{array}{l}-Y_{n-1}\\ -Y_n\\ -Y_{n+1}\end{array}\begin{array}{c}(\mathrm{6.3.11}\text{F})\\ (\mathrm{6.3.11}\text{C})\\ (\mathrm{6.3.11}\text{B})\end{array}$$

Figure 6.3 显示 $y(t)=\cos t$ 确实在 $t=2\pi$ 时完成一个圆，在时间步长 $\Delta t=2\pi /32$ 时，这三种差分方法在 $32$ 个时间步长并不能完成一个完美的圆。这些图片可以用特征值来解释：$$前向|\lambda |>1(螺旋向外\text{spiralout})中心|\lambda |=1(最优)后向|\lambda |<1(螺旋向内\text{spiralin})$$将 $2$ 步的方程（6.3.11）简化成 $1$ 步的系统 ${\mathbf{U}}_{n+1}=A{\mathbf{U}}_n$，由离散的未知数 ${\mathbf{U}}_n=(Y_n,Z_n)$ 代替 $\mathbf{u}=(y,y^{\prime })$，从 $\mathbf{U}(0)$ 开始，我们取 $n$ 个时间步长 $\Delta t$：

$$前向(\mathrm{6.3.11}\text{F})\boxed{\begin{array}{c}{Y}_{n+1}=Y_n+\Delta t{Z}_n\\ Z_{n+1}=Z_n-\Delta t{Y}_n\end{array}}变为\boxed{{\mathbf{U}}_{n+1}=\left[\begin{array}{cc}1& \Delta t\\ -\Delta t& 1\end{array}\right]\left[\begin{array}{c}{Y}_n\\ Z_n\end{array}\right]=A{\mathbf{U}}_n}(\mathrm{6.3.12})$$

这些和 $Y^{\prime }=Z、Z^{\prime }=-Y$ 很像，它们是包含时间 $n$ 和 $n+1$ 的一阶方程组，将 $Z$ 消去会还原到 “前向” 的二阶方程（6.3.11 $\text{F}$）。
问题：点 $(Y_n,Z_n)$ 是在圆 $Y^2+Z^2=1$ 上吗？答案是不在，它们会如 Figure 6.3 所示的无限长大。我们取幂 $A^n$ 而不是 $e^{At}$，所以我们检验的是 $|\lambda |$ 的大小而不是特征值的实部。

$$A的特征值\lambda =1\pm i\Delta t则|\lambda |>1且(Y_n,Z_n)螺旋向外$$

如果选择后向（6.3.11 $\text{B}$）将会得到相反的结果，如 Figure 6.4，注意新的 $A$：$$后向\begin{array}{c}{Y}_{n+1}=Y_n+\Delta t{Z}_{n+1}\\ Z_{n+1}=Z_n-\Delta t{Y}_{n+1}\end{array}就是\left[\begin{array}{cc}1& -\Delta t\\ \Delta t& 1\end{array}\right]\left[\begin{array}{c}{Y}_{n+1}\\ Z_{n+1}\end{array}\right]=\left[\begin{array}{c}{Y}_n\\ Z_n\end{array}\right]={\mathbf{U}}_n(\mathrm{6.3.13})$$这个矩阵的特征值是 $1\pm i\Delta t$，但是我们要取它的逆从 ${\mathbf{U}}_n$ 得到 ${\mathbf{U}}_{n+1}$，则 $|\lambda |<1$，这个就可以解释为什么后向差分的解螺旋向内到 $(0,0)。$
如果选择中心的方程，可以得到下面的差分方程：$$中心\begin{array}{l}{Y}_{n+1}=Y_n+\Delta t{Z}_n\\ Z_{n+1}=Z_n-\Delta t{Y}_{n+1}\end{array}即\left[\begin{array}{cc}1& 0\\ \Delta t& 1\end{array}\right]\left[\begin{array}{c}{Y}_{n+1}\\ Z_{n+1}\end{array}\right]=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right]\left[\begin{array}{c}{Y}_n\\ Z_n\end{array}\right]$$将上式左边的矩阵取逆，转换成 ${\mathbf{U}}_{n+1}=A{\mathbf{U}}_n$ 的模式：$$\left[\begin{array}{c}{Y}_{n+1}\\ Z_{n+1}\end{array}\right]=\left[\begin{array}{cc}1& \Delta t\\ -\Delta t& -(\Delta t{)}^2+1\end{array}\right]\left[\begin{array}{c}{Y}_n\\ Z_n\end{array}\right]=A{\mathbf{U}}_n$$Figure 6.4 的右侧是 $32$ 步的中心选择，如果 $\Delta t<2$ 时解很接近圆，这个是经常使用的蛙跳法（leapfrog method），第二差分 $Y_{n+1}-2Y_n+Y_{n-1}$ 跳过了式（6.3.11）的中心值 $Y_n$。

这是化学家研究分子运动的方法（分子动力学有海量的计算）；计算机科学一直生机勃发是因为一个微分方程可以被很多差分方程所代替 —— 有些不稳定，有些稳定，有些中性。
实际的工程和实际的物理处理的都是系统（不仅仅是某一点的单个质量），未知数 $\mathbf{y}$ 是一个向量，${\mathbf{y}}^{\prime \prime }$ 的系数是 $n$ 个质量的质量矩阵 $M$，$\mathbf{y}$ 的系数是一个刚度矩阵 $K$，不是一个数字 $k$，${\mathbf{y}}^{\prime }$ 的系数一个有可能是零的阻尼矩阵。
向量方程 $M{y}^{\prime \prime }+ky=\mathbf{f}$ 是计算力学的主要部分，它由 $K\mathbf{x}=\lambda M\mathbf{x}$ 中 $M^{-1}K$ 的特征值控制。

六、2 × 2 矩阵的稳定性

$d\mathbf{u}/dt=A\mathbf{u}$ 是一个基础问题，当 $t\to \infty$ 时，它的解趋近于 $\mathbf{u}=0$ 吗？这个问题通过消耗能量稳定吗？包含 $e^t$ 的解是不稳定的，解的稳定性与 $A$ 的特征值有关。
全解 $\mathbf{u}(t)$ 是由纯解 $e^{\lambda t}\mathbf{x}$ 得到的，如果特征值 $\lambda$ 是实数，则我们可以很明确的得知 $e^{\lambda t}$ 趋近于零的条件：数字 $\lambda$ 要为负数。如果特征值是复数 $\lambda =r+is$，实部的 $r$ 一定要是负数，此时 $e^{\lambda t}$ 可分解成 $e^{rt}{e}^{ist}$，因子 $e^{ist}$ 的绝对值固定为 $1$：$$e^{ist}=\cos st+i\sin st有|e^{ist}{|}^2={\cos }^{2}st+{\sin }^{2}st=1$$$\lambda$ 的实部控制着增长（ $r>0$）或衰减（$r<0$）。
问题是：哪些矩阵都是负特征值？ 更准确的说，何时 ${\lambda }^{\prime }s$ 的实部都是负数？$2\times 2$ 的矩阵这个答案很清楚：

稳定 \kern 6pt当所有的特征值 $\lambda$ 都有负实部时，$A$ 稳定且 $\mathbf{u}(t)\to 0$。则 $2\times 2$ 的矩阵 $A=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$ 一定满足下面两个条件：$$\begin{array}{ll}{\lambda }_1+{\lambda }_2<0& 迹T=a+d一定是负数\\ {\lambda }_1{\lambda }_2>0& 行列式D=ad-bc一定时正数\end{array}$$

原因： 如果 ${\lambda }^{\prime }s$ 都是实数且都为负数，则它们的和一定是负数，这个就是迹 $T$；它们的积是正数，这是行列式 $D$。这个命题反向也成立，如果 $D={\lambda }_1{\lambda }_2$ 是正数，则 ${\lambda }_1$ 和 ${\lambda }_2$ 同号；如果 $T={\lambda }_1+{\lambda }_2$ 是负数，则符号一定是负。我们可以验证 $T$ 和 $D$。
如果 ${\lambda }^{\prime }s$ 都是复数，它们的形式一定是 $r+is$ 和 $r-is$，否则 $T$ 和 $D$ 将不会是实数（ $A$ 是实矩阵）。行列式 $D$ 直接满足正数的条件，因为 $(r+is)(r-is)=r^2+s^2$；迹 $T$ 是 $r+is+r-is=2r$，所以迹 $T$ 是负数，这说明实部 $r$ 是负数且矩阵是稳定的。证毕！
Figure 6.5 是抛物线 $T^2=4D$ 将实数 ${\lambda }^{\prime }s$ 和复数 ${\lambda }^{\prime }s$ 分割开来，因为 ${\lambda }^2-T\lambda +D=0$ 的解包含有平方根 $\sqrt{T^2-4D}$，抛物线下面是实数，上面是虚数。稳定的区域是图形的左上方的四分之一区域 —— 这里迹 $T$ 是负数行列式 $D$ 是正数。