Chapter 6: Numerical Methods for Fractional PDEs

1. 作者

Christian Klein（通信作者），法国勃艮第大学（Université de Bourgogne）数学研究所，勃艮第数学研究所（IMB, UMR 5584），同时是法国大学研究院（Institut Universitaire de France）成员。Email: Christian.Klein@u-bourgogne.fr
Nikola Stoilov，勃艮第大学数学研究所，Email: Nikola.Stoilov@u-bourgogne.fr

Klein 是法国分数阶 PDE 数值分析的领军人物，特别在谱方法与Riemann 面方法的结合上独树一帜。本章是 Ch5（解析方法）的数值对照——同一类方程（分数阶 KdV、mKdV、NLS），但用完全不同的工具。本章由 EIPHI 项目（ANR-17-EURE-0002）和 EU Horizon 2020 IPaDEGAN 项目资助。

2. 内容概述

本章是全书数值方法的核心章节，专门讨论全实数轴上分数阶 Laplacian 的数值计算。三大部分：(1) 理论准备：分数阶导数的 Riemann-Liouville / Caputo / Riesz 定义，与Riemann 曲面的连接——分数阶积分的代数奇异性自然对应到 $\mathbb{Z}_q$ 曲线；(2) 数值方法综述：DFT/FFT、Chebyshev、多区域谱方法、Rational Mapping（$y = \tan(\theta/2)$）、Jacobi 展开、Runge-Kutta 卷积积分、有限差分（22 阶复平面 stencil）等；(3) 多区域谱方法的详细实现：把实数轴分为 3 个区域（$x < a$、$a \le x \le b$、$x > b$），用局部参数 $\xi_\pm = 1/(\pm x)^{1/q}$ 在 $\mathbb{Z}_q$ 曲面上处理积分；(4) 应用：用 Newton-GMRES 求解分数阶 KdV/mKdV/NLS 的孤子方程，比较 FFT 与多区域方法的精度——FFT 在 $\alpha$ 接近能量临界值 $1/3$ 时退化到 $10^{-3}$（plotting accuracy），多区域方法保持 $10^{-5}$。

3. 核心方程与概念

3.1 分数阶导数与积分（Riemann-Liouville / Caputo / Riesz）

Riemann-Liouville 分数阶导数（$\alpha > 0$，$n = [\alpha]$）： $$_0\partial_x^\alpha u(x) = \frac{1}{\Gamma(n-\alpha)}\partial_x^n \int_0^x \frac{u(y)}{(x-y)^{\alpha+1-n}}dy$$ 是分数阶积分的 $n$ 阶经典导数。

Caputo 分数阶导数： $$_x\partial^\alpha u(x) = \frac{1}{\Gamma(n-\alpha)}\int_0^x \frac{\partial_y^n u(y)}{(x-y)^{\alpha+1-n}}dy$$ 当 $n = 0$ 时两者一致。

Riesz 分数阶积分（全实数轴）： $$\tilde I_\alpha u(x) = \frac{1}{2\Gamma(\alpha)\sin(\alpha\pi/2)}\int_{-\infty}^\infty \frac{\text{sgn}(x-y)u(y)}{|x-y|^{1-\alpha}}dy$$ Fourier 变换：$\mathcal{F}(\tilde I_\alpha u) = -i\,\text{sgn}(k)|k|^{-\alpha}\hat u(k)$。

分数阶导数（关键）（全实数轴，$\alpha \in (0,1)$）： $$D^\alpha u(x) = \frac{1}{2\Gamma(1-\alpha)\sin((1-\alpha)\pi/2)}\partial_x \int_{-\infty}^\infty \frac{\text{sgn}(x-y)u(y)}{|x-y|^\alpha}dy$$ Fourier 形式：$\mathcal{F}(D^\alpha u) = |k|^\alpha \hat u(k)$。

等价的 Cauchy 主值形式： $$D^\alpha u = \frac{2^{2\alpha-1}\Gamma(1/2+\alpha/2)}{\sqrt\pi\Gamma(1-\alpha/2)}\int_\mathbb{R} \frac{u(x) - u(x+y)}{|y|^{1+\alpha}}dy$$ 缺点：$u(x)$ 单独存在使被积函数 $y \to 0$ 时奇异，数值上需要处理极限。

