Fokker-Planck方程

Posted 2020-11-24 immcrr

Fokker-Planck方程

对于一般随机过程有条件概率公式
[ P(y_n,t_n;y_{n-1},t_{n-1};cdots;y_1,t_1)=P(y_n,t_n|y_{n-1},t_{n-1};cdots;y_1,t_1)P(y_{n-1},t_{n-1};y_{n-2},t_{n-2};cdots;y_1,t_1) ]
Markov过程定义为随机变量未来的概率分布只和当前状态有关，和再往前的历史状态无关，即
[ P(y_n,t_n|y_{n-1},t_{n-1};cdots;y_1,t_1)=P(y_n,t_n|y_{n-1},t_{n-1}) ]
对于任意随机过程，满足Kolmogorov定理 ((t_1<t_2<t_3))，
[ P(y_3,t_3|y_1,t_1)=int ext{d}y_2P(y_2,t_2|y_1,t_1)P(y_3,t_3|y_2,t_2;y_1,t_1) ]
对于Markov过程上式化简为
[ P(y_3,t_3|y_1,t_1)=int ext{d}y_2P(y_2,t_2|y_1,t_1)P(y_3,t_3|y_2,t_2) ]
利用上式，可得 ((0<t<t+ au))，
[ P(y,t+ au|y_0,0)=int ext{d}xi P(y-xi,t|y_0,0)P(y,t+ au|y-xi,t) ]
对于小量( au)上式左边泰勒展开得到
[ frac{partial}{partial t}P(y,t|y_0,0)=frac{1}{ au}left[int ext{d}xi P(y-xi,t|y_0,0)P(y,t+ au|y-xi,t)-P(y,t|y_0,0) ight] ]
注意到上式括号内第二项可以写为
[ P(y,t|y_0,0)=int ext{d}xi P(y+xi,t+ au|y,t)P(y,t|y_0,0) ]
代回去得到
[ frac{partial}{partial t}P(y,t|y_0,0)=frac{1}{ au}left[int ext{d}xi P(y-xi,t|y_0,0)P(y,t+ au|y-xi,t)-int ext{d}xi P(y+xi,t+ au|y,t)P(y,t|y_0,0) ight] ]
上式是一个连续性方程：某区域内的概率变化量等于流入概率减流出概率。现在关注括号中第一项，令(z=y-xi)，则该项就是(z)的函数，写为
[ int ext{d}xi P(z,t|y_0,0)P(z+xi,t+ au|z,t) ]
上式在(z=y)处展开为泰勒级数，得到
[ int ext{d}xi sum_{n=0}^infty frac{(z-y)^n}{n!}left[frac{partial^n}{partial z^n}P(z,t|y_0,0)P(z+xi,t+ au|z,t) ight]Bigg|_{z=y} ]
注意上面求导最后在(z=y)点取值，因此可以直接改写为
[ int ext{d}xi sum_{n=0}^infty frac{(z-y)^n}{n!}frac{partial^n}{partial y^n}left[P(y,t|y_0,0)P(y+xi,t+ au|y,t) ight] ]
代回方括号中，得到
[ frac{partial}{partial t}P(y,t|y_0,0)=frac{1}{ au}left{int ext{d}xi sum_{n=0}^infty frac{(z-y)^n}{n!}frac{partial^n}{partial y^n}left[P(y,t|y_0,0)P(y+xi,t+ au|y,t) ight]-int ext{d}xi P(y+xi,t+ au|y,t)P(y,t|y_0,0) ight} ]
其中(n=0)项和括号第二项抵消，从而
[ frac{partial}{partial t}P(y,t|y_0,0)=frac{1}{ au}int ext{d}xi sum_{n=1}^infty frac{(-1)^n}{n!}xi^nfrac{partial^n}{partial y^n}left[P(y,t|y_0,0)P(y+xi,t+ au|y,t) ight] ]
取极限( au o0)，得到
[ frac{partial}{partial t}P(y,t|y_0,0)=sum_{n=1}^infty frac{(-1)^n}{n!}frac{partial^n}{partial y^n}left[P(y,t|y_0,0)int ext{d}xi lim_{ au o0}frac{xi^n}{ au}P(y+xi,t+ au|y,t) ight] ]
如果随机变量(y)满足在时间间隔( au)很小时变化同样很小，则上面的求和保留前几项即可。物理上的一维布朗运动就是一个例子，那里的(p(t))是Markov的，在演化时间( au)很小时，其变化也很小，所以只保留到二阶即可。

