CS236 補充課 B · 微分方程入門
CS236 補充課 B · 前置數學

微分方程入門：只學過微積分，從頭開始

這一頁假設你學過微積分，但完全沒學過微分方程。我們從「微分方程到底是什麼」開始，手算第一個解，再一路走到相圖、特徵值、數值解法，最後看懂 DDIM 為什麼能用很少的步數。每一步都會寫出來，不會跳過。

先修：微積分；第 7 節起用到補充課 A 的矩陣與特徵值 服務：第 7–8、14、16、17、18 講 預計閱讀 3 小時，可以分次讀

讀完這頁你應該能夠

  1. 說明微分方程在問什麼，並用「代回去」驗算一個解是不是對的。
  2. 用分離變數一步一步解出 \(\frac{dx}{dt} = -kx\)，解釋積分常數如何變成初始值。
  3. 看懂斜率場、向量場和相圖，並從圖上想像解的樣子。
  4. 說明為什麼軌跡不會交叉，以及為什麼這讓 ODE 可以倒著走。
  5. 推導為什麼從特徵向量 \(\mathbf v\) 出發的解是 \(e^{\lambda t}\mathbf v\)，以及為什麼兩個特徵值就決定了整張相圖。
  6. 寫出 Euler、Heun、RK4 的更新式，並用面積的圖解釋它們為什麼一個比一個準。
  7. 解釋 Euler 為什麼在 \(h\lambda > 2\) 時會爆掉，並一步一步推導指數積分器，看懂它和 DDIM 的關係。

00怎麼讀這一頁

只需要這三樣微積分工具：導數是變化率；積分是導數的反運算；連鎖律 \(\frac{d}{dt}f(g(t)) = f'(g(t))\,g'(t)\)。另外會反覆用到 \(\frac{d}{dt}e^{at} = ae^{at}\) 和 \(\frac{d}{dx}\ln x = \frac1x\)。

第 1–4 節只有一個變數，一定要慢慢讀。第 5–7 節進入多個變數，需要矩陣和特徵值；如果還不熟，先讀補充課 A 第 3–4 節。第 8–11 節是電腦怎麼解 ODE。綠色線標示的「手算」區塊，建議拿紙筆跟著算。

01微分方程是什麼：未知的是一個函數

國中解的方程式，未知的是一個數字。例如 \(2x + 1 = 7\)，答案是 \(x = 3\)。

微分方程未知的是一個函數，而且方程式描述的是這個函數的變化率。例如：

$$\frac{dx}{dt} = 2$$

意思是「有一個隨時間變化的量 \(x(t)\)，它每單位時間固定增加 2」。答案不是一個數字，而是一個函數，例如 \(x(t) = 2t\)。

三個一定要知道的詞

  • 解：一個函數，代入方程式後等號在所有時間都成立。
  • 初始條件：起點的值，例如 \(x(0) = 5\)。\(x = 2t\) 和 \(x = 2t + 5\) 都滿足 \(\frac{dx}{dt} = 2\)，只知道規則無法決定是哪一條，還要知道從哪裡出發。
  • ODE（ordinary differential equation，常微分方程）：只對一個變數（這裡是時間 \(t\)）微分的方程。對多個變數微分的叫 PDE，補充課 C 才會遇到。

這一頁的 ODE 都長這樣：

ODE 的一般形式 $$\frac{dx}{dt} = v(x, t)$$

把它讀成一條速度規則：如果現在時間是 \(t\)、位置是 \(x\)，那麼速度就是 \(v(x, t)\)。解 ODE 就是「照著這條規則，從起點一直走，看會走出什麼路線」。常見的另一種記法是在變數上加一點：\(\dot x\) 就是 \(\frac{dx}{dt}\)。

怎麼驗算一個解

解方程式很難，但驗算很簡單：把函數代回去，看左右兩邊是否相等。這個技巧整頁都會用。

手算

\(x(t) = 3e^{-2t}\) 是不是 \(\frac{dx}{dt} = -2x\) 的解？

  • 左邊：\(\frac{d}{dt}3e^{-2t} = -6e^{-2t}\)
  • 右邊：\(-2\cdot3e^{-2t} = -6e^{-2t}\)

兩邊相等，所以是解，而且 \(x(0) = 3\)。反過來，\(x(t) = t^2\) 代入後左邊是 \(2t\)、右邊是 \(-2t^2\)，不相等，所以不是解。

一句話帶走：ODE 是一條速度規則，解是順著規則走出來的函數；驗算時把函數代回去即可。
小練習：\(x(t) = \sin t\) 是 \(\frac{dx}{dt} = \cos t\) 的解嗎？\(x(t) = \sin t + 7\) 呢？
都是。兩者的導數都是 \(\cos t\)。差別只在初始條件：前者 \(x(0) = 0\)，後者 \(x(0) = 7\)。

02最簡單的情況：右邊只有 t，就是積分

如果速度只和時間有關，\(\frac{dx}{dt} = f(t)\)，這其實就是微積分：已知速度求位置，兩邊對 \(t\) 積分就好。

手算

\(\frac{dx}{dt} = 2t\)，\(x(0) = 1\)。

  1. 兩邊積分：\(x(t) = t^2 + C\)。\(C\) 是積分常數，因為常數微分後是 0，積分無法確定它。
  2. 用初始條件決定 \(C\)：代入 \(t = 0\)，得到 \(1 = 0 + C\)，所以 \(C = 1\)。
  3. 答案：\(x(t) = t^2 + 1\)。
  4. 驗算：\(\frac{dx}{dt} = 2t\) ✓，\(x(0) = 1\) ✓。
起點 x(0) = 1x = t² + 1tx
每一個 \(C\) 都是一條解，它們互相上下平移，合起來填滿整個平面。初始條件從中選出一條。

這裡已經出現了微分方程的基本結構：規則給出一整族解，初始條件選出其中一條。

一句話帶走：右邊只有 \(t\) 時，積分就解完了；積分常數由起點決定。

03第一個真正的 ODE：剩越多、掉越快

為什麼這題不能直接積分

放射性物質每單位時間衰變掉固定的比例，所以剩越多、衰變越快：

$$\frac{dx}{dt} = -k\,x\qquad(k > 0)$$

右邊是 \(x\)，不是 \(t\)。想直接積分會變成 \(x = -k\int x\,dt\)，但要積的正是我們還不知道的 \(x(t)\)，繞了一圈。這就是微分方程和普通積分真正不同的地方：速度取決於現在的位置，而位置又取決於過去的速度。

先不解，直接照規則走

就算不會解，也可以一步一步算。設 \(k = 0.1\)、\(x(0) = 100\)，每次走 \(h = 1\) 的時間。假設這一小段內速度不變，新值 ＝ 舊值 ＋ 速度 × 時間：

手算
\(t\)\(x\)速度 \(-0.1x\)下一個 \(x\)
0100−1090
190−981
281−8.172.9
372.9−7.2965.61

每一步都是乘上 \(0.9\)，所以 \(t = 10\) 時是 \(100\times0.9^{10} = 34.87\)。

但「這段時間內速度不變」並不準確，因為 \(x\) 一直在變小，速度也一直在變小。步長越小，這個假設越準：