3.2 分数阶积分与 Riemann 曲面

核心观察：$I_\pm(x) = \int_{-\infty}^x u(y)(x-y)^{-\alpha}dy$ 和 $\int_x^\infty u(y)(y-x)^{-\alpha}dy$ 是Abelian 积分——当 $\alpha$ 非有理时，定义在无穷亏格的 Riemann 曲面上。

有理 $\alpha = p/q$（$p, q$ 互素，$p < q$）时： - 曲面是紧致的，亏格有限； - 由代数曲线 $\mu^q = K - x$ 定义（$\mathbb{Z}_q$ 曲线），是复平面的 $q$ 叶覆盖，分支点在 $x$ 和 $\infty$； - 在 $x$ 附近：局部参数 $\lambda_x = (K-x)^{1/q}$； - 在 $\infty$ 附近：局部参数 $\xi_\pm = (\pm 1/x)^{1/q}$，双参数处理 $\pm x \to \infty$ 的不对称性。

孤子的衰减条件（Frank-Lenzmann [X1]）：$\lim_{|x|\to\infty} u(x)|x|^{1+\alpha} < \infty$——保证积分在 $\infty$ 收敛。

紧支集例子（特征函数 $u = \chi_{[a,b]}$）： $$D^\alpha \chi_{[a,b]}(x) \propto (x-a)^{-\alpha} + (b-x)^{-\alpha}$$ 奇异点在端点。用局部参数 $\lambda_a = (x-a)^{1/q}$，$\lambda_b = (b-x)^{1/q}$，积分可化为 $D^\alpha u \sim \lambda_a^{-p} + \lambda_b^{-p}$（Puiseux 展开）。

3.3 数值方法综述

DFT/FFT 方法：在 $[-L\pi, L\pi]$ 上离散 $N$ 个点，用 FFT 计算 $D^\alpha u = \mathcal{F}^{-1}[|k|^\alpha \hat u(k)]$，复杂度 $O(N\log N)$。优点：实现简单、平台优化好；缺点：被截断在 $[-L\pi, L\pi]$ 内，慢衰减函数（孤子）的尾部积分 $\int_{L\pi}^\infty$ 被忽略。精度对 $\alpha$ 接近临界值时退化到 $10^{-3}$。

Chebyshev 谱方法：$T_n(x) = \cos(n\arccos x)$，Chebyshev 配置点 $\xi_m = \cos(m\pi/N)$，系数通过快速余弦变换。导数用 Chebyshev 微分矩阵。

Rational Mapping（$y = \tan(\theta/2)$）—— [X2, X3] 用： - 把 $\mathbb{R}$ 映射到单位圆 $\theta \in [-\pi, \pi]$； - 展开基函数 $\phi_n(y) = (1+iy)^n/(1-iy)^{n+1}$ 代替三角函数； - 把 $\mathbb{R}$ 上的无界分数阶积分化为无截断的 FFT。 - 问题：$\phi_n, \chi_n$ 难以准确表示孤子的 $q$ 次根渐近行为。

Jacobi 多项式（分形分形多项式, GJF）—— [X4]：用 $t^\beta P_N(t)$ 作为基，Jacobi 多项式的分数阶积分可用 hypergeometric 性质解析表达。优点：稀疏矩阵；缺点：需要先验知道 $\beta$（奇性指数）。

Runge-Kutta 卷积积分（CoQ）：$O(x)$ 复杂度计算 Caputo 导数。

有限差分（Fornberg-Piret）：22 阶格式，$5\times 5$ 复平面 stencil，对解析函数高精度。

自适应有限元：边界附近加密网格（处理 $u \in [a,b]$ 端点奇异性）。

3.4 多区域谱方法（核心创新）

思路：把 $\mathbb{R}$ 分为 3 个区域： - I: $x < a < 0$（"远"负无穷） - II: $a \le x \le b$（"近"区域） - III: $x > b > 0$（"远"正无穷）

