M3.2 一千片葉子在亂流裡：Fokker–Planck 方程式

本篇重用M0.1閘門的進出：Divergence 與通量·M3.1水流加亂流：SDE 是加了噪聲的 ODE·M2.2站內人數變化 = 進 − 出：Continuity Equation·M0.2導航的三十秒：Taylor 展開與「假設直線」

一千片葉子，有亂流的河

M2.2 問過：一群葉子在沒有亂流的河裡，分佈怎麼演化？答案是 continuity equation——不追葉子，每個位置記帳「進 − 出」。M3.1 則把亂流加回河裡，但只追了一片葉子，看到它的路變成一把扇子。

現在兩件事合起來：一千片葉子，從橋下一小塊水面出發，河裡既有主流也有亂流。一分鐘後，這一千片葉子的分佈是什麼樣子？

先用直覺回答，並寫下你有多確定。多數人會說「整團被主流帶到下游，同時被亂流攤開、變寬變扁」。這個描述是對的。但它沒有回答：變寬的速度是多少？主流不均勻（有地方快有地方慢）時，「搬」與「攤」會怎麼互相影響？以及——能不能像 continuity equation 那樣，寫下一條分佈自己滿足的方程式，不用模擬一千片葉子就算出答案？

課堂提問Q1

翻成數學。「一千片葉子的分佈」是什麼物件？把 continuity equation 的「進 − 出」記帳法套到有亂流的河上，多出了什麼？

先想一想，再展開看整理後的答案

這一題可以直接對上前面的語言：一千片葉子是從 $p_0$ 抽出的樣本，每片各自照 SDE $dx=f\,dt+g\,dW$ 走（M3.1），一分鐘後它們是 $p_{60}$ 的樣本。我們要的是密度 $p(x,t)$ 隨時間的方程式。

記帳法還能用，但一塊地磚上的人數變化現在有兩種來源。第一種和以前一樣：drift $f$ 讓葉子以確定的速度穿過地磚的邊，淨流入是 $-\nabla\cdot(fp)$ 。第二種是新的：亂流讓葉子隨機跨過邊界——隔壁地磚葉子多，隨機跳過來的就多；這塊葉子多，隨機跳出去的也多。所以亂流造成的淨流入，正比於「隔壁比我多多少」，也就是密度的二階導數（M0.2：一點與它兩側平均值的差，就是二階導數）。密度凹下去的地方（局部低谷）會被填、凸起來的地方（局部高峰）會被削——這就是「攤平」。

所以要找的方程式長得像： $\partial_tp=(\text{drift 搬})+(\text{噪聲攤})$ ，第一項是 $-\nabla\cdot(fp)$ ，第二項與 $\Delta p=\nabla\cdot\nabla p$ （Laplacian，二階導數的總和）成正比，係數要算。

先看畫面：一坨顏料在流水裡

把一千片葉子想成一滴顏料滴進河裡。主流把整滴顏料搬到下游——形狀跟著水流變形，這是 M2.2 已經處理過的事。亂流則讓顏料滲開：邊緣越來越淡、範圍越來越大，中心的濃度越來越低，但顏料的總量不變。

兩件事同時發生：被搬走的同時也在滲開。關鍵是它們可以分開算再加起來——因為在一小步 $h$ 內，drift 的位移是 $fh$ 、噪聲的位移是 $g\sqrt h\,\varepsilon$ ，兩者對密度的影響到 $h$ 的一階是各自獨立的。

示意圖 · 規劃中搬與攤：Fokker–Planck 的兩項各做什麼

