随机分析、Langevin 方程与 Fokker-Planck 方程

1. Wiener 过程

首先考虑一个离散的一维随机漫步 (random walk)，记过程的起始点为 $x=0$，每一步行走距离为 $\Delta x$，向左向右走的概率均为 $1/2$，粒子每过 $\Delta t$ 的时间就移动一次。用 $X(t)$ 记录粒子在 $t$ 时刻的位置，则有

$$X(t) = \sum_{i=1}^{t/\Delta t} \eta_i,$$

其中

$$\eta_i = \begin{cases} +\Delta x, \text{ with prob. 1/2 }; \\ -\Delta x, \text{ with prob. 1/2 }. \end{cases}$$

由中心极限定理可知，当 $t$ 很大、或者说 $\Delta t$ 很小的时候，$X(t)$ 的分布趋近于高斯分布。此时，我们只需要知道 $X(t)$ 的均值和方差即可确定 $X(t)$ 的概率分布。均值很容易计算，即

$$\langle X(t) \rangle = \frac{t}{\Delta t}\langle\eta_i\rangle = 0;$$

其方差为

$$\sigma^2_{X(t)} = \langle X(t)^2\rangle = \frac{t}{\Delta t}\sigma^2_{\eta_i} = \frac{(\Delta x)^2}{\Delta t}t.$$

因此有 $X(t)\sim \mathcal N(0, \sqrt{(\Delta x)^2 t / \Delta t})$，由此我们可以计算很多与该过程有关的性质。
实际上，不严谨地说，Wiener 过程可以看作上述离散随机游走的连续极限，即取 $\Delta x \to 0$，$\Delta t \to 0$ 。为了使 $(\Delta x)^2 / \Delta t$ 为一个有限值，$\Delta x$ 需要与 $\sqrt{\Delta t}$ 同阶。记 $\Delta x^2 / \Delta t = \sigma^2$，由此得到的连续过程即为 Wiener 过程 $W(t) \sim \mathcal N(0, \sigma\sqrt{t})$ 。
从数学角度严格地说，Wiener 过程定义为满足下面三个性质的随机过程：

1. $W(t) \sim \mathcal N(0, \sigma\sqrt{t})$；
2. $W(t)$ 是马尔可夫过程；
3. 对于任意 $t' \ge 0$， $W(t+t') - W(t)$ 与 $W(t)$ 相互独立。

