扩散过程的停时问题

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

理论解

粒子的分布$\rho(t,x)$ρ(t,x)同样可以用Fokker-Planck方程：
$\frac{\partial \rho}{\partial t} = D \frac{\partial^2 \rho}{\partial x^2}$∂t∂ρ​=D∂x2∂2ρ​
初始条件可以写成$\rho(0,x)=\delta(x-1)$ρ(0,x)=δ(x−1)，边界条件可以写成$\rho(t,0)=0$ρ(t,0)=0，从而这个PDE的解为
$\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πt​1​[exp(−2t(x−1)2​)−exp(−2t(x+1)2​)]
现在考虑$T$T的分布，因为$T = \min_t \{t:x_t \le 0\}$T=mint​{t:xt​≤0}，所以
$F_T(t) = 1 - \int_0^{\infty} \rho(t,x)dx$FT​(t)=1−∫0∞​ρ(t,x)dx
这里用的思路其实是：如果$A^c,B^c$Ac,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])
所以$T$T的密度函数为
$f_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=0$x=0。带入PDE的解：
$f_T(t) = \frac{\exp(-t/2)}{\sqrt{2 \pi t^3}}$fT​(t)=2πt3​exp(−t/2)​

模拟解

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

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

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

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

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