通过惩罚最小残差迭代法求解Riesz分数阶反应-扩散方程的自适应谱求解器

《Mathematics and Computers in Simulation》：Adaptive spectral solver for Riesz fractional reaction–diffusion equations via penalized minimum residual iteration

高阶求解器被提出用于二维Riesz分数非线性反应扩散方程，采用中点启动和三阶后向差分公式实现二阶时间精度，结合加权雅可比谱空间离散获得近指数收敛。通过Newton线性化与惩罚Levenberg-Marquardt最小残差迭代（PLM-MRM），自适应调整阻尼和边界惩罚参数，无需边界适配基函数即可有效约束边界条件。数值实验表明该方法全局误差减少35%至2个数量级，时间步长可达现有四阶CN方法的四倍。

　　本文介绍了一种用于求解二维Riesz分数阶非线性反应扩散方程的高阶求解器。这种求解器结合了中点启动和三阶后向差分公式（BDF2），以实现时间上的二阶精度，同时采用加权Jacobi谱逼近方法，使解析解在空间方向上具有接近指数级的收敛性。通过牛顿线性化，每一步修正使用一种带有惩罚项的Levenberg-Marquardt最小残差方法（PLM-MRM）来获得。该迭代方法能够自适应地施加边界条件，无需使用边界适配的基础函数。我们建立了稳定性并证明了严格的先验误差界。广泛的数值实验验证了这些收敛速率，并且在几次PLM-MRM迭代中将残差降至机器精度水平。与传统的LM更新方法相比，全局误差减少了高达35%，并且相较于Galerkin-BDF或Crank-Nicolson（CN）基准方法，误差减少了约一个到两个数量级。对于给定的精度，该方案允许的时间步长是最近第四阶CN方法的四倍。

分数阶微积分研究非整数阶的导数和积分，从一种数学奇想发展为强大的建模工具。它能够捕捉记忆效应、遗传行为和多尺度传输，适用于许多现实世界系统。在过去二十年中，它被广泛用于建模粘弹性阻尼、多孔和生物介质中的异常扩散、长期金融相关性等现象。此外，当代发展还包括流行病学动力学和分布阶波型模型，以及分数变分问题的存在理论进展。

在这些分数模型中，二维Riesz空间分数阶非线性反应扩散方程引起了特别关注。这些方程用分数阶次代替标准拉普拉斯算子，引入了每个空间方向上的非局部Lévy型跳跃。反应项描述了局部生成、衰减或非线性相互作用。这些模型已被应用于跨血管药物输送、电化学存储、裂隙岩石水文学、生物粘性流、图像修复和微结构相演化等领域。

然而，Riesz算子的非局部特性带来了显著的数值挑战。它导致了密集的刚度矩阵，网格点之间存在强耦合，且需要在边界附近使用精细的空间分辨率。当存在非线性反应时，稳定性约束变得更加严格。这些挑战凸显了开发高阶、稳定且对边界鲁棒的数值方法的必要性，以解决二维Riesz分数阶反应扩散问题。

本文针对二维Riesz空间分数阶扩散和非线性反应扩散模型，提出了几种数值方法。Bu、Tang和Yang提出了一种基于后向差分离散化的全离散Galerkin有限元方案，并证明了其在分数Sobolev空间框架下的适定性、稳定性和收敛性。Zhu等人引入了一种用于Riesz空间分数阶Fisher方程的有限元方法，结合了CN时间步进和非结构化三角网格上的空间FEM，并建立了存在性、唯一性、稳定性和最优的L2误差估计。最近，Qu等人开发了一种用于非线性Riesz型反应扩散方程的第四阶有限差分方案，并设计了一种具有可证明网格独立收敛性和L2稳定性的τ矩阵预处理共轭梯度求解器。我们还注意到，最近的一些分布阶Riesz框架旨在处理相关设置中的紧密耦合时间-空间动力学。

尽管现有的方案已经提高了二维Riesz空间分数阶问题的数值处理能力，但它们在应用于光滑基准问题时，与当前求解器相比存在明显局限。Bu、Tang和Yang的Galerkin FEM依赖于一阶时间步进和低阶元素，导致仅具有代数收敛性，并且由于算子的非局部性而产生密集矩阵。Zhu的方案提高了时间精度，但其在非结构化网格上使用的一阶元素需要大量的自由度来达到固定的误差水平。Qu的方案引入了第四阶有限差分方案和高效的预处理方法，但其空间收敛性仍然是代数的，并且受限于规则网格。相比之下，我们的求解器在时间上具有二阶精度，在空间上实现了接近指数级的收敛，所需的基函数远少于其他方法。其基于残差的边界惩罚方法避免了边界适配元素的使用，并且有望自然扩展到混合或不规则边界，尽管这一特性将在未来的工作中进一步探讨。

