計算矩陣指數

為什麼我們從不加總那個級數

上一篇指南把矩陣指數 e^(At) 定義為與普通指數相同的那個無窮級數，只是把矩陣 A t 塞進每一個冪次：e^(At) = I + A t + (A t)^2 / 2! + (A t)^3 / 3! + 這樣永遠下去。那個定義既美麗又誠實——但你若真要對一個一般的 3×3 矩陣把它加起來，你會一路計算矩陣的冪、除以階乘，直到耐心耗盡，卻看不到任何封閉形式。級數告訴你這個對象「是」什麼；它幾乎從不告訴你這個對象「等於」什麼。

有一個值得放進口袋的快樂例外：若 A 是對角的，比如 A = [a, 0; 0, b]，那麼每個冪 A^k 不過是 a^k 與 b^k 排在對角線上，級數於是裂成兩個獨立的純量指數級數，你完全不費力就讀出 e^(At) = [e^(at), 0; 0, e^(bt)]。那個小小的奇蹟，正是第一個真正方法的種子——因為大多數矩陣並非對角，但許多矩陣可以藉由更換座標而「被弄成」對角。

方法一：對角化，乘上特徵基底這班車

假設 A 擁有一整組獨立的特徵向量——把它們收作矩陣 P 的各列，並把相應的特徵值放在 D 的對角線上。那麼 A = P D P^(-1)，而魔術就在這裡：每個冪都滿足 A^k = P D^k P^(-1)，因為中間那些 P^(-1) P 的配對會相互抵消。把它逐項代入級數，P 與 P^(-1) 會乾淨地分提到兩端，把 e^(Dt) 的級數困在中間。既然 D 是對角的，e^(Dt) 就是剛才那個容易的「指數排在對角線上」的矩陣。

這就給出以對角化計算 e^(At)的完整公式：e^(At) = P e^(Dt) P^(-1)。用白話說，換座標進入系統解耦成獨立純量片段的那個特徵基底、把每個純量取指數、再換回來。這與你已經見過、用於方程組的特徵值法——一條直線模態一條地手解——是完全相同的想法，只是重新打包成一個乾淨的矩陣。特徵值住在指數裡，特徵向量住在 P 裡，而矩陣指數只是把它們全部一次記好帳。

e^(At) = P e^(Dt) P^(-1),     where  A = P D P^(-1)

   D = [lambda1,   0   ]        e^(Dt) = [ e^(lambda1 t),     0       ]
       [   0,   lambda2]                 [     0,        e^(lambda2 t) ]

   columns of P  = eigenvectors
   diagonal of D = eigenvalues

對角化配方：唯一被取指數的是容易的對角塊 e^(Dt)；P 與 P^(-1) 把它送回原座標。

方法二：對角化失敗時，攀上 Jordan

對角化有一個真實的極限，而它正是糾纏重特徵值情形的那一個：有些矩陣根本沒有足夠的獨立特徵向量去填滿 P。一個經典的元兇是 [2, 1; 0, 2]——它唯一的特徵值是 2，卻只提供一個特徵向量方向，於是不存在由特徵向量構成的可逆 P，對角化這條路就卡住了。你無法憑空許願出第二個特徵向量；你必須拓寬工具箱。

補救之道是 Jordan 形式。任何矩陣都能寫成 A = P J P^(-1)，其中 J 是分塊對角的，每一塊在對角線上是一個特徵值 lambda，而緊鄰其上方有一串 1——P 的各列如今是特徵向量，再以廣義特徵向量補足以湊成完整的基底。同樣的分提技巧給出 e^(At) = P e^(Jt) P^(-1)，於是你再一次只需對塊 J 取指數。而這裡第二個小奇蹟救了你：Jordan 塊對角線上方那部分 N 是冪零的（它的某個冪是零矩陣），所以它的指數級數會在幾項之後「終止」，而不會永遠跑下去。

具體地說，對塊 [lambda, 1; 0, lambda]，Jordan 形式法給出 e^(Jt) = e^(lambda t) [1, t; 0, 1]。那個跑出來的因子 t，正是重根解中那個熟悉的 t e^(lambda t) 項的來源——你在二階方程裡見過的「多項式乘指數」的標記，其實就是某個 Jordan 塊的矩陣指數，換了一頂帽子戴。同樣的數學，終於統一了。

方法三與方法四：完全跳過特徵向量

求特徵向量（更糟的是求廣義特徵向量）既繁瑣又容易出錯，於是還有兩個方法繞開它們。第一個倚靠 Cayley-Hamilton 定理：一個矩陣滿足它自己的特徵方程。其深刻的後果是，一個 n×n 矩陣的「所有」冪——因而 e^(At) 本身——都能寫成 A 的、次數至多 n-1 的多項式。於是 e^(At) = c0(t) I + c1(t) A + ... + c_(n-1)(t) A^(n-1)，你只需找出那些純量係數函數，而不必求任何特徵向量。把這些係數有條理地磨出來的俐落辦法是 Putzer 演算法，它僅憑特徵值、透過一個小小的遞迴純量 ODE 系統把它們建構出來。

第四條路重新動用一個你已經信任的工具。因為 e^(At) 就是主基本矩陣——X' = A X 配 X(0) = I 的唯一解——你可以用拉普拉斯變換去解那個矩陣初值問題。對 X' = A X 變換會把導數變成乘法，給出 (s I - A) 乘上變換等於 I，於是 e^(At) 的變換就是矩陣逆 (s I - A)^(-1)。把每個元素逆變換回時域，你就得到答案。這就是以拉普拉斯變換計算 e^(At)：不需特徵向量，只需一個矩陣逆，以及對每個元素查一次表。

可對角化、特徵值相異或特徵向量夠多？用 e^(At) = P e^(Dt) P^(-1)——通常是手算最快的。
虧損（重特徵值、特徵向量不足）？用 Jordan 形式，e^(At) = P e^(Jt) P^(-1)，並記得每塊的級數會提早停止。
想避開特徵向量但已知特徵值？跑 Putzer／Cayley-Hamilton，把 e^(At) 寫成 A 的多項式。
對變換上手，或矩陣很雜亂？算出 (s I - A)^(-1)，再用拉普拉斯表逐元素逆變換。

同一個矩陣、四條路、誠實的但書

這四個方法計算出的是完全相同的 e^(At)——它們是通往同一個房間的不同門，哪扇門最容易，取決於你面前的矩陣，以及哪些工具在你手裡最自然。一旦有了 e^(At)，整個齊次系統就一筆解決：x(t) = e^(At) x(0)。而且因為 e^(At) 是那個錨定在單位矩陣的真正基本矩陣，它還服從流的定律 e^(A(t+s)) = e^(At) e^(As)——這個半群性質讓你能把相繼時間區間上的演化複合起來——接下來兩篇指南在引入外力時，正要利用這個事實。

在你慶祝之前，有兩個誠實的極限。其一，上面每個手算方法都暗地裡需要特徵值，這意味著要解特徵多項式——而對一個 5×5 或更大的矩陣，其根一般沒有公式可循，於是「精確計算 e^(At)」悄悄變得無法手算，工作便移交給電腦。其二，即使在電腦上，矩陣指數也以難搞著稱：當特徵向量近乎平行、或特徵值分佈極廣時，樸素的級數、乃至對角化都可能在數值上不穩定，這正是為什麼會有專門的演算法（縮放-平方法、Krylov 方法）存在。這裡的四條路，是給人類尺度問題用的精確算術工具；真正的大規模計算，自有其另一門手藝。