檔名：m3-stochastic-dynamics/imgs/m3-2-1.svg，16:9。
版面：三欄，各畫一維密度曲線 $p(x)$ 在 $t$ 與 $t+h$ 的疊圖（實線與虛線）。
左欄「只有 drift（ $g=0$ ）」：一個鐘形往右平移，形狀不變（ $f$ 常數）；下方標「 $\partial_tp=-\partial_x(fp)$ ：搬」。若 $f$ 隨 $x$ 變則附一個小圖顯示鐘形被壓窄或拉寬。
中欄「只有噪聲（ $f=0$ ）」：鐘形原地變寬變矮，兩側虛線高於實線、中央虛線低於實線；在三個點畫小箭頭標出 $\partial_{xx}p$ 的正負（中央負→下降，兩側正→上升）；下方標「 $\partial_tp=\tfrac{g^2}{2}\partial_{xx}p$ ：攤」。
右欄「兩者相加」：鐘形往右移並同時變寬；標「 $\partial_tp=-\partial_x(fp)+\tfrac{g^2}{2}\partial_{xx}p$ 」。
底部一行：三欄的曲線下面積都標 $=1$ 。
圖說：「drift 搬質量（一階導數項），噪聲攤平質量（二階導數項）；在一小步內兩者可以分開算再相加。總質量守恆。」

寫成數學：一步之內，搬與攤各貢獻多少

先算攤的係數。一維、 $f=0$ ：在一小步 $h$ 內每片葉子被隨機推 $g\sqrt h\,\varepsilon$ ，所以新密度是舊密度與 $\mathcal N(0,g^2h)$ 的 convolution：

p(x,t+h)=\mathbb E_\varepsilon\big[p(x-g\sqrt h\,\varepsilon,\,t)\big].

把 $p(x-y)$ 對 $y$ Taylor 展開（M0.2）到二階、再對 $\varepsilon$ 取期望： $\mathbb E[\varepsilon]=0$ 殺掉一階項， $\mathbb E[\varepsilon^2]=1$ 留下二階項：

p(x,t+h)=p-\,g\sqrt h\,\mathbb E[\varepsilon]\,\partial_xp+\tfrac{g^2h}{2}\,\mathbb E[\varepsilon^2]\,\partial_{xx}p+O(h^{3/2}) =p+\tfrac{g^2h}{2}\,\partial_{xx}p+O(h^{3/2}).

兩邊減 $p$ 、除以 $h$ 、令 $h\to0$ ： $\partial_tp=\tfrac{g^2}{2}\partial_{xx}p$ 。這就是「攤」的係數： $g^2/2$ 。 $\square$

搬的部分就是 continuity equation（M2.2）： $\partial_tp=-\partial_x(fp)$ 。一步之內兩者的貢獻都是 $O(h)$ ，交叉影響是更高階的小量，所以相加：

\boxed{\;\partial_tp=-\partial_x(f\,p)+\tfrac{g^2}{2}\,\partial_{xx}p\;}\qquad \text{（多維：}\ \partial_tp=-\nabla\cdot(f\,p)+\tfrac{g^2}{2}\,\Delta p\text{）}

這就是 Fokker–Planck 方程式（也叫 forward Kolmogorov equation）[1–4]。它是 SDE 的「守格子」版本：M3.1 追一片葉子，這裡直接演化整個分佈——正如 M2.2 之於 M2.0。

代入驗證的技巧。 解 PDE 很難，但驗證一個候選解很容易：把它代進兩邊，看是否相等。純 Brownian motion（ $f=0$ 、 $g=1$ ）從原點出發，M3.0 說 $p_t=\mathcal N(0,t)$ 。驗證：

p=\frac{1}{\sqrt{2\pi t}}e^{-x^2/2t},\qquad \partial_xp=-\frac{x}{t}\,p,\qquad \partial_{xx}p=\Big(-\frac1t+\frac{x^2}{t^2}\Big)p,\qquad \partial_tp=\Big(-\frac{1}{2t}+\frac{x^2}{2t^2}\Big)p .

右邊 $\tfrac12\partial_{xx}p=\big(-\tfrac{1}{2t}+\tfrac{x^2}{2t^2}\big)p=\partial_tp$ 。✓ 兩邊逐項相等， $\mathcal N(0,t)$ 確實是解。這個技巧會反覆用到：不推導、只驗證。

展開細節有 drift 的驗證：f = c 常數、g 常數，p_t = N(ct, g²t)