条件 (3) 称为独立增量，它描述了过程的空间均一性，即过程与起点无关。
由 Wiener 过程的定义即可计算 Wiener 过程的时间关联函数。由性质 (1)、(3) 可知 $\langle [W(t+t') - W(t)]W(t) \rangle = 0$，因此

$$\langle W(t+t')W(t)\rangle = \langle W(t)^2\rangle = \sigma^2t.$$

故 Weiner 过程在时间 $t_1$ 、$t_2$ 的时间关联函数为

$$\langle W(t_1)W(t_2)\rangle = \sigma^2\min\{ t_1, t_2 \}.$$

这一结论也可以从离散的随机游走看出来，即

$$\begin{aligned} \langle X(t_1)X(t_2) \rangle &= \sum_{i=1}^{t_1/\Delta t}\sum_{j=1}^{t_2/\Delta t}\langle \eta_i\eta_j\rangle \\ &= \Delta x^2\min\left\{ \frac{t_1}{\Delta t}, \frac{t_2}{\Delta t} \right\}. \end{aligned}$$

2. 白噪声

白噪声定义为 Weiner 过程的导数：

$$\xi(t) \equiv \frac{\mathrm d W(t)}{\mathrm d t}.$$

这个导数在数学分析层面无法良好定义，在此仅将其形式上写作这样。计算它的均值可得

$$\langle \xi(t)\rangle = \frac{\mathrm d }{\mathrm d t}\langle W(t)\rangle = 0,$$

时间关联函数为

$$\begin{aligned} \langle\xi(t_{1})\xi(t_{2})\rangle &= \left\langle\frac{\mathrm d}{\mathrm dt_{1}}W(t_{1})\frac{\mathrm d}{\mathrm dt_{2}}W(t_{2})\right\rangle \\ &= \frac{\mathrm d}{\mathrm dt_{1}}\frac{\mathrm d}{\mathrm dt_{2}}\sigma^{2} \min\{t_{1},t\} \\ &= \sigma^2\frac{\mathrm d }{\mathrm d t_1}\Theta(t_1 - t_2) \\ &= \sigma^2\delta(t_1 - t_2). \end{aligned}$$

其中 $\Theta$ 是阶梯函数，它的导数是 $\delta$ 函数。由关联函数可知，$\xi(t_1)$ 与 $\xi(t_2)$ 相互独立 ($t_1 \neq t_2$)，但当 $t_1 = t_2$ 时，$\langle \xi (t)^2 \rangle$ 发散。这样的 $\xi(t)$ 我们称为广义随机过程。
在离散情况下，我们可以算出

$$\frac{X(t') - X(t)}{t'-t} = \frac{1}{t'-t}\sum_{i=1}^{(t'-t)/\Delta t} \eta_i \sim \mathcal N\left(0, \frac{\sigma}{\sqrt{t'-t}}\right).$$

这一形式可以用于模拟白噪声。
我们将这一噪声称为白噪声，是因为其功率谱密度 $S_\xi(\omega)$ 为常数

$$S_\xi(\omega) = \int_{-\infty}^\infty\mathrm d t~e^{-\mathrm i\omega t}\langle\xi(t)\xi(0)\rangle = \sigma^2,$$

这与白光的谱密度形式相同。
我们可以将白噪声看作色噪声的特殊形式。若一个噪声 $Z(t)$ 的时间关联函数不是 $\delta$ 函数，那么该噪声就是色噪声。一般，我们将色噪声 $Z(t)$ 的时间关联函数写为

$$\langle Z(t)Z(t')\rangle = \frac{\sigma^2}{2\tau}e^{-\frac{|t-t'|}{\tau}}.$$

白噪声可以视为色噪声在 $\tau \to 0$ 的极限情况。

3. 随机积分

考虑随机微分方程

$$\dot{x}(t) = f(x(t)) + g(x(t))\xi(t),$$

其中 $\xi(t)$ 为高斯白噪声。由于噪声项的出现，这一方程为随机微分方程。这一方程有微分形式的表达式，即

$$\mathrm d x(t) = f(x(t))\mathrm d t + g(x(t))\mathrm d W(t),$$

其中 $W(t)$ 是 Weiner 过程。
对上述方程积分，形式上可以得到

$$x(t) = x(0) + \int_0^t\mathrm d t'~f(x(t')) + \int_0^t\mathrm d t'~g(x(t'))\mathrm d W(t').$$

这一形式告诉我们，要解随机微分方程，就要先定义随机积分。然而，随机积分并不容易定义。它的定义取决于不同表示方式的选取。
考虑一个一般的随机积分

$$\int_0^\tau \psi(t)\mathrm d W(t).$$

黎曼积分定义为如下求和的极限

$$\sum_{k=1}^{\tau/\Delta t}\psi(t_k^*)[W(t_k) - W(t_{k_1})],$$

其中 $t_k = k\Delta t$，$t_k^* \in [t_{k-1}, t_k]$ 。黎曼积分的良好性质在于，积分的值与 $t_k^*$ 的选取无关，我们可以选区间内的任意一点。然而，对于随机积分，选取不同的 $t_k^*$ 会使得最终的积分结果不同，因此随机积分的结果取决于表示形式。
有两种经常被选用的表示，一种是 Ito 表示，一种是 Stratonovich 表示：

$$\begin{cases} t_k^* = t_{k-1}: \text{ Ito}, \\ t_k^* = \frac{1}{2}(t_{k-1}+t_k) = \left(k-\displaystyle\frac{1}{2}\right)\Delta t: \text{Stratonovich}. \end{cases}$$

我们用一个例子来展示两种表示的不同积分结果。考虑随机积分

$$\int_0^\tau \mathrm d t~W(t)\mathrm d W(t).$$

对于 Ito 表示，我们可以计算积分结果的均值为

$$\begin{aligned} & \left\langle \sum_{k=1}^{\tau/\Delta t}W(t_{k-1})[W(t_{k})-W(t_{k-1})]\right\rangle \\ ={}& \sum_{k=1}^{\tau/\Delta t}\left\langle W(t_{k-1})\right\rangle\left\langle [W(t_{k})-W(t_{k-1})]\right\rangle \\ ={}& 0. \end{aligned}$$

对于 Stratonovich 表示，可得

$$\begin{aligned} & \left\langle \sum_{k=1}^{\tau/\Delta t}W(t_k^*)[W(t_{k})-W(t_{k-1})]\right\rangle \\ ={}& \sum_{k=1}^{\tau/\Delta t}\left\langle W(t_k^*)W(t_{k})\right\rangle - \left\langle W(t_k^*)W(t_{k-1})\right\rangle \\ ={}& \sum_{k=1}^{\tau/\Delta t}(\sigma^2 t_k^* - \sigma^2 t_{k-1}) \\ ={}& \sigma^2\sum_{k=1}^{\tau/\Delta t}(t_k^* - t_{k-1}) \\ ={}& \frac{\sigma^2}{2}\tau. \end{aligned}$$

其中第二个等号用了 $W(t)$ 的时间关联函数表达式。容易看出，这两种表示给出的计算结果是不同的。
现在可以与微积分中的计算方式做一个对比。对于确定函数的微积分，积分结果应当为

$$\begin{gathered} \int_0^\tau W(t)\mathrm d W(t) = \left.\frac{W(t)^2}{2}\right|_{0}^\tau = \frac{W(\tau)^2}{2} \\ \implies \left\langle \frac{W(t)^2}{2} \right\rangle = \frac{\sigma^2}{2}\tau. \end{gathered}$$

实际上，可以证明，Stratonovich 表示与确定函数的微积分给出的结果相同，而 Ito 表示则有所不同。计算 Stratonovich 表示的随机微积分时，我们可以直接使用链式法则。

4. Ito 微积分

为什么传统的微积分对伊藤积分无效？从前面对于离散情况的分析我们已经可以略知一二，即 $(\Delta x)^2$ 与 $\Delta t$ 同阶。
我们可以从离散情况出发推导 Ito 表示下的 $W(t)\mathrm d W(t)$ 是什么。

$$\begin{aligned} {}&W(t) [W(t+\Delta t) - W(t)] \\ ={}&\frac{1}{2}[W(t+\Delta t)^2 - W(t)^2] - \frac{1}{2}[W(t+\Delta t) - W(t)]^2, \end{aligned}$$

而 $\mathrm d W(t)^2 / 2$ 的离散表达式为

$$\frac{\mathrm d W(t)^2}{2} = \frac{1}{2}[W(t+\Delta t)^2 - W(t)^2],$$

若二阶项不能忽略，则它们不一样。我们已经计算出了二阶项的均值为

$$\langle [W(t+\Delta t) - W(t)]^2 \rangle = \sigma^2\Delta t,$$

可以发现，Weiner 过程的二阶项与时间的一阶项同阶，这就是随机分析与普通微积分的不同之处。
伊藤微积分的思想便是，保持式中的二阶项 $(\mathrm d W)^2$，并将其记为一阶项 $\sigma^2\mathrm d t$ 。
由此，我们可以得到随机微分方程的变量代换公式，该公式由伊藤引理给出。令 $x(t)$ 为随机微分方程

$$\mathrm d x = f(x)\mathrm d t + g(x)\mathrm d W$$

在 Ito 表示下的解，那么 $y(t) = y(x(t))$ 是如下方程的解：

$$\begin{aligned} \mathrm d y &= \frac{\mathrm d y}{\mathrm d x}\mathrm d x + \frac{1}{2}\frac{\mathrm d^2y}{\mathrm d x^2}(\mathrm d x)^2 \\ &= \frac{\mathrm d y}{\mathrm d x}[f(x)\mathrm d t + g(x)\mathrm d W] + \frac{1}{2}\frac{\mathrm d^2 y}{\mathrm d x^2}[f(x)\mathrm d t + g(x)\mathrm d W]^2 \\ &= \left[\frac{\mathrm d y}{\mathrm d x}f(x)+\frac{\sigma^2}{2}\frac{\mathrm d^2y}{\mathrm d x^2}g(x)^2\right]\mathrm d t + \frac{\mathrm d y}{\mathrm d x}g(x)\mathrm d W. \end{aligned}$$

上述公式可以简记为

$$\mathrm d y = \left(y'f+\frac{\sigma^2}{2}y''g^2\right)\mathrm d t + y'g\mathrm d W.$$

我们称 $\frac{\sigma^2}{2}y''g^2\mathrm d t$ 为 Ito 项。

5. 随机微分方程

在解随机微分方程之前，我们需要指明选用的表示形式，因为不同的表示会给出不同的积分结果。对于一个一般的随机微分方程

$$\dot x = f(x) + g(x)\xi,$$

为了区分不同的表示方法，通常我们记 Ito 表示为 $g(x)\xi$，记 Stratonovich 表示为 $g(x)\circ \xi$ 。对应的微分表达为

$$\begin{cases} \mathrm d x = f(x)\mathrm d t + g(x)\mathrm d W: \text{Ito interpertation}; \\ \mathrm d x = f(x)\mathrm d t + g(x)\circ \mathrm d W: \text{Stratonovich interpertation.} \end{cases}$$

若我们考虑色噪声，再令 $\tau \to 0$，那么将会得到 Stratonovich 表示下的白噪声。
若 $g(x)$ 为常数，则 Ito 表示和 Stratonovich 表示给出的解相同。Langevin 方程就是 $g(x)$ 为常数的情况，因此以往的处理并不需要区分这两种不同的表示形式。然而，若我们想要在 Langevin 系统上定义热力学量，那么我们需要选择 Stratonovich 表示。
对于 Ito 积分，我们需要使用伊藤引理，它非常适合描述离散时间系统，例如股票等；对于 Stratonovich 积分，我们可以直接使用与普通微积分相同的形式，它适用于描述连续时间系统，例如物理和化学系统。
对于速度为 $v(t)$ 的布朗粒子，其运动方程为（欠阻尼） Langevin 方程

$$m\dot v(t) = -\gamma v(t) + \xi(t).$$

其中 $\xi(t)$ 为高斯白噪声。此时两种表述给出的解相同。
由于没有外力，系统的能量即为粒子的动能

$$E(t) = \frac{1}{2}mv(t)^2.$$

在 Stratonovich 表示下能量满足的运动方程为

$$\dot E(t) = mv(t)\dot v(t) = -\gamma v(t)^2 + v(t)\circ\xi(t).$$

在 Ito 表示下，使用 Ito 引理，运动方程可以写为

$$\dot E(t) = mv\dot v + \frac{\sigma^2}{2}m\left( \frac{1}{m} \right)^2 = -\gamma v^2+\frac{\sigma^2}{2m} + v\xi.$$

可以看出，Stratonovich 表示给出的能量变化率的均值为零，而原 Langevin 方程对于这两种表示会给出相同的结果，因此 Ito 表示下的能量运动方程的均值也应该为零。故有

$$\langle \dot E(t)\rangle = -\frac{2\gamma}{m}\langle E(t)\rangle + \frac{\sigma^2}{2m} = 0.$$

由此即可得到

$$\frac{\sigma^2}{2m} = \frac{2\gamma}{m}\frac{kT}{2} \implies \sigma^2 = 2\gamma kT.$$

这就是涨落耗散定理的最简单形式。
下面我们给出两种表示下随机微分方程之间的关系。若 $x(t)$ 是方程

$$\dot x(t) = f(x)+g(x)\xi$$

的解，那么 $x(t)$ 也是方程

$$\dot{x}(t) = f(x) - \frac{\sigma^2}{2}g(x)g'(x) + g(x)\circ \xi$$

的解。

证明：
令

$$y(x) = \int\frac{\mathrm d x'}{g(x')}, \quad \frac{\mathrm d y}{\mathrm dx} = \frac{1}{g(x)}; \quad \frac{\mathrm d^2 y}{\mathrm d x^2} = -\frac{g'(x)}{g(x)^2}.$$

由伊藤引理可得

$$\begin{aligned} \dot y &= \frac{1}{g}(f+g\xi) + \frac{\sigma^2}{2}\left( -\frac{g'}{g^2} \right)g^2 \\ &= \frac{f}{g} - \frac{\sigma^2}{2}g' + \xi, \end{aligned}$$

再由普通的微积分可得

$$\begin{aligned} \dot x(t) &= g(x)\dot y \\ &= g\left( \frac{f}{g} - \frac{\sigma^2}{2}g' + \xi \right) \\ &= f - \frac{\sigma^2}{2}gg' + g\circ \xi. \end{aligned}$$

反过来也有类似的结论，若 $x(t)$ 是方程

$$\dot x(t) = f(x)+g(x)\circ\xi$$

的解，那么 $x(t)$ 也是方程

$$\dot{x}(t) = f(x) + \frac{\sigma^2}{2}g(x)g'(x) + g(x)\xi$$

的解。注意此定理与上述定理的正负号不同。

6. Fokker-Planck 方程

Fokker-Planck 方程是概率密度 $p(x, t)$ 的演化方程。
首先处理 Ito 表示。对于一般的随机微分方程

$$\dot x(t) = f(x(t)) + g(x(t))\xi(t),$$

考虑一个任意函数 $A(t) = A(x(t))$，使用伊藤引理，$A(t)$ 的运动方程为

$$\dot A(t) = A'(x(t))[f(x(t)) + g(x(t))\xi(t)] + \frac{\sigma^2}{2}A''(x(t))g(x(t))^2.$$

对上述方程求期望可得

$$\langle \dot A(t)\rangle = \langle A' f\rangle + \frac{\sigma^2}{2}\langle A'' g^2\rangle.$$

而 $A(t)$ 的期望为

$$\langle A(x(t)) \rangle = \int\mathrm d x~\rho(x, t)A(x).$$

对上式求时间导数，再对照先前的式子可得

$$\int\frac{\partial \rho}{\partial t}A = \int\mathrm d x~\rho\left( fA' + \frac{\sigma^2}{2}g^2 A'' \right).$$

对上式右边的两项做分部积分。第一项分部积分得到

$$\int_a^b \mathrm d x~\rho(x, t)f(x)A'(x) = -\int\mathrm d x~[\rho(x, t)f(x)]'A(x),$$

其中分部积分得到的第一项由于边界条件而为零。对第二项做两次分部积分，同理可得

$$\int\mathrm d x~\rho(x, t)g^2(x) A''(x) = \int\mathrm d x~[\rho(x, t)g^2(x)]''A(x).$$

将这两项放在一起可得

$$\int \mathrm d x~A(x)\left[ \frac{\partial \rho}{\partial t} + \frac{\partial f\rho}{\partial x} - \frac{\sigma^2}{2}\frac{\partial^2g^2\rho}{\partial x^2} \right] = 0.$$

此式对任意 $A(x)$ 都成立，因此有关于 $\rho(x, t)$ 的方程

$$\frac{\partial \rho(x, t)}{\partial t} = -\frac{\partial}{\partial x}[f(x)\rho(x, t)] + \frac{\sigma^2}{2}\frac{\partial^2}{\partial x^2}[g(x)^2 \rho(x, t)].$$

此式即为伊藤表示下的一维 Fokker-Planck 方程。Stratonovich 表示下的 Fokker-Planck 方程可以由上面的方程变换得到，其结果为

$$\begin{aligned} \frac{\partial \rho(x, t)}{\partial t} =&{} -\frac{\partial}{\partial x}[f(x)\rho(x, t)] + \frac{\sigma^2}{2}\frac{\partial}{\partial x}[g(x)g'(x)\rho(x, t)] \\ &{}+ \frac{\sigma^2}{2}\frac{\partial^2}{\partial x^2}[g(x)^2\rho(x, t)] \\ =&{} -\frac{\partial}{\partial x}[f(x)\rho(x, t)] + \frac{\sigma^2}{2}\frac{\partial}{\partial x}\left[ g(x)\frac{\partial}{\partial x}g(x)\rho(x, t) \right]. \end{aligned}$$

若将 $\partial_x g$ 记为算符，那么上式可以写为

$$\partial_t \rho = -\partial_x(f\rho) + \frac{\sigma^2}{2}(\partial_x g)^2\rho.$$