Fokker-Planck方程

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\text{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\text{d}y_2P(y_2,t_2|y_1,t_1)P(y_3,t_3|y_2,t_2)
\]
利用上式,可得 (\(0<t<t+\tau\)),
\[
P(y,t+\tau|y_0,0)=\int\text{d}\xi P(y-\xi,t|y_0,0)P(y,t+\tau|y-\xi,t)
\]
对于小量\(\tau\)上式左边泰勒展开得到
\[
\frac{\partial}{\partial t}P(y,t|y_0,0)=\frac{1}{\tau}\left[\int\text{d}\xi P(y-\xi,t|y_0,0)P(y,t+\tau|y-\xi,t)-P(y,t|y_0,0)\right]
\]
注意到上式括号内第二项可以写为
\[
P(y,t|y_0,0)=\int\text{d}\xi P(y+\xi,t+\tau|y,t)P(y,t|y_0,0)
\]
代回去得到
\[
\frac{\partial}{\partial t}P(y,t|y_0,0)=\frac{1}{\tau}\left[\int\text{d}\xi P(y-\xi,t|y_0,0)P(y,t+\tau|y-\xi,t)-\int\text{d}\xi P(y+\xi,t+\tau|y,t)P(y,t|y_0,0)\right]
\]
上式是一个连续性方程:某区域内的概率变化量等于流入概率减流出概率。现在关注括号中第一项,令\(z=y-\xi\),则该项就是\(z\)的函数,写为
\[
\int\text{d}\xi P(z,t|y_0,0)P(z+\xi,t+\tau|z,t)
\]
上式在\(z=y\)处展开为泰勒级数,得到
\[
\int\text{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+\tau|z,t)\right]\Bigg|_{z=y}
\]
注意上面求导最后在\(z=y\)点取值,因此可以直接改写为
\[
\int\text{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+\tau|y,t)\right]
\]
代回方括号中,得到
\[
\frac{\partial}{\partial t}P(y,t|y_0,0)=\frac{1}{\tau}\left\{\int\text{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+\tau|y,t)\right]-\int\text{d}\xi P(y+\xi,t+\tau|y,t)P(y,t|y_0,0)\right\}
\]
其中\(n=0\)项和括号第二项抵消,从而
\[
\frac{\partial}{\partial t}P(y,t|y_0,0)=\frac{1}{\tau}\int\text{d}\xi \sum_{n=1}^\infty \frac{(-1)^n}{n!}\xi^n\frac{\partial^n}{\partial y^n}\left[P(y,t|y_0,0)P(y+\xi,t+\tau|y,t)\right]
\]
取极限\(\tau\to0\),得到
\[
\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\text{d}\xi \lim_{\tau\to0}\frac{\xi^n}{\tau}P(y+\xi,t+\tau|y,t)\right]
\]
如果随机变量\(y\)满足在时间间隔\(\tau\)很小时变化同样很小,则上面的求和保留前几项即可。物理上的一维布朗运动就是一个例子,那里的\(p(t)\)是Markov的,在演化时间\(\tau\)很小时,其变化也很小,所以只保留到二阶即可。

而有势能场的一维布朗运动是扩展版,一个粒子的一维坐标、动量组成的向量\(\{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\text{d}\xi_1\int\text{d}\xi_2 \lim_{\tau\to0}\frac{1}{\tau}\left(\xi_1\frac{\partial}{\partial x_1}+\xi_2\frac{\partial}{\partial x_2}\right)^n\left[P(x_1,x_2,t|x_{10},x_{20},0)P(x_1+\xi_1,x_2+\xi_2,t+\tau|x_1,x_2,t)\right]
\]
应用到布朗运动上,其中\(x_1\)对应于坐标\(x\),\(x_2\)对应于动量\(p\),\(\tau\)对应于\(\delta t\),\(\xi_1=x(t+\tau)-x(t)=x(t+\delta t)-x(t)\)记为\(\delta x\),同样的\(\xi_2=\delta y\),而且对于时间起点\(t=0\)记为\(t_0\),则泰勒展开保留到二阶给出 (略去极限\(\lim_{\delta t\to 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)\right]-\frac{\partial}{\partial p}\left[\frac{\overline{\delta p}}{\delta t}P(x,p,t|x_0,p_0,t_0)\right]\\+\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)\right]+\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)\right]\\+\frac{\partial^2}{\partial x\partial p}\left[\frac{\overline{(\delta x)(\delta p)}}{\delta t}P(x,p,t|x_0,p_0,t_0)\right]
\]

