$x_t$ 和 $B$ 是 1 乘 n 矩陣，$A$ 和 $M$ 是 n 乘 n 矩陣，0 是 n 乘 n 的零矩陣或 1 乘 n 的零矩陣，   $I$ 是 n 乘 n 的單位矩陣，

$\begin{aligned}\alpha(M) \equiv e^M \equiv I + M + \frac{M^2}{2!} + \frac{M^3}{3 !} + \cdots \\\beta(M) \equiv I + \frac{M}{2 !} + \frac{M^2}{3 !} + \frac{M 3}{4 !} + \cdots = M^{-1} (e^M - I) = (e^M - I) M^{-1}\end{aligned}$

$\alpha$ 函數和 $\beta$ 函數可以和任何係數是數字的也是$M$的多項式函數相乘順序交換。$\alpha$ 函 數 有個性質： 如果 $M_1 M_2 = M_2 M_1$ 則 $\alpha(M_1 + M_2) = \alpha(M_1) \alpha(M_2) = \alpha(M_2) \alpha(M_1)$ 。如果 $M_1 M_2 \ne M_2 M_1$ 則 $\alpha(M_1 + M_2)$ 和 $\alpha(M_1) \alpha(M_2)$ 沒有關係。

線性常微分方程 $\frac{d x_t}{d t} = x_t A + B$ ，初始條件是 $x_0$ 的解答是

$x_t = B A^{-1} \left(e^{t A} - I\right) + x_0 e^{t A}= t B(t A)^{-1} \left(e^{t A} - I \right) + x_0 e^{t A}= t B \beta(t A) + x_0 \alpha(t A)$

$t A$ 是沒有物理衡量單位的矩陣， $B$ 的物理衡量單位是 $x_t$ 的物理衡量單位除以時間。最終值相對於初始值仍是線性，所以會發散還是收斂和初始值無關，只和 $A$ 有關。故意改寫成這樣的理由是因為 $A$ 未必是可逆矩陣，寫成這樣電腦程式才不會程式當掉。例如當 $A = 0$ ，則因為 $\beta(t A) = I ,\alpha(t A) = I$ ，則 $x_t = t B + x_0$

當 $t$ 趨於無限大時如果 $e^{t A}$ 趨於 0，則解答最終會穩定在 $- B A^{-1}$ ，此時 $\frac{d x_t}{d t} = 0$ ，所以最終狀態的方程式也是 $0 = x A + B$

$\alpha(t A)$ 和 $\beta(t A)$ 要展開泰勒幾個項精確度才會夠的問題只和 $t A$ 這個矩陣有關，和測量單位無關。要做 symbolic 分析之前一定會有一個前提：要設定哪些對角線矩陣 entry 是 super positive或 super negative，然後針對它做以下 $e^{t A} = P(t) e^{t D} Q(t)$ 的方法的動作。

片段式的線性常微分方程式

週期是 $T$ 但是每個周期細分有數個片段，每個片段都是一個普通的線性微分方程。以 $x_0$ 表示時間點 0 的狀態， $x_1$ 表示時間點 $t_1$ 的狀態， $x_2$ 表示時間點 $t_1 + t_2$ 的狀態，.... ，$x_{K - 1}$ 表示時間點 $t_1 + t_2 + \cdots + t_{K - 1}$ 的狀態， $x_K$ 表示時間點 $t_1 + t_2 + \cdots + t_{K - 1} + t_K = T$ 的狀態。最後狀態趨於在一個周而復始的穩定狀態，也就是 $x_0 \rightarrow x_1 \rightarrow x_2 \rightarrow \cdots \rightarrow x_{K -1} \rightarrow x_K = x_0$

以 $A_k$ 和 $B_k$ 表示第 $k$ 個片段的 $A$ 和 $B$。因為

$\begin{aligned}x_1 &= t_1 B_1 \beta(t_1 A_1) + x_0 \alpha(t_1 A_1) \\x_2 &=t_2 B_2 \beta(t_2 A_2) + x_1 \alpha(t_2 A_2)\\&\cdots\\x_K &=t_K B_K \beta(t_K A_K) + x_{K -1} \alpha(t_K A_K)\end{aligned}$

所以 (為方便，定義 $\alpha(empty) = I$ )

$x_0 = x_K = x_0 \prod_{j = 1}^K \alpha(t_j A_j) + \sum_{k = 1}^K {t_k B_k \beta(t_k A_k) \prod_{j = k + 1}^K \alpha(t_j A_j)}$

所以可以解出 $x_0$

定義每個片段的比例 $d_j \equiv \frac{t_j}{T},A \equiv d_1 A_1 + d_2 A_2 + \cdots + d_K A_K,B \equiv d_1 B_1 + d_2 B_2 + \cdots + d_K B_K$

