一种求解奇异分数阶 Lane–Emden 方程的高精度长程混合方法

《Communications in Nonlinear Science and Numerical Simulation》：A High-Precision Long-Range Hybrid Method for Singular Fractional Lane–Emden Equations

　　

摘要 含Caputo导数的奇异分数阶Lane-Emden方程出现在具有记忆效应的径向对称模型中，但由于原点处的奇异系数以及局部级数近似收敛范围的有限性，其数值处理非常困难。本文发展了一种三阶段混合计算方案，该方案结合了用于在奇异点附近进行正则初始化的偶数索引分数阶幂级数展开、用于长程延拓的分数阶Adams-Bashforth-Moulton预报-校正方法，以及用于精度提升的Richardson外推法。所提出的方法避免了从奇异原点进行直接数值积分，并提供了一个实用框架，在该框架中局部解析结构与在扩展域上的稳定数值延拓相结合。我们推导了以级数启动误差和分数阶ABM离散化误差表示的全局误差估计。数值实验证实了该方法在经典和分数阶Lane-Emden问题中的精度和鲁棒性。在经典极限下，该方法能够准确重现近临界多边形指数的基准第一零点值，包括n=4.99。在分数阶情况下，计算结果表明在(α, n)平面中存在一个过渡边界，将有限半径剖面与单调无穷延拓解分离开来。这些结果表明，所提出的混合方案是高效可靠的工具，可用于对长计算区间上的奇异非线性分数阶模型进行高精度模拟。

引言 Lane-Emden方程

y″(x) + 2x y′(x) + y^n = 0, x > 0,

其中 y(0) = 1, y′(0) = 0，多边形指数 n ≥ 0，在天体物理学中作为恒星结构和自引力气体球的基本模型出现[1], [2]。在物理上，自变量 x 表示从恒星中心出发的无量纲径向距离，而 y(x) 表示气体的无量纲密度剖面。参数 n 表征状态方程，决定了恒星物质的热力学性质。天体物理学的主要目标是求解此方程以揭示恒星的内部结构，具体来说，即追踪密度如何从致密核心（y(0) = 1）衰减到恒星表面，恒星表面在数学上由密度消失（y(x) = 0）的第一个零点定义。

由于奇异项 2/x 的存在，仅对 n 的少数特殊值可获得精确解，大量数值和半解析技术已被提出用于近似求解，包括Adomian分解方法[3], [4], [5]、变分迭代方法[6]、同伦分析法[7]、同伦摄动法[8], [9]、带Padé近似的同伦摄动法[10]、差分变换方法[11]和幂级数递推方法[12]。特别是，文献[12]中证明(1)的解是偶函数，可写为 y(x) = Σ_{j≥0} a_j x^{2j}，从而得到更简单的递推关系，与标准幂级数展开相比系数数量减少了大约一半。此外，文献[12]中的数值比较表明，偶数索引级数可以比几种经典半解析方法（包括Adomian分解方法、同伦摄动法和差分变换方法）以较低的计算成本实现更高的精度。

近年来，经典微分模型的分数阶扩展已被广泛用于捕捉复杂系统中的记忆和非局部效应[13], [14], [15], [16]。对于经典Lane-Emden方程，通过将二阶导数替换为阶数为 2α 的Caputo导数，将一阶导数替换为阶数为 α 的Caputo导数，可得分数阶推广：

D_x^{2α} y + k x^α D_x^α y + y^n = 0, x > 0,

初始条件为

y(0) = y_0, (D_x^α y)(0) = 0,

其中 0 < α ≤ 1, k > 0，D_x^α 表示Caputo分数阶导数。此处 y_0 > 0 表示恒星中心的归一化无量纲密度，在标准Lane-Emden公式中 y_0 = 1。条件 (D_x^α y)(0) = 0 被采用作为与奇异项 x^{-α} D_x^α y 相关联的分数阶正则性条件。对于 1/2 < α ≤ 1，正则解满足

y(x) = y_0 + O(x^{2α}),

因此 y′(0) = 0，恢复了原点处的经典正则性条件。

对于 0 < α ≤ 1/2，y(0) = y_0 是所需的初始值，而 (D_x^α y)(0) = 0 作为与奇异结构相关联的附加正则性条件。当 α = 1 时，该条件退化为经典Lane-Emden条件 y′(0) = 0。

