求解薛定谔-亥姆霍兹方程的两种高效、高阶精度、质量与能量守恒的数值方法

《Communications in Nonlinear Science and Numerical Simulation》：Two efficient, high-order accurate, mass- and energy-preserving numerical methods for solving the Schr?dinger–Helmholtz Equations

　　

# 摘要

本文提出了一种线性标量辅助变量（LSAV）方法，其中辅助变量直接定义为能量泛函中的非线性势能部分，无需任何经验参数。应用LSAV方法，我们将Schr?dinger–Helmholtz方程组改写为等价的方程组系统，并给出一个近似系统，其中能量被修改为线性-二次结构。采用Gauss型Runge–Kutta方法在时间方向离散、伪谱方法在空间方向离散该近似系统，我们推导出了一种完全离散的数值方法，该方法保持质量和能量守恒，易于实现，并且在空间和时间上均具有高阶精度。为进一步改善所提方案的计算效率和精度，引入预测-校正技术，开发了一种高效且保守的替代方案，该方案对光滑解可实现任意高阶时间精度。广泛的数值实验表明，所提方案具有高阶精度、保守性质和计算效率。

# 引言

本文考虑如下d维（d=1,2,3）Schr?dinger–Helmholtz（SH）方程组：

$$\begin{cases} i u_t + \Delta u + \phi u = 0, & \mathbf{x} \in \Omega, \, t > 0, \\ \alpha \phi - \beta^2 \Delta \phi = |u|^2, & \mathbf{x} \in \Omega, \, t > 0, \\ u(\mathbf{x}, 0) = u_0(\mathbf{x}), & \mathbf{x} \in \overline{\Omega}, \end{cases}$$

其中 $i$ 为虚数单位，$\mathbf{x}$ 表示空间坐标，$t > 0$ 为时间。该系统在有限区域 $\Omega$ 上考虑周期性边界条件或齐次Dirichlet/Neumann边界条件。未知量 $u(\mathbf{x}, t)$ 是具有给定初始状态 $u_0(\mathbf{x})$ 的复值波函数，而 $\phi(\mathbf{x}, t)$ 作为实值势。非负常数 $\alpha$ 和 $\beta$ 满足 $\alpha + \beta \neq 0$。作为一个重要的耦合非线性模型，系统(1.1)被广泛应用于非线性光学、等离子体物理和量子力学[1]。它是一个广义框架，当 $\alpha = 0$ 时可退化为Schr?dinger–Poisson（SP）方程[2],[3]，当 $\beta = 0$ 时退化为标准三次非线性Schr?dinger（NLS）方程[4]。在整个工作中，我们关注参数 $\alpha > 0$，$\beta > 0$ 的SH系统。值得注意的是，本文开发的数值方案可直接应用于求解SP方程和NLS方程而无需修改。

如文献[5]所建立的那样，SH系统(1.1)保持总质量

$$\mathcal{M}(t) = \int_\Omega |u(\mathbf{x}, t)|^2 d\mathbf{x} \equiv \mathcal{M}(0), \quad t > 0,$$

和能量

$$\mathcal{E}(t) = \int_\Omega \left(|\nabla u(\mathbf{x}, t)|^2 - \frac{1}{2}\phi(\mathbf{x}, t)|u(\mathbf{x}, t)|^2\right) d\mathbf{x} \equiv \mathcal{E}(0), \quad t > 0.$$

对SH方程已进行了广泛的数学和数值研究。从数学角度，Cao等人[1]将SH系统作为广义NLS方程的数值正则化引入，并建立了其全局适定性。从数值观点来看，已开发了各种数值方案来求解SH系统，主要关注其时空解耦。例如，Yang等人[6]为线性化的一阶精度后向Euler Galerkin有限元方法（FEM）建立了无条件误差估计。为实现更高精度，提出了若干具有二阶时间精度的数值方案。Wang[7]发展了线性化Crank–Nicolson型Galerkin FEM，而Shi等人[5]分析了线性化且二阶精度的BDF2（两步后向差分公式）Galerkin FEM。最近，Guo等人[8]提出了二阶精度且保持质量和能量守恒的虚拟元方法（VEM），Qin等人[9]提出了线性隐式蛙跳方案，同样保持了二阶时间精度。尽管这些现有数值方法实现了高空间精度，但值得注意的是，现有针对SH系统的守恒和/或线性隐式方案主要局限于一阶或二阶时间精度。这种低阶时间精度通常需要极小的时间步长来维持长期模拟中所需的分辨率，从而大大降低计算效率。据我们所知，关于同时保持SH系统质量和能量性质的高阶时间方法的研究仍十分有限。

