扩散过程的停时问题

简单一点，考虑一维的随机扩散问题，假设扩散系数D=12D=\frac{1}{2}D=21，粒子源的初始位置为x0=1x_0=1x0=1，0≤xt≤10 \le x_t \le 10≤xt≤1。假设TTT是粒子第一次扩散到x=0x=0x=0的时间，求TTT的分布。

理论解

粒子的分布ρ(t,x)\rho(t,x)ρ(t,x)同样可以用Fokker-Planck方程：
∂ρ∂t=D∂2ρ∂x2\frac{\partial \rho}{\partial t} = D \frac{\partial^2 \rho}{\partial x^2}∂t∂ρ=D∂x2∂2ρ
初始条件可以写成ρ(0,x)=δ(x−1)\rho(0,x)=\delta(x-1)ρ(0,x)=δ(x−1)，边界条件可以写成ρ(t,0)=0\rho(t,0)=0ρ(t,0)=0，从而这个PDE的解为
ρ(t,x)=12πt[exp⁡(−(x−1)22t)−exp⁡(−(x+1)22t)]\rho(t,x) = \frac{1}{\sqrt{2 \pi t}} \left[ \exp \left( -\frac{(x-1)^2}{2t}\right) - \exp \left( -\frac{(x+1)^2}{2t}\right) \right]ρ(t,x)=2πt1[exp(−2t(x−1)2)−exp(−2t(x+1)2)]
现在考虑TTT的分布，因为T=min⁡t{t:xt≤0}T = \min_t \{t:x_t \le 0\}T=mint{t:xt≤0}，所以
FT(t)=1−∫0∞ρ(t,x)dxF_T(t) = 1 - \int_0^{\infty} \rho(t,x)dxFT(t)=1−∫0∞ρ(t,x)dx
这里用的思路其实是：如果Ac,BcA^c,B^cAc,Bc独立，则
P[A∩B]=1−P[(A∩B)c]=1−P[Ac⊔Bc]=1−(P[Ac]+P[Bc])P[A \cap B] = 1-P[(A \cap B)^c] = 1-P[A^c \sqcup B^c] = 1 - (P[A^c]+P[B^c])P[A∩B]=1−P[(A∩B)c]=1−P[Ac⊔Bc]=1−(P[Ac]+P[Bc])
所以TTT的密度函数为
fT(t)=∂FT(t)∂t=−∫01∂ρ(t,x)∂tdx=−D∫01∂2ρ(t,x)∂x2dx=D∂ρ(t,0)∂xf_T(t) = \frac{\partial F_T(t)}{\partial t} = -\int_0^{1} \frac{\partial \rho(t,x)}{\partial t} dx = -D\int_0^{1} \frac{\partial^2 \rho(t,x)}{\partial x^2} dx = D\frac{\partial \rho(t,0)}{\partial x} fT(t)=∂t∂FT(t)=−∫01∂t∂ρ(t,x)dx=−D∫01∂x2∂2ρ(t,x)dx=D∂x∂ρ(t,0)
这个量又叫diffusive flux to the wall x=0x=0x=0。带入PDE的解：
fT(t)=exp⁡(−t/2)2πt3f_T(t) = \frac{\exp(-t/2)}{\sqrt{2 \pi t^3}}fT(t)=2πt3exp(−t/2)

模拟解

一个比较naive的想法是，模拟nnn次粒子第一次打到x=0x=0x=0的过程，这样就有nnn个TTT的样本了，根据这nnn个样本可以画出TTT的经验密度函数。贴一张我老师slides上的code示例：

第二行iii是for循环第iii次模拟，这段code模拟了一万次粒子第一次打到x=0x=0x=0的过程。第三行选取的τ=0.0001\tau=0.0001τ=0.0001并设置初始值x0=1x_0 = 1x0=1。第四行用while做判断，x>0x>0x>0的样本是还没打到x=0x=0x=0的，我们不需要，当x≤0x \le 0x≤0时会退出这个while循环，这时的t就是第一次打到x=0x=0x=0的一个样本。3是设置的一个时间上限，超过这个上限还没打到x=0x=0x=0就认为粒子可能移动到x>1x>1x>1那段了。While循环内部做更新：
xi+1=xi+τζx_{i+1} = x_i + \sqrt{\tau} \zetaxi+1=xi+τζ
这个就是RK2导出的递推式。

这个是我老师的输出，实线是解析解，空心点是模拟一万个样本的结果，显然这一万个样本方差有点大；实心点是模拟一百万个样本的结果，倒是与解析解接近了很多，但模拟样本数量太大了，耗时会很久，所以这种naive的算法不推荐。

我们再来审视一下这个naive算法，如果粒子不会打到x=0x=0x=0上，那么while循环要做3000次直到t≥3t\ge3t≥3才会退出，这会浪费很多时间。为了提高计算速度，可以考虑更灵活一点的步长。

