分数阶扩散方程四阶差分逼近格式的分析与稳定性研究

《The ANZIAM Journal》：ANALYSIS OF FOURTH-ORDER DIFFERENCE APPROXIMATION SCHEMES FOR FRACTIONAL DIFFUSION EQUATIONS

本文针对分数阶扩散方程数值求解中高阶格式稳定性难题，提出基于Lubich型二阶逼近的两种四阶差分格式。研究结合Crank-Nicolson时间离散，通过加权移位Grünwald差分(WSGD)和准紧致逼近(QCA)方法构建新算法，严格证明其在α∈[1,2]范围内的无条件稳定性与收敛性。数值实验验证了四阶收敛精度，为分数阶偏微分方程高精度计算提供重要理论工具。

　　

分数阶微积分作为描述反常扩散、非局部效应等复杂物理过程的核心工具，近年来在物理、工程和金融等领域展现出强大应用价值。然而，分数阶微分方程的数值求解始终面临严峻挑战：经典整数阶微分方程的稳定数值格式直接移植到分数阶情形时会出现难以预料的不稳定性。特别是Grünwald差分(GD)逼近虽具有一阶精度，但与隐式Euler法或Crank-Nicolson(CN)格式结合时会产生发散解，这严重制约了分数阶模型的实际应用。

为解决这一难题，研究者们先后发展了移位Grünwald差分(SGD)、加权移位Grünwald差分(WSGD)等改进格式。但现有高阶方法（如三阶WSGD）在临界分数阶阶数α∈(1,2]范围内仍存在稳定性缺陷，亟需构建兼具高精度和强稳定性的新数值框架。本文正是在此背景下，通过创新性地融合Lubich型生成函数与准紧致逼近技术，开发出两种新型四阶数值格式，为分数阶扩散问题的高效计算开辟了新途径。

本研究主要采用以下关键技术方法：基于Lubich型二阶生成函数W_2,r(z)的移位逼近理论，通过不同移位组合(r₁,r₂)=(1,0)和(1,-1)构建加权移位格式；结合准紧致差分(QCA)技术消除主导误差项，将精度提升至四阶；采用Crank-Nicolson时间离散格式保持时间方向二阶精度；应用Toeplitz矩阵生成函数理论进行稳定性分析，并通过模型问题验证算法性能。

4. 分数阶扩散方程的逼近

将构建的四阶格式应用于空间分数阶扩散方程?u/?t=K₁D_x-^αu+K₂D_x+^αu+f(x,t)，其中扩散系数K₁,K₂为非负常数。通过引入算子P_h=I+h²?₂δ_h^{2对空间导数进行预处理，结合Crank-Nicolson格式得到全离散格式：(P_h-τ/2·L_h)Um+1=(P_h+τ/2·L_h)Um+τP_hFm+1/2，截断误差为O(τ3+τh4)。}

5. 稳定性和收敛性

通过建立离散内积和范数体系，证明了算子P_h的自伴性和正定性，其满足1/5‖v‖≤(P_hv,v)≤9/5‖v‖的能量估计。关键创新在于通过生成函数G_r(x)=W_2,r(e^-ix)e^irx的极坐标分解，解析推导出CN-QCA1和CN-QCA2格式的实部恒为负值。

图1显示H_c,(1,0)(α,x)在α∈[1,2]区间内恒为负值，结合径向函数R_r(x)=(2sin(x/2))^αS_r^α(x)的正定性，证实了生成函数实部在全域的负定性。类似地，CN-QCA2格式的稳定性也通过图2的等高线分析得到验证。

收敛性分析通过误差能量估计完成，证明数值解与精确解的误差满足‖e^m‖≤25τ∑_l=0^m-1‖S^l‖的有界性，确保当网格尺寸h,τ→0时数值解一致收敛于精确解。

6. 数值结果

选取测试问题u(x,t)=x^m(1-x)^me^-t进行算法验证。表1-4分别展示了m=4,6,2,3时CN-QCA1和CN-QCA2的收敛性能。当解足够光滑(m≥4)时，两种格式均表现出清晰的四阶收敛特性；而当解正则性较低(m=2,3)时，收敛阶降至m阶，这与边界层效应的理论预期一致。DTTS迭代法在σ=1时仅需2次迭代即可达到10^-7容差，证实算法的高效性。

7. 结论

本研究成功构建了两种基于Lubich型二阶生成函数的四阶差分格式，通过严格的数学分析证明了其在α∈[1,2]范围内的无条件稳定性。创新性地将生成函数理论与矩阵分析工具相结合，解决了长期困扰分数阶数值计算的高阶格式稳定性问题。数值实验验证了理论结果，为分数阶偏微分方程的高精度计算提供了可靠工具。特别值得注意的是，该方法框架可扩展至更复杂的分数阶模型，为后续研究奠定了坚实基础。

这项工作的重要意义在于首次实现了分数阶扩散方程四阶格式的全范围稳定性突破，弥补了传统WSGD方法在临界参数区的理论缺陷。所建立的稳定性分析范式为分数阶数值方法的设计提供了新思路，对推动分数阶微积分在工程实际中的应用具有深远影响。