一种用于动力学与物理学中的非线性时间分数阶对流-反应-扩散-波动方程的精确局部化无网格技术及快速算法

《Communications in Nonlinear Science and Numerical Simulation》：An accurate localized meshfree technique with fast algorithm for nonlinear time-fractional convection-reaction-diffusion-wave equations arising in dynamics and physics

　　

摘要 本文提出并分析了一种用于求解二维 (2D) 非线性对流 - 反应 - 扩散 - 波方程的快速局部化无网格方法，该方程由阶数为 α ∈ (1, 2) 的 Caputo 分数阶导数控制，并表现出初始时间奇异性。核心策略涉及通过引入辅助变量，将原始的 α 阶方程简化为两个 (α/2) 阶分数阶方程的等效耦合系统。建立了两组基于一般非均匀时间网格的时间半离散格式。第一组利用直接 L1 和 Alikhanov 公式来近似 Caputo 导数，并结合 Newton 线性化处理非线性项。第二组采用通过指数和 (SOEs) 近似开发的快速 L1 和快速 Alikhanov 公式来处理分数阶导数。空间离散随后由局部径向基函数单位分解 (LRBF-PU) 技术处理，总共产生四种全离散算法。严格的能量分析证明了所有半离散格式在 L2- 和 H1-范数下的稳定性和收敛性。收敛结果表明，最优网格参数选择使得 L1 型算法具有 (2?α/2) 阶时间精度，而 Alikhanov 型算法具有 2 阶时间精度。涉及正则和不规则区域上 Sine-Gordon 型、Klein-Gordon 型和 Fisher 型非线性的数值实验确认了理论发现。此外，快速 L1 和快速 Alikhanov LRBF-PU 无网格算法表现出比直接对应方法显著更高的计算效率，同时保持相当的精度，突显了所提出的快速方法的实际优势。

引言 时间分数阶扩散 - 波方程 (TF-DWEs) 自然捕捉了经典整数阶偏微分方程 (PDEs) 中不存在的非局部记忆效应，使其非常适合建模反常传输过程，如多孔介质中的流体渗流和反常热传导 [1]。与仅描述瞬时状态演变的整数阶模型不同，分数阶公式考虑了依赖历史的动态。为了容纳包括流体流动和污染物迁移在内的多物理传输行为 [2]，引入对流和反应项来构建时间分数阶对流 - 反应 - 扩散 - 波方程 (TF-CRDWEs)。进一步引入非线性源项产生了重要的模型，如分数阶 Sine-Gordon 和分数阶 Klein-Gordon 方程，这些模型表征了与长程记忆耦合的非线性波现象，如晶格振荡和量子场动力学 [3]。然而，非局部 Caputo 分数阶算子与非线性源项之间的强耦合使得精确解析解几乎难以获得。因此，构建准确、高效的数值离散格式构成了分数阶 PDE 数值分析中的核心研究主题，并为物理和工程领域的实际模拟奠定了必要基础。

本工作涉及以下非线性 TF-CRDWE 的数值逼近：
0Dtαu=μΔu+ν·?u+κu+f(x,t)+g(u),x∈Ω,t∈(0,T],
受初始条件约束
u(x,0)=?(x),ut(x,0)=φ(x),x∈Ω∪?Ω,
以及 Dirichlet 边界条件
u(x,t)=b(x,t),x∈?Ω,t∈(0,T].
此处 α ∈ (1, 2), T > 0 表示终止时间，ut=?u/?t，且 x=(x,y)∈Ω?R2，其中 Ω 是边界为 ?Ω 的有界区域。常数 μ > 0 和 κ 分别表示扩散和反应系数，ν=(ν1,ν2) 是恒定对流速度向量。函数 ?, φ 和 b 是给定的连续初始 - 边界数据；f 和 g(u) 分别表示线性和非线性源项。符号 ? 和 Δ 分别表示梯度算子和 Laplace 算子。阶数为 α 的 Caputo 时间分数阶导数 0Dtαu [4] 定义为
0Dtαu(x,t)=∫0tω2?α(t?ζ)?2u(x,ζ)?ζ2dζ,
其中 Abel 核函数为 ωα(ζ)=ζα?1/Γ(α)，且 Γ( · ) 表示 Gamma 函数。

