sangkuriang：用于 Korteweg-de Vries 孤立子模拟的伪谱 Python 库 sangkuriang: A pseudo-spectral Python library for Korteweg-de Vries soliton simulation
开源 Python 伪谱求解器，高效模拟浅水 KdV 方程孤立子动力学并验证守恒性。
前置知识
Korteweg-de Vries (KdV) 方程
KdV 方程是一个经典的非线性偏微分方程，形式为 $\partial_t u + \varepsilon u \partial_x u + \mu \partial_x^3 u = 0$。它描述了弱非线性、弱色散条件下长波在浅水中的传播行为。方程中 $\varepsilon u \partial_x u$ 是非线性项（导致波形变陡），$\mu \partial_x^3 u$ 是色散项（导致不同频率的波以不同速度传播）。当这两者达到平衡时，就会产生一种形状保持不变、碰撞后仅发生相移的特殊波——孤立子（soliton）。KdV 方程是完全可积的 Hamilton 系统，拥有无穷多守恒量，其中质量 $M = \int u \, dx$、动量 $P = \int u^2 \, dx$ 和能量 $E = \int \left(\frac{\varepsilon}{2}u^3 - \frac{3\mu}{2}(\partial_x u)^2\right) dx$ 是前三个物理意义明确的守恒律。
KdV 方程是本文求解器的核心对象，理解它的数学结构（非线性与色散的平衡、孤立子解、守恒律）是理解整个验证框架的基础。
伪谱方法（Pseudo-spectral method）
伪谱方法是一种求解偏微分离散化的高精度数值技术。核心思想是：在空间域上用离散傅里叶变换（DFT）将物理空间的函数变换到频率（波数）空间，在频率空间中微分运算变为简单的乘法——例如 $\partial_x u$ 对应 $ik_n \hat{u}_k$，$\partial_x^3 u$ 对应 $(ik_n)^3 \hat{u}_k$。然后再通过逆变换回到物理空间。对于光滑周期函数，这种方法能达到指数级收敛精度（spectral accuracy），远优于有限差分方法的多项式精度。计算复杂度为 $O(N \log N)$，借助 FFT 算法实现。
伪谱方法是本文求解器的核心数值技术，直接决定了空间离散化的精度和效率。理解它才能理解为什么该方法能精确保持孤立子结构和守恒律。
DOP853 自适应时间积分
DOP853（Dormand-Prince 8(5,3)）是一种八阶显式 Runge-Kutta 方法，嵌入五阶和三阶误差估计器用于自适应步长控制。每步需要 12 次函数求值，通过比较八阶和五阶解的差异来估计局部误差 $\text{err}_{\text{local}} = \|u^{(8)}_{n+1} - u^{(5)}_{n+1}\|$。若误差满足容差条件 $\text{err}_{\text{local}} \leq \text{rtol} \cdot \|u_{n+1}\| + \text{atol}$，则接受该步并按公式调整下一步步长；否则缩小步长重新计算。
时间积分器的选择直接影响长时间模拟的精度和稳定性。DOP853 虽然不是辛积分器，但八阶精度使其在本文考虑的积分时间内能保持良好的守恒性。
Numba JIT 编译与并行化
Numba 是 Python 的即时编译器，能将标注了 @njit 装饰器的数值函数通过 LLVM 编译为优化的机器码。设置 parallel=True 后，使用 prange 替代标准 range 可自动将循环分配到多个 CPU 核心（类似 OpenMP 线程并行）。cache=True 选项将编译后的机器码缓存到磁盘，避免重复编译开销。对于 $N \geq 512$ 的网格规模，典型加速比可达 5x 到 50x 相比纯 NumPy 实现。
Numba 加速是本文在普通笔记本上实现研究级计算效率的关键，理解其工作原理才能评估求解器的实际性能和可复现性。
递归量化分析（Recurrence Quantification Analysis, RQA）
RQA 是一种非线性时间序列分析方法，用于刻画动力系统的相空间结构。基本步骤：(1) 将系统状态投影到低维空间（如归一化守恒量 $(M/M_0, P/P_0, E/E_0)$）；(2) 计算状态对间的欧氏距离矩阵 $D_{ij}$；(3) 设定阈值 $\varepsilon$，构建二元递归矩阵 $R_{ij} = \Theta(\varepsilon - D_{ij})$；(4) 从中提取标量指标——递归率 RR（递归点对占总对数的比例）和确定性 DET（形成对角线结构的递归点比例）。DET 接近 1 表示确定性动力学（规则运动在不变流形上），DET 接近 0 表示混沌动力学。
RQA 是本文独特的诊断工具之一，用于独立验证数值解是否保持了 KdV 方程的可积结构。它提供了一种不依赖解析解的模型无关评估方法。
研究动机
Korteweg-de Vries 方程在地球物理流体力学中具有基础性地位，是描述浅水和层化海洋环境中长内孤立波和表面重力波传播的原型模型。尽管 KdV 方程的解析理论（逆散射变换）已经非常成熟，但在实际研究中，数值模拟仍然不可或缺：需要探索解析方法难以覆盖的参数空间、验证近似初始条件的有效性、以及刻画多孤立子相互作用的瞬态动力学。然而，现有的数值工具存在几个具体问题：首先，许多求解器缺乏系统性的守恒律诊断，无法定量评估长时间模拟中能量、动量等物理量的数值耗散；其次，大多数验证仅依赖目视检查波形是否保持形状，缺少信息论和相空间分析等更严格的定量指标；第三，现有的 Python 科学计算生态虽然强大（NumPy、SciPy），但直接用 Python 实现的 KdV 求解器往往计算效率不足，难以在普通硬件上完成研究级模拟；最后，输出格式不统一，难以与海洋学标准分析工具（如 NetCDF/CF 惯例）无缝对接。
本文的目标是本文的直接目标是开发一个名为 sangkuriang 的开源 Python 库，提供一个兼具研究级精度和计算效率的 KdV 方程伪谱求解器。具体可量化的目标包括：(1) 通过四个从简到繁的理想化海洋波浪场景（孤立子传播、对称相互作用、超越碰撞、三体相互作用）系统验证求解器；(2) 确保所有守恒律（质量、动量、能量）的相对误差低于 $O(10^{-4})$；(3) 验证测量的孤立子速度与理论振幅-速度关系 $v = \varepsilon A/3$ 的偏差在 5% 以内；(4) 使用信息论诊断（谱熵、LMC 复杂度、Fisher 信息）和递归量化分析证明数值解保持了 KdV 方程的可积 Hamilton 结构；(5) 实现在普通笔记本工作站上数分钟内完成所有测试用例的计算效率。
与已有工作不同的是，本文的独特切入点在于将「诊断框架」提升到与「求解器本身」同等重要的地位。大多数 KdV 数值研究止步于波形可视化和守恒量时间序列的检查，而本文引入了四个层面的互补诊断：守恒律误差分析（含漂移率检测）、孤立子轨迹追踪与速度验证、谱信息论特征化（谱熵、LMC 统计复杂度、Fisher 信息）、以及相空间递归量化分析。这些诊断工具是模型无关的，不仅适用于理想化 KdV 系统，也为未来研究更复杂的变系数 KdV 扩展方程（用于真实海洋层化和地形变化的建模）提供了评估框架。此外，本文将所有输出标准化为 NetCDF/CF 惯例格式，降低了与现有海洋学工作流集成的门槛。
核心方法
sangkuriang 的方法架构可以用一个烹饪的比喻来理解：将 KdV 方程的求解过程想象为「烹饪」一道非线性波动菜肴。食材是初始波形（$\text{sech}^2$ 孤立子轮廓），烹饪分两步：第一步「切菜」——用傅里叶伪谱方法将空间导数的计算从物理域搬到频率域（在那里微分变成简单乘法），就像用料理机高效地处理食材；第二步「控火」——用 DOP853 自适应时间积分器推进时间演化，八阶精度保证每一步都足够精确，自适应步长机制像一个智能温控系统，在动力学变化剧烈时自动加密时间步。整个烹饪过程用 Numba JIT 编译加速，相当于给厨房装了涡轮增压。最后，成品（解的时间序列）打包为 NetCDF 格式，附带完整的诊断报告（守恒量、轨迹、谱特征、相空间结构）。
本文的核心创新不在于发明新的数值方法——伪谱空间离散化和高阶 Runge-Kutta 时间积分都是成熟技术——而在于系统性地将「诊断」提升为核心产出。具体来说，sangkuriang 的独特之处是四层验证框架：第一层是传统的守恒律误差监测，但不仅看最大误差，还用线性回归检测是否存在系统性漂移（漂移率 $\gamma_Q$）；第二层是孤立子轨迹追踪，通过峰值检测和振幅排序自动识别和跟踪每个孤立子，验证振幅-速度关系；第三层是信息论分析——谱熵 $S_k$ 衡量谱能量分布的宽度，LMC 复杂度 $C_{\text{LMC}} = H \times D$ 区分结构化非平衡态和无特征随机态，Fisher 信息 $F$ 检测空间局域化程度；第四层是递归量化分析，通过递归率 RR 和确定性 DET 判断相空间轨迹是否保持规则（非混沌）结构。这种多层次、多视角的诊断框架使得验证不限于「看起来对不对」，而是能定量回答「数值解在多大程度上保持了原方程的几何结构」。
方法步骤详情
求解器的完整算法流程如下：(1) 初始化阶段：给定初始条件 $u_0(x)$，在 $N$ 个等距网格点 $x_j = x_{\min} + j\Delta x$ 上采样得到初始状态向量 $\mathbf{u}(0) \in \mathbb{R}^N$；预计算波数向量 $k_n$ 并存储；通过 Numba 的 set_num_threads(n) 初始化并行线程。(2) 右端算子计算（每个时间步重复）：对当前状态 $\mathbf{u}(t)$ 做 FFT 得到 Fourier 系数 $\hat{u}_k$；乘以 $ik_n$ 和 $(ik_n)^3$ 分别得到一阶和三阶导数的 Fourier 系数；通过 IFFT 变回物理空间得到 $\partial_x u$ 和 $\partial_x^3 u$；用 JIT 编译的并行函数组装右端项 $F_j(u) = -\varepsilon u_j (\partial_x u)_j - \mu (\partial_x^3 u)_j$。(3) 时间推进：DOP853 积分器调用右端算子，每步 12 次函数求值，比较八阶和五阶解估计局部误差；若满足容差条件则接受并调整步长，否则缩小步长重算。(4) 诊断输出：在指定输出时刻存储解状态，用梯形法则计算守恒量 $M(t)$、$P(t)$、$E(t)$ 并记录相对误差；输出 NetCDF 格式数据，生成动画可视化和诊断报告。对于四个测试用例，参数设置分别为：Case 1-3 使用 $N=512$，$\text{rtol}=10^{-10}$，$\text{atol}=10^{-12}$；Case 4 使用 $N=1024$，$\text{rtol}=10^{-11}$，$\text{atol}=10^{-13}$。
技术新颖性
从技术新颖性角度分析，sangkuriang 的贡献不在于单一技术创新，而在于工程整合和诊断深度两个维度的组合创新。在工程整合方面，它将伪谱方法的指数收敛精度、DOP853 的八阶自适应精度、Numba 的 JIT 并行加速、以及 NetCDF/CF 的标准化输出整合为一个开箱即用的 Python 库，通过 PyPI 一键安装——这在已有工具中是少见的。大多数 KdV 求解器要么是 Fortran/MATLAB 的研究代码（不便于 Python 生态集成），要么是纯 NumPy 实现（计算效率不足）。在诊断深度方面，本文引入的谱熵、LMC 复杂度、Fisher 信息和 RQA 四种诊断工具的组合使用在 KdV 数值研究中是新颖的——特别是 RQA 的引入，它提供了一种不依赖解析解、不依赖特定方程结构的模型无关验证方法。此外，本文强调的「漂移率检测」（通过线性回归分离系统性漂移和有界振荡）也是一个实用的诊断创新，对判断长时间模拟的可靠性尤为关键。
实验结果
sangkuriang 求解器在四个逐步复杂的测试用例中表现出色，全面验证了其数值精度和物理保真度。在守恒律方面，所有测试用例的守恒误差均远低于实际模拟可接受的 $10^{-2}$ 阈值：单孤立子（Case 1）的质量、动量、能量最大相对误差分别为 $7.64 \times 10^{-5}$、$2.98 \times 10^{-7}$、$1.73 \times 10^{-7}$；等幅双孤立子（Case 2）为 $1.57 \times 10^{-4}$、$1.90 \times 10^{-6}$、$1.45 \times 10^{-6}$；超越碰撞（Case 3）为 $3.33 \times 10^{-4}$、$6.49 \times 10^{-6}$、$2.60 \times 10^{-6}$；三孤立子（Case 4）为 $2.56 \times 10^{-4}$、$8.86 \times 10^{-6}$、$6.20 \times 10^{-6}$。一个值得注意的层次结构是：质量误差比动量误差大约两个数量级，动量误差又比能量误差大——这反映了不同守恒量对解的高频分量的敏感性差异，能量泛函中的空间导数项更侧重于光滑的孤立子核心。线性回归分析表明漂移率在 $10^{-9}\,\text{s}^{-1}$ 量级，说明 DOP853 积分器引入的主要是有界振荡误差而非系统性偏差。在孤立子速度验证方面，所有追踪轨迹的线性拟合 $R^2 > 0.99$，测量速度与理论关系 $v = \varepsilon A/3$ 的偏差在 5% 以内。例如，Case 1 的测量速度为 0.349 m/s，理论预测为 0.337 m/s（偏差 3.6%）；Case 3 的快孤立子测量速度 0.508 m/s，慢孤立子 0.185 m/s。信息论诊断方面，谱熵 $S_k$ 在 0.46-0.59 范围内保持有界，碰撞事件期间出现瞬态增加后恢复；LMC 复杂度在 0.022-0.034 范围内；Fisher 信息在碰撞期间出现显著尖峰（Case 3 达到 1.8 m$^{-2}$），反映瞬态梯度变陡。RQA 分析显示确定性 DET > 0.99，证实相空间轨迹保持规则（非混沌）结构。计算效率方面，在配备 Intel Core i7-8550U 的普通笔记本上，Case 1-3 的吞吐量约 1300-1370 步/秒，Case 4（$N=1024$）约 1030 步/秒，最耗时的 Case 4 在 534 秒内完成。
查看结构化数据
| 任务 | 指标 | 本文 | 基线 | 提升 |
|---|---|---|---|---|
| 单孤立子传播 (Case 1) | 质量守恒误差 | 7.64×10⁻⁵ | 实用阈值 10⁻² | 低于阈值约 2.5 个数量级 |
| 单孤立子传播 (Case 1) | 动量守恒误差 | 2.98×10⁻⁷ | 实用阈值 10⁻² | 低于阈值约 4.5 个数量级 |
| 单孤立子传播 (Case 1) | 能量守恒误差 | 1.73×10⁻⁷ | 实用阈值 10⁻² | 低于阈值约 4.8 个数量级 |
| 超越碰撞 (Case 3) | 质量守恒误差 | 3.33×10⁻⁴ | 实用阈值 10⁻² | 低于阈值约 1.5 个数量级 |
| 三孤立子 (Case 4) | 能量守恒误差 | 6.20×10⁻⁶ | 实用阈值 10⁻² | 低于阈值约 3.2 个数量级 |
| 单孤立子速度验证 | 速度偏差 | 3.6% | 理论 v=εA/3 | 在 5% 以内 |
| 轨迹线性度 | R² | >0.99 | 完美线性 R²=1 | 极接近理论预测 |
| 递归量化分析 | 确定性 DET | 0.999 | 完全可积理想值 1 | 极接近理想可积结构 |
局限与改进
作者明确承认了四个主要局限：(1) 测试用例使用近似 $\text{sech}^2$ 初始条件而非逆散射理论给出的精确 $N$-孤立子解，这意味着初始阶段存在调整瞬态，严格来说不如直接用精确解验证来得严格；(2) 周期性边界条件排除了孤立子与边界、非周期波背景的相互作用，而这些在真实海岸环境中很重要（波浪遇到海岸线、岬角、河流羽流）；(3) 标准 KdV 方程仅捕捉了非线性与色散的领头阶平衡，扩展到修正 KdV、Kawahara 方程等高阶模型以及考虑缓慢变化地形和层化的变系数方程才能覆盖更广泛的海洋环境；(4) 仅考虑了表面波类比，扩展到连续层化流体中的内波动力学需要引入模态结构和密度依赖系数。从独立观察来看，还有几个值得补充的局限：DOP853 不是辛积分器，虽然在本文考虑的积分时间内误差主要表现为有界振荡，但对需要数小时甚至数天的跨大陆架孤立子传播模拟，误差累积可能成为问题；$N=512/1024$ 的网格分辨率对于真实海洋尺度的模拟（千米级域、米级孤立子）是否足够尚未验证；此外，所有测试都在一维空间中进行，而真实海洋内波具有显著的二维/三维特征（如 Kadomtsev-Petviashvili 方程描述的斜相互作用）。
独立分析的弱点
从独立分析的角度，sangkuriang 存在以下几个可改进的弱点。时间积分器的结构保持性：DOP853 虽然精度高（八阶），但不是辛积分器，不保证相空间体积守恒。在本文的短时间积分中（50-80 秒），这不构成问题，但对于海洋学家关心的跨大陆架孤立子传播模拟（可能需要积分数小时），辛积分器（如隐式中点法或分裂步方法）可能更合适。改进方向：实现辛积分器作为可选后端，或至少提供能量漂移的长期监控指标。非线性项的混叠控制：非线性项 $-u \partial_x u$ 在 Fourier 空间中引入二次耦合，可能产生混叠误差。论文提到在 $N=512/1024$ 下这些效应相对于指定容差可以忽略，但没有实施标准的 2/3 去混叠规则。改进方向：添加可选的去混叠策略（如 Orszag 的 2/3 规则），这将允许在更低的网格分辨率下获得同等精度。二维扩展的缺失：当前仅支持一维 KdV 方程，而真实海洋内波具有显著的横向结构。改进方向：将伪谱框架扩展到 KP（Kadomtsev-Petviashvili）方程，支持二维孤立子动力学。初始条件的局限：仅提供 $\text{sech}^2$、tanh、Gaussian 等预设轮廓类，缺少逆散射变换的自动化接口。改进方向：集成标准的逆散射数值实现（如 Marchenko 方程的离散化），使用户能从任意初始条件自动生成精确的 $N$-孤立子解。
未来方向
作者在论文中提出了几个明确的未来方向：(1) 集成逆散射变换例程以生成精确初始条件；(2) 开发变系数 KdV 扩展以模拟波浪在变化地形上的传播——这是从理想化走向真实海岸海洋应用的关键一步；(3) 通过 Ostrovsky 方程引入旋转效应，适用于地球自转尺度下重要的内波；(4) 扩展到二维 KP 方程以研究斜相互作用和横向不稳定性。基于本文的成果还可以延伸出更多方向：将诊断框架（特别是 RQA 和信息论指标）应用于已有的非可积 KdV 变体（如阻尼 KdV、随机 KdV），定量评估扰动对可积结构的破坏程度；开发 GPU 加速版本（通过 CuPy 替代 NumPy 的 FFT），支持大规模参数扫描和集合模拟；将 sangkuriang 作为教育工具，为地球物理流体力学课程提供交互式孤立子实验平台；以及将求解器与观测数据同化框架耦合，用于从卫星 SAR 图像反演海洋内波参数。
复现评估
sangkuriang 在可复现性方面做得非常好：源代码在 GitHub 以 MIT 许可证开源，可通过 pip install sangkuriang-ideal-solver 从 PyPI 一键安装；数据分析脚本单独开源；所有模拟输出（动画、NetCDF 数据、图形、计算日志、统计分析）存档在 Open Science Framework。算力需求极低——所有测试在配备 Intel Core i7-8550U（8 核，4.0 GHz）的普通笔记本上完成，最耗时的 Case 4 仅需 534 秒（约 9 分钟）。依赖链简洁（NumPy、SciPy、Numba、Matplotlib、netCDF4），均为 Python 科学计算的标准库。复现难度很低：安装依赖后，直接运行 sangkuriang 命令即可执行预定义测试用例，或通过 Python API 编程调用。唯一的门槛是需要理解 KdV 方程的物理背景来解读输出结果。
论文图表
汇总了四个测试用例的网格点数 N、空间分辨率、最终时间、自适应时间步数和 CPU 时间。Case 4（三孤立子，N=1024）最耗时，534 秒完成，吞吐量约 1030 步/秒。
这张表直接回答「在普通硬件上完成这些模拟需要多长时间」，是评估求解器实际可用性的关键数据。