关键技巧：用 $\mathbb{Z}_q$ 曲面上的局部参数 $$\xi_- = (-x)^{-1/q}, \quad \xi_+ = x^{-1/q}$$ 使被积函数在 $\infty$ 处光滑。具体地，定义 $$u_I(\xi_-) := u(x)(-x)^{1+\alpha}, \quad x < a$$ $$u_{II}(x) := u(x), \quad a \le x \le b$$ $$u_{III}(\xi_+) := u(x)x^{1+\alpha}, \quad x > b$$

把 $I_-$ 分解为两个子积分（在 $\xi_-^q/2$ 处切分），用 $s = -1/y$ 变换到 $\xi_-$ 局部参数。每个子积分的被积函数光滑，用 Clenshaw-Curtis 算法计算。

Clenshaw-Curtis 积分：等价于 Gauss 积分，展开为 Chebyshev 多项式 $$\int_{-1}^1 h(l)dl \approx \sum_m w_m h(\xi_m)$$ 加权 Chebyshev 配置点。配合重心插值（barycentric interpolation）实现 $u$ 在任意点的查询。

精度对比（来自 [X5]）： - FFT 方法：$\alpha$ 远离临界时 $10^{-6}$，$\alpha \to 1/3$ 时退化到 $10^{-3}$（plotting accuracy）； - 多区域方法：$10^{-5}$ 在 $\alpha = 0.4$ 仍保持。

3.5 应用：分数阶 KdV / mKdV / NLS 孤子

行波方程（共同形式）： $$D^\alpha Q_c + cQ_c - \kappa Q^n = 0$$ $n = 2$（KdV）或 $3$（mKdV/NLS）。

衰减行为（Frank-Lenzmann [X1]）：$|Q(x)| \sim |x|^{-(1+\alpha)}$，比 KdV 的指数衰减慢——是数值困难的核心原因。

能量（守恒量）： - FKdV: $H(u) = \int_\mathbb{R} \left[\frac{1}{2}|D^{\alpha/2}u|^2 - \frac{1}{6}u^3\right]dx$ - mKdV: $\frac{1}{2}|D^{\alpha/2}u|^2 \mp \frac{1}{6}u^4$ - FNLS: $\frac{1}{2}|D^{\alpha/2}u|^2 \mp \frac{1}{4}|u|^4$

临界指数（标度不变性 $u_\lambda(x,t) = \lambda^\alpha u(\lambda x, \lambda^{\alpha+1}t)$）： - FKdV：$L^2$ 临界 $\alpha = 1/2$，能量临界 $\alpha = 1/3$； - mKdV：$L^2$ 临界 $\alpha = 1$，能量临界 $\alpha = 1/2$； - FNLS：$L^2$ 临界 $\alpha = 1$，能量临界 $\alpha = 1/2$。

数值方法（构造孤子）：用 Newton-GMRES 求解非线性系统 $$Q^{(m+1)} = Q^{(m)} - \text{Jac}(F)^{-1}|_{Q^{(m)}} F(Q^{(m)})$$ $\alpha$ 从 $0.9$ 开始，用 BO 孤子（$\alpha = 1$ 解析解 $Q = 4/(1+z^2)$）作初值，continuation 到 $\alpha = 0.4$（多区域方法可达 $0.4$，FFT 难以接近 $1/3$ 临界值）。

孤子的形状随 $\alpha$ 演化（图 4, 9, 10, 11）： - $\alpha \downarrow$：中心峰更尖锐，但尾部更平（$|x|^{-(1+\alpha)}$ 衰减变慢）； - $\alpha \to 1/3$（FKdV 临界）：孤子趋向"$\delta$ 函数样"——数值上极难； - 多区域方法的 $\xi_\pm$ 表示（图 8）展示Puiseux 级数结构（无穷远附近是 $1/q$ 次根）。

4. 关键结论