而有势能场的一维布朗运动是扩展版，一个粒子的一维坐标、动量组成的向量({x,p})是Markov过程，有对应的Fokker-Planck方程。唯一不同的是，有势能场的布朗运动那里，满足Markov的是向量({x,p})，这涉及到上面方程的二维版本，这只需把泰勒展开换成二元函数的泰勒展开即可。对于二维的情形，对应于(y)的是(x_1,x_2)，对应于(xi)的是(xi_1,xi_2)，从而有
[ frac{partial}{partial t}P(x_1,x_2,t|x_{10},x_{20},0)=sum_{n=1}^infty frac{(-1)^n}{n!}int ext{d}xi_1int ext{d}xi_2 lim_{ au o0}frac{1}{ au}left(xi_1frac{partial}{partial x_1}+xi_2frac{partial}{partial x_2} ight)^nleft[P(x_1,x_2,t|x_{10},x_{20},0)P(x_1+xi_1,x_2+xi_2,t+ au|x_1,x_2,t) ight] ]
应用到布朗运动上，其中(x_1)对应于坐标(x)，(x_2)对应于动量(p)，( au)对应于(delta t)，(xi_1=x(t+ au)-x(t)=x(t+delta t)-x(t))记为(delta x)，同样的(xi_2=delta y)，而且对于时间起点(t=0)记为(t_0)，则泰勒展开保留到二阶给出 (略去极限(lim_{delta t o 0})不写)
[ frac{partial}{partial t}P(x,p,t|x_0,p_0,t_0)=-frac{partial}{partial x}left[frac{overline{delta x}}{delta t}P(x,p,t|x_0,p_0,t_0) ight]-frac{partial}{partial p}left[frac{overline{delta p}}{delta t}P(x,p,t|x_0,p_0,t_0) ight]\+frac{1}{2}frac{partial^2}{partial x^2}left[frac{overline{(delta x)^2}}{delta t}P(x,p,t|x_0,p_0,t_0) ight]+frac{1}{2}frac{partial^2}{partial p^2}left[frac{overline{(delta p)^2}}{delta t}P(x,p,t|x_0,p_0,t_0) ight]\+frac{partial^2}{partial xpartial p}left[frac{overline{(delta x)(delta p)}}{delta t}P(x,p,t|x_0,p_0,t_0) ight] ]

其中

[ overline{(delta x)^2}=int ext{d}(delta x)int ext{d}(delta p)(delta x)^nP(x+delta x,p+delta p,t+delta t|x,p,t) ]

两边同时乘以(P(x_0,p_0,t_0))并对(x_0)和(p_0)积分，则(P(x,p,t|x_0,p_0,t_0))可以替换为(P(x,p,t))。

文献中常见的简记形式为
[ frac{partial}{partial t}P(x_1,x_2,t)=-frac{partial}{partial x_1}mathscr{M}_1P(x_1,x_2,t)-frac{partial}{partial x_2}mathscr{M}_2P(x_1,x_2,t)+sum_{i,j=1,2}frac{partial^2}{partial x_1partial x_2}mathscr{D}_{ij}P(x_1,x_2,t) ]
偏导对后面整体作用，其中(x_1,x_2)分别为(x,p)，且有
[ mathscr{M}_i=frac{overline{delta x_i}}{delta t},mathscr{D}_{ij}=frac{1}{2!}frac{overline{delta x_idelta x_j}}{delta t} ]
上面就是一维布朗运动中常见的Fokker-Planck方程的原型，其中的各个系数需要根据运动方程来确定。