关于初边值问题 (IBVP) (1.1)–(1.3) 的适定性，[5] 为有界 Lipschitz 区域上带有 Caputo 导数的时间分数阶波方程建立了线性模型的弱解框架，并在 [6] 中扩展到半线性分数阶波方程，假设非线性源项具有 C1-正则性。在我们的基本假设下，特别是 g∈C2(R)，非线性项满足所有所需的正则性条件。对于包含扩散、传输和反应项的模型，相关的弱适定性理论见 [7]。在恒定对流速度下，那项工作中采用的散度形式传输算子简化为此处考虑的非散度对流项。因此，本工作重点关注全耦合非线性 TF-CRDWE 系统 (1.1)–(1.3) 的数值离散以及稳定性和误差分析。

根据源项 g(u) 的非线性，TF-CRDWE (1.1) 可分为线性 (g=0) 和非线性 (g ≠ 0) 情况。线性情况作为基本基准，而非线性情况表现出丰富的动力学，并包括几个著名的模型，如 Fisher 型 (g(u)=c1u(1?u))、Sine-Gordon 型 (g(u)=?c12sin(u)) 和 Klein-Gordon 型 (g(u)=c1u+c2u2+c3u3) 方程，其中 ci 为常数系数。指数、双曲和对数非线性也常用于建模各种物理过程。现有的 TF-DWEs 数值方法通常根据对时间分数阶导数的处理方式分为两类。第一类直接使用基于插值的公式对 0Dtαu 进行离散，如 L1 [8]、L2 [9] 和 H2N2 [10]。第二类通过 Riemann–Liouville 积分将问题重构为等效的偏积分 - 微分方程，然后采用后向差分或 Crank–Nicolson 方法结合卷积积分 [11], [12], [13]。在足够的时间正则性假设下，大多数建立的方法实现了典型的收敛阶 O(N?1)、O(N?(3?α)) 或 O(N?2)，其中 N 表示时间子区间的数量。

对于阶数为 α ∈ (0, 1) 的 Caputo 导数，经典的基于插值的格式，包括 L1、L2-1σ 和 L2 公式 [14], [15], [16]，提供了可靠的直接离散化。然而，分数阶导数的非局部性质使得这些格式计算昂贵且内存密集，严重限制了其在长时间和大规模模拟中的应用。为了减轻这种计算负担，快速核近似技术引起了广泛关注，包括基于多极展开的方法 [17], [18] 和通过高斯积分实现的指数和 (SOEs) 近似 [19], [20]。通过将 SOEs 技术与 L1 [21] 和 L2-1σ 公式 [22] 耦合，构建了各种快速时间格式。其他框架内分数阶导数的额外快速近似策略见 [23], [24]。

然而，一个关键挑战依然存在：Caputo 时间分数阶 PDEs 固有的初始时间奇异性。对于 α ∈ (0, 1) 阶扩散方程，Stynes 等人 [25] 严格证明了在 t=0 处解 u 保持连续，而其一阶时间导数 ut 在 t→0+ 时发散。这种降级的时间正则性导致均匀网格上的基于插值的格式无法保持其最优精度。为了解决这个问题，开发了分级和一般非均匀时间网格结合 L1 和 L2-1σ 公式以恢复最优收敛速率 [25], [26], [27], [28]，相应的快速变体也已在最近建立 [29], [30]。