这一段代码替换的是naive算法中for循环的内容。第二行定义初值，第三行做while循环的判断，如果x>0x>0x>0，粒子还没有到x=0x=0x=0处，就做第四行开始的循环体。循环体第一句做时间的更新：
ti+1=ti+dt,dt=xi2t_{i+1} = t_i + dt,dt=x_i^2ti+1=ti+dt,dt=xi2
这里的时间步长就不是固定的，这句code也不是对所有问题都能用的，因为dt=xi2dt=x_i^2dt=xi2，说明粒子越靠近x=0x=0x=0的时候，步长会变得越来越小，抽取的样本会变多；粒子越远离x=0x=0x=0的时候，步长会按偏离x=0x=0x=0的程度的二次方扩大，减少退出while循环的计算次数。循环体第二句做位置的更新，
xi+1=xi+dtζx_{i+1} = x_i + \sqrt{dt} \zetaxi+1=xi+dtζ
接下来的一个当ttt不少于一个元素时开始执行的while循环比较有意思，先介绍一下这个while循环背后的直觉。因为dt=xi2dt=x_i^2dt=xi2，再加上x0=1x_0=1x0=1，所以一开始的时候步长是相当大的，因为终点x=0x=0x=0也就离初始值一个单位。第一步就跨过x=0x=0x=0的概率正好等于1−Φ(1)1-\Phi(1)1−Φ(1)，为了让这些点不浪费，我们定义一个更宽松的标准：在循环中记录下粒子的轨迹，如果轨迹上相邻两点的中点满足P(x=0)P(x=0)P(x=0)比较小，我们就不再考虑这条轨迹了。内层的while循环只有一个if else结构，if后面的条件就是在执行这个标准，只是这里把这条标准更细化了一下：and前半句要求时间步长不能太短，and后半句的含义是相邻两个点偏离x=0x=0x=0不太远。如果不满足这条标准，就不再考虑这条轨迹，执行else后面的删除命令。如果满足这条轨迹，就取相邻两点的中点，从中点出发重新抽样。这段代码貌似能解决naive算法的耗时长的问题，但实际上如果遇到的是厚尾分布的话同样会被耽误很久的时间。事实上对于这种没有形成封闭区域的扩散问题计算上都会有这个难点，模拟一个粒子的运动的时候很难说粒子未来会不会飘回来打到边界上，所以某个维度上没有边界的扩散问题的停时都会是厚尾分布。

UA MATH575B 数值分析下VI 统计物理的随机模拟方法2相关推荐

UA MATH575B 数值分析下VI 统计物理的随机模拟方法1
UA MATH575B 数值分析下VI 统计物理的随机模拟方法1 模拟SDE的解用蒙特卡罗方法估计SDE解的自相关函数模拟SDE的解以最简单的随机微分方程为例,考虑随机游走 dxt=ξtdtdx ...
UA MATH575B 数值分析下计算统计物理例题2
UA MATH575B 数值分析下计算统计物理例题2 理论解法 C-K方程法特征值法(近似解) 模拟解法 Rejection Sampling Importance Sampling 一个位于原点 ...
UA MATH575B 数值分析下计算统计物理例题1
UA MATH575B 数值分析下计算统计物理例题1 统计物理方法的解析解 Markov链理论解数值解 Monte Carlo模拟. 一道有趣的统计物理的题目.下面这个简单的迷宫中,一只老鼠一开 ...
UA MATH575B 数值分析下统计物理的随机模拟方法5
UA MATH575B 数值分析下统计物理的随机模拟方法5 Ising Model Gibbs Sampling Glauber Dynamics 这一讲介绍Ising Model,它是MCMC与G ...
UA MATH575B 数值分析下统计物理的随机模拟方法4
UA MATH575B 数值分析下统计物理的随机模拟方法4 这一讲介绍MCMC方法,这个方法最早出现在Metropolis在1953年发在J Chem Phys上的Equation of state ...
UA MATH575B 数值分析下IV 带约束的优化
UA MATH575B 数值分析下IV 带约束的优化问题带等式约束的优化问题带不等式约束的优化问题同时带等式约束与不等式约束的优化问题今天不想敲公式,就不写理论了,反正方法也就是前面的Newt ...
UA MATH575B 数值分析下III 图像恢复
UA MATH575B 数值分析下III 图像恢复图像去噪的优化模型算法实现图像去噪的优化模型假设VVV代表原图,观测到的图像是WWW,WWW是原图与噪声的叠加W=V+ξW=V+\xiW=V+ ...
UA MATH575B 数值分析下I 梯度下降
UA MATH575B 数值分析下I 梯度下降梯度下降最速下降下降算法实现举例梯度下降的应用对于凸优化问题min⁡xf(x)\min_x f(x)minxf(x),最主流的数值计算方法是下 ...
UA MATH575B 数值分析下II 牛顿算法
UA MATH575B 数值分析下II 牛顿算法 Pure Newton算法 Damped Newton算法 Levenberg-Marquardt算法 Quasi-Newton算法割线法 Broy ...

UA MATH575B 数值分析下VI 统计物理的随机模拟方法2

UA MATH575B 数值分析下VI 统计物理的随机模拟方法2

扩散过程的停时问题

理论解

模拟解

UA MATH575B 数值分析下VI 统计物理的随机模拟方法2相关推荐

最新文章

热门文章