$\begin{aligned}x_0 = -T \left(\sum_{k = 1}^K {d_k B_k \beta(T d_k A_k) \prod_{j = k + 1}^K \alpha(T d_j A_j)}\right) \left(\prod_{j = 1}^K \alpha(T d_j A_j) - I\right)^{-1}\\= -\left(\sum_{k = 1}^K {d_k B_k \beta(T d_k A_k) \prod_{j = k + 1}^K \alpha(T d_j A_j)} \right) \left(\frac{1}{T} \left(\prod_{j = 1}^K \alpha(T d_j A_j) - I\right)\right)^{-1}\end{aligned}$

如果系統是用 column 來表示 vector，上述數學式就 Transpose。這是數學上的解。當 $T$ 很小的時候，因為 $\beta$ 和 $\alpha$   幾乎是 $I$，所以 $x_0 \approx - B A^{-1}$ ，和單純的線性微分方程的解的樣子很像。如果 $T$ 不小或者想根據 $T$ 稍微作泰勒展開也可以看想展開到 $T$ 的幾次方項，$\beta(T d_k A_k)$ 和 $\alpha(T d_k A_k)$ 都是 $T$ 的的多項式或多項式的倒數，比如說要看到 $T$ 的三次方，就把展開之後把三次方以上的項砍掉。例如

$\begin{aligned}(\beta(T A))^{-1} \approx \left( I + \frac{T}{2} A + \frac{T^2}{6} A^2 + \frac{T^3}{24} A^3 \right)^{-1} \\ \approx I - \left(\frac{T}{2}A + \frac{T^2}{6} A^2 + \frac{T^3}{24} A^3 \right) + \left(\frac{T}{2} A + \frac{T^2}{6} A^2 + \frac{T^3}{24} A^3 \right)^2 - \left(\frac{T}{2}A + \frac{T^2}{6} A^2 + \frac{T^3}{24} A^3 \right)^3 \\\approx I - \frac{T}{2} A - \frac{T^2}{6} A^2 - \frac{T^3}{24} A^3 + \frac{T^2}{4} A^2 + \frac{T^3}{6} A^3 - \frac{T^3}{8} A^3 \\\approx I - \frac{T}{2} A + \frac{T^2}{12} A^2\end{aligned}$

$\sum_{k = 1}^K {d_k B_k \beta(T d_k A_k) \prod_{j = k + 1}^K \alpha(T d_j A_j)}$ 這一項看似複雜，程式上有比較方便的寫法：

$\begin{aligned}ans \leftarrow d_1 B_1 \beta(T d_1 A_1)\\ans \leftarrow d_2 B_2 \beta(T d_2 A_2) + ans \cdot \alpha(T d_2 A_2)\\ans \leftarrow d_3 B_3 \beta(T d_3 A_3) + ans \cdot \alpha(T d_3 A_3)\\ \cdots \\ans \leftarrow d_K B_K \beta(T d_K A_K) + ans \cdot \alpha(T d_K A_K)\end{aligned}$

到底需要多少個 term 的 polynomial 來描述，端看 $T d_j A_j$ 這些 metric-free 的矩陣性質而定。去逼近 $\alpha(T d_j A_j)$ 或 $\beta(T d_k A_k)$ 也不一定只能用指數的 power series ，也可以混用 exp 和 power series像是 $\sum_i e^{\omega_i t} \left(a_{i 0} + a_{i 1} t + a_{i 2} t^2 + \cdots \right)$

事實上，越有效率的逼近就是項數越少卻還能維持精確度的逼近，所謂 eigensystem 就是選最準的 exp 來描述此時 $\omega_i$ 是 eigen value，只需要 n 個 exp其中 n 是矩陣的行數但是也很難有 symbolic formula。

exp(t A) = P(t) exp(t D) Q(t) 的方法

給定一個方矩陣 $A$，以 $E$ 表示 eigen column vectors 構成的矩陣，$\lambda_i$ 代表 eigen values ，則

$A = E \begin{bmatrix}\lambda_1&&\\&\ddots&\\&&\lambda_n\end{bmatrix} E^{-1}$

而且

$exp(t A) = E \begin{bmatrix}e^{t \lambda_1}&&\\& \ddots &\\&&e^{t \lambda_n}\end{bmatrix} E^{-1}$

根據這個啟示，所以這裡的逼近 idea 設定是

$\begin{aligned}exp(t\frac{A}{2})P(0)=P(t)exp(t\frac{D}{2})\\Q(0)exp(t\frac{A}{2})=exp(t\frac{D}{2})Q(t)\end{aligned}$