相同的奇异性行为扩展到 α ∈ (1, 2) 阶 TF-DWEs [31]：u 和 ut 在 t=0 处保持连续，但 utt 在 t→0+ 时发散。为了将奇异性适应技术扩展到这一区域，采用了降阶 (OR) 策略。对于线性情况 (g=0, ν=0, κ=0)，[32] 设 v=ut，将方程重构为整数阶方程和 (α?1) 阶时间分数阶扩散方程 (TF-DE) 的耦合系统，并在分级网格上构建了具有 SOEs 加速的基于 L1 的有限差分格式；然而，得到的时间收敛阶 O(N?min{γ α/2,3?α}) 是非最优的。为了纳入非线性效应，[33] 开发了一种对称分数阶降阶 (SFOR) 方法，将控制方程解耦为两个 (α/2) 阶 TF-DEs。结构的对称性使得能够在非均匀网格上进行统一离散化，并大幅简化了基于能量的稳定性和收敛性分析。SFOR 方法此后已应用于几个简化模型：具有分级网格上 L1 和 L2-1σ 格式的线性 TF-DWEs (g=0, ν=0, κ=0) [34], [35]，具有分级网格上 L1 格式的线性时间分数阶反应 - 扩散 - 波方程 (g=0, ν=0, κ ≤ 0) [36]，以及具有一般非均匀网格上快速 L2-1σ 格式的非线性 Sine-Gordon 型方程 (g=?sin(u), ν=0, κ ≤ 0) [37]。

上述所有研究 [31], [32], [33], [34], [35], [36], [37] 共有的一个显著限制是缺少对流项 ν · ?u，这将现有的高阶格式限制为 (1.1) 的简化特殊情况。从这一观察中，我们确定了以下未解决的研究空白。首先，在 g∈C2(R) 的假设下，一般非均匀网格上使用 L1、L2-1σ 公式及其快速对应方法对 TF-CRDWEs 进行严格的时间离散和数值分析仍然缺失。其次，现有的空间离散主要采用有限差分、有限元或紧差分方法处理正则区域，缺乏处理实际应用中常见的不规则区域的几何灵活性。第三，据我们所知，几乎所有的现有分数阶波方程无网格方法都采用均匀时间网格并忽略了初始时间奇异性。迄今为止，尚未为同时包含对流、反应、扩散和非线性源项的全耦合 TF-CRDWEs 构建严格的无网格数值框架。

上述研究空白促使本工作采用的两个核心方法论选择：SFOR 时间重构和 LRBF-PU 空间离散。SFOR 方法 [33] 将原始 α 阶方程 (1.1) 重构为两个 (α/2) 阶 TF-DEs 的等效耦合系统。我们采用这种重构有两个原因。在分析上，其对称结构促进了基于一般非均匀网格上 L1 和 L2-1σ 离散的基于能量的稳定性和收敛性统一分析。在计算上，将高阶分数阶导数解耦为两个低阶分量使得所得半离散系统能够在顺序 (级联) 时间递推框架内求解。这种方法避免了构建大型单体代数系统，并减少了每个时间步的计算成本。

为了克服基于网格方法固有的几何不灵活性，我们诉诸径向基函数 (RBF) 近似。在各种 RBF 技术中，全局 RBF (GRBF) 方法 [38], [39] 可以实现谱收敛，但它们生成稠密、病态的线性系统，因此限于大规模模拟 [40]。局部 RBF (LRBF) 方法通过局部近似减轻了这种困难，导致稀疏且条件良好的离散系统。在本工作中，我们采用 LRBF-PU 框架。该方法将计算区域分解为重叠的补丁，在每个补丁上构建局部 RBF 近似，并使用紧支撑单位分解 (PU) 函数组装全局数值解 [41], [42], [43], [44]。LRBF-PU 技术已应用于对流 - 扩散模型 [45]、椭圆边值问题 [46] 和时间依赖 PDEs [47], [48], [49], [50], [51], [52]。然而，其在涉及初始时间奇异性的时间分数阶波方程中的应用仍未被探索。

本工作的主要目标是为定义在正则和不规则区域上的非线性 TF-CRDWEs (1.1)–(1.3) 开发高效、高精度且理论严格的 LRBF-PU 无网格算法。所提出的算法旨在同时解决多个困难：初始时间奇异性、非局部记忆效应、非线性和耦合的多物理传输项。为了实现这一目标，我们的方法论策略如下。SFOR 变换 [33] 将原始高阶分数阶方程重构为包含两个 α/2 阶子方程的等效耦合系统，这构成了我们离散策略的基础。为了解决奇异行为并实现最优近似精度，分数阶导数由一般非均匀时间网格上的直接 L1 和 L2-1σ 公式近似。进一步引入其 SOEs 加速变体以降低计算成本和内存消耗，同时保持数值精度。非线性源项通过 Newton 线性化处理，以避免每个时间层级内迭代计算引起的额外计算开销。最后，LRBF-PU 无网格空间近似使得所提出的算法在正则和不规则区域上的模拟中具有灵活性和鲁棒性。此后，L2-1σ 和快速 L2-1σ 公式分别称为 Alikhanov 和快速 Alikhanov 公式。主要贡献总结如下。

