使用三次样条和无网格近似方法对分布阶Caputo–Prabhakar时间与Riesz空间分数阶模型的数值模拟

《Engineering Analysis with Boundary Elements》：Numerical simulation of distributed-order Caputo–Prabhakar time and Riesz space fractional models using cubic spline and meshless approximations

　　

## 摘要

本文提出了一种新的分布阶次分数阶偏微分方程模型，该模型采用Caputo-Prabhakar时间导数与变阶次Riesz空间导数，能够捕捉复杂的记忆效应和非局部效应。本文开发了一种完全离散化的数值方法，在时间方向采用三次样条近似，在空间方向采用无网格移动最小二乘法进行离散化。收敛性与稳定性分析证实，该方法在空间方向具有二阶精度，在时间方向具有三阶精度。文中给出两个示例：其中一个具有精确解以验证精度，另一个没有精确解以展示方法的鲁棒性。数值结果表明，所提出的方法高效、精确，适用于求解科学与工程应用中复杂的分数阶偏微分方程。

## 引言

分数阶模型[1]、[2]、[3]、[4]、[5]、[6]、[7]、[8]已在现代工业和工程中得到了广泛应用，因为它们能够为具有记忆效应、遗传特性和空间非局部性的系统提供更为精确的描述。在聚合物制造、电池技术和石油开采等工业过程中，经典的整数阶模型无法捕捉长期记忆效应和异常输运现象。因此，分数阶模型被应用于热工程、信号处理、机器人技术和控制系统中，以提高预测精度和稳定性分析能力。在机械工程领域，分数阶模型被用于描述粘弹性材料、阻尼装置和复杂结构动力学；而在电气工程中，分数阶模型则出现在超级电容器和具有记忆效应的传输线建模中。Riesz分数阶导数在物理学和工程学中起着至关重要的作用，因为它代表了对称形式的空间非局部性，且与Lévy飞行和异常扩散密切相关。Riesz分数阶导数已成功应用于等离子体物理、量子力学以及异质介质中的热传导。在工程科学中，Riesz导数被用于建模多孔材料中的扩散、非局部梁的振动以及复杂网络中的输运。其通过左右Riemann-Liouville导数进行表述的数学形式，能够捕捉双向长程相互作用，这对于描述真实的物理系统至关重要。

Caputo-Prabhakar分数阶导数[9]、[10]近年来引起了极大的关注，这得益于其能够以Mittag-Leffler型和Prabhakar型行为描述广义记忆核的能力。在物理学中，它被用于描述复杂介质中的弛豫过程、介质极化和异常热传递。与经典的Caputo导数相比，Caputo-Prabhakar导数提供了更大的灵活性，因为它涉及多个参数来控制记忆核的衰减速率和形状。这种额外的自由度使其更适合拟合实验数据和描述不同扩散机制之间的交叉现象。因此，本文选择Caputo-Prabhakar导数以更好地表征标准Caputo导数无法捕捉的时变记忆效应。

研究空间变阶次和分布阶次变量的重要性在于，物理系统在不同空间尺度上往往表现出不同的扩散机制。在异质介质和分形介质中，有效扩散阶次可能随位置或相互作用范围的变化而变化。分布阶次公式允许将多个分数阶次组合为单一模型，从而提供更为现实和灵活的数学描述。此类模型可以描述正常扩散与异常扩散之间的转换，以及亚扩散和超扩散行为的共存。因此，研究空间中的变分数阶和分布分数阶对于理解复杂输运过程以及开发具有强物理和工程应用相关性的精确数学模型至关重要。

设$v(\mathbf{r},t)$，其中$\mathbf{r}=(x,y)$，为定义在域$(\mathbf{r},t)\in(0,1)\times(0,1)\times(0,1]$上的未知函数。我们考虑以下非线性分布阶次空间分数阶与Caputo-Prabhakar型时间分数阶模型：

