一种用于近似空间分数阶淬火-燃烧方程的全离散非线性方法

《Journal of Computational and Applied Mathematics》：A fully discrete nonlinear method for approximating spatially fractional-order quenching-combustion equations

　　

摘要
在本研究中，我们提出了一种全离散非线性有限差分方法，用于求解空间分数阶熄灭型方程。利用标准 Grünwald-Letnikov 公式与移位 Grünwald-Letnikov 公式的加权平均来近似 Riemann-Liouville 分数阶导数，同时采用自适应 Crank-Nicolson 积分器进行时间离散。所得全离散格式在时间和空间方向均达到二阶精度。此外，数值解能够保持单调性和非负性等关键物理性质，并在适当网格条件下保持稳定。数值实验表明，临界长度、熄灭时间和熄灭位置等重要熄灭特征被精确捕捉。我们还进行了系统性分析，考察分数阶次对非线性熄灭行为的影响。计算结果进一步验证了所提全离散方法的准确性、有效性和鲁棒性。

引言
反应-扩散偏微分方程在建模化学、生物医学和制造业中的现象方面起着关键作用。其中，熄灭-燃烧非线性模型近年来因其在内燃优化等医学和工程过程中的独特应用而受到广泛关注 [1]-[4]。当某些环境参数超过临界阈值时，熄灭解及其时间导数可能在极短的时间段内表现出强烈的奇异性 [2], [5], [6]。当熄灭发生时，偏微分方程中的非线性源项倾向于迅速发散 [2], [7]-[9]。Sel 等人 [10] 研究了 Neumann 边界条件下的一类多维半线性抛物方程的熄灭解，观察到精确的有限时间熄灭，并确定了解的熄灭集合。Chan [11] 采用上下解方法处理了球形区域中多维半线性抛物方程的熄灭集合，证明解仅在单点处熄灭，而时间导数的发散则扩散到整个物理域。Xu 等人 [12] 采用并行时间有限差分方法求解反应-扩散方程。Derakhshan 等人 [13]-[15] 提出了分数阶方程的高精度差分格式和切比雪夫矩阵方法，并比较了 Weber 与 Riemann-Liouville 分数阶系统的渐近稳定性，以丰富分数阶稳定性理论。Dar 等人 [16]-[19] 将 Grünwald-Letnikov 离散格式与虚拟元方法相结合，为各类时间分数阶偏微分方程构造了相应的全离散数值格式。

人们还注意到，传统的整数阶熄灭型微分方程在描述熄灭现象方面存在严重的理论局限，其基于理想化假设的公式无法准确建模实际应用中复杂的物理情况。相比之下，分数阶微分算子具有独特的非局部记忆性质。空间分数阶导数能够反映粒子扩散行为的全局相关性和历史依赖性，非常适合描述工程燃烧和生物医学输运问题中的反常扩散过程和复杂质量传递现象，因此在建模实际复杂物理过程方面表现出明显优势。Padgett [20] 系统研究了具有 Caputo 时间分数阶导数和空间分数阶拉普拉斯算子的 Kawarada 问题，并定量分析了熄灭发生对区域大小的依赖性。另一方面，Zhou 等人 [21] 研究了具有记忆效应的非线性源函数，并通过有限差分近似成功捕捉了熄灭点处时间导数的熄灭行为。Beauregard 等人 [22] 将这些研究进一步推广到多维情形。基于分数阶模型与熄灭现象之间定量关系发现的启发，Xu [23] 考察了具有退化源项的分数阶扩散方程，观察到在适当左右侧分数阶导数下扩散对称时熄灭的存在性。Cao 和 Li [24], [25] 继续利用 Caputo-Hadamard 时间导数探索分数阶 Kawarada 问题，分析了解的存在性、唯一性和熄灭特征。此外，Zhu 等人率先对一系列一维单侧和双侧 Riemann-Liouville 空间分数阶对流-扩散模型方程进行了系统研究，深入考察了分数阶导数对临界长度、熄灭时间和熄灭位置的影响 [26]-[30]。

