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

物理·数学·科研 · 2023-08-21 · 2607 人浏览

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 x + 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 积分,我们需要使用伊藤引理。它非常适合描述“无法预测未来”的系统,即在这样的系统中我们只能用过去事件的信息对未来做出推断,但无法准确预测,例如股票等。出于此特点,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(x)^2A''(x) = \int\mathrm d x~[\rho(x, t)g(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. $$

随机过程 随机分析
Theme Jasmine by Kent Liao