其中 $I = P(0)Q(0)$ 。左邊乘左邊，右邊乘右邊：

$exp(t \frac{A}{2}) P(0) Q(0) exp(t \frac{A}{2}) = P(t) exp(t \frac{D}{2}) exp(t \frac{D}{2}) Q(t)$ 也就是 $exp(t A) = P(t) exp(t D) Q(t)$

現在問題是去決定 $P(t)$ 和 $Q(t)$ 這兩個級數的係數。以下的 $A$ 是上面的 $\frac{A}{2}$ ，$D$ 是上面的 $\frac{D}{2}$ ，

$P_0(t) \equiv P(t),Q_0(t) \equiv Q(t),P_i(t) \equiv \frac{d^i}{d t^i} P(t)$

不斷對 $t$ 微分：

$\begin{aligned}e^{t A} P(0) = P(t) e^{t D}\\e^{t A} A P(0) = \left(P_ 1(t) + P_0(t) D \right) e^{t D}\\e^{t A} A^2 P(0) = \left(P_2(t) + P_1(t) D \right) e^{t D} + (P_1(t) + P_0(t) D) D e^{t D} = \left(P_2(t) + 2 P_1(t) D + P_0(t) D^2 \right) e^{t D}\\e^{t A} A^3 P(0) = \left(P_3(t) + 2 P_2(t) D + P_1(t) D^2 \right) e^{t D} + \left(P_2(t) + 2 P_1(t)D + P_0(t) D^2 \right) D e^{t D} \\= \left(P_3(t) + 3 P_2(t) D + 3 P_1(t) D^2 + P_0(t) D^3 \right) e^{t D}\\ \cdots\end{aligned}$

在 $t=0$ 取值：

$\begin{aligned}A P(0) = P_1(0) + P_0(0) D\\A^2 P(0) = P_2(0) + 2 P_1(0) D + P_0(0) D^2\\A^3 P(0) = P_3(0) + 3 P_2(0) D + 3 P_1(0) D^2 + P_0(0) D^3\\\cdots\end{aligned}$

所以，${}_j^i C$ 代表二項式係數：

$\begin{aligned}P_1(0) = A P(0) - P_0(0) D\\P_2(0) = A^2 P(0) - \left(2 P_1(0) D + P_0(0) D^2 \right)\\P_3(0) = A^3 P(0) - \left(3 P_2(0) D + 3 P_1(0) D^2 + P_0(0)D^3 \right)\\P_i(0) = A^i P(0) - \left( {}_1^i C P_{i-1}(0) D + {}_2^i C P_{i-2}(0) D^2 + \cdots + {}_{i-2}^i C P_2(0) D^{i-2} + {}_{i-1}^i C P_1(0) D^{i-1} + P_0(0) D^i \right)\end{aligned}$

依序就得到 $P_i(0)$ 。類似的，

$Q_i(0) = Q(0)A^i - \left( {}_1^i C D Q_{i-1}(0) + {}_2^i C D^2 Q_{i-2}(0)+ \cdots + {}_{i-2}^i C D^{i-2}Q_2(0) + {}_{i-1}^i C D^{i-1}Q_1(0) + D^i Q_0(0) \right)$

$P_i(0)$ 和 $Q_i(0)$ 這些方陣都知道後，

$\begin{aligned}P(t) = P(0) + \sum_{i = 1} \frac{P_i(0)}{i!} t^i\\Q(t) = Q(0) + \sum_{i = 1} \frac{Q_i(0)}{i!} t^i\end{aligned}$

當 $D$ 就是 eigen values 對角線矩陣而且 $P(0)$ 就是 eigen column vectors 的時候， $any(X)$ 是任何 $X$ 的級數，就有： $any(t A) \cdot P(0) = P(0) \cdot any(t D)$ 也就是

現在$any(t A) \cdot P(0) (any(t D))^{-1} = P(0)e^{t \frac{A}{2}} P(0) = P(t) e^{t \frac{D}{2}}$ 也就是 $e^{t \frac{A}{2}} P(0) (e^{t \frac{D}{2}})^{-1} = P(t)$

當 $D$ 猜得很接近 eigen values 對角線矩陣而且 $P(0)$ 猜得很接近 eigen column vectors ，也就是 $P(t)$ 也會很接近 $P(0)$ ，所以 $P(t)$ 這個級數的次數就不用太多，這個式子 $P(t) e^{t D} Q(t)$ 就會是一個計算 $e^{t A}$ 很經濟或是可以方便 symbolic 分析問題的式子。除非使用數值方法，eigen values 一般只可遠觀不可褻玩，Abel's impossibility theorem 說五次以上方程式是沒有公式的。