分数阶 Laplacian 的数值计算有两个等价路径：(1) Fourier 域（FFT），实现简单但被有界域截断影响慢衰减函数；(2) 物理空间（多区域 + Riemann 曲面），实现复杂但精度高。对孤子问题，多区域方法不可替代。
分数阶积分自然定义在 Riemann 曲面上——$\alpha = p/q$ 时曲面是代数 $\mathbb{Z}_q$ 曲线，亏格有限。局部参数 $\xi_\pm = 1/(\pm x)^{1/q}$ 解决了 $\infty$ 处的 $q$ 次根奇异性。
FFT 局限：精度上限约 $10^{-3}$（plotting accuracy），对时间演化或Newton 迭代会累积。多区域方法可保 $10^{-5}$，但实现复杂。
孤子的 $|x|^{-(1+\alpha)}$ 衰减是分数阶 PDE 的本征特征——比经典 KdV 的指数衰减慢得多，使数值计算本质困难。
临界 $\alpha$ 值（FKdV $\alpha_c = 1/3$，FNLS $\alpha_c = 1/2$）对应能量临界——$\alpha < \alpha_c$ 无光滑孤子，是分数阶 PDE 的"零孤子极限"。
BO 极限（$\alpha = 1$）的精确解 $Q = 4/(1+z^2)$ 是 Newton continuation 的基准。
mKdV 存在"反相孤子"（$\sigma = \pm 1$）和"呼吸子"——可积系统的可加性结构。

5. 挑战和开放性问题

能量临界极限（$\alpha \to 1/3$ 或 $1/2$）——孤子趋向 $\delta$-分布，当前所有数值方法都失效。多区域方法可能突破，但需动态重缩放配合 Newton 迭代。
时间演化的精度问题——FFT 在 $\alpha$ 接近临界时累积误差。多区域 + 时间 FFT 的组合（如 [X6]）需要陡降路径识别。
Blow-up 数值描述——对 $L^2$-临界 gKdV（$\alpha = 1/2$），自相似 blow-up 已知；分数阶情形的 blow-up 几何与 blow-up profile 数值上未解决。
多维分数阶 Laplacian——$d = 2, 3$ 情形的多区域方法实现困难，Riemann 曲面在多维下的推广仍开放。
$\alpha$ 无理数的数值处理——多区域方法假设 $\alpha = p/q$，无理 $\alpha$ 对应无穷亏格 Riemann 曲面。
长时间演化——孤子解的分辨率猜想（soliton resolution）需要在大时间尺度上保持精度，这对任何数值方法都是挑战。
Breathers（mKdV 的呼吸子解）——Klein-Saut-Wang (2022) 给出存在性，但数值构造没有系统化方法。
FNLS 在 1/2 < $\alpha$ < 1 范围内——孤子解形状的细节不明确。

6. 个人反思与批判性分析

优点： 1. 数值实验与解析理论完美呼应——多区域方法的精度（$10^{-5}$）与 Ch5 的小振幅展开（$O(a^4)$ 精度）匹配，数值验证了解析定理。 2. Riemann 曲面的引入是真正的创新——把"奇异性在 $\infty$"的物理问题转化为"代数曲线在 $\mathbb{Z}_q$ 上"的纯数学对象，给出严格的拓扑/代数结构。 3. 方法比较客观——FFT 优点（$O(N\log N)$）和缺点（$10^{-3}$ 精度极限）都讲清楚，不夸张多区域方法。 4. 多区域 + Chebyshev + 重心插值的组合在工程上完整可行——不是空中楼阁。 5. 孤子的 Puiseux 展开（图 8）是分数阶 PDE 数值方法的独特视角——经典 PDE 的孤子在 $\infty$ 是指数衰减，没有这种代数奇性。 6. 从 FKdV 推到 mKdV 推到 FNLS 的应用面覆盖广，展示了方法的普适性。

缺点： 1. Riemann 曲面概念对非代数几何背景的读者理解门槛高——许多读者可能跳过 §2.2。 2. 多区域方法的实现细节（如 Clenshaw-Curtis 切分点选择、$\xi_-^q/2$ 的具体计算）不充分——读者难以复现。 3. FFT 与多区域方法的实际效率比较缺失——只给了精度，没给计算时间对比。在 $\alpha = 0.8$ 时 FFT 是否已经"够好"？多区域方法额外开销多少？ 4. Blow-up 数值只占 1 段，多维分数阶 PDE（$d \ge 2$）完全没有数值例子——这些是当代计算数学的活跃方向。 5. 机器学习/神经网络方法（如 Physics-Informed Neural Networks, PINN）近年在分数阶 PDE 数值上有显著进展，本章完全未涉及。 6. 与 Ch3 的关系——Ch3 §7.1-7.2 已详细讨论时间分数阶和空间分数阶的数值方法，但未涉及本章的"全实数轴 + 孤子"问题。两章应互相引用但都未做。 7. $\alpha$ 接近 $1/3$ 临界值时多区域方法的实际表现未充分测试——作者承认这是"future research"。