記 m = ct、s = g²t，p = (2πs)^(−1/2) exp(−(x−m)²/(2s))。 ∂ₓp = −((x−m)/s)·p；∂ₓₓp = (−1/s + (x−m)²/s²)·p。 ∂ₜp：對 log p 微分，m′ = c、s′ = g²，得 ∂ₜp = [−g²/(2s) + c(x−m)/s + g²(x−m)²/(2s²)]·p。右邊：−∂ₓ(cp) = c(x−m)/s·p；(g²/2)∂ₓₓp = [−g²/(2s) + g²(x−m)²/(2s²)]·p。相加恰為 ∂ₜp。✓ 這也是 SDE 那一篇「x_T ~ N(x₀+cT, g²T)」的 PDE 版證明。

順帶一提：(g²/2)Δp = ∇·((g²/2)∇p) = ∇·(p·(g²/2)∇log p)，所以 Fokker–Planck 也可以寫成 continuity equation 的形式 ∂ₜp + ∇·(p·u) = 0，u = f − (g²/2)∇log p。同一條 p_t，可以由「SDE 的一群葉子」產生，也可以由「ODE 的一群葉子」產生——這是 continuity equation 那一篇「同一分佈路徑有多個速度場」的又一個例子，只是這次其中一個速度場藏在噪聲裡。

回到河面：直覺對了，多了三個數字

起點問題的答案：分佈的中心被 $f$ 帶著走，寬度的 variance 以 $g^2$ 的速率線性增加（標準差按 $g\sqrt t$ 長），總質量守恆。直覺的「搬過去同時攤開」完全正確；它補不上的是攤開的速率 $g^2/2$ 、以及「攤開」的精確意思——密度的 Laplacian 為正的地方上升、為負的地方下降，所以高峰被削、低谷被填，多峰的分佈會被抹成單峰。

主流不均勻時兩項會互動：收攏的 drift（ $\partial_xf<0$ ）把分佈壓窄，噪聲把它攤寬，最後可能達到平衡——下一篇（M3.3）就是這件事。

依賴的假設：每一步的噪聲零均值、variance 有限，而且更高階的矩不會太大（Taylor 展開裡 $\mathbb E[\varepsilon^3]$ 以上的項才會消失得夠快。Gaussian 只是最常見的例子，不是必要條件——下面的 quiz c 會問這件事；真正會壞的是重尾噪聲，偶爾一跳很遠、 $\mathbb E[\varepsilon^2]$ 無限，二階項就不夠描述），以及 $g$ 不依賴 $x$ （否則二階項變成 $\tfrac12\partial_{xx}(g^2p)$ ， $g$ 要進到微分裡面，見 M3.1 展開框的 Itô 討論）。

回到河面：直方圖對上 PDE，再換一種噪聲

import numpy as np
rng = np.random.default_rng(0)
c, g, T, h = 0.5, 0.3, 60.0, 0.1
x = np.zeros(20000)
for _ in range(int(T/h)):                      # 一千片（這裡兩萬片）葉子各走 SDE
    x = x + h*c + g*np.sqrt(h)*rng.standard_normal(x.size)
hist, edges = np.histogram(x, bins=60, density=True)
mid = 0.5*(edges[1:] + edges[:-1])
p_theory = np.exp(-(mid - c*T)**2/(2*g**2*T))/np.sqrt(2*np.pi*g**2*T)   # Fokker–Planck 的解
print(np.abs(hist - p_theory).max())          # 小，且隨葉子數增加而縮小
# 違反假設：把 standard_normal 換成 standard_cauchy，直方圖尾巴變厚、方差不收斂

直方圖與 Fokker–Planck 的解 $\mathcal N(cT,g^2T)$ 重合。把噪聲換成 Cauchy（重尾），直方圖的尾巴遠比 Gaussian 厚、中央反而更尖，而且加更多葉子也不會收斂到任何 Gaussian——「攤平」不再由二階導數描述。

互動 demo · 規劃中 Fokker–Planck：葉子的直方圖 vs. 方程式的解