步長 \(h\)每步乘上\(t = 10\) 時的 \(x\)
20.832.77
10.934.87
0.50.9535.85
0.10.9936.60
極限 \(h\to0\)36.79

數字收斂到某個值。微積分學過 \((1 - kh)^{t/h}\to e^{-kt}\)（\(h\to0\)），所以極限是 \(100e^{-1} = 36.79\)。指數函數就是這樣出現的：每一小步都乘上一個接近 1 的數，乘無限多次就是指數。

每步 h = 5，沿切線走精確解 100e^(−0.1t)t100
橙色折線每走一段就重新看一次速度。步長越小，它越貼近黑色的精確解。

Demo B1 · 照著速度規則一步一步走

怎麼看：橙色是上面表格的算法，黑色是精確解。把 \(h\) 拉小，折線會貼上曲線。再把 \(h\) 拉大到 \(kh > 1\)：一步就衝過 0，變成負的，這個現象第 10 節會解釋。
精確解 100·e^(−kt)每步用當下的速度走 h

精確解法：分離變數

現在正式解出來。想法是：把和 \(x\) 有關的都搬到左邊，和 \(t\) 有關的都搬到右邊，兩邊就都能積分了。

  1. 兩邊除以 \(x\)（先假設 \(x\neq0\)）： $$\frac{1}{x}\frac{dx}{dt} = -k$$
  2. 左邊正是某個函數的導數。由連鎖律，\(\frac{d}{dt}\ln x(t) = \frac{1}{x}\frac{dx}{dt}\)。所以方程式變成： $$\frac{d}{dt}\ln x(t) = -k$$ 這就是 \(\ln x\) 版本的第 2 節！右邊只是常數。
  3. 兩邊對 \(t\) 積分： $$\ln x(t) = -kt + C$$
  4. 兩邊取指數，把 \(\ln\) 去掉： $$x(t) = e^{-kt + C} = e^{C}\,e^{-kt}$$
  5. \(e^C\) 只是一個還沒決定的常數，記成 \(A\)。代入 \(t = 0\)：\(x(0) = A\cdot e^0 = A\)。所以積分常數就是初始值。
結果 $$\frac{dx}{dt} = -kx\ \Rightarrow\ x(t) = x(0)\,e^{-kt}$$

驗算：\(\frac{d}{dt}x(0)e^{-kt} = -k\,x(0)e^{-kt} = -k\,x(t)\) ✓。

你在其他書上可能會看到更簡短的寫法：把 \(dx\) 和 \(dt\) 當成分數拆開，得到 \(\frac{dx}{x} = -k\,dt\)，再兩邊分別積分。它和上面是同一件事，只是省略了連鎖律那一步。

\(k\) 的正負決定一切：\(\dot x = ax\) 的解是 \(x(0)e^{at}\)。\(a \lt 0\) 衰減到 0，\(a > 0\) 爆炸，\(a = 0\) 不動。第 7 節會反覆用到這個結果。

手算：另一個可以分離的例子

\(\frac{dx}{dt} = x^2\)，\(x(0) = 1\)。除以 \(x^2\)：\(x^{-2}\frac{dx}{dt} = 1\)。左邊是 \(\frac{d}{dt}(-x^{-1})\)，積分得 \(-\frac1x = t + C\)。代入 \(t = 0, x = 1\)，得 \(C = -1\)，所以 \(x(t) = \frac{1}{1-t}\)。

驗算：\(\frac{dx}{dt} = \frac{1}{(1-t)^2} = x^2\) ✓。注意它在 \(t = 1\) 衝到無窮大：ODE 的解不一定對所有時間都存在。

大多數 ODE 無法這樣分離，例如 \(\frac{dx}{dt} = \sin(2t) - x\)（第 11 節會解）；速度場是神經網路時，更不可能有公式解。所以後面才需要圖形和數值方法。

一句話帶走：\(\dot x = -kx\) 的解是 \(x(0)e^{-kt}\)：除以 \(x\) 後左邊是 \(\ln x\) 的導數，積分再取指數，積分常數就是起點。
小練習：\(\frac{dx}{dt} = 0.5x\)，\(x(0) = 2\)。\(x(t)\) 是什麼？驗算一次。
\(x(t) = 2e^{0.5t}\)。驗算：\(\frac{dx}{dt} = 2\cdot0.5\,e^{0.5t} = e^{0.5t}\)，而 \(0.5x = e^{0.5t}\) ✓。

04斜率場：不解方程也能看出解的樣子

公式解很少見，但有一個辦法永遠可用：把規則畫出來。在 \((t, x)\) 平面的每一點，規則都告訴我們解在那裡的斜率 \(\frac{dx}{dt}\)。畫一小段那個斜率的線段，就是斜率場。

例如 \(\frac{dx}{dt} = -x\)。在 \(x = 2\) 的地方斜率是 \(-2\)（陡陡地往下）；\(x = 0\) 的地方斜率是 0（水平）；\(x = -1\) 的地方斜率是 \(+1\)（往上）。

tx
灰色線段是斜率場。解曲線處處和線段相切，就像順著水流漂的落葉。

平衡點

讓速度為 0 的位置叫平衡點：粒子放在那裡就不會動。平衡點分兩種，看附近的粒子是被吸回來還是被推開：

手算

族群成長 \(\frac{dx}{dt} = x(1-x)\)。令右邊等於 0，得到平衡點 \(x = 0\) 和 \(x = 1\)。

  • \(x = 0.1\)：速度 \(0.09 > 0\)，往上走，離開 0。所以 \(x = 0\) 是不穩定的。
  • \(x = 0.9\)：速度 \(0.09 > 0\)，往上走，靠近 1。\(x = 1.2\)：速度 \(-0.24\)，往下走，也靠近 1。所以 \(x = 1\) 是穩定的。

Demo B2 · 斜率場與解曲線（點一下畫布加一條解）

怎麼看：點畫布任一點，就會畫出通過那一點的解。多點幾條，你會發現它們都不會交叉。切換方程：
  • 衰減：所有解都往 0 靠近，越來越貼近，但不會碰在一起。
  • 族群成長：解被 \(x = 1\) 吸過去，被 \(x = 0\) 推開。
  • 被推動的衰減：不論從哪裡出發，過一段時間都被拉到同一條擺動的曲線上，起點被「忘掉」了。第 16 講說「雜訊夠多時，從哪裡出發都一樣」，是同一種現象。
  • 爆炸：\(x(0) = 1\) 的解在 \(t = 1\) 衝到無窮大。
斜率場預設的幾條解你點的解
一句話帶走：每一點的規則就是一個斜率，解是和這些斜率處處相切的曲線。
小練習：\(\frac{dx}{dt} = x - x^3\) 的平衡點有哪些？哪些是穩定的？
\(x(1-x)(1+x) = 0\)，平衡點是 \(0, 1, -1\)。\(x = 0.5\) 速度是正的（往 1 走），\(x = 2\) 速度是負的（也往 1 走），所以 \(x = 1\) 穩定；對稱地 \(x = -1\) 也穩定。\(x = 0\) 兩側都被推開，不穩定。

05二維：向量場與相圖

為什麼需要多個變數

一張圖片有上百萬個像素，是一個上百萬維的點。Diffusion model 就是在這個空間裡推動粒子。維度太高畫不出來，但道理在二維時就看得到。