为实现高阶时间精度同时在离散意义上保持总质量和能量，已制定了各种高阶能量保持策略，包括Hamiltonian边值方法[10]、时间有限元方法[11]、连续阶段Runge–Kutta（RK）方法[12],[13]和平均向量场方法[14],[15]。特别是，在过去十年中，不变能量正交化（IEQ）[16],[17]和标量辅助变量（SAV）[18],[19],[20]方法取得了巨大成功。通过引入辅助变量，这些方法将非线性能量转化为二次形式，从而能够高效地半隐式处理非线性项[21]。众所周知，标准RK方法保持线性不变量，而辛RK方法保持二次不变量[22]。基于能量正交化概念，研究人员将IEQ或SAV框架与RK方法结合，构建了高阶能量保持方案[23],[24],[25]。虽然这些方案以最优阶数实现了高精度，但它们是完全非线性的，在每个时间步都需要耗时的计算迭代。为改善计算效率并避免非线性迭代，最近已提出并分析了若干基于SAV和IEQ方法的线性化高阶方法用于梯度流模型[20],[26],[27],[28],[29]。通过采用显式外推技术[26]或后向差分公式（BDF）[20]，这些线性化SAV/IEQ方案显著降低了计算成本。随后，此类线性化高效方法被扩展到守恒系统，发展了特定的高阶质量保持方案[23]和高阶能量保持方案[24]。尽管具有计算效率，但当应用于如SH系统这样的复杂守恒模型时，这些现有线性化高阶方法往往无法同时保持离散质量和能量。此外，由于缺乏严格的结构保持性，它们在长期模拟中经常遭受辅助变量漂移问题。

为构建高阶且高效的结构保持方案，需要精心处理非线性项。标准SAV方法因其线性且半隐式的特性已取得巨大成功。然而，其平方根变换严格要求非线性能量有下界，并且需要一个人为稳定常数 $c_0$。由于 $c_0$ 的选择往往是经验的，过大的值会严重降低数值精度和计算稳定性。为解决这些问题，已提出了各种高级变体。例如，指数SAV（ESAV）方法[30]利用指数映射消除能量的下界限制并避免 $c_0$ 的选择。此外，广义SAV（GSAV）方法[31]引入任意可逆函数，扩展了仅总自由能有界的复杂系统的适用性。然而，现有高级变体专注于为耗散流量身定制的非线性映射。当应用于守恒系统时，这些变换破坏了能量的线性-二次结构，严重降低了辛Runge–Kutta（SRK）方法的结构保持性质。因此，所得的完全离散方案需要计算昂贵的非线性迭代求解器，意味着守恒律的精确保持严重依赖于迭代容差足够小。因此，开发一种新方法以保持线性-二次能量结构而无需此类昂贵迭代是非常令人期待的。

受此启发，我们提出了一种简化的变换方法，称为线性标量辅助变量（LSAV）方法。与使用复杂的非线性变换不同，我们直接将辅助变量定义为关于原始非线性势能泛函本身的线性函数。因此，与引入复杂指数和对数计算的Exponential SAV（ESAV）方法[30]相比，所提出的LSAV方法具有更简单的结构。这种直接代入将能量泛函改写为线性-二次形式，使SRK方法能够在每个时间步无需非线性迭代的情况下精确保持离散意义上的总质量和能量。此外，与容易受到零分母奇点和符号反转问题困扰的Generalized SAV（GSAV）方法[31]不同，我们的LSAV方法采用了一种新颖的二次有理乘数，其分母严格被1从下方有界。这种公式本质上避免了奇点，保证了非线性项的符号一致性，并有效地稳定了时间步进，防止长期模拟中的数值爆炸。