$$PD^{\mu,\nu,\vartheta}_{0+}{}^{\omega}v(\mathbf{r},t) = \theta_1\int_{\frac{1}{2}}^{1}\chi_1(\eta)\,D_{\mathbf{r}}^{\eta}v(\mathbf{r},t)\,d\eta + \theta_2\int_{0}^{1}\chi_2(\xi)\,D_{\mathbf{r}}^{\xi}v(\mathbf{r},t)\,d\xi + \mathcal{R}(v(\mathbf{r},t)) + \int_{\Omega}\mathcal{H}(\mathbf{r},\mathbf{s},v(\mathbf{s},t))\,d\mathbf{s},$$

其中$\Omega=[0,1]\times[0,1]$。该模型满足齐次Dirichlet边界条件：

$$v(0,y,t)=v(1,y,t)=v(x,0,t)=v(x,1,t)=0,\quad t\in(0,1],$$

以及初始条件：

$$v(\mathbf{r},0)=\Phi(\mathbf{r}),\quad \mathbf{r}\in[0,1]\times[0,1].$$

式(1.1)中出现的时间分数阶算子是一种广义Caputo-Prabhakar型导数，其定义为[11]、[12]：

$$PD^{\mu,\nu,\vartheta}_{0+}{}^{\omega}\psi(\mathbf{r},t) = \int_{0}^{t}\frac{(t-\tau)^{\nu-1}}{\Gamma(\nu)}\,M_{\mu,1-\nu,\vartheta-\omega}^{\vartheta}\!\left((t-\tau)^{\mu}\right)\frac{\partial\psi(\mathbf{r},\tau)}{\partial\tau}\,d\tau,$$

其可以写为卷积形式：

$$PD^{\mu,\nu,\vartheta}_{0+}{}^{\omega}\psi = (W_{\mu,\nu,\vartheta}^{\omega}\ast\partial_t\psi)(\mathbf{r},t),$$

其中

$$W_{\mu,\nu,\vartheta}^{\omega}(s) = s^{-\nu}\,M_{\mu,1-\nu,\vartheta-\omega}^{\vartheta}\!\left(\vartheta s^{\mu}\right),$$

而

$$M_{\mu,1-\nu,\vartheta-\omega}^{\vartheta}\!\left(\vartheta s^{\mu}\right) = \frac{1}{\Gamma(-\omega)}\,{}_{1}\Theta_{1}\!\left(\begin{matrix}(\mu,1-\nu)\\(1,-\omega)\end{matrix}\right);\vartheta s^{\mu}.$$

Fox-Wright函数定义为[13]、[14]、[15]：

$${}_{1}\Theta_{1}\!\left(\begin{matrix}(\mu,1-\nu)\\(1,-\omega)\end{matrix}\right);s = \sum_{k=0}^{\infty}\frac{s^k}{k!\,\Gamma(k-\omega)\,\Gamma(\mu k+1-\nu)}.$$

式(1.1)中的空间算子$D_{\mathbf{r}}^{\alpha}$表示Riesz分数阶导数[16]、[17]，其通过左右Riemann-Liouville分数阶导数表示为：

$$D_{\mathbf{r}}^{\alpha}v(\mathbf{r},t) = -\frac{1}{2}\cos\!\left(\frac{\pi\alpha}{2}\right)\left[\,{}_0^{RL}D_x^{\alpha}v(\mathbf{r},t) + {}_x^{RL}D_1^{\alpha}v(\mathbf{r},t) + {}_0^{RL}D_y^{\alpha}v(\mathbf{r},t) + {}_y^{RL}D_1^{\alpha}v(\mathbf{r},t)\right],$$

其中左Riemann-Liouville分数阶导数定义为：

$${}_0^{RL}D_x^{\alpha}v(x,y,t) = \frac{1}{\Gamma(n-\alpha)}\frac{\partial^n}{\partial x^n}\int_{0}^{x}(x-\xi)^{n-\alpha-1}v(\xi,y,t)\,d\xi,$$

右Riemann-Liouville分数阶导数定义为：

