具有回火α-Stable等待时间的分数阶Fokker-Planck方程

The Fractional Fokker-Planck Equation with Tempered α-Stable Waiting Times

摘要: 本文介绍了一种具有回火稳定等待时间(tempered stable waiting times)的朗之万型次扩散(Langevin-type subdiffusion)模型。我们考虑了空间相关外力场的情况。该模型在小时间尺度上表现出次扩散行为，并在大时间尺度上收敛到正常扩散。我们从回火稳定过程(tempered stable processes)的理论推导出回火反常扩散(tempered anomalous diffusion)的一般性质，特别是推导了对应于回火次扩散的分数阶福克–普朗克方程(fractional Fokker-Planck equation)的形式。此外，我们构建了一种算法来模拟所引入过程的样本路径。我们应用该算法来近似求解分数阶福克–普朗克方程，并使用蒙特卡洛方法(Monte Carlo methods)研究回火次扩散的统计特性。

Abstract: In this paper, we introduce a Langevin-type subdiffusion model with tempered stable waiting times. We consider the case of a spatially dependent external force field. The model exhibits subdiffusion behavior at small time scales and converges to normal diffusion at large time scales. We derive the general properties of tempered anomalous diffusion from the theory of tempered stable processes, particularly identifying the form of the fractional Fokker-Planck equation corresponding to tempered subdiffusion. Additionally, we construct an algorithm to simulate sample paths of the introduced process. We apply this algorithm to approximate solutions of the fractional Fokker-Planck equation and investigate the statistical properties of tempered subdiffusion using Monte Carlo methods.

1. 引言

在过去的几年中，反常扩散引起了越来越多的关注，在物理和相关科学的各个领域都有观察到。反常扩散的特征是均方位移(MSD)的幂律形式，它不同于众所周知的布朗运动的扩散特性。根据反常指数$\alpha$ ，将反常扩散分为次扩散($0<\alpha <1$ )和超扩散($\alpha >1$ ) [1]。受外力场($\alpha =1$ )影响的布朗运动可以用福克–普朗克方程来描述。

Lévy飞行不具有有限的均方位移，这导致了关于它们的物理意义的问题，因为有限质量的粒子不应该能够进行瞬时跳跃。然而，近年来，使用Lévy飞行来描述物理模型变得越来越流行。它们被用来模拟各种过程，如多孔玻璃和透镜材料的表面扩散，胶束系统或非均质岩石的输运，反应动力学中的特殊问题，单分子光谱中的特殊问题，以及等待和开关弛豫[2] [3]。

著名的分数阶福克–普朗克方程(FFPE)描述了在外部势能($V\left(x\right)$ )的影响下，次扩散和Lévy飞行之间的竞争。分数阶福克–普朗克方程的一般形式为[4]：

$\frac{\partial \omega \left(x,t\right)}{\partial t}={}_{0}D{}_t^{1-\alpha }\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+K{\nabla }^n\right]\omega \left(x,t\right)$ (1)

这里，算子定义为：

${}_0{D}_t^{1-\alpha }f\left(t\right)=\frac{1}{\text{Γ}\left(\alpha \right)}\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^t{\left(t-s\right)}^{\alpha -1}f\left(s\right)\text{d}s}$ (2)

其中$0<\alpha <1$ ，${}_0{D}_t^{1-\alpha }$ 为Riemann-Liouville型的分数阶导数，${\nabla }_{\mu }$ 是Riesz分数阶导数，具有傅里叶变换$ℱ\left\{{\nabla }^{\mu }f\left(x\right)\right\}=-{\left|k\right|}^{\mu }\tilde{f}\left(k\right)$ ，算子${}_0{D}_t^{1-\alpha }$ 的出现是由两者之间的重尾等待时间引起的，而${\nabla }_{\mu }$ 则与底层连续时间随机游走(CTRW)方案中跳跃的重尾分布有关[5]。方程(1)首先从一个广义主方程中推导出来。常数$$ \mathrm{K}$$为反常扩散系数，而$\eta$ 为广义摩擦常数。当$\mu =2$ 时，我们得到了根据均方位移描述次扩散的FFPE，而当$\alpha =1$ 时，方程(1)简化为马尔可夫Lévy飞行。当$\mu =2$ ，$\alpha =1$ 时对应于标准的福克–普朗克方程，分数阶福克–普朗克方程适用于描述一类有限的反常扩散过程的概率密度函数(PDF)的时间演化，其中次扩散指数α不随时间变化。