•算法贡献 — 直接和快速无网格格式。我们为 2D 非线性 TF-CRDWEs (1.1)–(1.3) 构建了两族直接和快速 LRBF-PU 无网格算法。直接格式采用 L1 和 Alikhanov 公式结合 Newton 线性化，而快速格式采用其 SOEs 加速对应方法。所有格式共享统一的 LRBF-PU 空间离散框架。

•算法贡献 — 级联时间递推框架。每个全离散格式生成三个相互依赖的源自 SFOR 解耦变换的线性子系统。不同于 [52] 中的单体组装策略，我们开发了一种新颖的顺序 (级联) 求解框架。在每个时间层级，三个子系统被依次求解，其中前者子系统的数值结果作为后者子系统的必要输入。这种策略避免了组装大规模全局矩阵，减少了每个时间步的矩阵维度，并显著提高了长时间大规模模拟的计算效率。

•理论贡献 — 扩展的误差分析。我们将时间截断误差估计扩展到正则性区域 σ ∈ (2, 3)，而现有分析 [26], [28], [29], [30], [33], [37] 仅覆盖 σ ∈ (0, 1) ∪ (1, 2)。利用 L2 和 H1 能量论证结合数学归纳法，我们严格建立了所有四种半离散格式的稳定性、锐先验误差估计和数值解的一致有界性。这些结果为所提出框架的可靠性提供了理论基础。

•数值验证 — 精度和效率验证。

涉及各种非线性源项（正弦 - 戈登型、克莱因 - 戈登型和费雪型非线性）在规则和不规则区域上的数值实验表明，所开发的 L1 型和 Alikhanov 型无网格算法分别实现了最优的时间收敛阶 (2?α/2) 和 2。此外，快速 LRBF-PU 无网格算法保持了相当的数值精度，并在长时间大规模模拟中比其直接对应算法表现出更优越的计算效率。与四阶紧格式有限差分法的基准比较进一步证实了所提出无网格框架的准确性和效率优势。

这项工作的几个构建模块——SFOR 变换、L1/Alikhanov 离散化、SOEs 加速和 LRBF-PU 近似——各自都是成熟的。本研究的创新性不在于组件本身，而在于它们的协同集成以及由此产生的理论和算法进展。理论上，贡献 3 将误差分析扩展到更广泛的正则性范围，为 L1 型和 Alikhanov 型格式在非均匀网格上的行为提供了新的见解。算法上，贡献 2 引入了级联求解框架，这与整体式方法有根本区别，并为耦合分数阶系统实现了可观的计算节省。集成上，核心创新是针对完全耦合非线性 TF-CRDWEs 的第一个 LRBF-PU 无网格方法——一个同时解决初始时间奇异性、不规则几何、非局部记忆、非线性和耦合传输项的框架。据我们所知，文献中尚不存在能实现这种组合的方法。贡献 1 和 4 通过系统的构建和全面的数值实验验证了这一框架。总的来说，这些理论见解和算法改进构成了一个完整且新颖的针对不规则区域上非线性 TF-CRDWEs 的数值框架。

本文的其余部分组织如下。第 2 节介绍初步材料，包括函数空间、SFOR 变换、时间网格设置、正则性假设和时间截断误差估计。第 3 节构建 L1 和 Alikhanov 半离散时间格式，并给出严格的稳定性和收敛性证明。第 4 节开发使用 SOEs 近似的相应快速变体，并建立其理论性质。第 5 节详述完全离散的 LRBF-PU 算法及其矩阵形式。第 6 节展示验证理论预测并证明计算效率的数值实验。结论和未来研究方向在第 7 节给出。