本文的贡献基于两个核心创新点。首先，每个牛顿更新通过一个PLM最小残差步骤获得，该步骤动态调整阻尼和边界惩罚参数，从而增强鲁棒性和收敛性。其次，狄利克雷条件通过在残差中最小化边界不匹配项来施加，这消除了对边界对齐基函数的需求，并为复杂几何形状提供了更大的灵活性。

本文的其余部分组织如下。第二部分介绍了支撑分析的加权Sobolev空间、迹结果和辅助引理。第三部分构建了二维Jacobi张量基，并证明了其关键的逼近估计。第四部分描述了中点/三阶时间积分器、Jacobi空间投影以及由此产生的牛顿线性化。第五部分汇总了残差公式和提供每个牛顿修正的自适应PLM最小残差迭代。第六部分分析了局部一致性并建立了时间、空间和迭代误差的完全分割先验误差界。第七部分建立了在轻微步长条件下的加权L2能量稳定性定理。第八部分通过一系列二维Riesz分数阶反应扩散测试验证了该方法，并将其与有限差分和有限元方案进行了基准比较。第九部分总结了全文，并概述了未来工作的方向。

在本部分中，我们首先介绍了一些基本概念。定义2.1给出了分数阶导数的基本定义。假设α在1和2之间，且u属于L1空间，左和右的Riemann-Liouville分数阶导数在区间（a，b）内定义。定义2.2给出了Riesz分数阶导数的定义，其中假设α在1和2之间，且u属于L1空间，左和右的Riemann-Liouville导数存在。定义2.3则涉及一些参数，这些参数定义了空间上的加权Sobolev空间。

接下来，我们讨论了二维一般Jacobi多项式基。假设μ_a和μ_b为大于-1的实参数，经典Jacobi多项式在区间[-1,1]上关于权重函数正交。这些多项式被定义为一组正交基，并且它们的正交性被证明。此外，我们还介绍了如何通过正交化方法得到一组归一化的基函数。

在时间离散化部分，我们采用了BDF2方法来离散时间导数。时间步长τ被定义为T除以时间网格点数N_t，并且时间层次被定义为t_i = iτ，其中i从0到N_t。对于任何时间依赖函数v，其在时间层次t_i的值被表示为v_i(x, y)。特别地，初始条件给出u_0(x, y) = u(x, y, 0) = ?(x, y)。对于i ≥ 2，我们采用BDF2公式进行时间更新，从而确保时间精度达到二阶。

在最小残差解部分，我们讨论了如何求解牛顿修正。我们限制搜索空间为有限维子空间，并在该子空间中寻找一个最小残差解。该解通过最小化合适的残差泛函来确定。我们证明了该解的存在性，并给出了其收敛性的相关分析。

在误差分析部分，我们讨论了如何评估求解器的精度。我们定义了全局最大点误差和边界最大误差，以衡量数值解与解析解之间的差异。这些误差指标为验证求解器的收敛性和稳定性提供了依据。我们还分析了误差的来源，并探讨了如何通过不同的方法进行误差控制。

在稳定性分析部分，我们讨论了如何确保数值方案的稳定性。我们将非线性算子分解为扩散部分和反应部分，并假设扩散部分在加权内积下严格耗散。同时，我们假设反应部分满足单侧增长条件。这些假设为证明方案的稳定性提供了理论基础。

在数值实验部分，我们展示了所提出方案的数值结果。通过不同的测试案例，我们验证了PLM-MRM方法的准确性。我们计算了全局最大点误差和边界最大误差，并比较了不同网格尺寸和时间步长下的误差变化。这些实验结果表明，该求解器在各种情况下都能保持高精度，并且在边界处理上具有良好的鲁棒性。

在结论部分，我们总结了本文的研究成果。我们提出了一种用于二维Riesz分数阶非线性反应扩散方程的高阶求解器，该求解器在时间上采用中点启动和BDF2更新，空间上使用加权Jacobi谱展开。每个牛顿步骤通过自适应的惩罚项Levenberg-Marquardt最小残差迭代进行修正，其中阻尼和边界惩罚参数被联合更新，从而正则化雅可比矩阵并弱施加狄利克雷条件。这些方法在提高计算效率和数值稳定性方面表现出色，适用于复杂和高维问题。

本文的作者贡献声明指出，Chaoyue Guan负责撰写初稿、软件开发和概念设计，而Jian Zhang负责撰写修订稿和概念设计。作者声明他们没有已知的财务利益或个人关系可能影响本文报告的研究。本文得到了湛江科技计划（项目编号：2025B01056）的支持。