[论文解读] Rational approximation to the fractional Laplacian operator in reaction-diffusion problems

Lidia Aceto, Paolo Novati|arXiv (Cornell University)|Jul 14, 2016

Fractional Differential Equations Solutions参考文献 13被引用 62

一句话总结

本文提出了一种有理逼近方法，以高效求解在有界区域上具有空间分数阶拉普拉斯算子的分数阶反应-扩散方程。通过将密集的分数阶矩阵幂 $ L^{rac{eta}{2}} $ 近似为两个带状矩阵 $ M^{-1}K $ 的乘积，该方法将半线性系统转化为稀疏且计算高效的格式，从而在一维和二维空间中实现精确且可扩展的求解，相较于标准矩阵转移技术具有显著的计算优势。

ABSTRACT

This paper provides a new numerical strategy to solve fractional in space reaction-diffusion equations on bounded domains under homogeneous Dirichlet boundary conditions. Using the matrix transform method the fractional Laplacian operator is replaced by a matrix which, in general, is dense. The approach here presented is based on the approximation of this matrix by the product of two suitable banded matrices. This leads to a semi-linear initial value problem in which the matrices involved are sparse. Numerical results are presented to verify the effectiveness of the proposed solution strategy.

研究动机与目标

• 解决涉及离散拉普拉斯矩阵的密集矩阵幂的分数阶反应-扩散方程求解所带来的高计算成本问题。
• 克服标准矩阵转移技术效率低下的问题，该技术会产生需要昂贵求逆或分解的密集矩阵。
• 开发一种与维度无关的数值策略，在保持精度的同时减少内存占用和计算负载。
• 在不改变核心求解方法的前提下，将该方法推广至二维和三维空间。
• 为分数阶PDE提供一种实用且可扩展的替代方案，替代通过轮廓积分或Krylov方法计算矩阵函数。

提出的方法

• 将分数阶矩阵幂 $ L^{\alpha/2} $ 近似为 $ M^{-1}K $，其中 $ M $ 和 $ K $ 是通过 $ z^{\alpha/2 - 1} $ 的有理逼近得到的带状矩阵。
• 使用高斯-雅克比求积规则对矩阵函数的积分表示进行离散化，以实现稳定且精确的有理逼近。
• 通过基于高斯-雅克比规则选择最优极点来构建有理逼近，其中参数 $ \tau $ 经调优以实现最佳收敛性和精度。
• 将原始的半线性ODE系统 $ \frac{du}{dt} = -\kappa_\alpha h^{-\alpha} L^{\alpha/2} u + f $ 转化为 $ M \frac{du}{dt} = -\kappa_\alpha h^{-\alpha} K u + M f $，以保持稀疏性。
• 利用 $ M $ 和 $ K $ 的带状结构，通过标准ODE求解器（如 `ode15s`）实现高效的时域积分。
• 在一维和二维问题中应用该方法，分别采用标准的五点和三点有限差分格式（3点和5点）。

实验结果

研究问题

• RQ1 能否构造出对矩阵函数 $ L^{\alpha/2} $ 的有理逼近，使得所得的矩阵分解 $ M^{-1}K $ 保持稀疏性，并支持高效的时域积分？
• RQ2 高斯-雅克比求积参数 $ \tau $ 的选择如何影响有理逼近的精度和收敛速率？
• RQ3 与标准矩阵转移技术相比，该方法在内存使用和求解时间方面具有哪些计算优势？
• RQ4 该有理逼近策略能否在不损失效率或改变求解框架的前提下，推广至二维和三维空间？
• RQ5 逼近误差如何随求积阶数 $ k $ 变化？数值实验中观察到的收敛行为如何？

主要发现

• 该有理逼近方法在精度上可与矩阵转移技术相媲美，对于 $ k = 3 $ 和 $ k = 5 $，最大误差分别处于 $ 10^{-6} $ 和 $ 10^{-5} $ 量级。
• 对于 $ N = 40 $ 的二维问题，为达到相同精度，矩阵转移技术的计算成本是该有理逼近方法的三倍。
• 有理逼近的误差随求积阶数 $ k $ 增加而减小，$ k = 5 $ 时的误差显著低于 $ k = 1 $，该结果在示例3中 $ t = 0.5 $ 处得到验证。
• 该方法成功捕捉了扩散对分数阶阶数 $ \alpha $ 的依赖关系，在示例2中，$ \alpha = 1.1 $ 和 $ \alpha = 1.9 $ 时观察到不同的解分布。
• 该方法在不同初始条件和源项下均保持稳定和高精度，包括示例1、3和4中的非平凡解析解。
• 该方法具有鲁棒性和可扩展性，从一维扩展到二维问题时无需改变算法结构，充分体现了其与维度无关的特性。