本文结构如下：第二章求出了具有回火α-stable等待时间的分数阶福克–普朗克方程及其稳态解，第三章构建了一种算法来模拟所引入过程的样本路径，并利用蒙特卡洛方法研究其统计特性，第四章描述了其相关结论及应用。

2. 随机表示

分数阶福克–普朗克方程是一个非常方便的工具，用于研究外力场存在下的反常扩散。根据从属过程给出了方程(1)的朗之万方程[6]，

$Z\left(t\right)=X\left(S\left(t\right)\right)$ (3)

这里，过程$X\left(\tau \right)$ 由Itô随机微分方程(SDE)给出[7]，

$\text{d}X\left(\tau \right)=-V^{\prime }\left(X\left(\tau \right)\right){\eta }^{-1}\text{d}\tau +K^{1/u}\text{d}{L}_u\left(\tau \right)$ (4)

由傅里叶变换$\langle {\text{e}}^{ik{L}_{\mu }\left(\tau \right)}\rangle ={\text{e}}^{-\tau {\left|k\right|}^{\mu }}$ 的对称α-stable Lévy运动$L_u\left(\tau \right)$ 驱动，$S\left(t\right)$ 是定义为$S\left(t\right)=\inf \left\{\tau :U\left(\tau \right)>t\right\}$ 的逆α-stable从属子，其中$U\left(\tau \right)$ 是严格递增的α-stable Lévy运动。扩散过程$X\left(\tau \right)$ 支配着运动的空间特性，而逆从属子$S\left(t\right)$ 引入了重尾阱(粒子保持静止的周期)的机制。捕获事件显著减缓了整体运动，导致测试粒子的时间均方位移呈亚线性。因此，为了实现从次扩散到正常扩散的过渡(即最终消除捕获周期对粒子运动的影响)，需要适当修改等待时间是必要的。因此，我们修改了方程(1)中等待时间的重尾分布，将其替换为回火α-stable定律。

回火α-stable随机变量$T_{\alpha ,\lambda }$ 通过其拉普拉斯变换$E\left({\text{e}}^{-u{T}_{\alpha ,\lambda }\left(\tau \right)}\right)={\text{e}}^{-\tau \left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]}$ 方便地定义，其中常数$\lambda >0$ 为回火参数，$0<\alpha <1$ 为稳定性参数。注意，当$\lambda \to 0$ 时，我们恢复了单侧稳定分布的拉普拉斯变换。值得一提的是，$T_{\alpha ,\lambda }$ 的PDF形式为$c{\text{e}}^{-\lambda x}{f}_{\alpha }\left(x\right)$ ，其中$f_{\alpha }\left(x\right)$ 为单侧稳定分布的概率密度函数，$c>0$ 为归一化常数。因此，$T_{\alpha ,\lambda }$ 的所有阶矩都是有限的，这使得回火α-stable分布对物理应用特别有吸引力[8]。

给定一个无限可分的回火α-stable随机变量$T_{\alpha ,\lambda }$ ，我们可以通过考虑时间来扩展我们的定义，并通过其拉普拉斯变换引入相应的回火α-stable Lévy过程$T_{\alpha ,\lambda }\left(\tau \right)$ ，它满足拉普拉斯变换$E\left({\text{e}}^{-u{T}_{\alpha ,\lambda }\left(\tau \right)}\right)={\text{e}}^{-\tau \left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]}$ [9]。因此，我们定义逆回火α-stable从属变量$S_{\alpha ,\lambda }\left(t\right)=\inf \left\{\tau >0:T_{\alpha ,\lambda }\left(\tau \right)>t\right\}$ ，$t\ge 0$ ，该过程$S_{\alpha ,\lambda }\left(t\right)$ 是系统的一个新的运行时间，将时间尺度从内部时间$\tau$ 改变为物理时间t。最后，将回火次扩散模型定义为隶属朗之万过程[10]，

