上一讲我们主要讲了差分方程 (difference equation) 和矩阵的幂 (powers of matrix) 之间的联系。主要的 insight 是把差分方程的每次递归, i.e., 从 { a k ,   a k − 1 ,   . . . } \{a_k,~a_{k-1},~...\} {ak​, ak−1​, ...} 到 a k + 1 a_{k+1} ak+1​, 表示成一种矩阵关系。这样一来，若干次递归之后的结果相当于初始条件乘以一个矩阵的 K K K 次方。而差分方程最后的解是否 stable 取决于矩阵的特征值是否在复平面的单元圆中。

这一讲 Prof. Strang 着重讲矩阵和其特征值的在微分方程 (differential equations) 中的应用。

矩阵的指数函数 e A e^{\bm{A}} eA

学过高等数学都知道泰勒展开，比如
1 1 − x = 1 + x + x 2 + x 3 + . . . ,     ( if  ∣ x ∣ < 1 ) \frac{1}{1-x} = 1+x+x^2+x^3+...,~~~(\text{if}~|x|<1) 1−x1​=1+x+x2+x3+...,   (if ∣x∣<1)

e x = 1 + x + 1 2 ! x 2 + 1 3 ! x 3 + . . . e^x = 1+x+\frac{1}{2!} x^2+\frac{1}{3!}x^3+... ex=1+x+2!1​x2+3!1​x3+...

同样的，矩阵也可以进行同样的操作

( I − A ) − 1 = 1 + A + A 2 + A 3 + . . . ,     ( if  A ∞ → 0 ) (\bm{I}-\bm{A})^{-1}=1+\bm{A}+\bm{A}^2+\bm{A}^3+...,~~~(\text{if}~\bm{A}^\infty\to 0) (I−A)−1=1+A+A2+A3+...,   (if A∞→0)

e A = I + A + 1 2 ! A 2 + 1 3 ! A 3 + . . .        ( 1 ) e^\bm{A} = \bm{I}+\bm{A}+\frac{1}{2!} \bm{A}^2+\frac{1}{3!}\bm{A}^3+...~~~~~~(1) eA=I+A+2!1​A2+3!1​A3+...      (1)

这样就引出了矩阵的指数函数的定义。特别的，当 A \bm{A} A 有 n n n 个线性无关的特征向量时，我们可以把 A \bm{A} A 相似对角化
A = S Λ S − 1 \bm{A=S\Lambda S^{-1}} A=SΛS−1

e A = I + S Λ S − 1 + 1 2 ! S Λ 2 S − 1 + 1 3 ! S Λ 3 S − 1 + . . . = S e Λ S − 1      ( 2 ) e^\bm{A} = \bm{I}+\bm{S\Lambda S^{-1}}+\frac{1}{2!} \bm{S\Lambda^2 S^{-1}}+\frac{1}{3!}\bm{S\Lambda^3 S^{-1}}+...=\bm{Se^\Lambda S^{-1}}~~~~(2) eA=I+SΛS−1+2!1​SΛ2S−1+3!1​SΛ3S−1+...=SeΛS−1    (2)

看到没，只要 A \bm{A} A 可以对角化，我们便能把 e A e^\bm{A} eA 转换为 S e Λ S − 1 \bm{Se^\Lambda S^{-1}} SeΛS−1. 这样做有什么好处尼？当然是 e Λ e^\Lambda eΛ 比较简单啦，请看
e Λ = I + Λ + 1 2 ! Λ 2 + 1 3 ! Λ 3 + . . . = [ e λ 1 0 0 ⋯ 0 e λ 2 0 ⋯ 0 0 e λ 3 ⋯ ⋮ ⋮ ⋮ ⋱ ]      ( 3 ) e^\Lambda=\bm{I}+\Lambda+\frac{1}{2!} \Lambda^2+\frac{1}{3!}\Lambda^3+... =\begin{bmatrix} e^{\lambda_1} & 0 & 0 & \cdots \\ 0 & e^{\lambda_2} & 0 & \cdots \\ 0 & 0 & e^{\lambda_3} & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix}~~~~(3) eΛ=I+Λ+2!1​Λ2+3!1​Λ3+...=⎣⎢⎢⎢⎡​eλ1​00⋮​0eλ2​0⋮​00eλ3​⋮​⋯⋯⋯⋱​⎦⎥⎥⎥⎤​    (3)

也就是说，对角阵的指数函数直接就是各个对角线元素的分别求指数后得到的矩阵。