從一條方程變成一組方程

粒子的位置現在有兩個座標 \((x, y)\)，規則也變成兩條：

$$\frac{dx}{dt} = v_1(x, y),\qquad\frac{dy}{dt} = v_2(x, y)$$

寫成向量就是 \(\frac{d\mathbf x}{dt} = \mathbf v(\mathbf x)\)，和一維長得一模一樣。平面上每一點都有一支速度箭頭 \(\mathbf v(\mathbf x)\)，合起來叫向量場。把粒子放進去，它就順著箭頭走，走出來的路線叫軌跡。

手算

\(\frac{dx}{dt} = -y,\ \frac{dy}{dt} = x\)，從 \((1, 0)\) 出發。在 \((1, 0)\) 的速度是 \((0, 1)\)，朝上；在 \((0, 1)\) 的速度是 \((-1, 0)\)，朝左。看起來在繞圈。

猜 \(x = \cos t,\ y = \sin t\) 並驗算：\(\frac{dx}{dt} = -\sin t = -y\) ✓，\(\frac{dy}{dt} = \cos t = x\) ✓。解確實是單位圓。

相圖是什麼

同一個解有兩種畫法。第一種是分別畫出 \(x, y\) 對 \(t\) 的曲線。第二種是不畫時間，只在 \((x, y)\) 平面上畫出粒子走過的路線。第二種就叫相圖（phase portrait）。它通常會畫很多條軌跡，每條從不同的起點出發，讓你一眼看到所有可能的未來。

x(t)y(t)t對時間畫t = 0xy相圖：時間藏起來
左邊是 \(x(t) = \cos t,\ y(t) = \sin t\) 對時間的圖，右邊是同一個解的相圖，也就是一個圓。相圖上看不到走多快，但能看到整體的形狀；箭頭標示前進的方向。

Demo B3 · 向量場與流動的粒子

怎麼看：灰色箭頭是速度場，藍色是相圖，按按鈕放入一群粒子。在score 場裡，粒子被吸到幾個峰頂，只靠確定性的流動，它們會擠成幾個點，而不是散成一個分佈。Langevin 取樣多加一點雜訊（補充課 F）才能讓粒子散開；probability flow ODE 則用一個隨時間改變的向量場，讓粒子剛好在終點散成資料分佈。
速度方向流線（從格點出發的軌跡）粒子

這就是 CS236 裡 probability flow ODE 和 flow matching 的畫面：學一個向量場，把一團雜訊粒子推成資料的形狀。

一句話帶走：向量場在每一點給一支速度箭頭，軌跡是順著箭頭走的路線，相圖把許多軌跡畫在同一張圖上。

06軌跡不交叉，所以可以倒著走

解是唯一的

只要速度場夠平滑，從每個起點出發的解都存在而且唯一。「夠平滑」的嚴格說法是 Lipschitz 連續：存在常數 \(L\)，對任意兩點都有

$$\Vert\mathbf v(\mathbf x) - \mathbf v(\mathbf y)\Vert\le L\Vert\mathbf x - \mathbf y\Vert$$

白話就是：兩點靠近時，它們的速度也要靠近，速度場不能有無限陡的斜率。

不交叉

唯一性直接帶來一個看得見的結果：兩條軌跡不可能交叉。

未來 1？未來 2？在交點只有一支速度箭頭，不可能同時往兩邊走
如果兩條軌跡在紅點相遇，從紅點出發就有兩個不同的未來，違反唯一性。

精確地說：速度場隨時間改變時，是「同一時刻不會有兩顆粒子在同一位置」。相圖上的路線，有可能在不同時刻經過同一點，但粒子不會同時在那裡相遇。

可以倒著走

如果把時間倒過來，令 \(\mathbf y(s) = \mathbf x(T - s)\)，由連鎖律：

$$\frac{d\mathbf y}{ds} = -\frac{d\mathbf x}{dt}\Big|_{t = T - s} = -\mathbf v(\mathbf y(s))$$

所以把速度反過來走，就會原路回去。把 ODE 從 0 解到 \(T\) 是一個可逆的映射：不同的起點一定走到不同的終點（否則就交叉了），而且可以用 \(-\mathbf v\) 走回來。第 14、17 講說「ODE 的解唯一，所以 diffusion model 可以變成 normalizing flow」，指的就是這件事。

Demo B4 · 往前流，再倒著流回來

怎麼看：向量場是 Demo B3 的 score 場。中間的圖案被擠成幾團，看起來資訊都「丟了」，但右圖倒流回來又還原了：資訊其實還在，只是被壓得很密。把解法換成 Euler、步數拉少，還原得就歪歪的，因為每一步的數值誤差不會自動抵銷（Euler 的定義見第 8 節）。
t = 0：起點的圖案
流到 t = T
用 −v 倒流回 t = 0

反例：速度場不平滑時，解可以不唯一

手算

\(\frac{dx}{dt} = \sqrt{|x|}\)，\(x(0) = 0\)。

  • \(x(t) = 0\) 是一個解：左邊 0，右邊 0。
  • \(x(t) = t^2/4\)（\(t\ge0\)）也是：左邊 \(t/2\)，右邊 \(\sqrt{t^2/4} = t/2\)。

粒子可以在原點停任意久再出發，未來不唯一。問題出在 \(\sqrt{|x|}\) 在 0 的斜率是無窮大，不是 Lipschitz 連續。神經網路用平滑的激活函數、權重有限，所以 diffusion 的速度場通常沒有這個問題。

出處：Lai et al. (2026) 附錄 A.1.2，局部與全域的存在唯一性定理，以及「Uniqueness Implies Non-Intersection of Solutions」一段。

一句話帶走：平滑的速度場解唯一，所以軌跡不交叉，流動是可逆的。
在 CS236 哪裡用到：第 14、17 講把 diffusion 的 ODE 當成 normalizing flow，算精確的 likelihood（補充課 C）。

07線性系統：為什麼特徵值決定了一切

為什麼需要它

一般的 ODE 沒有公式解，但有一大類可以完全解出來：線性系統，也就是速度是位置的矩陣乘法。

$$\frac{d\mathbf x}{dt} = A\mathbf x$$

它在 CS236 很重要：diffusion 的前向過程、第 18 講的離散 diffusion 都是線性的。而且任何平滑的 ODE，放大到一點附近看，都近似線性（Jacobian 就是那個矩陣），所以這裡學到的東西到處都用得上。

這一節會用到特徵值和特徵向量。如果忘了，先看補充課 A 第 4 節：\(A\mathbf v = \lambda\mathbf v\) 的 \(\mathbf v\) 是不會轉彎的方向，\(\lambda\) 是沿那個方向拉長的倍數。

第一步：沒有互相影響的情況

如果 \(A\) 是對角矩陣，兩個變數互不影響，就是兩個獨立的第 3 節：

$$\frac{d}{dt}\begin{bmatrix}x\\y\end{bmatrix} = \begin{bmatrix}-1&0\\0&-3\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}$$

也就是 \(\dot x = -x,\ \dot y = -3y\)，解是 \(x(t) = x(0)e^{-t},\ y(t) = y(0)e^{-3t}\)。很容易。