已开发了多种半解析技术用于求解分数阶Lane-Emden方程和相关的Lane-Emden型模型，包括修正差分变换方法[17]、同伦摄动法[18]、结合Jacobi多项式的幂级数方法[19]、带Laplace变换的残差幂级数方法[20]、结合Laplace变换的q-同伦分析法[21]以及采用Elzaki和Laplace变换的混合同伦摄动方案[22]。尽管这些方法在解决奇异点附近的解行为方面是有效的，但由于级数解收敛半径有限，这些方法通常局限于短计算区间。

从数值角度来看，也提出了众多方法用于分数阶Lane-Emden型方程。其中，分数阶Adams-Bashforth-Moulton (ABM) 预报-校正方案已成为Caputo型分数阶微分方程广泛接受的基准，并与其他数值求解器进行了广泛比较[23], [24]。其他代表性数值方法包括均匀和非均匀Haar小波方法[25]、结合Adomian分解方法的Haar小波[26]、基于Chebyshev正交函数的配置方法[27]和正交Bernstein多项式运算矩阵[28]、结合有限差分策略的Chebyshev小波[29]、基于Legendre多项式的数值方案[30]、结合拟线性化的第二类Chebyshev小波[31]、Chebyshev小波方法[32]以及带拟线性化的Legendre小波[33]。

尽管这些数值方案能够在固定有界域上近似解，但许多方法依赖于全局离散化，需要在整个计算区间上预先选择基函数或空间网格。因此，它们在奇异点 x = 0 附近的精度对离散化参数高度敏感，通常需要进行大量细化才能解析局部行为。此外，扩展计算域通常需要成比例地增加配置节点或基函数，导致更高的计算成本和降低的效率。另外，纯数值方案可能由于奇异系数而在原点附近遭受不稳定性或精度损失。这个问题在分数阶情况下更为显著，因为记忆效应导致在扩展域上累积误差，且局部解析结构没有明确地纳入数值构造中。

尽管已提出了广泛的半解析和数值技术用于分数阶Lane-Emden方程，但大多数报道的结果仍局限于短计算区间，通常为 x ∈ [0, 1]。据我们所知，目前尚无统一且鲁棒的框架，能够在保留奇异点附近解析结构的同时，在统一的数值方法内实现稳定和高精度的解延拓。

主要贡献 相比之下，本工作采用了一种根本不同的策略来同时解决初始奇异性和长程积分的挑战。我们不依赖全局基函数离散化或针对性网格细化，而是提出了一种鲁棒的三阶段混合框架。本文的主要贡献总结如下：

• 我们推导了奇异分数阶Lane-Emden方程的偶数索引分数阶幂级数表示。所得递推消除了所有奇数索引系数，并提供了解析正则化的奇异点附近局部近似。

• 我们构建了三阶段混合计算框架，将偶数索引局部级数与分数阶Adams-Bashforth-Moulton预报-校正延拓相结合。这种耦合避免了从奇异原点进行直接数值积分，并实现了稳定的长程模拟。

• 我们将Richardson外推法纳入混合框架以减少分数阶ABM延拓的导数离散化误差。这在不需要禁止性精细网格的情况下显著提高了精度。

• 我们建立了全局误差估计，表明数值误差由级数启动误差和分数阶ABM离散化误差控制。

• 我们将方法与经典基准值进行了验证，并用于在扩展域上计算近临界Lane-Emden剖面。此外，我们识别了分离有限半径和无穷延拓解区域的分数阶过渡边界。

本文结构如下。在第2节中，我们介绍关于Caputo导数的必要预备知识。第3节提出了我们策略的第一阶段，在其中我们严格分析了分数阶Lane-Emden方程并推导了用于初始化的偶数索引递推关系。第4节概述了第二阶段：将解析级数与分数阶Adams-Bashforth-Moulton (ABM) 方案相结合的混合计算框架。第5节详细阐述了第三阶段，重点是通过Richardson外推法进行精度提升。在第6节中，我们通过数值实验展示了方法的性能，包括近临界区域（n = 4.99）中的稳定性分析和临界边界 α_crit 的识别。本文在第7节中总结我们的发现。

Section 片段