本文的主要贡献总结如下：

- 通过引入新颖的LSAV方法，其中辅助变量被定义为与Schr?dinger方程非线性项相关的能量泛函部分的线性函数，我们将SH方程组改写为能量泛函呈现线性-二次结构的等价方程组系统。该方法完全避免了标准SAV方法所需的稳定常数 $c_0$ 的人工选择。

- 所提出的两种数值方案可实现任意高阶时间精度和空间谱精度（前提是解足够光滑），并在离散意义上精确保持总质量和能量。得益于高阶精度特性，数值模拟中可以采用相对较粗的网格而不影响目标数值精度。

- 两种提出的线性数值方案都易于实现，并可通过快速傅里叶或正弦变换高效计算。这是在LSAV框架下，通过在时间方向耦合时间外推技术与SRK方法、在空间方向使用伪谱方法来实现的。

总之，所提出的两种线性方案可在时间和空间上实现任意高阶精度，在离散意义上精确保持总质量和能量，并保持高计算效率。因此，它们可作为SH系统长期模拟的高度可靠的数值工具。

本文其余部分组织如下。第2节，我们引入LSAV技术，推导SH方程的等价模型改写，该模型继承了质量和能量守恒性质。第3节，我们采用SRK方法，提出两种在离散意义上同时保持质量和能量的半离散高阶线性方案。第4节，采用傅里叶或正弦伪谱方法进行空间离散化，得到两种完全离散方案。第5节，进行了广泛的数值实验以验证收敛性、精度、守恒律和计算效率。最后，第6节得出简洁的结论。

# 使用LSAV方法的模型改写

在本节中，我们利用LSAV方法将SH方程组改写为等价方程组系统。通过引入修正的标量辅助变量，能量泛函被改写为线性-二次形式，且改写后的系统严格保持质量和能量守恒律。

为方便后续推导，我们首先引入一些标准数学记号。设 $(\cdot, \cdot)$ 表示区域 $\Omega$ 上的复 $L^2$ 内积，定义为...

# 线性化、高阶精度且保守的半离散方案

在本节中，我们利用外推技术为改写模型(2.11)–(2.13)构造线性化LSAV系统，该系统严格继承了质量和能量守恒律。为实现时间上的高阶精度，引入Gauss型Runge–Kutta方法离散该线性系统，得到线性半离散LSAV-SRK方案。此外，我们应用预测-校正策略构造另一个线性化半离散LSAV-SRK-PC方案，该方案仍然...

# 空间离散化

SH方程通常施加周期性或齐次Dirichlet边界条件。对于后者，我们采用正弦伪谱方法进行空间离散化。该离散方法的综合讨论详见[24]，我们将在第5节中提供相应的数值验证。本节关注周期性边界值问题，对此采用傅里叶伪谱方法。

为简便起见，我们...

# 数值结果

在本节中，我们展示一系列数值实验以验证所提LSAV-SRK和LSAV-SRK-PC方案的时间收敛阶、空间谱精度和守恒性质。此外，我们评估LSAV-SRK-PC方案的计算效率并研究耦合SH方程的复杂动态行为。在所有数值测试中，我们假设精确解足够光滑，使得空间谱精度可以...

# 结论

本文开发了两种高效、高阶精度且结构保持的Schr?dinger–Helmholtz方程数值方案。通过引入新颖的LSAV方法，我们直接将辅助变量定义为能量泛函中对应于Schr?dinger方程非线性项的部分。这明确避免了标准SAV方法所需的经验稳定常数 $c_0$。结合外推策略、SRK方法、伪谱方法和... 未引用参考文献[34]。

CRediT作者贡献声明

Chuchu Sun：撰写 – 初稿、可视化、调查、概念化。Jialing Wang：撰写 – 审稿与编辑、形式分析。

王庭春：论文写作——审阅与修订、方法设计、调查研究、形式化分析、概念构思。利益冲突声明：作者声明，他们不存在已知的竞争经济利益或人际关系，这些因素可能会影响本文所报告的工作。孙初初 | 王嘉凌 | 王庭春 南京信息工程大学数学与统计学院，中国南京，210044。