章节片段
函数空间与不等式
对于整数 m > 0，设 Ω?Rm 为具有边界?Ω 的有界开区域。对于任何 x=(x1,x2,?,xm)∈Ω，Ck(Ω) 表示 Ω 上 k 次连续可微函数的空间，其中 k∈N∪{0}。对于 1 ≤ q ≤ ∞，勒贝格空间 Lq(Ω) 构成一个巴拿赫空间，其范数为
∥w∥L∞(Ω)=esssupx∈Ω|w(x)|,∥w∥Lq(Ω)q=∫Ω|w(x)|qdx,1≤q<∞。
特别地，L2(Ω) 是一个希尔伯特空间，其内积为〈w,v〉=∫Ωw(x)v(x)dx，由此诱导的范数为∥w∥L2(Ω)=〈w,w〉。对于所有 w,v∈L2(Ω)，该范数满足赫尔德的

非均匀 L1 和 Alikhanov 近似
设 N1,jξ 表示 ξ 在时间节点 tj?1,tj 处的线性牛顿插值，对于 1 ≤ j ≤ n，以及 N2,jξ 表示 ξ 在 tj?1,tj,tj+1 处的二次牛顿插值，对于 1≤j≤n?1。它们的一阶导数分别为
(N1,jξ)′(t)=?τξjτj,(N2,jξ)′(t)=?τξjτj+2(t?tj?1/2)τj(τj+τj+1)(ρj?τξj+1??τξj)。
对于 0 < β < 1，Caputo 分数阶导数 (0Dtβξ)(tn) 由经典 L1 公式近似为
(Dτβξ)n:=∑j=1n∫tj?1tjω1?β(tn??)(N1,jξ)′(?)d?=∑j=1nan?j(n,0)?τξj,1≤n≤N,
其中
an?j(n,0)=1τj∫tj?1tjω1?β(tn??)d?,1≤j≤n。

基于非均匀 L1 和 Alikhanov 的快速时间离散化：公式与分析
引理 4.1 [19], [21], SOEs 近似方法
给定参数?∈(0, 2) 和终止时间 T > 0，对于选定的截断时间τ¨和绝对误差容差ε?1，存在一个整数 Nex > 0，以及两组正实数序列{??l}l=1Nex 和{Sl}l=1Nex，分别为权重和求积节点，使得
|t???∑l=1Nex??le?Slt|≤ε,?t∈[τ¨,T],
其中 Nex 满足 Nex=O(logε?1(loglogε?1+log(T/ε))+logτ¨?1(loglogε?1+logτ¨?1))。

注 4.1
实际上，容差ε通常选择在机器精度水平，使得

LRBF-PU 无网格空间离散化
本节为问题 (1.1)–(1.3) 构建四个完全离散算法。这些格式是通过 LRBF-PU 配点法对半离散格式 (3.13) 和 (4.8) 进行空间离散化获得的。

数值实验
本节展示了所提出的四种算法应用于方程 (1.1) 的三种变体的数值结果：非线性正弦 - 戈登型、克莱因 - 戈登型和费雪型模型。这些数值测试旨在验证所提出算法的准确性和效率，并验证前面的理论分析。考虑了三个数值例子：第一个利用制造出的精确解，而第二个和第三个对应于没有

结论
这项工作开发了一个稳健且高效的 LRBF-PU 无网格框架，带有基于快速 SOEs 的加速，用于求解阶数α∈(1, 2) 的 2D 非线性 Caputo 型 TF-CRDWEs。所提出的方法系统地解决了分数阶 PDE 的关键数值挑战：初始时间奇异性、非局部记忆效应、非线性和耦合多物理传输项。一种分数阶降阶变换将原始高阶方程转换为耦合低阶

利益冲突
所有作者披露没有已知的竞争性财务利益或个人隶属关系可能影响这项工作的结果。

CRediT 作者贡献声明
Yang Cao: 写作 – 初稿，可视化，验证，软件，方法论，调查，形式分析，概念化。Zhijun Tan: 写作 – 审查 & 编辑，监督，方法论，资金获取，概念化。

竞争利益声明
作者声明他们没有已知的竞争性财务利益或个人关系可能看起来影响本论文报告的工作。

Yang Cao|Zhijun Tan
中山大学计算机科学与工程学院，中国广州，510006