其中

\[
\overline{(\delta x)^2}=\int\text{d}(\delta x)\int\text{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_1\partial 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_i\delta x_j}}{\delta t}
\]
上面就是一维布朗运动中常见的Fokker-Planck方程的原型,其中的各个系数需要根据运动方程来确定。

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

继续考虑外势不为零的Langevin方程,
\[
\dot{x}(t)=\frac{p(t)}{M}\\dot{p}(t)=-\gamma p(t)-\frac{\partial U(x)}{\partial x}+F(t)
\]
现在在大于\(\tau_c\)尺度上,向量\(\{x(t),p(t)\}\)是Markov的。现在利用上面方程来求前面提到的\(\mathscr{M}_i\)和\(\mathscr{D}_{ij}\)。需要特别注意,例如,\(\mathscr{M}_2\)的定义为
\[
\mathscr{M}_2=\lim_{\delta t\to 0}\frac{\overline{\delta p}}{\delta t}=\lim_{\delta t\to 0}\frac{1}{\delta t}\int\text{d}(\delta p)\int\text{d}(\delta p) \delta xP(x+\delta x,p+\delta p,t+\delta t|x,p,t)
\]
它并不等于
\[
\lim_{\delta t\to0}\frac{\langle p(t+\delta t)\rangle-\langle p(t)\rangle}{\delta t}=\frac{\partial \langle p(t)\rangle}{\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')\text{d}t'
\]
两边对系综求平均,利用\(\langle F(t)\rangle=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\to 0}\frac{1}{\delta t}\int_t^{t+\delta t}\text{d}t'\int_t^{t+\delta t}\text{d}t'' \langle F(t')F(t'')\rangle
\]
计算这个二重积分要小心。首先注意到积分区域是正方形,于是可以分为两个直角三角形区域来积分,有
\[
\int_t^{t+\delta t}\text{d}t'\int_t^{t+\delta t}\text{d}t'' \langle F(t')F(t'')\rangle=2\int_t^{t+\delta t}\text{d}t'\int_{t'}^{\infty}\text{d}t'' \langle F(t')F(t'')\rangle
\]
因为考虑到\(\langle F(t')F(t'')\rangle=2D\delta (t'-t'')\),所以第二个积分上限改为正无穷。因为$\delta $函数是偶函数,于是
\[
\int_{t'}^\infty\text{d}t''\delta(t'-t'')=\frac{1}{2}
\]
所以上面二重积分等于\(2D\delta 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}\right)P(x,p,t)\right]+D\frac{\partial^2}{\partial p^2}P(x,p,t)
\]

原文地址:https://www.cnblogs.com/immcrr/p/12004136.html

时间: 2024-08-30 17:06:00

Fokker-Planck方程的相关文章

利用图形窗口分割法将极坐标方程:r=cos(θ/3)+1/9用四种绘图方式画在不同的窗口中

利用图形窗口分割法将极坐标方程:r=cos(θ/3)+1/9用四种绘图方式画在不同的窗口中. 解:MATLAB指令: theta=0:0.1:6*pi;rho=cos(theta/3)+1/9; >> polar(theta,rho) >> >> plot(theta,rho) >> semilogx(theta,rho) >> grid >> hist(rho,15) 结果分别如下图: 图1 图2 图3 图4

BZOJ 1041: [HAOI2008]圆上的整点【数论,解方程】

1041: [HAOI2008]圆上的整点 Time Limit: 10 Sec  Memory Limit: 162 MBSubmit: 4210  Solved: 1908[Submit][Status][Discuss] Description 求一个给定的圆(x^2+y^2=r^2),在圆周上有多少个点的坐标是整数. Input 只有一个正整数n,n<=2000 000 000 Output 整点个数 Sample Input 4 Sample Output 4 HINT 科普视频 So

洛谷 P2312 解方程

题目描述 已知多项式方程: $a_0+a_1x+a_2x^2+..+a_nx^n=0$//用LaTex好看多了 求这个方程在[1, m ] 内的整数解(n 和m 均为正整数) 输入输出格式 输入格式: 输入文件名为equation .in. 输入共n + 2 行. 第一行包含2 个整数n .m ,每两个整数之间用一个空格隔开. 接下来的n+1 行每行包含一个整数,依次为a0,a1,a2..an 输出格式: 输出文件名为equation .out . 第一行输出方程在[1, m ] 内的整数解的个