$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ (5)

其中，过程$X\left(\tau \right)$ 是Itô随机微分方程(SDE)的解。与纯次扩散情况类似，$X\left(\tau \right)$ 负责粒子的跳跃，而$S_{\alpha ,\lambda }\left(t\right)$ 控制捕获事件，这些事件现在遵循回火α-stable定律。事实证明，通过将$S_{\alpha ,\lambda }\left(t\right)$ 作为系统的新时间，我们能够恢复从反常扩散到正常扩散的所需转变。的确，我们可以证明$S_{\alpha ,\lambda }\left(t\right)\propto t$ 为$t\to \infty$ [11]。因此，对于大时间尺度，过程Y表现为标准扩散X。此外，对于无力情况，回火次扩散的二阶矩满足$\langle Y^2\left(t\right)\rangle \propto t^{\alpha }$ 对于足够小的t，这是典型的次扩散动力学。

接下来我们证明回火次扩散$Y\left(t\right)$ 的PDF满足以下广义分数阶福克–普朗克方程：

$\frac{\partial \omega \left(x,t\right)}{\partial t}={\Phi }_t\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+k{\nabla }^n\right]\omega \left(x,t\right)$ (6)

其中$\omega \left(x,0\right)=\delta \left(x\right)$ ，这里${\Phi }_t$ 定义为积分微分算子，定义为：

${\Phi }_{t}f\left(t\right)=\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^{t}M}\left(t-y\right)f\left(y\right)\text{d}yv$ (7)

其中内存核$M\left(t\right)$ 通过它的拉普拉斯变换定义，

$\widehat{M}\left(u\right)={\displaystyle {\int }_0^{\infty }{\text{e}}^{-ut}}M\left(t\right)\text{d}t=\frac{1}{{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}$(8)

我们可以发现，当$\lambda \to 0$ 时，微积分算子${\Phi }_t$ 与分数Riemann-Liouville导数成正比，我们恢复了分数阶福克–普朗克方程(1)。此外，当$M\left(t\right)=1$ 时，式(6)简化为经典的福克–普朗克方程，与微积分算子${\Phi }_t$ 回火分数阶导数密切相关。

由于过程$X\left(\tau \right)$ 由Itô随机微分方程(SDE) (4)给出，因此它的PDF $f\left(x,\tau \right)$ 满足普通的福克–普朗克方程：

$\frac{\partial f\left(x,\tau \right)}{\partial \tau }=\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+k{\nabla }^n\right]f\left(x,\tau \right)$ (9)

对(9)两边同时在Laplace变换可以得到：

$k\widehat{f}\left(x,k\right)-f\left(x,0\right)=\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }\widehat{f}\left(x,k\right)+K{\nabla }^{\mu }\widehat{f}\left(x,k\right)$ (10)

同样，拉普拉斯空间中的(6)有：