differential equation

好，我们先来看看这样一个一阶微分方程
d   u d   t = A u        ( 4 ) \frac{d\,\bm{u}}{d\,t}=\bm{Au}~~~~~~(4) dtdu​=Au      (4)

其中 A \bm{A} A 可以相似对角化； u \bm{u} u 是一个列向量，它实际上有 n n n 个随时间变化的entries，但是每个entry的变化率都被 matrix A \bm{A} A 搅在了一起，相互关联。

为了更具体些，我们可以看一个例子，
A = [ − 1 2 1 − 2 ] \bm{A}=\begin{bmatrix} -1 & 2 \\ 1 & -2 \\ \end{bmatrix} A=[−11​2−2​]

展开微分方程
{ d   u 1 d   t = − u 1 + 2 u 2 d   u 2 d   t = u 1 − 2 u 2 \begin{cases} \frac{d\,{u}_1}{d\,t}=-{u_1} + 2{u_2}\\ \frac{d\,{u}_2}{d\,t}={u_1} - 2 {u_2} \\ \end{cases} {dtdu1​​=−u1​+2u2​dtdu2​​=u1​−2u2​​

可以看出这是一个线性系统，有两个变量 u 1 {u}_1 u1​ 和 u 2 u_2 u2​ 且两者的变化率相互 couple 在一起。稍后，我们将解出这个例子，但是在此之前，我们先看看我们解微分方程的思路。

---

Train of thought: 为了解 (4)，我们可以尝试把 u \bm{u} u 用 A \bm{A} A 一组线性无关的特征向量构成的矩阵表示出来，i.e.,
u = S v \bm{u}=\bm{Sv} u=Sv

因此 (4) 是可以重新写为
S d   v d   t = A S v \bm{S}\frac{d\,\bm{v}}{d\,t}=\bm{ASv} Sdtdv​=ASv

d   v d   t = S − 1 A S v = Λ v      ( 5 ) \frac{d\,\bm{v}}{d\,t}=\bm{S}^{-1}\bm{ASv}=\bm{\Lambda v}~~~~(5) dtdv​=S−1ASv=Λv    (5)

可以看出, (4) 式被顺利地转换为 (5) 式，其中，(5) 式的系数矩阵是对角阵 – 这就意味着两个变量的变化被分割开了 – 微分方程因此可解。

---

Solution: 对于每一个decoupled 方程
d   v i ( t ) d   t = λ i v i ( t ) \frac{d\,v_i(t)}{d\,t}=\lambda_i v_i(t) dtdvi​(t)​=λi​vi​(t)

我们可以解出
v i ( t ) = c e λ i t v_i(t)=ce^{\lambda_i t} vi​(t)=ceλi​t

其中常数 c = v i ( 0 ) c=v_i(0) c=vi​(0) 由初始条件
v ( 0 ) = S − 1 u ( 0 ) \bm{v}(0)=\bm{S}^{-1}\bm{u}(0) v(0)=S−1u(0)

各项给出。

To summarize，
v ( t ) = [ v 1 ( 0 ) e λ 1 t v 2 ( 0 ) e λ 2 t ⋯ v n ( 0 ) e λ n t ] = [ e λ 1 t 0 0 ⋯ 0 e λ 2 t 0 ⋯ 0 0 e λ 3 t ⋯ ⋮ ⋮ ⋮ ⋱ ] [ v 1 ( 0 ) v 2 ( 0 ) ⋯ v n ( 0 ) ] = e Λ t v ( 0 ) \bm{v}(t)=\begin{bmatrix} v_1(0)e^{\lambda_1 t} \\ v_2(0)e^{\lambda_2 t} \\ \cdots \\ v_n(0)e^{\lambda_n t} \\ \end{bmatrix}=\begin{bmatrix} e^{\lambda_1 t} & 0 & 0 & \cdots \\ 0 & e^{\lambda_2 t} & 0 & \cdots \\ 0 & 0 & e^{\lambda_3 t} & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} \begin{bmatrix} v_1(0) \\ v_2(0) \\ \cdots \\ v_n(0) \\ \end{bmatrix}=e^{\Lambda t} \bm{v}(0) v(t)=⎣⎢⎢⎡​v1​(0)eλ1​tv2​(0)eλ2​t⋯vn​(0)eλn​t​⎦⎥⎥⎤​=⎣⎢⎢⎢⎡​eλ1​t00⋮​0eλ2​t0⋮​00eλ3​t⋮​⋯⋯⋯⋱​⎦⎥⎥⎥⎤​⎣⎢⎢⎡​v1​(0)v2​(0)⋯vn​(0)​⎦⎥⎥⎤​=eΛtv(0)

u = S v = S e Λ t S − 1 u ( 0 ) = e A t u ( 0 )      ( 6 ) \bm{u}=\bm{Sv}=\bm{S}e^{\Lambda t} \bm{S}^{-1}\bm{u}(0)=e^{\bm{A}t}\bm{u}(0)~~~~(6) u=Sv=SeΛtS−1u(0)=eAtu(0)    (6)

这便是微分方程的解。

---

好，让我们回到本节开始时候的例子
d   u d   t = A u \frac{d\,\bm{u}}{d\,t}=\bm{Au} dtdu​=Au

A = [ − 1 2 1 − 2 ] ,    u ( 0 ) = [ 1 0 ] \bm{A}=\begin{bmatrix} -1 & 2 \\ 1 & -2 \\ \end{bmatrix},~~ \bm{u}(0)=\begin{bmatrix} 1 \\ 0 \\ \end{bmatrix} A=[−11​2−2​],  u(0)=[10​]

因为我们已经知道了微分方程的解由 (6) 式给出，因此我们应该分析 A \bm{A} A 的对角化矩阵 S \bm{S} S. 首先， A \bm{A} A 不满秩，因此一定有一个 λ 1 = 0 \lambda_1=0 λ1​=0, 由矩阵的trace知道另一个 λ 2 = − 3 \lambda_2=-3 λ2​=−3. 我们可以进一步把他们对应的特征向量求出
λ 1 = 0 ,    x 1 = [ 2 , 1 ] ⊤ \lambda_1=0,~~\bm{x_1}=[2,1]^\top λ1​=0,  x1​=[2,1]⊤

λ 2 = − 3 ,    x 2 = [ 1 , − 1 ] ⊤ \lambda_2=-3,~~\bm{x_2}=[ 1,-1]^\top λ2​=−3,  x2​=[1,−1]⊤

因此 A \bm{A} A 可以被相似对角化为
A = S Λ S − 1 = [ 2 1 1 − 1 ] [ 0 0 0 − 3 ] [ 1 / 3 − 1 − 1 − 2 / 3 ] \bm{A}=\bm{S\Lambda S^{-1}}=\begin{bmatrix} 2 & 1\\ 1 & -1\\ \end{bmatrix} \begin{bmatrix} 0 & 0\\ 0 & -3\\ \end{bmatrix} \begin{bmatrix} 1/3 & -1\\ -1 & -2/3\\ \end{bmatrix} A=SΛS−1=[21​1−1​][00​0−3​][1/3−1​−1−2/3​]

代入 (6), 我们有
u = S e Λ t S − 1 u ( 0 ) = S Λ S − 1 = [ 2 1 1 − 1 ] [ e 0 t 0 0 e − 3 t ] [ 1 / 3 1 / 3 1 / 3 − 2 / 3 ] [ 1 0 ] \bm{u}=\bm{S}e^{\Lambda t} \bm{S}^{-1}\bm{u}(0)=\bm{S\Lambda S^{-1}}=\begin{bmatrix} 2 & 1\\ 1 & -1\\ \end{bmatrix} \begin{bmatrix} e^{0t} & 0\\ 0 & e^{-3t}\\ \end{bmatrix} \begin{bmatrix} 1/3 & 1/3\\ 1/3 & -2/3\\ \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\ \end{bmatrix} u=SeΛtS−1u(0)=SΛS−1=[21​1−1​][e0t0​0e−3t​][1/31/3​1/3−2/3​][10​]

= 1 3 [ 2 1 1 − 1 ] [ 1 0 0 e − 3 t ] [ 1 1 ] = 1 3 [ 2 1 1 − 1 ] [ 1 e − 3 t ] =\frac{1}{3}\begin{bmatrix} 2 & 1\\ 1 & -1\\ \end{bmatrix} \begin{bmatrix} 1 & 0\\ 0 & e^{-3t}\\ \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ \end{bmatrix}=\frac{1}{3}\begin{bmatrix} 2 & 1\\ 1 & -1\\ \end{bmatrix} \begin{bmatrix} 1\\ e^{-3t}\\ \end{bmatrix} =31​[21​1−1​][10​0e−3t​][11​]=31​[21​1−1​][1e−3t​]

= 1 3 [ 2 1 ] + 1 3 e − 3 t [ 1 − 1 ] =\frac{1}{3}\begin{bmatrix} 2 \\ 1 \\ \end{bmatrix}+\frac{1}{3}e^{-3t}\begin{bmatrix} 1\\ -1\\ \end{bmatrix} =31​[21​]+31​e−3t[1−1​]

---

对于上面这个例子，Prof. Strang是直接给出一个通解 (两个特解的linear combination)
u = e A t u ( 0 ) = c 1 e λ 1 t x 1 + c 2 e λ 2 t x 2       ( 7 ) \bm{u}=e^{\bm{A} t} \bm{u}(0)=c_1e^{\lambda_1 t}\bm{x_1}+c_2e^{\lambda_2 t}\bm{x_2}~~~~~(7) u=eAtu(0)=c1​eλ1​tx1​+c2​eλ2​tx2​     (7)

其中常数 c 1 c_1 c1​ 和 c 2 c_2 c2​ 由初始条件 u ( 0 ) \bm{u}(0) u(0) 给出。

他主要想对比上节课给出的difference equation的解
F k = c 1 λ 1 k v 1 + c 2 λ 2 k v 2 F_k= c_1\lambda_1^k\bm{v_1}+c_2\lambda_2^k\bm{v_2} Fk​=c1​λ1k​v1​+c2​λ2k​v2​

A k → 0 \bm{A}^k\to 0 Ak→0 的条件是特征值在单位圆内，而 e A t → 0 e^{\bm{A t}}\to 0 eAt→0 的条件是 A \bm{A} A 所有特征值的实部都小于 0 (如下图所示)。 从 (7) 可以看出， λ 1 = 0 \lambda_1=0 λ1​=0 因此那一项保留下来了， λ 2 < 0 \lambda_2<0 λ2​<0 因此那一项随着 t t t 的增大变成了 0 0 0. 而如果有复数特征值 λ = a + b j \lambda=a+bj λ=a+bj,
∣ e a + b j ∣ = ∣ e a ∣ ∣ e j b ∣ = ∣ e a ∣ |e^{a+bj}|=|e^a| |e^{jb}|=|e^a| ∣ea+bj∣=∣ea∣∣ejb∣=∣ea∣

换句话说， e j b e^{jb} ejb 并不改变模值，只是让特征向量 x \bm{x} x 的每个 entry 忽大忽小而已 (最终模值不变)。 因此，正如下入所示，如果想要 e A t → 0 e^{\bm{A t}}\to 0 eAt→0， A \bm{A} A 所有特征值的实部都必须小于 0。如果有一个 λ = 0 \lambda=0 λ=0, 则那一项的系数会保留下来，系统会慢慢进入steady state (上面那个例子就是如此)。

高阶微分方程

在这节课的最后， Prof. Strang 给出了高阶微分方程的解题思路。

比如我们有一个二阶微分方程
y ′ ′ + b y ′ + k y = 0 y''+by'+ky=0 y′′+by′+ky=0

解决技巧在于构造新函数 (就像上节中的斐波那契数列一样)。构造
u = [ y ′ y ] \bm{u}=\begin{bmatrix} y'\\ y \end{bmatrix} u=[y′y​]

我们便能得到
u ′ = [ y ′ ′ y ′ ] = [ − b − k 1 0 ] [ y ′ y ] = [ − b − k 1 0 ] u \bm{u'}=\begin{bmatrix} y''\\ y' \end{bmatrix}=\begin{bmatrix} -b & -k\\ 1 & 0 \end{bmatrix}\begin{bmatrix} y'\\ y \end{bmatrix}=\begin{bmatrix} -b & -k\\ 1 & 0 \end{bmatrix}\bm{u} u′=[y′′y′​]=[−b1​−k0​][y′y​]=[−b1​−k0​]u

这样便把一个二阶微分方程转换为了一阶微分方程，用这节的通解即可以解决。

如果更高阶尼？比如五阶？方法还是一样的，构造
u = [ y ′ ′ ′ ′ y ′ ′ ′ y ′ ′ y ′ y ] \bm{u}=\begin{bmatrix} y''''\\ y''' \\ y''\\ y'\\ y \end{bmatrix} u=⎣⎢⎢⎢⎢⎡​y′′′′y′′′y′′y′y​⎦⎥⎥⎥⎥⎤​

即可不断地降维.

That’s the end of today’s lecture~