$${}_x^{RL}D_1^{\alpha}v(x,y,t) = \frac{(-1)^n}{\Gamma(n-\alpha)}\frac{\partial^n}{\partial x^n}\int_{x}^{1}(\xi-x)^{n-\alpha-1}v(\xi,y,t)\,d\xi,$$

其中$n-1<\alpha0$，$\Re(\nu)>0$。

常数$\lambda,\sigma,\rho$描述非线性反应项和非线性Fredholm相互作用的强度，而$q>1$决定非线性的程度。函数$\Phi(\mathbf{r})$表示给定的初始分布。

式(1.1)中引入的模型可以看作是一个非线性分布阶次空间分数阶与Caputo-Prabhakar型时间分数阶积分微分方程。选择该模型的原因在于，它同时通过广义Caputo-Prabhakar导数纳入了时间方向的记忆效应，通过分布阶次Riesz分数阶导数纳入了空间非局部相互作用，为描述异常输运和扩散过程提供了现实的框架。非线性反应项与非线性Fredholm积分算子的存在显著增加了动力学的复杂性，反映了经典模型无法捕捉的长程相互作用和反馈机制。由于其多参数分数阶结构，该模型表现出丰富的行为谱，包括亚扩散、超扩散以及它们之间的过渡机制。这种复杂的结构使得数学分析具有挑战性，并促使对解的存在性、唯一性和稳定性进行严格的研究。从物理角度来看，该模型非常适合表征异质多孔介质、粘弹性材料和复杂流体中的输运过程。在工程科学中，它应用于具有记忆效应的热传导、复杂网络中的信号传播以及具有遗传特性的控制系统。因此，所提出的分数阶模型不仅具有重要的理论意义，而且在物理学和工程学中具有相当大的实际应用价值。

Caputo-Prabhakar型和Riesz型分数阶导数近年来引起了广泛关注，这得益于其描述复杂记忆效应和异常输运现象的能力。目前已开发了多种Caputo-Prabhakar导数的数值近似方法，包括针对线性和非线性分数阶模型的高阶格式和预测-校正方法[11]、[18]、[19]。这些方法已成功应用于分数阶对流-扩散方程和混沌动力系统，证明了其高效性和准确性[11]、[20]。此外，基于小波和配点法的技术也被提出用于变阶次和积分微分分数阶问题，提供了稳定而灵活的计算框架[21]。基于这些进展，本工作采用高阶时间离散化策略对Caputo-Prabhakar分数阶导数进行数值近似。由此构造了半离散时间格式，并通过改编分数阶扩散和演化方程的标准技术来考察其收敛行为[18]、[19]。该分析确保了时间近似在保持数值稳定性的同时，保留了Caputo-Prabhakar算子的本质记忆特征。

为了离散化空间分数阶算子，Riesz空间分数阶导数通过高阶差分公式进行近似。文献中已提出并分析了若干用于Riesz导数的高效算法[22]、[23]、[24]、[25]，这些方法已被广泛应用于分数阶扩散和对流-弥散模型[26]、[27]。在本框架中，这些空间近似与时间Caputo-Prabhakar离散化相结合，从而获得完全离散的数值格式。通过采用离散能量技术和相容性论证，按照空间-时间分数阶模型的标准方法[23]、[24]、[27]，对完全离散格式的收敛性进行了研究。理论结果表明，在适当的正则性假设和网格约束条件下，所提出的数值解收敛于精确解。此外，还对完全离散格式的稳定性进行了分析，表明该方法在长时间积分中保持稳定，且对分数阶次变化具有鲁棒性。

最后，进行了若干数值实验以验证理论发现。考虑了涉及时间分数阶对流-扩散方程和空间分数阶扩散方程的测试问题，数值结果证实了所提出格式的有效性、准确性和稳定性。这些模拟表明，Caputo-Prabhakar时间离散化与Riesz空间分数阶近似的组合为求解复杂分数阶偏微分方程提供了一个可靠的计算工具[18]、[19]、[28]。