外势场为零的一维Langevin方程为
[ dot{p}(t)=-gamma p(t)+F(t) ]
其中(p(t))是动量，(F(t))是随机力。有阻尼项的存在，动量(p(t))有衰减，(T_ ext{R}=1/gamma)是描述(p(t))衰减的特征时间，即在该时间尺度上，(p(t))有显著衰减。假设这里的随机力由粒子和周围流体微粒碰撞而造成，碰撞发生在时间尺度( au_c)上，满足( au_cll T_ ext{R})。在大于( au_c)的时间尺度上，这里面(p(t))是Markov的，下一时刻(t+ ext{d}t)时的(p)仅和目前的(p)值相关，而和过去的(p)值无关。首先随机力(F(t))是零均值的，即(langle F(t) angle=0)，其中均值是对系综求均值，或者说是对各可能的随机力函数求均值。并且认为随机力在大于( au_c)尺度上是无记忆的，即如果(t'-tgg au_c)，则相关函数(langle F(t)F(t') angle=2Ddelta(t-t'))，这一点可以用来化简问题。直接求解Langevin方程给出
[ p(t)=p_0 ext{e}^{-gamma(t-t_0)}+int_{t_0}^t ext{d}t'F(t') ext{e}^{-gamma(t-t')} ]
由上式可以得到
[ langle p(t) angle=p_0 ext{e}^{-gamma(t-t_0)} ]
由上面两式可以得到
[ sigma_p^2(t)=langle[p(t)-langle p(t) angle]^2 angle=int^t_{t_0} ext{d}t'int_{t_0}^t ext{d}t''langle F(t)F(t') angle ext{e}^{-gamma(t-t')} ext{e}^{-gamma(t-t'')} ]
利用(langle F(t)F(t') angle=2Ddelta(t-t'))，得到
[ sigma_p^2(t)=frac{D}{gamma}left[1- ext{e}^{-2gamma(t-t_0)} ight] ]
对于( au_cll t-t_0ll T_ ext{R})有(sigma_p^2(t)=2D(t-t_0))，对于(t-t_0gg T_ ext{R})有(sigma_p^2(t)=langle p^2 angle-langle p angle^2=langle p^2 angle=D/gamma)，另一方面，经过长时间后粒子和周围流体达到稳态，满足
[ frac{langle p^2 angle}{2M}=frac{1}{2}k_ ext{B}T ]
从而有(D=Mgamma k_ ext{B}T)，称为爱因斯坦关系。

继续考虑外势不为零的Langevin方程，
[ dot{x}(t)=frac{p(t)}{M}\dot{p}(t)=-gamma p(t)-frac{partial U(x)}{partial x}+F(t) ]
现在在大于( au_c)尺度上，向量({x(t),p(t)})是Markov的。现在利用上面方程来求前面提到的(mathscr{M}_i)和(mathscr{D}_{ij})。需要特别注意，例如，(mathscr{M}_2)的定义为
[ mathscr{M}_2=lim_{delta t o 0}frac{overline{delta p}}{delta t}=lim_{delta t o 0}frac{1}{delta t}int ext{d}(delta p)int ext{d}(delta p) delta xP(x+delta x,p+delta p,t+delta t|x,p,t) ]
它并不等于
[ lim_{delta t o0}frac{langle p(t+delta t) angle-langle p(t) angle}{delta t}=frac{partial langle p(t) angle}{partial t} ]
根据(mathscr{M}_2)的定义，从(t)时刻一个确定的({x(t),p(t)})值开始，对(dot{p}(t)=-gamma p(t)-frac{partial U(x)}{partial x}+F(t))两边在(delta t)短时间内积分，得到
[ delta p=-gamma p(t)delta t-frac{partial U(x)}{partial x}delta t+int_t^{t+delta t}F(t') ext{d}t' ]
两边对系综求平均，利用(langle F(t) angle=0)得到
[ frac{overline{delta p}}{delta t}=-gamma p(t)-frac{partial U(x)}{partial x} ]
同样计算得到
[ frac{overline{delta x}}{delta t}=frac{p(t)}{M} ]
对于(frac{overline{(delta p)^2}}{delta t})，首先把上面的(delta p)平方，注意(delta t)的二次项贡献为零，因此只剩下最后一项有贡献，因此
[ frac{overline{(delta p)^2}}{delta t}=lim_{delta t o 0}frac{1}{delta t}int_t^{t+delta t} ext{d}t'int_t^{t+delta t} ext{d}t'' langle F(t')F(t'') angle ]
计算这个二重积分要小心。首先注意到积分区域是正方形，于是可以分为两个直角三角形区域来积分，有
[ int_t^{t+delta t} ext{d}t'int_t^{t+delta t} ext{d}t'' langle F(t')F(t'') angle=2int_t^{t+delta t} ext{d}t'int_{t'}^{infty} ext{d}t'' langle F(t')F(t'') angle ]
因为考虑到(langle F(t')F(t'') angle=2Ddelta (t'-t''))，所以第二个积分上限改为正无穷。因为$delta $函数是偶函数，于是
[ int_{t'}^infty ext{d}t''delta(t'-t'')=frac{1}{2} ]
所以上面二重积分等于(2Ddelta t)，于是
[ frac{overline{(delta p)^2}}{delta t}=2D ]
因为(delta t)的二次项没有贡献，因此
[ frac{overline{(delta x)^2}}{delta t}=frac{overline{(delta x)(delta p)}}{delta t}=0 ]
至此，(mathscr{M}_i)和(mathscr{D}_{ij})都计算完毕，得到Fokker-Planck方程为
[ frac{partial P(x,p,t)}{partial t}=-frac{p}{M}frac{partial}{partial x}P(x,p,t)+frac{partial}{partial p}left[left(gamma p+frac{partial U}{partial x} ight)P(x,p,t) ight]+Dfrac{partial^2}{partial p^2}P(x,p,t) ]