$u\stackrel{⌢}{\omega }\left(x,u\right)-\omega \left(x,0\right)=\frac{u}{{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+K{\nabla }^n\right]\stackrel{⌢}{\omega }\left(x,u\right)$ (11)

接下来，让我们分别用$h\left(t,\tau \right)$ 和$g\left(\tau ,t\right)$ 表示$T_{\alpha ,\lambda }\left(\tau \right)$ 和$S_{\alpha ,\lambda }\left(t\right)$ 的概率密度函数，利用$\mathbb{P}\left(S_{\alpha ,\lambda }\left(t\right)\le \tau \right)=\mathbb{P}\left(T_{\alpha ,\lambda }\left(\tau \right)\ge t\right)$ 的性质，我们可以得到：

$g\left(\tau ,t\right)=-\frac{\partial }{\partial \tau }{\displaystyle {\int }_{-\infty }^{t}h}\left(t^{\prime },\tau \right)\text{d}{t}^{\prime }$ (12)

其中，(12)中$g\left(\tau ,t\right)$ 对t的Laplace变换为：

$\widehat{g}\left(\tau ,u\right)=\frac{{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}{u}{\text{e}}^{-\tau \left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]}$(13)

利用总概率公式和$X\left(\tau \right)$ 与$S_{\alpha ,\lambda }\left(t\right)$ 的独立性，我们得到$X\left(S_{\alpha ,\lambda }\left(t\right)\right)$ 的概率密度函数$p\left(x,t\right)$ 由

$p\left(x,t\right)={\displaystyle {\int }_0^{\infty }f}\left(x,\tau \right)g\left(\tau ,t\right)\text{d}\tau$ (14)

将(14)作Laplace变换可以得到：

$\widehat{p}\left(x,u\right)={\displaystyle {\int }_0^{\infty }{\text{e}}^{-ut}}p\left(x,t\right)\text{d}t={\displaystyle {\int }_0^{\infty }f}\left(x,\tau \right)\widehat{g}\left(\tau ,u\right)\text{d}\tau$(15)

将(13)代入到(15)可以得到：

$\begin{array}{c}\widehat{p}\left(x,u\right)={\displaystyle {\int }_0^{\infty }{\text{e}}^{-ut}p\left(x,t\right)\text{d}t}={\displaystyle {\int }_0^{\infty }f}\left(x,\tau \right)\widehat{g}\left(\tau ,u\right)\text{d}\tau \\ ={\displaystyle {\int }_0^{\infty }f}\left(x,\tau \right)\frac{{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}{u}{\text{e}}^{-\tau \left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]}\text{d}\tau \\ =\frac{{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}{u}\widehat{f}\left(x,{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right).\end{array}$(16)

现在通过(10)中的变量$u\to \left[\begin{array}{c}{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\end{array}\right]$ 的变化，我们可以得到：

$\left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]\widehat{f}\left(x,{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right)-f\left(x,0\right)=\left(\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+K{\nabla }^n\right)\widehat{f}\left(x,{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right)$ (17)

最后，由(15)和$f\left(x,0\right)=p\left(x,0\right)$ ，我们可以得出$\widehat{p}\left(x,u\right)$ 在Laplace空间中满足：

$u\widehat{p}\left(x,u\right)-p\left(x,0\right)=\frac{u}{{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+K{\nabla }^n\right]\widehat{p}\left(x,u\right)$ (18)

对(18)作逆Laplace变换可以得到：

$\frac{\partial p\left(x,t\right)}{\partial t}={\Phi }_t\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+k{\nabla }^n\right]p\left(x,t\right)$ (19)

将上述与(6)比较，我们有$\omega \left(x,t\right)=p\left(x,t\right)$ ，因此，我们已经证明了广义分数阶福克–普朗克方程(6)的解$\omega \left(x,t\right)$ 描述了从属过程$X\left(S_{\alpha ,\lambda }\left(t\right)\right)$ 的PDF动力学。

现在我们可以证明(6)的稳态解满足玻尔兹曼–吉布斯分布[12]。我们把(6)的右边写成：

$\frac{\partial \omega \left(x,t\right)}{\partial t}={\Phi }_t\frac{\partial S\left(x,t\right)}{\partial x}$ (20)

其中，

$S\left(x,t\right)=\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+K{\nabla }^n\right]\omega \left(x,t\right)$ (21)

为概率电流。用$V\left(x\right)$ 表示外部电势。因此，$F\left(x\right)=-V^{\prime }\left(x\right)$ 。如果$\omega \left(x,t\right)$ 达到系统达到稳态，其性质变得与时间无关，则$S\left(x,t\right)$ 必须是常数。因此，如果$S_{st}\left(x\right)$ 在任意点满足$S_{st}\left(x\right)=0$ ，那么它在任何地方都必须等于零。因此，式(6)的平稳解满足：

$\left[\frac{\partial }{\partial x}\frac{V^{\prime }\left(x\right)}{\eta }+K{\nabla }^n\right]{\omega }_{st}\left(x\right)=0$ (22)

从中我们可以很容易地得到玻尔兹曼–吉布斯分布：

${\omega }_{st}\left(x\right)=K{\nabla }^{n}c\exp \left(\frac{-V\left(x\right)}{\eta K}\right)$ (23)

其中c为归一化常数。

3. 样本路径的数值逼近

现在演示如何模拟回火扩散过程的样本路径$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 。利用了$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 是一个从属过程的事实，并利用其固有属性来定义回火扩散过程。现在，假设我们想要在区间$\left[0,t\right]$ 上模拟$Y\left(t\right)$ ，其中t是时间范围。所提出的算法由两个步骤组成：

在第一步中，我们近似逆从属子$S_{\alpha },\lambda \left(t\right)$ 的轨迹，我们提出以下近似方案来模拟复合过程$S_{\alpha },\lambda \left(t\right)$ ，它结合了计算效率和控制精度：

$S_{\alpha ,\lambda ,\Delta t}\left(t\right)=\left[\min \left\{n\in N:T_{\alpha ,\lambda }\left(n\Delta t\right)>t\right\}-1\right]\Delta t$ (24)

其中，$\Delta t$ 为步长，$\Delta t\in \left[0,t\right]$ ，可以证明近似过程$S_{\alpha ,\lambda ,\Delta t}\left(t\right)$ 满足：

$\underset{t\in \left[0,T\right]}{\sup }\left|S_{\alpha ,\lambda ,\Delta \text{}\text{}t}\left(t\right)-S_{\alpha ,\lambda }\left(t\right)\right|\le \Delta t$ (25)

因此，我们选择的$\Delta t$ 越小，我们得到的近似就越准确。鉴于(24)要数值模拟过程$S_{\alpha ,\lambda ,\Delta \text{}\text{}t}\left(t\right)$ 只需要生成值$T_{\alpha ,\lambda }\left(n\Delta t\right)$ ，其中$n=1,2,\cdots$ 这可以通过将Lévy过程的独立和平稳增量求和的标准方法$T_{\alpha ,\lambda }\left(t\right)$ 来完成，

$T_{\alpha ,\lambda }\left(0\right)=0$ (26)

$T_{\alpha ,\lambda }\left(n\Delta t\right)=T_{\alpha ,\lambda }\left(\left[n-1\right]\Delta t\right)+Z_i$ (27)

其中$Z_i$ 为独立的回火α-stable随机变量，每个随机变量具有相同的Laplace变换$E\left({\text{e}}^{-u{Z}_i}\right)={\text{e}}^{-\Delta t\left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]}$ [13]。下面，我们详细介绍计算$Z_i$ 的具体算法。假设我们想要生成具有拉普拉斯变换$E\left({\text{e}}^{-uZ}\right)={\text{e}}^{-\Delta t\left[{\left(u+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }\right]}$ 的回火随机变量$Z>0$ 。其方法如下[14]：

1) 生成均值为${\lambda }^{-1}$ 的指数随机变量E。

2) 使用以下公式生成完全偏斜的α-stable随机变量S：

$S=\text{Δ}{t}^{1/\alpha }\frac{\sin \left(\alpha \left(U+\frac{\pi }{2}\right)\right)}{\cos {\left(U\right)}^{1/\alpha }}{\left[\frac{\cos \left(U-\alpha \left(v+\frac{\pi }{2}\right)\right)}{w}\right]}^{\left(1-\alpha \right)/\alpha }$ (28)

这里，U在$\left[\begin{array}{c}-\pi /2,\pi /2\end{array}\right]$ 上均匀分布，w服从均值为1的指数分布。

3) 如果$E>S$ 则$Z=S$ ，否侧转到第1步。

在第二步中，我们的目标是逼近由SDE(4)给出的扩散过程$X\left(\tau \right)$ 。请注意，上一步中考虑的近似过程$S_{\alpha ,\lambda ,\text{}\text{}\Delta \text{}\text{}t}\left(t\right)$ 只取形式$k\Delta t$ (其中$k=0,1,\cdots ,N$ )，其中N是满足$N\text{Δ}t=S_{\alpha ,\lambda ,\text{Δ}t}\left(T\right)$ 的适当的最后一个指标。因此，我们必须仅在时间点${\tau }_k=k\Delta t$ ，其中$k=0,1,\cdots ,N$ 。这是由经典的欧拉格式[15]完成的。

$X\left(0\right)=0$ (29)

$X\left({\tau }_k\right)=X\left({\tau }_{k-1}\right)-\frac{V^{\prime }\left(X\left({\tau }_{k-1}\right)\right)}{\eta }\text{Δ}\tau +K^{1/u}{\left(\text{Δ}\tau \right)}^{1/u}{\xi }_k$ (30)

其中$k=0,1,\cdots ,N$ 。这里的${\xi }_k$ 是具有标准对称$\mu$ -stable分布的随机变量。随机变量${\xi }_k$ 的计算机模拟实现方法如下：

${\xi }_k=\frac{\sin \left(\mu \overline{V}\right)}{{\left[\cos \left(\overline{V}\right)\right]}^{\frac{1}{\mu }}}{\left(\frac{\cos \left(\overline{V}-\mu \overline{V}\right)}{\overline{W}}\right)}^{\frac{1-\mu }{\mu }}$ (31)

这里，随机变量$\overline{V}$ 均匀分布在$\left[\begin{array}{c}-\pi /2,\pi /2\end{array}\right]$ 上，$\overline{W}$ 呈指数分布，均值为1。从算法的第一步开始，假设$\overline{V}$ 和$\overline{W}$ 是相互独立的，并且是独立随机变量。

最后，我们在模拟算法的初始步骤推导出$S_{\alpha },\lambda \left(t\right)$ 的近似轨迹和第二步中推导出$X\left(\tau \right)$ 的近似轨迹的基础上，将前两个步骤进行整合，生成回火扩散过程的样本模拟路径$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 。本文采用的数值算法基于从属过程分解和欧拉离散化方法。算法的收敛性可以从以下两个方面进行分析：首先，对于逆从属子的近似，由式(24)可知，近似误差与步长成正比；其次，对于扩散过程的欧拉离散化，其弱收敛阶为1。综合考虑，整个算法的整体弱收敛阶为1。

在上述算法中，不管外力$F\left(x\right)$ 怎么设置，过程$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 的样本路径依赖于参数$\alpha \in \left(0,1\right)$ ，$\mu \in \left(0,2\right)$ ，并且$\lambda >1$ 。在图1中，我们看到了在无力情况下逆从属子和回火次扩散的典型轨迹，根据图中三个轨迹之间的恒定间距，可以推断该过程的捕获时间遵循由回火α-stable定律支配的分布。在图2中显示了过程$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 的分位数线，提供了关于轨迹时间演化的基本统计信息，通过10,000个粒子的蒙特卡洛模拟，我们得到了分位数线的估计。随机过程$Y\left(t\right)$ 的p-分位数线$p\in \left(0,1\right)$ 是由关系$P\left(Y\left(t\right)\le q_p\left(t\right)\right)=p$ 给出的函数$q_p\left(t\right)$ 。

Figure 1. An exemplary sample path of the inverse tempered α-stable subordinator $S_{\alpha }$ , $\lambda \left(t\right)$ (for the first figure), the primary process $X\left(\tau \right)$ exhibits anomalous diffusion (for the second figure), the tempered subdiffusion process (for the third figure). The parameters are $\alpha =0.95$ , $\lambda =0.01$ and $\mu =1.0$

图1. 逆回火α-stable从属子$S_{\alpha }$ ，$\lambda \left(t\right)$ (第一张图)的示例示例路径，主过程$X\left(\tau \right)$ 表现出反常扩散(第二张图)，回火次扩散过程(第三张图)。参数为$\alpha =0.95$ ，$\lambda =0.01$ ，$\mu =1.0$

Figure 2. Estimated quantile lines with the two sample paths of the anomalous diffusion $Y\left(t\right)$ , with parameters as shown in Figure 1

图2. 反常扩散$Y\left(t\right)$ 的两条样本路径估计的分位数线，参数如图1所示

Figure 3. Evolution in time of the PDF of the anomalous diffusion $Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ . The parameters are $\alpha =0.95$ ，$\lambda =0.01$ ，$\mu =1.0$ and $F\left(x\right)=-V^{\prime }\left(x\right)=-\left(x^3-16x\right)/20$ . These results were obtained via Monte Carlo methods using 10,000 simulated trajectories

图3. 反常扩散$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 的PDF在不同时间的演化，参数为$\alpha =0.95$ ，$\lambda =0.01$ ，$\mu =1.0$ 和$F\left(x\right)=-V^{\prime }\left(x\right)=-\left(x^3-16x\right)/20$ ，这些结果是通过蒙特卡罗方法获得的，使用了10,000个模拟轨迹

Figure 4. Evolution in time of the PDF of the anomalous diffusion $Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ The parameters are $\alpha =0.95$ ，$\lambda =0.01$ ，$\mu =1.5$ and $F\left(x\right)=-V^{\prime }\left(x\right)=-\left(x^3-16x\right)/20$ . These results were obtained via Monte Carlo methods using 10,000 simulated trajectories

图4. 双势阱中回火次扩散$Y\left(t\right)=X\left(S_{\alpha },\lambda \left(t\right)\right)$ 的分析与估计稳态pdf的比较。我们观察到这两条路线几乎完全一致。参数为$\alpha =0.95$ ，$\lambda =0.01$ ，$\mu =1.5$ 和$F\left(x\right)=-V^{\prime }\left(x\right)=-\left(x^3-16x\right)/20$ ，这些结果是通过蒙特卡罗方法获得的，使用了10,000个模拟轨迹

我们采用蒙特卡洛方法数值逼近广义分数阶福克–普朗克方程(6)的解，图3给出了双势阱$V\left(x\right)=\left(x^4/4-8x^2\right)/20$ 在不同时间点t下过程$Y\left(t\right)$ 的概率密度函数。这些概率密度函数也是$F\left(x\right)=-V^{\prime }\left(x\right)=-\left(x^3-16x\right)/20$ 时分数阶福克–普朗克方程(6)的解。在引入了回火参数λ之后，可以修正了纯幂律模型的非物理特性，使模型既能描述短时反常扩散，又能反映长时正常化行为。相比之下，这种截断特性使得模型更符合物理实际，因为真实系统中的粒子不可能具有无限跳跃或无限等待时间。如图4所示，理论稳态概率密度函数与Rosenblatt-Parzen核密度估计值(从10,000条模拟轨迹中获得)非常吻合。

回火参数λ在布朗运动和Lévy飞行驱使的分数阶福克–普朗克方程中均起到关键作用，但其影响机制和效果存在显著差异。对于布朗运动驱使的分数阶福克–普朗克方程，回火参数主要通过修正时间分数阶导数，将短期的次扩散行为过渡到长期的正常扩散，同时保持稳态的玻尔兹曼分布形式。

4. 结论

本文将回火次扩散模型推广到空间依赖力场的情况。模型结合了两个独立的随机机制：一个是经典的扩散过程$X\left(\tau \right)$ ，它控制着测试粒子的空间特征(跳跃)——这里，$X\left(\tau \right)$ 是Lévy飞行；另一种是由逆服从过程$S_{\alpha },\lambda \left(t\right)$ (回火α-stable等待时间)表征的捕获事件机制。

本文导出了一个广义分数阶福克–普朗克方程(FFPE)，该方程描述了从属过程$X\left(S_{\alpha },\lambda \left(t\right)\right)$ 的概率密度函数(PDF)的演化。当$\lambda =0$ 时，恢复经典的次扩散方程。路径模拟算法基于朗之万方法，构建了一个数值格式来模拟这一过程的样本路径。该算法允许通过蒙特卡罗方法研究系统的许多相关属性。它可以用来近似具有任意空间相关力($F\left(x\right)$ )的广义FFPE的解，并验证回火次扩散的各种物理性质。