第二步：互相影響時，找不會互相影響的方向

現在考慮

$$A = \begin{bmatrix}-2&1\\1&-2\end{bmatrix}$$

也就是 \(\dot x = -2x + y,\ \dot y = x - 2y\)。\(x\) 的速度取決於 \(y\)，\(y\) 的速度取決於 \(x\)，兩條方程綁在一起，不能分開解。

換個問法：有沒有哪個方向，粒子從那裡出發後會一直留在同一條直線上？如果有，沿那條線就只剩一維的問題了。

設粒子一直留在方向 \(\mathbf v\) 上，也就是 \(\mathbf x(t) = f(t)\,\mathbf v\)，其中 \(f(t)\) 是一個純量函數，表示它在線上的位置。代入 ODE：

  • 左邊：\(\frac{d}{dt}\big(f(t)\mathbf v\big) = f'(t)\,\mathbf v\)（\(\mathbf v\) 是固定的向量）
  • 右邊：\(A\big(f(t)\mathbf v\big) = f(t)\,A\mathbf v\)

要讓 \(f'(t)\,\mathbf v = f(t)\,A\mathbf v\) 成立，\(A\mathbf v\) 必須和 \(\mathbf v\) 同方向，也就是 \(A\mathbf v = \lambda\mathbf v\)：\(\mathbf v\) 必須是特徵向量。代入後兩邊的 \(\mathbf v\) 消掉，只剩：

$$f'(t) = \lambda f(t)$$

這正是第 3 節解過的方程，答案是 \(f(t) = f(0)e^{\lambda t}\)。所以：

直線解

如果 \(A\mathbf v = \lambda\mathbf v\)，從 \(\mathbf v\) 出發的解就是

$$\mathbf x(t) = e^{\lambda t}\mathbf v$$

也可以直接驗算：左邊 \(\frac{d}{dt}e^{\lambda t}\mathbf v = \lambda e^{\lambda t}\mathbf v\)，右邊 \(A(e^{\lambda t}\mathbf v) = e^{\lambda t}A\mathbf v = e^{\lambda t}\lambda\mathbf v\) ✓。

手算

上面的 \(A\)：特徵方程 \(\lambda^2 + 4\lambda + 3 = (\lambda+1)(\lambda+3) = 0\)，所以 \(\lambda_1 = -1\)、\(\lambda_2 = -3\)。

\(\lambda_1 = -1\)：\(A + I = \begin{bmatrix}-1&1\\1&-1\end{bmatrix}\)，要求 \(v_1 = v_2\)，取 \(\mathbf v_1 = (1, 1)\)。

\(\lambda_2 = -3\)：\(A + 3I = \begin{bmatrix}1&1\\1&1\end{bmatrix}\)，要求 \(v_1 = -v_2\)，取 \(\mathbf v_2 = (1, -1)\)。

沿 \((1,1)\) 方向，\(\dot x = -2x + x = -x\)，確實以 \(e^{-t}\) 衰減；沿 \((1,-1)\) 方向（\(y = -x\)），\(\dot x = -2x - x = -3x\)，以 \(e^{-3t}\) 衰減。

沿 (1,1)：∝ e^(−t)沿 (1,−1)：∝ e^(−3t)
綠色虛線是兩個特徵方向。從線上出發的粒子永遠留在線上，只是以各自的速率靠近原點。

第三步：其他起點怎麼辦？拆開再加起來

起點不一定在特徵方向上，但線性系統有一個關鍵性質：兩個解相加還是解。如果 \(\dot{\mathbf a} = A\mathbf a\)、\(\dot{\mathbf b} = A\mathbf b\)，則

$$\frac{d}{dt}(\mathbf a + \mathbf b) = A\mathbf a + A\mathbf b = A(\mathbf a + \mathbf b)$$

所以做法是：

  1. 把起點拆成兩個特徵向量的組合：\(\mathbf x(0) = c_1\mathbf v_1 + c_2\mathbf v_2\)。
  2. 每一部分各自依第二步演化。
  3. 加起來：
一般解 $$\mathbf x(t) = c_1e^{\lambda_1t}\mathbf v_1 + c_2e^{\lambda_2t}\mathbf v_2$$
手算

同一個 \(A\)，起點 \((1, 0)\)。拆解：\((1, 0) = \tfrac12(1, 1) + \tfrac12(1, -1)\)，所以 \(c_1 = c_2 = \tfrac12\)。

$$\mathbf x(t) = \tfrac12e^{-t}\begin{bmatrix}1\\1\end{bmatrix} + \tfrac12e^{-3t}\begin{bmatrix}1\\-1\end{bmatrix}$$

也就是 \(x(t) = \tfrac12(e^{-t} + e^{-3t})\)、\(y(t) = \tfrac12(e^{-t} - e^{-3t})\)。

驗算 \(\dot x = -2x + y\)。左邊：\(\tfrac12(-e^{-t} - 3e^{-3t})\)。右邊：\(-(e^{-t} + e^{-3t}) + \tfrac12(e^{-t} - e^{-3t})\)，合併後是 \(-\tfrac12e^{-t} - \tfrac32e^{-3t}\)，和左邊相等 ✓。

t = 0t = 0.5t = 1.5
黑點是粒子，藍色和橙色箭頭是它在兩個特徵方向上的分量。橙色（\(e^{-3t}\)）縮得比藍色（\(e^{-t}\)）快，所以粒子最後會貼著藍色的方向進入原點。

Demo B5 · 把起點拆成兩個特徵方向

怎麼看：點畫布選起點，拖動時間 \(t\)。藍色和橙色箭頭是兩個分量，各自只沿自己的綠線伸縮；黑點是它們的和。讀數比較「公式算出來的」和「直接用數值方法解 ODE」的位置，兩者相同。切換到鞍點，把起點點在靠近 \((1,-1)\) 方向的地方：橙色分量縮小、藍色分量放大，所以粒子先靠近原點，再沿藍色方向離開。
c₁ e^(λ₁t) v₁，v₁ = (1, 1)c₂ e^(λ₂t) v₂，v₂ = (1, −1)粒子軌跡

為什麼只要看兩個特徵值，就知道整張相圖

現在可以回答這個問題了。一般解告訴我們，任何一條軌跡都只是 \(e^{\lambda_1t}\) 和 \(e^{\lambda_2t}\) 這兩個數字的組合，不同的起點只是改變 \(c_1, c_2\)。所以兩個特徵值的正負就決定了所有軌跡的共同命運：

穩定節點兩個負實數鞍點一正一負螺旋複數，實部負中心純虛數
四種典型的相圖。特徵向量決定直線的方向，特徵值決定沿直線進來還是出去、有多快。
特徵值相圖為什麼例子
兩個負實數穩定節點\(e^{\lambda t}\) 都衰減，所有軌跡進入原點DDPM 前向 SDE 的 drift \(-\tfrac12\beta\mathbf x\)
兩個正實數不穩定節點都放大，所有軌跡離開原點Demo B5 第三個矩陣
一正一負鞍點一個方向縮小、一個方向放大Demo B3 的鞍點
複數，實部負螺旋收斂邊轉邊衰減（見下）阻尼振盪
複數，實部正螺旋發散邊轉邊放大
純虛數中心只轉不縮無阻尼振盪

這也解釋了長期行為：時間久了，\(e^{\lambda t}\) 中 \(\lambda\) 最大的那一項會壓過其他項，粒子的方向會貼近對應的特徵向量。

特徵值是複數時：為什麼是旋轉

第 5 節的旋轉 \(A = \begin{bmatrix}0&-1\\1&0\end{bmatrix}\) 特徵方程是 \(\lambda^2 + 1 = 0\)，特徵值 \(\pm i\)，也沒有實數的特徵方向，符合「每個方向都會轉彎」。更一般地看：

$$A = \begin{bmatrix}\sigma&-\omega\\\omega&\sigma\end{bmatrix},\qquad\lambda = \sigma\pm i\omega$$

它的解是（可以代回去驗算）：

$$\begin{aligned}x(t) &= e^{\sigma t}(x_0\cos\omega t - y_0\sin\omega t)\\y(t) &= e^{\sigma t}(x_0\sin\omega t + y_0\cos\omega t)\end{aligned}$$

實部 \(\sigma\) 決定半徑怎麼變，虛部 \(\omega\) 決定旋轉速度。這和指數函數一致：由 Euler 公式 \(e^{i\omega t} = \cos\omega t + i\sin\omega t\)，所以 \(e^{(\sigma + i\omega)t} = e^{\sigma t}(\cos\omega t + i\sin\omega t)\)，複數指數本來就同時包含了放大和旋轉。

Demo B6 · 線性系統的相圖

怎麼看：拖動 \(A\) 的四個數字，左圖的相圖會變。右圖把所有 2×2 矩陣攤在「trace–determinant 平面」上，橙點是目前的矩陣：
  • 由 \(\lambda = \frac12(\operatorname{tr}\pm\sqrt{\operatorname{tr}^2 - 4\det})\)，拋物線 \(\det = \operatorname{tr}^2/4\) 以上根號裡是負的，特徵值是複數（旋轉）；以下是實數。
  • \(\det \lt 0\) 時兩個特徵值一正一負（鞍點），因為行列式是特徵值的積。
  • \(\det > 0\) 時，\(\operatorname{tr} \lt 0\) 穩定、\(\operatorname{tr} > 0\) 不穩定，因為 trace 是特徵值的和。
軌跡特徵向量方向右圖：目前的 (tr A, det A)

矩陣指數：把解寫成一個式子

一維的解是 \(x(t) = e^{at}x(0)\)。多維能不能也寫成 \(\mathbf x(t) = e^{tA}\mathbf x(0)\)？可以，只要定義「矩陣的指數」。普通指數的級數是 \(e^{y} = 1 + y + \frac{y^2}{2!} + \cdots\)，把 \(y\) 換成 \(tA\)，1 換成單位矩陣：

$$e^{tA} = I + tA + \frac{(tA)^2}{2!} + \frac{(tA)^3}{3!} + \cdots$$

逐項對 \(t\) 微分，會得到 \(\frac{d}{dt}e^{tA} = A\,e^{tA}\)，和一維完全一樣，所以 \(e^{tA}\mathbf x(0)\) 確實是解。

它和上面的一般解是同一件事。如果 \(\mathbf v\) 是特徵向量，\(A^k\mathbf v = \lambda^k\mathbf v\)，代入級數：

$$\begin{aligned}e^{tA}\mathbf v &= \Big(1 + t\lambda + \frac{(t\lambda)^2}{2!} + \cdots\Big)\mathbf v\\&= e^{\lambda t}\mathbf v\end{aligned}$$

實際上，矩陣指數最常用的算法就是先拆成特徵方向，每個方向乘上自己的 \(e^{\lambda t}\)，再加回來。

一句話帶走：沿特徵向量方向，系統退化成一維的 \(\dot f = \lambda f\)，解是 \(e^{\lambda t}\mathbf v\)；任何起點都是這些方向的組合，所以特徵值決定整張相圖。
小練習：\(A = \begin{bmatrix}1&0\\0&-2\end{bmatrix}\)，起點 \((1, 1)\)。寫出解。這是哪一種相圖？
對角矩陣，所以 \(x(t) = e^{t}\)、\(y(t) = e^{-2t}\)。特徵值 \(1, -2\)，一正一負，是鞍點：粒子先靠近 \(x\) 軸，再沿 \(x\) 軸往外跑。
在 CS236 哪裡用到：第 18 講的離散 diffusion 是線性 ODE \(\frac{d\mathbf p}{dt} = Q\mathbf p\)，解就是 \(\mathbf p_t = e^{tQ}\mathbf p_0\)（補充課 G）；補充課 E 的 OU 過程與 VP SDE，drift 都是線性的，所以能寫出閉式解。

08Euler 法：讓電腦解 ODE

為什麼需要它

速度場是神經網路時沒有公式解，只能用電腦一步一步走。CS236 的取樣，本質上就是在數值解一個 ODE 或 SDE。

方法

第 3 節的表格其實已經用了這個方法：把時間切成長度 \(h\) 的小段，每段都假設速度不變。

Euler 法 $$\mathbf x_{n+1} = \mathbf x_n + h\,\mathbf v(\mathbf x_n, t_n)$$

它就是一階 Taylor 展開（補充課 A 第 1 節）：\(\mathbf x(t + h)\approx\mathbf x(t) + h\,\dot{\mathbf x}(t)\)，丟掉了 \(O(h^2)\) 的項。

誤差有多大

  • 每一步丟掉 \(O(h^2)\) 的誤差，這叫局部誤差。
  • 走到固定時間 \(T\) 要 \(T/h\) 步。誤差大致累積起來：\(\frac{T}{h}\times O(h^2) = O(h)\)，這叫全域誤差。

所以步數加倍（\(h\) 減半），最後的誤差大約只減半。全域誤差是 \(O(h)\) 的方法叫一階方法。第 3 節的表格可以驗證：\(t = 10\) 的誤差在 \(h = 1\) 時是 \(36.79 - 34.87 = 1.92\)，\(h = 0.5\) 時是 \(0.94\)，大約減半。

一個看得到誤差的例子：繞圈

第 5 節的 \(\dot x = -y,\ \dot y = x\)，精確解是單位圓，半徑永遠是 1。Euler 的一步是：

$$\begin{aligned}x_{n+1} &= x_n - hy_n\\y_{n+1} &= y_n + hx_n\end{aligned}$$

算新的半徑平方，交叉項會互相抵銷：

$$\begin{aligned}&x_{n+1}^2 + y_{n+1}^2\\={}&(x_n - hy_n)^2 + (y_n + hx_n)^2\\={}&(1 + h^2)(x_n^2 + y_n^2)\end{aligned}$$

半徑每步乘上 \(\sqrt{1 + h^2} > 1\)。切線永遠指向圓外，所以每一步都走到圓外一點，誤差只會累積，不會抵銷。繞一圈用 \(N = 12\) 步時，\(h = 2\pi/12\)，半徑變成 \((1 + h^2)^{6}\approx4.28\) 倍。

Demo B7 · 繞圈的 ODE：Euler 會越繞越外面

怎麼看：灰色圓是精確解，拖動步數，比較 Euler 的半徑和讀數裡的理論值。勾選後會同時畫出下一節的 Heun 與 RK4，它們幾乎貼著圓。
精確解（單位圓）EulerHeunRK4
一句話帶走：Euler 法每步假設速度不變，局部誤差 \(O(h^2)\)，累積後全域誤差 \(O(h)\)。
小練習：\(\dot x = -2x\)，\(x(0) = 3\)，步長 \(0.25\)。用 Euler 算到 \(t = 0.5\)，和精確值比較。
每步乘上 \(1 - 0.5 = 0.5\)，兩步得 \(3\times0.25 = 0.75\)。精確值 \(3e^{-1}\approx1.104\)。

09Heun 與 RK4：多看幾眼速度場

換個角度看：解 ODE 就是估計面積

把 ODE 兩邊從 \(t\) 積分到 \(t + h\)，就得到完全精確的式子：

$$\mathbf x(t + h) = \mathbf x(t) + \int_t^{t+h}\mathbf v(\mathbf x(s), s)\,ds$$

積分就是速度曲線下的面積，所以每種解法的差別，只在於怎麼估計這塊面積。這和微積分裡的數值積分是同一回事：

Euler：左端點Heun：梯形RK4：拋物線
黑色曲線是一步之內的速度。Euler 用左端點的矩形，Heun 用兩端點的梯形，RK4 用頭、中、尾三處的拋物線。

但有一個困難：要知道終點的速度，得先知道終點的位置，而那正是我們要算的東西。解決的辦法是：先用便宜的方法試走一步，用試走的位置去看速度。

Heun 法（改良的 Euler，二階）

梯形的面積是兩端速度的平均乘上 \(h\)：

$$\begin{aligned}\mathbf k_1 &= \mathbf v(\mathbf x_n, t_n)\\\mathbf k_2 &= \mathbf v(\mathbf x_n + h\mathbf k_1,\ t_n + h)\\\mathbf x_{n+1} &= \mathbf x_n + \frac h2(\mathbf k_1 + \mathbf k_2)\end{aligned}$$

\(\mathbf k_1\) 是起點的速度。\(\mathbf k_2\) 是用 Euler 試走一步之後，在那裡看到的速度。最後用兩者的平均，從起點重新走一步。

為什麼 Heun 的全域誤差是 \(O(h^2)\)

只看 \(\dot x = v(x)\) 的一維情況。精確解的二階 Taylor 展開是 \(x(t+h) = x + hv + \tfrac12h^2\ddot x + O(h^3)\)，而由連鎖律 \(\ddot x = \frac{d}{dt}v(x(t)) = v'v\)。

Heun 的 \(k_2 = v(x + hv)\approx v + hv'v\)，所以一步是 \(x + \frac h2(2v + hv'v) = x + hv + \tfrac12h^2v'v\)。和精確解相同，一直到 \(h^2\) 項都對，局部誤差是 \(O(h^3)\)，累積後全域誤差是 \(O(h^2)\)。

RK4（四階 Runge–Kutta）

多試幾個位置，在頭、中、中、尾各看一次速度，再加權平均：

$$\begin{aligned}\mathbf k_1 &= \mathbf v(\mathbf x_n,\ t_n)\\\mathbf k_2 &= \mathbf v(\mathbf x_n + \tfrac h2\mathbf k_1,\ t_n + \tfrac h2)\\\mathbf k_3 &= \mathbf v(\mathbf x_n + \tfrac h2\mathbf k_2,\ t_n + \tfrac h2)\\\mathbf k_4 &= \mathbf v(\mathbf x_n + h\mathbf k_3,\ t_n + h)\\\mathbf x_{n+1} &= \mathbf x_n + \frac h6(\mathbf k_1 + 2\mathbf k_2 + 2\mathbf k_3 + \mathbf k_4)\end{aligned}$$

權重 \(1, 2, 2, 1\) 從哪裡來？如果速度只和時間有關，\(\mathbf k_2 = \mathbf k_3\)，公式就變成 \(\frac h6(\mathbf k_1 + 4\mathbf k_2 + \mathbf k_4)\)，這正是微積分學過的 Simpson 法則（用拋物線估計面積）。一般情況下全域誤差是 \(O(h^4)\)，證明要展開很多項，這裡不推導，下面的 demo 會直接量出來。

手算：同一步，三種方法

\(\dot x = -x\)，\(x_0 = 1\)，\(h = 0.5\)。精確值 \(e^{-0.5} = 0.60653\)。

Euler：\(1 + 0.5(-1) = 0.5\)，誤差 \(0.107\)。

Heun：\(k_1 = -1\)。試走到 \(1 + 0.5(-1) = 0.5\)，所以 \(k_2 = -0.5\)。\(x_1 = 1 + \frac{0.5}{2}(-1.5) = 0.625\)，誤差 \(0.018\)。

RK4：

  • \(k_1 = -1\)
  • \(k_2 = -(1 - 0.25\cdot1) = -0.75\)
  • \(k_3 = -(1 - 0.25\cdot0.75) = -0.8125\)
  • \(k_4 = -(1 - 0.5\cdot0.8125) = -0.59375\)

括號裡是 \(-1 - 1.5 - 1.625 - 0.59375 = -4.71875\)，所以 \(x_1 = 1 - \frac{0.5}{6}\times4.71875 = 0.60677\)，誤差 \(0.00024\)。

比較時要看 NFE

Heun 每步算 2 次速度場，RK4 算 4 次。在 diffusion model 裡，算一次速度場就是跑一次神經網路，這是最花時間的部分，叫做一次 NFE（number of function evaluations）。所以比較解法時要看「同樣的 NFE 下誰比較準」，而不是同樣的步數。

方法每步 NFE全域誤差步數加倍，誤差變成
Euler1\(O(h)\)\(1/2\)
Heun2\(O(h^2)\)\(1/4\)
RK44\(O(h^4)\)\(1/16\)

Demo B8 · 誤差隨步數下降的速度（log–log 圖）

怎麼看：兩軸都取 log，所以誤差 \(\propto h^p\) 會變成斜率 \(-p\) 的直線。讀數裡的實測斜率大約是 \(-1, -2, -4\)，就是上面的階數。把橫軸改成 NFE，RK4 的優勢會縮小，但在要求高精確度時仍然最好。
EulerHeunRK4
  • Karras et al. (2022) 的 EDM 用 Heun 解 diffusion 的 ODE，在 35 次左右的 NFE 達到很好的品質（出處見延伸閱讀）。
  • 實務上還會自動調整步長：在速度場變化快的地方走小步、平緩的地方走大步（例如 RK45 / dopri5）。第 14 講用 ODE 算 likelihood，Song et al. 用的就是這種求解器。
Picard 迭代：同時猜整段軌跡再修正

本節開頭的積分式還能這樣用：先猜整段軌跡（例如全部等於起點），代入右邊的積分，得到一條新的軌跡，再代入，反覆做。這叫 Picard 迭代。因為每次都可以同時算所有時間點的速度，所以能平行計算。第 17 講的 ParaDiGMS 就是這個想法。

一句話帶走：解 ODE 就是估計速度曲線下的面積；多看幾個位置的速度，誤差階數就更高，但每步要多跑幾次網路。
小練習：\(\dot x = x\)，\(x_0 = 1\)，\(h = 0.1\)。用 Heun 走一步，和 \(e^{0.1}\) 比較。
\(k_1 = 1\)，\(k_2 = 1 + 0.1 = 1.1\)，\(x_1 = 1 + 0.05(2.1) = 1.105\)。\(e^{0.1}\approx1.10517\)，誤差約 \(0.00017\)。

10Stiffness：為什麼 Euler 步長太大會爆掉

測試方程

回到最簡單的 \(\dot x = -\lambda x\)（\(\lambda > 0\)）。精確解 \(e^{-\lambda t}\) 永遠平穩地衰減到 0，這個方程本身非常穩定。但 Euler 的一步是：

$$x_{n+1} = x_n + h(-\lambda x_n) = (1 - h\lambda)\,x_n$$

每步乘上 \(1 - h\lambda\)，所以 \(x_n = (1 - h\lambda)^nx_0\)。會發生什麼，完全看這個係數：

\(h\lambda\)係數 \(1 - h\lambda\)\(x_0, x_1, x_2, x_3, \dots\)行為
0.50.51, 0.5, 0.25, 0.125平順衰減
101, 0, 0, 0一步到 0
1.5−0.51, −0.5, 0.25, −0.125上下振盪，但會衰減
2.5−1.51, −1.5, 2.25, −3.375越振越大，爆掉
0 < hλ < 1平順衰減1 < hλ < 2振盪但衰減hλ > 2越振越大
\(|1 - h\lambda| > 1\) 時每步都放大。\(\lambda > 0\) 時，這個條件就是 \(h\lambda > 2\)。
Euler 的穩定條件 $$h\lambda\le2,\quad\text{也就是}\quad h\le\frac{2}{\lambda}$$

為什麼會這樣？速度正比於 \(x\) 且 \(\lambda\) 很大時，一步會走太遠，直接越過 0 衝到另一邊，而且衝得比原來還遠。精確解是一路減速地靠近 0，但 Euler 假設速度不變，所以煞車不及。

什麼是 stiff

上面的例子解本身也變化很快，用小步長好像理所當然。真正麻煩的是下面這種情況：

$$\frac{dx}{dt} = -\lambda\,(x - \cos t),\qquad\lambda\text{ 很大}$$

規則是「不論現在在哪裡，都被很快地拉向 \(\cos t\)」。解很快就貼上 \(\cos t\)，之後變化非常平順，看起來用大步長就夠了。但 Euler 的一步是：

$$x_{n+1} = (1 - h\lambda)x_n + h\lambda\cos t_n$$

係數還是 \(1 - h\lambda\)。只要 \(h\lambda > 2\)，和 \(\cos t\) 的誤差就會每步放大。解明明很平順，步長卻被一個早就結束的快速過程綁住，這種情況就叫 stiff。多維系統也一樣：特徵值中最大的 \(|\lambda|\) 決定 Euler 能用多大的步長，就算那個方向的分量早就衰減完了。

第 3 節的 Demo B1 已經看過一半的現象：\(kh > 1\) 時一步衝過 0，變成負的。補充課 F 的 Langevin 取樣也有同樣的步長限制。

一句話帶走：Euler 每步把誤差乘上 \(1 - h\lambda\)，所以 \(h\lambda > 2\) 就爆掉；解很平順卻被迫用小步長，就叫 stiff。
小練習：\(\dot x = -100x\)。Euler 最大能用多大的步長？
\(h\le2/100 = 0.02\)。想要平順不振盪，要 \(h\le0.01\)。

11指數積分器：把好解的部分精確解掉，以及 DDIM

想法

Stiff 的麻煩來自 \(-\lambda x\) 這一項，但它正好是我們會精確解的（第 3 節）。那能不能把這部分精確解掉，只對剩下的部分做近似？這類方程叫半線性（semilinear）：

$$\frac{dx}{dt} = -\lambda x + g(t)$$

一個線性項，加上一個其他的東西 \(g\)。上一節的例子就是 \(g(t) = \lambda\cos t\)。

積分因子：一步一步解

直接分離變數行不通，因為右邊有兩項。關鍵的技巧是換一個變數，讓線性項消失。定義

$$u(t) = e^{\lambda t}x(t)$$

也就是用「精確解的衰減速度」把 \(x\) 放大回來。對 \(u\) 微分，用乘積法則：

$$\frac{du}{dt} = \lambda e^{\lambda t}x + e^{\lambda t}\frac{dx}{dt}$$

把 ODE 代入 \(\frac{dx}{dt}\)：

$$\frac{du}{dt} = \lambda e^{\lambda t}x + e^{\lambda t}\big(-\lambda x + g(t)\big)$$

前後兩個 \(\lambda e^{\lambda t}x\) 剛好抵銷：

$$\frac{du}{dt} = e^{\lambda t}g(t)$$

\(x\) 不見了！右邊只有 \(t\)，所以這是第 2 節的情況，直接積分就好。\(e^{\lambda t}\) 叫積分因子。從 \(t_n\) 積到 \(t_n + h\)：

$$\begin{aligned}&e^{\lambda(t_n + h)}x(t_n + h) - e^{\lambda t_n}x(t_n)\\={}&\int_{t_n}^{t_n + h}e^{\lambda s}g(s)\,ds\end{aligned}$$

兩邊除以 \(e^{\lambda(t_n + h)}\)：

精確的一步（沒有任何近似） $$\begin{aligned}x(t_n + h) ={}& e^{-\lambda h}x(t_n)\\&+ \int_{t_n}^{t_n + h}e^{-\lambda(t_n + h - s)}g(s)\,ds\end{aligned}$$

第一項是線性部分的精確解，第二項是 \(g\) 累積的影響，而且越早的貢獻衰減得越多。

唯一的近似

\(g\) 一般沒辦法積分（它可能是神經網路），所以假設它在這一小段內是常數 \(g(t_n)\)。積分就能算了：

$$\int_{t_n}^{t_n + h}e^{-\lambda(t_n + h - s)}\,ds = \frac{1 - e^{-\lambda h}}{\lambda}$$
指數積分器 $$x_{n+1} = e^{-\lambda h}x_n + \big(1 - e^{-\lambda h}\big)\frac{g(t_n)}{\lambda}$$

這個式子有很好的解讀。\(g\) 固定時，\(x\) 的平衡點（速度為 0 的位置）是 \(g/\lambda\)，也就是「目標」。令 \(w = e^{-\lambda h}\)，更新式就是

$$x_{n+1} = w\,x_n + (1 - w)\,\frac{g}{\lambda}$$

\(w\) 永遠在 0 和 1 之間，所以 \(x_{n+1}\) 永遠是現在位置和目標的加權平均，一定落在兩者之間，不可能衝過頭。不論 \(h\) 多大都不會爆。把 \(g(t) = \lambda\cos t\) 代入，就是 \(x_{n+1} = e^{-\lambda h}x_n + (1 - e^{-\lambda h})\cos t_n\)。

x目標 g/λx目標 g/λEuler 走過頭了（hλ = 2）指數積分器停在 x 和目標之間
同一步 \(h\lambda = 2\)。Euler 走到目標的另一邊，再大一點就會爆掉；指數積分器永遠落在中間。
手算

\(\dot x = -10x + 10\)，\(x_0 = 0\)，步長 \(h = 0.25\)（\(h\lambda = 2.5\)）。目標是 \(g/\lambda = 1\)，精確解是 \(1 - e^{-10t}\)。

步精確Euler指數積分器
10.9182.50.918
20.993−1.250.993
30.9994.3750.999
41.000−4.06251.000

Euler 每步：\(x_1 = 0 + 0.25(10) = 2.5\)，接著 \(x_2 = 2.5 + 0.25(-25 + 10) = -1.25\)。指數積分器每步：\(w = e^{-2.5} = 0.0821\)，\(x_1 = 0.0821\cdot0 + 0.9179\cdot1 = 0.918\)。這裡 \(g\) 真的是常數，所以它和精確解完全相同。

Demo B9 · 同一個 stiff ODE：Euler vs. 指數積分器

怎麼看：方程是 \(\dot x = -\lambda(x - \cos t)\)。把步數往下拉，看讀數裡的 \(h\lambda\)。在 1 和 2 之間，Euler 上下鋸齒狀地振盪；一超過 2，Euler 就爆掉。指數積分器始終貼著精確解。
精確解Euler指數積分器cos t

把 \(g\) 用更高次的多項式近似（不只是常數），就得到高階的指數積分器。

和 DDIM 的關係

現在可以看懂第 17 講投影片 33「exponential integrator (semi-linear ODE)」的意思了。VP diffusion 的 probability flow ODE 正好是半線性的：

$$\frac{d\mathbf x}{dt} = -\tfrac12\beta(t)\mathbf x + (\text{含神經網路的部分})$$
本節diffusion
線性項 \(-\lambda x\)\(-\tfrac12\beta(t)\mathbf x\)，可以精確解
其他部分 \(g(t)\)含神經網路的部分，只能近似
積分因子 \(e^{\lambda t}\)因為 \(\beta\) 隨時間改變，換成 \(\exp\big(\tfrac12\int\beta\,dt\big)\)，步驟相同

在 \(t\) 接近 0 時 \(\beta\) 變化很劇烈，直接用 Euler 需要很小的步長。DDIM 等於把線性部分精確解掉、只對神經網路的部分做近似；後來的 DPM-Solver、DEIS 把它推廣成高階的指數積分器。這是 DDIM 只要 10–50 步的主要原因之一。另外，第 17 講 Demo A 把步數集中在雜訊小的尾段（二次間隔），也是在處理同一件事：在變化快的地方多走幾步。

出處：Lai et al. (2026) 附錄 A.1.3「Semilinear ODEs and Exponential Integration Factors」；第 17 講頁面 §3「DDIM 為什麼是一個好的 ODE 解法」。

一句話帶走：乘上積分因子讓線性項消失，再只對剩下的部分做近似；更新式是現在位置和目標的加權平均，所以永遠不會爆。
小練習：第 3 節提到 \(\frac{dx}{dt} = \sin(2t) - x\) 無法分離變數。用積分因子解它的第一步：\(u = ?\)，\(\frac{du}{dt} = ?\)
這裡 \(\lambda = 1\)、\(g = \sin 2t\)。\(u = e^{t}x\)，\(\frac{du}{dt} = e^{t}\sin 2t\)，右邊只有 \(t\)，可以直接積分。

12在 CS236 哪裡用到

概念出現的地方
速度場推動粒子第 14 講的 probability flow ODE、第 16 講「ODE 把 VAE 變成 flow」、flow matching
解唯一、軌跡不交叉 → 可逆第 14、17 講：diffusion 的 ODE 是一個 normalizing flow，可以精確算 likelihood（補充課 C）
線性 ODE、矩陣指數第 18 講的 \(e^{tQ}\)；VP / VE SDE 的閉式 \(q(\mathbf x_t\mid\mathbf x_0)\)
Euler、Heun、RK、NFE第 16–17 講：DDPM 是反向 SDE 的 Euler–Maruyama、DDIM、Heun（EDM）、自適應步長的 likelihood 計算
指數積分器、stiffness第 17 講的 DDIM、DPM-Solver、二次時間間隔
Picard 迭代第 17 講的平行求解 ParaDiGMS

常見誤解

誤解：「步數加倍，誤差一定減半」。那只對一階方法（Euler）成立。Heun 會變成 1/4，RK4 會變成 1/16（Demo B8）。
誤解：「RK4 一定比 Euler 好」。要用同樣的 NFE 比較。精確度要求很低、或速度場很不平滑時，高階方法的優勢會縮小。
誤解：「流到 \(T\) 之後粒子擠在一起，資訊就丟了，所以不可逆」。精確解永遠可逆，只是被壓得很密（Demo B4）。實務上的不可逆來自數值誤差和浮點精度。
誤解：「Euler 爆掉是因為 ODE 本身不穩定」。Demo B9 的 ODE 本身非常穩定（被快速拉回 \(\cos t\)），爆掉的是數值方法：步長乘上 \(\lambda\) 超過 2。
誤解：「ODE 取樣沒有隨機性，所以不能產生多樣的樣本」。隨機性全部在起點：從不同的雜訊出發，就得到不同的樣本（第 17 講「像 reparameterization」那段）。
誤解：「特徵值只有兩個，所以線性系統只有兩條軌跡」。軌跡有無限多條，每個起點一條。只是它們都由同樣的兩個 \(e^{\lambda t}\) 組合而成，所以形狀有共同的規律。

自我測驗

1. \(\dot x = -2x\)，\(x(0) = 3\)。\(x(1)\) 是多少？用 Euler、步長 0.5 算兩步得到多少？
精確 \(3e^{-2}\approx0.406\)。Euler：\(x_1 = 3(1 - 1) = 0\)，\(x_2 = 0\)。\(h\lambda = 1\)，一步就把解打到 0。
2. 為什麼兩條 ODE 軌跡不能交叉？這對 diffusion model 有什麼用？
若相交，從交點出發會有兩個不同的未來，違反解的唯一性。因此從時間 0 到 \(T\) 的映射是一對一、可逆的，diffusion 的 probability flow ODE 可以當作 normalizing flow，算精確的 likelihood。
3. \(\dot{\mathbf x} = A\mathbf x\)，\(A = \begin{bmatrix}1&2\\0&-1\end{bmatrix}\)。特徵值是什麼？相圖是哪一種？
上三角矩陣的特徵值就是對角線：\(1\) 和 \(-1\)。一正一負，是鞍點。也可以用 \(\det A = -1 \lt 0\) 直接判斷。
4. 上一題從 \(\mathbf v = (1, 0)\) 出發的解是什麼？
\(A(1, 0) = (1, 0)\)，所以 \((1, 0)\) 是特徵值 1 的特徵向量，解是 \(e^{t}(1, 0)\)：沿 \(x\) 軸往外跑。
5. 某解法把步數加倍後誤差變成 1/4。它是幾階？RK4 呢？
\(2^{-p} = 1/4\)，所以 \(p = 2\)，二階。RK4 是四階，誤差變成 1/16。
6. 為什麼 DDIM 這類方法可以用很大的步長？
VP diffusion 的 ODE 是半線性的。DDIM 把線性部分精確解掉，只對神經網路的部分做近似，所以不受線性項的 stiffness 限制。

延伸閱讀