位置：本節之後。這是本 class 的核心互動。檔名：m3-stochastic-dynamics/m3-2-1.html。
版面：上半是一維的密度圖，兩條曲線疊在一起—— $N$ 片葉子照 Euler–Maruyama 走出來的直方圖（柱狀）、以及 Fokker–Planck 用有限差分往前解的 $p(x,t)$ （實線）。下半是俯視的河，畫出葉子的位置點雲。共用時間滑桿 $t\in[0,60]$ 。
控制項：滑桿 $N$ （ $100$ – $20000$ ，對數）；滑桿「主流 $c$ 」；滑桿「亂流 $g$ 」（可調到 $0$ ，退回 continuity equation）；下拉「drift 型」：等速、收攏 $f=-0.1(x-x^*)$ 、雙井（兩個匯合點）；下拉「噪聲分佈」：Gaussian、均勻（同 variance）、Cauchy（重尾）；下拉「初始分佈 $p_0$ 」：窄 Gaussian、雙峰。
雙面板開關「分開看兩項」：把上半拆成三條細曲線——只搬（ $g=0$ 的解）、只攤（ $c=0$ 的解）、兩者相加的解，讀者拖 $t$ 看三條線的關係。
Hover 上半任一 $x$ ：顯示該處的 $p$ 、 $-\partial_x(fp)$ 、 $\tfrac{g^2}{2}\partial_{xx}p$ 三個數值與它們的和 $\partial_tp$ ，並在點上畫一支小箭頭表示密度此刻在上升或下降。
行為：Gaussian noise 下， $N$ 越大直方圖越貼實線； $g=0$ 時實線退回 continuity equation 的解、直方圖不攤開；選雙峰 $p_0$ 並拖 $t$ ，看兩峰被抹成一峰、右上角的「峰數」從 2 變 1；選均勻噪聲時直方圖仍貼實線（只要 variance 相同，二階項一樣——Gaussian 不是必要的，有限 variance 才是）；選 Cauchy 時直方圖尾巴脫離實線並亮警示「 $\mathbb E[\varepsilon^2]=\infty$ ，二階項不夠描述」。
底部觀察句：「觀察：drift 搬、噪聲攤，兩項相加就是分佈的演化，不必追任何一片葉子。攤的速率是 $g^2/2$ 倍的二階導數——高峰被削、低谷被填。噪聲只要 variance 有限就行；重尾噪聲讓方程式失效。」
匯出 PNG：上下兩半併成 16:9，含當前 $t$ 、 $N$ 、 $c$ 、 $g$ 與噪聲型。

先消化一下

想一想

一杯靜水裡滴一滴墨水（沒有攪動）。三十秒後墨水的濃度分佈，該用哪個工具描述？

想一想

你把 $p_t=\mathcal N(0,t)$ 代進 $\partial_tp=\tfrac12\partial_{xx}p$ 驗證成功。若改用 $p_t=\mathcal N(0,t^2)$ （標準差按 $t$ 長）代進去：

想一想

把 SDE 的 Gaussian noise 換成同樣 variance 的均勻分佈噪聲（每步 $g\sqrt h\,\varepsilon$ ， $\varepsilon$ 均勻於 $[-\sqrt3,\sqrt3]$ ），Fokker–Planck 方程式：

參考文獻

Särkkä, S., Solin, A. Applied Stochastic Differential Equations. Cambridge University Press, 2019.（第 5 章：Fokker–Planck 方程式的推導與 Gaussian 情形的 closed form。）
Pavliotis, G. A. Stochastic Processes and Applications: Diffusion Processes, the Fokker–Planck and Langevin Equations. Springer, 2014.（第 4 章：forward Kolmogorov 方程式與其性質。）
Risken, H. The Fokker–Planck Equation: Methods of Solution and Applications. Springer, 2nd ed., 1989.（專書；第 4 章推導，第 5 章一維解法。）
Gardiner, C. Stochastic Methods: A Handbook for the Natural and Social Sciences. Springer, 4th ed., 2009.（第 3、4 章：從 Kramers–Moyal 展開到 Fokker–Planck；本篇 Taylor 展開的來源。）