对本人研究（血管生物力学 / SMC 介导的血管重塑）的关联： - 多区域谱方法的思想可应用于血管壁的层状结构建模——内弹力板、中膜、外膜用不同区域，每区域用不同本构模型（Ogden / Neo-Hookean / Fung）。这种"区域自适应"在血管数值中已有先例（如 Goriely 2017 的分层管壁模型），但与 Riemann 曲面无关。 - $\mathbb{Z}_q$ 曲线的代数奇异性分析对血管壁粘弹性的"分数阶 Prony 级数"有启发——长时记忆对应"代数衰减"而非指数衰减。 - FKdV 的 $\alpha = 1/3$ 临界值与"血管平滑肌细胞体外培养"实验中观察到的"粘弹性指数 = 1/3" 现象可对比。 - Blow-up 的自相似性可类比动脉瘤破裂前兆的临界动力学——血管壁的弱化可能是自相似 blow-up 过程。 - BO 方程的 Hilbert 变换与医学成像中 X 射线 CT 的 Radon 变换有数学同构——分数阶算子与 Radon 变换都涉及"非局部积分"。

7. 重要参考文献（按出现顺序编号）

[X1] R.L. Frank, E. Lenzmann. Acta Math. 210, 261 (2013). —— FKdV 孤子 $|x|^{-(1+\alpha)}$ 衰减的精确证明。
[X2] A. Iserles, S.P. Nørsett. —— Rational mapping 处理无界域。
[X3] C. Klein, C. Sparber, P. Markowich. Proc. R. Soc. A 470, 20140364 (2014). —— 分数阶 NLS 的数值解。
[X4] (Jacobi 多项式 / GJF 的原始工作) —— [X4] 出现在文中引用。
[X5] C. Klein, N. Stoilov. —— 多区域谱方法的奠基论文。
[X6] (B. Fornberg, C. Piret 复平面 stencil 22 阶方法) —— [X6] 出现在文中引用。
[X7] C. Klein, J.-C. Saut, Y. Wang. Nonlinearity 35, 1170 (2022). —— 分数阶 mKdV 的精确解。
[X8] R. Hilfer. Threefold Introduction to Fractional Derivatives (Wiley-VCH, 2008). —— 分数阶导数综述。
[X9] D. Lannes. The Water Waves Problem (AMS, 2013). —— 水波的非局部方程。
[X10] C. Klein, F. Linares, D. Pilod, J.C. Saut. Stud. Appl. Math. 140, 133 (2018). —— Whitham 方程。
[X11] J.-P. Berland, C. Klein, C. Stoilov. —— 分数阶 PDE 的多区域数值（与本章作者合作）。
[X12] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang. Spectral Methods in Fluid Dynamics (Springer). —— 谱方法经典教材。
[X13] A. Iserles. —— Chebyshev 展开与分数阶算子。
[X14] C.E. Kenig, G. Ponce, L. Vega. —— KdV 的适定性（Ch5 也引）。
[X15] Y. Saad. Iterative Methods for Sparse Linear Systems (SIAM, 2003). —— GMRES 算法的标准参考。
[X16] J. Fröhlich, B.L.G. Jonsson, E. Lenzmann. Commun. Math. Phys. 274, 1 (2007). —— Boson-star 模型。
[X17] V.E. Tarasov. Fractional Dynamics (Springer, 2010). —— 分数阶动力学综述。
[X18] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo. —— 分数阶积分和导数。
[X19] D. Mumford. Tata Lectures on Theta I, II (Birkhäuser, 1983). —— $\theta$ 函数与 Riemann 曲面。
[X20] P.G. Grinevich, P.M. Santini. —— 分数阶 KdV 孤子的小振幅展开（与 Ch5 [X1] 同）。

本章由主 agent 亲自从 PDF 提取并撰写，无子任务委托。
撰写日期：2026-06-01