预备知识 我们回顾Caputo分数阶导数的定义和基本性质[13], [14], [15], [16]。我们还给出了广义幂级数中非线性项 y^n 的J.C.P. Miller公式[34], [35], [36]。

定义 2.1 Riemann-Liouville分数阶积分 对于 μ > 0 和 f ∈ L^1[0, b]，阶数为 μ 的Riemann-Liouville分数阶积分定义为

(I^μ f)(x) = 1/Γ(μ) ∫_0^x (x-t)^{μ-1} f(t) dt, x ∈ [0, b]。

定义 2.2 Caputo导数 设 μ ∈ (0, 2] 且令 m = ?μ? ∈ {1, 2}。阶数为 μ 的Caputo分数阶导数定义为

D_x^μ y(x) = 1/Γ(m-μ) ∫_0^x (x-t)^{m-μ-1} y^{(m)}(t) dt,

前提是……

分数阶幂级数解 定理 3.1 设 0 < α ≤ 1, k > 0, n > 0 和 y_0 > 0。则问题(2)-(3)存在幂级数解

y(x) = Σ_{m=0}^{∞} c_m x^{mα}, c_0 = y_0, c_1 = 0,

且对于 m ≥ 0，系数满足

c_{m+2} = -F_m · Γ(mα+1)Γ((m+1)α+1) / [Γ((m+2)α+1)(Γ((m+1)α+1) + kΓ(mα+1))],

其中 F_m 由引理 2.2 给出。

证明 将

y(x) = Σ_{m=0}^{∞} c_m x^{mα}

代入(2)并对级数逐项应用引理 2.1，得

D_x^α y(x) = c_1/Γ(α+1) + Σ_{m=0}^{∞} c_{m+2} · Γ((m+2)α+1)/Γ((m+1)α+1) · x^{(m+1)α},

且

D_x^{2α} y(x) = Σ_{m=0}^{∞} c_{m+2} · Γ((m+2)α+1)/Γ(mα+1) · x^{mα}。

因此，

k x^α D_x^α y(x) = k c_1 Γ(α+1) x^{-α} + k Σ_{m=0}^{∞} c_{m+2} Γ((m+2)α+1)/Γ((m+1)α+1) …

Phase 2: 扩展域的混合方法 定理 3.4 中建立的偶数索引分数阶幂级数在奇异点 x = 0 附近提供了高度精确的解。然而，对于大的计算域，级数的收敛半径构成了根本性限制。为了将解稳健地扩展到 L 任意大的 [0, L] 区间，我们开发了一种混合策略，将第一阶段中的偶数索引级数初始化与分数阶Adams-Bashforth-Moulton (ABM) 预报-校正延拓[23], [24]相结合。

Phase 3: 通过Richardson外推法提升精度 尽管第4节中的混合方法在扩展域上提供了稳定的解，但其精度最终受到分数阶Adams-Bashforth-Moulton (ABM) 预报-校正方案离散化误差的影响。对于 0 < α < 1，在适当的正则性假设下，分数阶Adams方法具有收敛阶 O(h^{1+α}) [24]。对分数阶Adams方法的更详细分析表明，误差结构通常涉及整数和分数阶……

长程渐近行为的考察 所提出混合方案的一个显著优势是其卓越的稳定性，使其能够可靠地积分远超文献中通常探索的域。在本节中，我们将计算域扩展到 L ≥ 100，以研究解的渐近行为，重点关注分数阶数 α 与Lane-Emden方程振荡特性之间的相互作用。所有结果均使用混合方法结合两级……计算得出。

结论 本文发展了一种针对含Caputo导数的奇异分数阶Lane-Emden方程的三阶段混合计算方案。所提方法的主要新颖性在于三个互补组件的解析正则化耦合：奇异原点附近的偶数索引分数阶幂级数初始化、扩展域上的分数阶Adams-Bashforth-Moulton延拓，以及用于提高实际数值……的两阶段外推细化。

CRediT作者贡献声明 Sita Charkrit：撰写—审阅与编辑，撰写—原始草稿，可视化，验证，监督，软件，资源，项目管理，方法论，调查，经费获取，形式分析，数据整理，概念化。

利益冲突声明 作者声明不存在利益冲突。

Sita Charkrit 泰国清迈，清莱大学理学院数学系