## 无网格数值方法

近年来，无网格数值方法因其处理复杂几何形状的灵活性和无需生成网格的特点而备受关注。在分数阶偏微分方程的领域，移动最小二乘法（MLS）及相关无网格技术得到了广泛应用。文献[29]中开发了一种用于时间相关分数阶对流-扩散方程的隐式MLS无网格格式，并在文献[30]中将其推广至二维时间相关分数阶扩散-波动方程。文献[31]中提出了用于分布阶次时间分数阶扩散-波动方程的改进无网格策略，并推导了误差估计。文献[32]中构造了用于双侧空间分数阶波动方程的高阶无网格格式，并在文献[33]中通过改进的移动最小二乘近似进一步优化。文献[34]中对移动克里金插值与无单元Galerkin方法在分数阶偏微分方程中的应用进行了对比研究。

无网格方法也被成功应用于非线性和变阶次分数阶模型。文献[35]中引入了一种用于涉及Mittag-Leffler核的非线性变阶次时间分数阶反应-扩散方程的无网格方法。文献[36]中提出了一种用于求解非线性多阶次时间分数阶微分方程的配点离散最小二乘无网格格式。此外，文献[37]、[38]中利用无网格局部Petrov-Galerkin格式处理了由时间分数阶偏微分方程控制的反演和不适定问题。

结合无网格方法与经典离散化的混合技术也得到了研究。文献[39]中将耦合有限差分与径向基函数无网格方法应用于分数阶随机对流-扩散方程。文献[40]中提出了一种用于时间分数阶非线性扩散-波动方程的局部无网格方法，文献[41]中进行了分数阶对流-扩散-反应方程的数值模拟。文献[42]、[43]中提出了校正离散最小二乘无网格方法和基于Voronoi的MLS格式用于自由表面流动问题。最近，无网格技术被进一步推广至复杂物理模型。文献[44]中开发了一种基于移动克里金的无网格方法用于磁流体力学分数阶Maxwell流体方程。Sun等人的研究开发了一种改进的无单元Galerkin格式，采用维度分割MLS策略，能够准确模拟复杂不规则几何形状中的二维势问题。类似地，Lv等人提出了一种用于时间分数阶Black-Scholes方程的混合插值Galerkin无网格格式，提供了严格的误差估计，增强了金融建模中的数值可靠性。这些进展表明，无网格方法为求解广泛类别的分数阶偏微分方程——包括扩散方程、波动方程、反应-扩散方程和反演问题——提供了一个强大而灵活的框架。

无网格方法已被广泛研究并成功应用于广泛的线性和非线性偏微分方程，包括分数阶和积分微分模型。在这些方法中，移动最小二乘法（MLS）及相关无网格格式因其在无需生成网格的情况下处理离散节点和复杂几何形状的灵活性而备受关注。特别是，基于MLS的无网格局部Petrov-Galerkin（MLPG）方法已被有效地用于非线性和分数阶偏微分方程。Dehghan和Mirzaei将MLPG框架应用于广义二维非线性Schr?dinger方程，展示了其在非线性波动传播问题上的准确性和稳定性[45]。这项工作突出了无网格格式在避免网格构建的同时保持高阶空间精度的效率。

最近，MLS技术被推广至分数阶和积分微分方程。Elbostani和El Jid提出了一种基于MLS的时间分数阶偏微分积分方程方法，表明MLS为涉及记忆效应和弱奇点的问题提供了准确的数值近似[46]。类似地，Solhi等人开发了一种用于Hammerstein型随机分数阶Volterra积分微分方程的增强MLS方法，强调了随机环境下改进的稳定性以及收敛性质[47]。