大多数现有的熄灭问题有限差分格式局限于半离散公式。由于在时间方向采用数值积分，会引入额外的时间方向误差。研究已证明时间方向的精度无法达到二阶 [26], [27], [29], [30]。相比之下，全离散有限差分格式能够同时离散时间和空间，直接导出迭代代数方程，从而实现空间和时间方向精度的统一设计，有利于构造快速算法，是工程计算的优选方案。文献中鲜有关于全离散格式的报告。受此空缺的启发，本研究为半线性空间分数阶反应-扩散方程的熄灭问题构造了一种具有时间空间二阶精度的全离散有限差分格式。严格建立了数值解的正性、单调性和稳定性，并通过数值实验准确捕捉并展示了数值解熄灭行为的演化。

本文其余部分结构如下。第 2 节介绍空间分数阶模型以及时间和空间域的离散，并提出全离散有限差分格式。第 3 节给出所提离散方法的稳定性和收敛性分析，并在时间和空间步长的适当约束下建立数值解序列的正性和单调性。第 4 节通过两个数值模拟实验讨论临界长度、熄灭时间和熄灭位置。最后，第 5 节总结全文并展望未来的研究方向。

分数阶模型与差分近似
设 λ > 0。我们考虑如下熄灭型空间分数阶反应-扩散问题：
?u(x,t)/?t = μ(x)·?^β u(x,t)/?x^β + ?(u), x ∈ [0, λ], t ∈ (0, T);
u(λ, t) = u(0, t) = 0, t ∈ (0, T);
u(x, 0) = ?(x), x ∈ [0, λ]，
其中 β ∈ (1, 2]，变扩散系数 μ(x) > 0，非线性源函数 ?(u) = 1/(1-u)。分数阶导数取 Riemann-Liouville 意义 [26]。0 < ?(x) ? 1 为初值，T < ∞。

取 h = λ/N，τ = T/K，N, K ∈ ?。引入均匀网格区域 {(x_i, t_k)}_{N,K}^{0,0}。

稳定性
定义 3.1[34] 在数值代数中，若矩阵 A = [a_ij] ∈ ?^{N×N} 满足以下性质，则称为 M-矩阵：(i) 对所有 i ≠ j，a_ij ≤ 0，且对所有 i，a_ii > 0；(ii) 严格对角占优；(iii) 矩阵每行元素之和为正，则弱行和准则成立。

定义 3.2[34] 设 A = [a_ij] ∈ ?^{N×N}。若 A 为 M-矩阵且 A^{-1} 的所有元素为正，则称矩阵 A 为逆正且单调的。

定理 3.1 当 (√17 - 1)/2 ≤ β ≤ 2 时，矩阵 I - τ_{k+1}/(2)·B 为 M-矩阵。
证明 设 L = m_sub

模拟实验
本节设计若干数值实验，以验证在不同源项下熄灭型空间分数阶反应-扩散问题 (2.1) 中熄灭现象所预期的临界长度、熄灭时间和熄灭位置。设 μ(x) = 1.0，h = λ/(N+1)，τ = 5×10^{-5}，N = 500。最后，通过广义 Milne 方法验证数值解及其时间导数的收敛阶。

结论与展望
首先，我们为具有一维空间分数阶的反应-扩散熄灭方程构造了一种时间空间均为二阶精度的全离散有限差分格式。本工作弥补了现有半离散格式仅在空间方向达到二阶精度的不足。其次，我们在适当网格约束下严格证明了所提格式的稳定性与收敛性，并验证了数值解保持非负性。

利益冲突声明
作者声明不存在可能影响本文所报告工作的已知竞争性经济利益或个人关系。

杨兰｜朱林｜魏建英｜盛琴
宁夏大学数学与统计学院，宁夏大学力学与科学计算跨学科重点实验室，宁夏基础科学研究中心数学研究所，宁夏银川，中国，750021