用python解方程和微积分

用python解方程: from sympy import * x = Symbol('x')  y = Symbol('y') print solve([2* x - y -3,3* x + y -7],[x, y]) 2. 求极限: 代码中的oo就代表无穷. from sympy import * n = Symbol('n') s = ((n+3)/(n+2))**n print limit(s, n, oo) 3. 求定积分: integrate函数用于积分问题. from sympy 

[BZOJ 3751][NOIP2014]解方程(哈希)

Description 已知多项式方程: a0+a1*x+a2*x^2+...+an*x^n=0 求这个方程在[1,m]内的整数解(n和m均为正整数). Solution 一道很久很久以前就应该做的noip的题 一定要放上来是要见证我人品崩坏的一下午 生无可恋…QAQ 题解其实也很简单啦 随便找几个素数取模验证是不是等于0就好了 随便 随便 随便 随便找几个…素数 在WA\TLE\OLE间切换,最后还是抄了别人的几个素数 怀疑人生[望天 (BZOJ上的数据是加强了的,如果是ccf的数据那当然就随

NOIP201410解方程(C++)

NOIP201410解方程 难度级别:A: 运行时间限制:1000ms: 运行空间限制:51200KB: 代码长度限制:2000000B 试题描述 已知多项式方程: a0+a1*x+a2*x^2+a3*x^3+-+an*x^n=0 求这个方程在[1, m]内的整数解(n 和 m 均为正整数). 输入 输入共 n+2 行.第一行包含 2 个整数 n.m,每两个整数之间用一个空格隔开.接下来的 n+1 行每行包含一个整数,依次为a0,a1,a2-an.  输出 第一行输出方程在[1, m]内的整数解

解方程

传送门 题目描述 已知多项式方程: a0+a1x+a2x^2+..+anx^n=0 求这个方程在[1, m ] 内的整数解(n 和m 均为正整数) 输入输出格式 输入格式: 输入文件名为equation .in. 输入共n + 2 行. 第一行包含2 个整数n .m ,每两个整数之间用一个空格隔开. 接下来的n+1 行每行包含一个整数,依次为a0,a1,a2..an 输出格式: 输出文件名为equation .out . 第一行输出方程在[1, m ] 内的整数解的个数. 接下来每行一个整数,按

codevs3732==洛谷 解方程P2312 解方程

P2312 解方程 195通过 1.6K提交 题目提供者该用户不存在 标签数论(数学相关)高精2014NOIp提高组 难度提高+/省选- 提交该题 讨论 题解 记录 题目描述 已知多项式方程: a0+a1x+a2x^2+..+anx^n=0 求这个方程在[1, m ] 内的整数解(n 和m 均为正整数) 输入输出格式 输入格式: 输入文件名为equation .in. 输入共n + 2 行. 第一行包含2 个整数n .m ,每两个整数之间用一个空格隔开. 接下来的n+1 行每行包含一个整数,依次

二元一次不定方程

定义: a,b,c是整数,且ab!=0,那么形如ax+by=c的方程称为二元一次不定方程 定理1: 设a,b是整数且 d = gcd(a,b) 如果d|c,那么返程存在无穷多个整数解,否则不存在整数解. 这里大概可以发现二元一次不定方程和同余方程可以相互转化,如在a>0且b>0的条件下,求二元一次方程ax+by=c整数解等价于求一元线性同余方程ax≡c(mod b)的整数解. 给出poj 2142这个题练练手 The Balance Time Limit: 5000MS   Memory Li

《有限元分析基础教程》(曾攀)笔记二-梁单元方程推导

上图是<有限元分析基础教程>中的图. 这是<材料力学>(孙训方)里面给出的图. 之所以给出这两幅图,是因为在推导公式的时候,第一幅图让我误解了:红箭头标注的微端中的外荷载$\bar{p(x)}$,看起来像是面荷载,实际推导公式的时候,$\bar{p(x)}$是个线荷载.由于两幅图中的一些符号不一致,我下面的推导均采用曾攀老师书中的符号,有的时候也将$\bar{p(x)}$简写成$\bar{p}$. 我们取梁的挠度方程$v(x)$作为整个推导的基本量,同时也是整个推导的核心,所有的推