除偏微分方程应用外，MLS及相关最小二乘格式也被用于数据逼近和计算建模。Gu等人提出了一种用于测量数据拟合的鲁棒移动全最小二乘法，改善了噪声和规则采样条件下数值稳定性[48]。此外，Wang和Liu开发了一种基于径向基函数的点插值无网格方法，为现代无网格离散化技术奠定了重要基础，并进一步影响了基于MLS的方法[49]。这些研究共同表明，基于MLS的无网格方法对于求解分数阶、非局部和积分微分方程非常有效，兼具理论灵活性和实践鲁棒性。

在本研究中，这些优势被扩展到了具有Riesz空间算子的分布式阶Prabhakar型分数阶模型。所提出的分数阶模型的选取原因在于它将分布式阶Caputo–Prabhakar时间导数与分布式阶Riesz空间分数阶算子相结合，从而能够对复杂记忆效应和非局部空间相互作用进行精确描述。与经典整数阶模型相比，该框架通过Prabhakar核捕捉多尺度时间记忆效应，并通过变阶Riesz分数阶算子捕捉非均质异常扩散。这种结构特别适用于表现出长程依赖、非局部输运和多速率弛豫现象的物理与工程系统。对于时间离散化，采用基于三次样条的近似方法对一阶时间导数进行近似，然后将其嵌入Caputo–Prabhakar分数阶算子中。这一选择提供了更高阶的时间精度，并确保了解历史的光滑重构，这对于具有记忆效应的分数阶导数至关重要。与标准有限差分或L1型格式相比，基于样条的方法在保持Prabhakar核非局部结构的同时，具有更高的精度和更好的稳定性。对于空间离散化，采用无网格移动最小二乘（MLS）近似方法，通过Riemann–Liouville算子近似分布式阶Riesz分数阶导数。该方法避免了结构化网格的构建，能够灵活处理复杂几何形状和散射节点。与有限元或有限差分方法不同，基于MLS的近似方法不需要单元连接，并且提供了数值解的更高光滑性，这对于涉及高阶空间微分和弱奇异核的分数阶导数求值至关重要。所提出的近似方法与传统格式之间的主要区别在于时间样条离散化与空间无网格近似的结合使用。经典方法通常依赖于时间上的低阶有限差分公式和基于网格的空间离散化，这可能导致精度降低且难以处理非局部算子。然而，本方法在时间上实现高阶精度，空间上实现二阶精度，同时保持了模型的固有非局部特征。所提出数值方法的效率和优势总结如下：它为时间和空间分数阶算子提供了高精度；它对不规则空间域具有灵活性；它自然地纳入了分布式阶效应；并且表现出良好的稳定性和收敛性。这些特征使该方法非常适合模拟由分布式阶分数阶偏微分方程控制的复杂物理过程。本文其余部分组织如下。第2节，我们给出了Prabhakar型分数阶导数的数值近似。第3节介绍了模型的半离散时间离散化。第4节，开发了一种无网格移动最小二乘（MLS）方法来近似分布式阶Riesz分数阶导数。第5节构建了分布式阶Prabhakar型分数阶偏微分方程的全离散格式。第6节分析了全离散格式的收敛性，随后在第7节进行稳定性分析。第8节提供了数值模拟，其中考虑了两个说明性例子：一个具有精确解以验证精度，另一个没有精确解以展示所提方法的鲁棒性。最后，第9节对研究发现和潜在应用进行了总结。第2节数值近似Prabhakar型分数阶导数在这一节中，我们首先通过三次样条插值近似一阶时间导数，然后将该近似代入分数阶积分表示，从而构建Prabhakar型分数阶导数PDμ,ν,?,0+ωv(r,t)的数值近似。设tn=nΔt，n=0,1,…,N，为步长为Δt>0的均匀时间网格，设v(t)为充分光滑的函数。v(t)在[0,T]上的三次样条插值是一个分段三次多项式S(t)，使得半离散时间离散化模型在这一节中，我们考虑一种半离散方法来近似模型PDμ,ν,?,0+ωu(x,y,t)=?[u](x,y,t)的时间演化（参见公式(1.1)），同时保持空间变量x和y连续。设tn=nΔt，n=0,1,…,N，为步长为Δt>0的均匀时间网格，设un(x,y)表示u(x,y,tn)的近似。利用基于三次样条的一阶时间导数近似，Prabhakar型分数阶导数在t=tn处的半离散化为PDμ,ν,?,0+ω无网格MLS近似分布式阶Riesz分数阶导数近年来，在Riesz空间分数阶算子和分布式阶分数阶模型的高阶数值格式发展方面取得了显著进展。特别是，谱配置方法已成功扩展到非线性Riesz分布式阶分数阶微分方程，为光滑解提供了高精度，但依赖于全局基函数和结构化离散化[50]。与此同时，高阶有限差分全离散格式分布式阶Prabhakar型分数阶偏微分方程从时间上的半离散模型（参见公式(3.29)）出发：∑j=0n?1[ωn,j(0)bj(x,y)+ωn,j(1)cj(x,y)+ωn,j(2)dj(x,y)]=?[un](x,y)，n=1,2,…,N，我们现在用分布式阶Riesz导数的全离散MLS近似（参见公式(4.41)）来替代空间算子?[un](x,y)。利用记号vh(r,tn)≈un(x,y)，我们有?[un](x,y)≈∑l=1Nηwl(η)χ1(ηl)∑i∈N(r)ci,l(ηl)un(ri)+∑m=1Nξw?m(ξ)χ2(ξm)∑i∈N(r)c?i,m(ξm)un(ri)+?h[un]+∑i∈NΩ?h(r,ri,un(ri))，其中•ci,l(ηl)和c?i,m(ξm)是全离散格式的收敛性分析在本节中，我们针对分布式阶Prabhakar型分数阶偏微分方程PDμ,ν,?,0+ωu(x,y,t)=?[u](x,y,t)的全离散数值格式进行详细的收敛性分析，其中(x,y)∈Ω，t∈(0,T]，初始条件为u(x,y,0)=u0(x,y)，(x,y)∈Ω，并配备适当的边界条件。设uhn(ri)表示空间节点ri和时间tn=nΔt处的全离散解，其中i=1,…,Ns，n=0,1,…,N。定义节点ri处的离散误差为ein=u(xi,yi,tn)?uhn(ri)，设en={ein}i=1Ns为全离散格式的稳定性分析在本节中，我们严格分析分布式阶Prabhakar型分数阶偏微分方程PDμ,ν,?,0+ωu(x,y,t)=?[u](x,y,t)的全离散数值格式的稳定性，其中(x,y)∈Ω，t∈(0,T]，初始条件为u(x,y,0)=u0(x,y)，(x,y)∈Ω，并配备适当的边界条件。设uhn(ri)表示空间节点ri和时间tn=nΔt处的全离散解，i=1,…,Ns，n=0,1,…,N。定义扰动?in为?in=?hn(ri)?uhn(ri)，其中?hn为扰动后的数值解。记误差向量?n=数值模拟在本节中，我们进行数值模拟以验证分布式阶Prabhakar型分数阶偏微分方程全离散格式的精度和收敛性。所有计算均使用MATLAB完成。设uexact(x,y,t)表示解析解（或充分精细的参考解），uhn(xi,yj)表示空间节点(xi,yj)和时间tn=nΔt处的全离散近似。我们定义离散绝对误差为E∞n=maxi,j|uexact(xi,yj,tn)?uhn(xi,yj)|。时间上的结论本文提出了一种针对分布式阶Caputo–Prabhakar型分数阶偏微分方程的统一数值和解析框架，其创新之处既体现在全离散格式的结构上，也体现在相应的理论结果中。主要创新在于构建了一种离散时间近似，将Prabhakar核的完整记忆效应嵌入到高阶三次样条表示中。通过构建样条函数致谢本基金没有资金来源声明。利益冲突声明作者声明不存在可能影响本文所报告工作的已知竞争性经济利益或个人关系。致谢所有作者均批准了待发表的稿件版本。MohammadHossein Derakhshan