适用于二维无界域中含波动算子的非线性薛定谔方程的自适应时间步进格式：收敛性及守恒性研究

《Mathematics and Computers in Simulation》：Adaptive time-stepping scheme for the nonlinear Schr?dinger equation with wave operator in a 2D unbounded domain: Convergence and conservation

　　姜慧玲|胡东东
江西理工大学理学院，南昌330022，中华人民共和国

摘要
本文提出了一种新型的线性隐式、能量守恒的时间步进方案，用于二维无界域中的非线性薛定谔方程及其波算子问题。为此，我们首先将原方程重新表述为一个等价的修改系统，并结合Crank–Nicolson方法和线性插值，为修改系统提出了一种二阶变步长时间步进方案。经证明，所提出的方案具有能量守恒性、唯一解性和收敛性。同时，我们选择了映射切比雪夫谱Galerkin方法来处理无界域问题，并详细提供了完全离散化方案的实施细节。大量的数值比较表明，理论结果正确，所提出的方案在长时间计算中具有高效性。

引言
非线性薛定谔（NLS）方程是量子力学的基本方程，在物理、化学和工程等多个领域都有应用[1]。有趣的是，带有波算子的NLS方程是NLS方程的一种推广形式，它来源于不同的物理应用，例如Klein–Gordon方程的非相对论极限[2]和等离子体中的Langmuir波包络近似[3]。在本文中，我们考虑了二维（2D）无界域中的带有波算子的NLS方程：
i?tu + ?ttu ? κ?2u + β|u|2 = 0, (x,y,t)∈R2×(0,T), u(x,y,0) = u?(x,y), ut(x,y,0) = v?(x,y), (x,y)∈R2, u(x,y,t)→0 当 |x|→∞ 或 |y|→∞，
其中 R=(?∞,∞)，?=(?x,?y)? 表示梯度算子，i 是虚数单位，κ 是色散系数，β 是非线性强度，u 是一个足够光滑的复值波函数。那么，NLS方程满足以下质量、能量和动量守恒定律：
dM/dt = 0，其中 M(t)=∫R2[|u|2 + 2Im(utu?)]dxdy = 0，
dE/dt = 0，其中 E(t)=∫R2[|ut|2 + κ|?u|2 + β2|u|?]dxdy = 0，
dM?/dt = 0，其中 M?(t)=∫R2[1/2Im(u?u?) ? Re(ut?u?)]dxdy = 0，
这里 u? 是 u 的复共轭，Im(•) 和 Re(•) 分别表示 ‘•’ 的虚部和实部。

如何为非线性守恒系统构建尽可能多不变量的数值方案是一项具有挑战性的任务。在过去的几十年中，学者们系统地开发了结构保持算法[4][5]。同时，有效维护原始问题的物理不变量也成为评估数值算法优劣的重要指标之一。在这些年里，逐渐提出了一些高效的守恒系统算法，例如辛Runge–Kutta（sRK）方法[6]、哈密顿边界值（HBV）方法[7]、平均矢量场方法[8]、离散梯度方法[9]、指数积分器方法[10]等。值得注意的是，上述方法大多是完全隐式的，需要在每个时间步进行非线性迭代。

最近，李等人[11][12]提出了一类保持系统原始能量的显式方法。其核心思想基于松弛技术和指数Rosenbrock方法，显著提高了能量守恒算法的计算效率。特别是，研究人员提出了不变能量二次化方法[13][14]和标量辅助变量（SAV）方法[15]来构建能量稳定的方案。它们的主要思想是通过引入辅助变量将原始能量转化为等价的二次型。相关方法包括指数SAV方法[16]和广义SAV（GSAV）方法[17]。上述新型方案仅在每个时间步解决一类线性代数系统。这激发了我们为二维无界域中的NLSW方程设计线性隐式能量守恒时间步进方案的动机。

据我们所知，无界域中拉普拉斯算子的数值方法主要可以分为以下三类：(i) 人工边界条件方法[18][19][20][21]，它引入人工边界将无界域截断为有界计算域，并构建接近原始问题的边界条件；(ii) 完美匹配层方法[22][23]，它在计算域外添加吸收层通过坐标拉伸来抑制反射；(iii) 无限元素和谱方法[24][25][26][27][28][29][30][31]，它们使用特殊基函数（如Laguerre函数、Hermite函数、映射Jacobi函数）来处理无界域。

据我们所知，对于有界域中的NLSW方程，存在一些使用均匀时间步进的能量守恒时间步进方案，例如HBV方法[32]、Crank-Nicolson方法[33][34][35]、sRK方法[36]、辅助变量方法[37][38]等。最近，郭等人[40]提出了一种适用于多维无界域中NLSW方程的质量、能量和动量守恒方案，并采用自适应时间步进策略。然而，由于方案的复杂性，无法进行收敛性分析。受参考文献[40]的启发，我们为二维无界域中的NLSW方程设计了一种高效的自适应时间步进能量守恒方案，并对所提出的方案进行了严格的收敛性分析。这项工作的主要贡献可以总结为：
- 通过结合变步长Crank–Nicolson方法和GSAV方法，提出了一种线性隐式自适应时间步进能量守恒方案。所提出的时间步进方案被证明对于二维无界域中的NLSW方程是唯一可解且收敛的。
- 我们为二维无界域中的拉普拉斯算子提出了一种快速傅里叶类切比雪夫谱Galerkin方法，并构建了所提出时间步进方案的完全离散化实现。
- 我们通过使用均匀和时间自适应步进策略解决了所提出的方案，结果表明自适应时间步进策略在计算效率上优于均匀时间步进策略。

本文的其余部分安排如下：第2节中，我们为二维无界域中的NLSW方程提出了一种自适应时间步进能量守恒方案。第3节中，我们证明了所提出时间步进方案的收敛性。第4节中，我们应用映射切比雪夫谱Galerkin作为空间离散化，并详细提供了完全离散化方案的实现。第5节中，我们报告了一些数值结果以验证理论结果。最后，在第6节中，我们得出了一些结论。

模型重构
遵循GSAV方法的思想，我们引入了一个辅助变量 s(t)=∫R2β2|u|?dxdy + C?，其中C?是给定的非负常数。对s(t)求时间导数，我们得到 st=∫R2β|u|2(uut? + u?ut)dxdy = 2Re∫R2β|u|2uut?dxdy = 2Re∫R2Υ(u,s)ut?dxdy，其中Υ(u,s)=β|u|2us∫R2β2|u|?dxdy + C?。因此，通过引入辅助变量 v=ut，我们将原始哈密顿能量重新表述为以下修改形式：
E[u,v,s](t)=∫R2[|v|2 + κ|?u|2]dxdy + s(t) ? C?。

辅助引理
在证明收敛性之前，我们介绍一些有用的引理。
引理3.1 [41, Proposition 1.3] 假设 u∈H2(R2)，那么有 ‖u‖∞≤C‖u‖2。
引理3.2 对于时间序列 {η?,η?,…,η?} 和 {r?,r?,…,r?}，有 |2∑?=??tkrkδtη??1/?| ≤ ∑?=??tk|η??1|2 + ∑?=??tk|δtrk??/?|2 + 1/2|η?|2 + 2|r?|2 + |η?|2 + |r?|2。

证明
通过一些简单的计算，我们有 2∑?=??tkrkδtη??1/? = 2∑?=??rkη??? + ∑?=??rkη??? = 2∑?=????rkη??? + 2rnη???r?η? = ?2∑?=??tkδtrk??/?η??? + 2rnη???r?η? ≤ ∑?=??tk|η??1|2 + ∑?=??tk|δtrk??/?|2 + 1/4??|η?|2 + ??|2rn|2 + 1/4??|η?|2 + ??|2r?|2。证明完毕。

映射切比雪夫谱Galerkin方法
为了处理无界域，我们通过代数映射函数 ξ=χ?+χ?, χ=ξ??ξ?, χ∈(?1,1), χ∈R 来提出谱Galerkin方法。对于 χ∈R 和整数 n≥0，我们定义映射切比雪夫函数 [25] 为 T?(χ)=1/cnπ21?ξ2, T?(ξ)=1/cnπ21+χ2, T?(ξ)，其中 c?=2 且 cn=1（对于 n≥1），T? 表示n阶切比雪夫函数[24]，它由 T?(ξ)=cos(narccos(ξ)) 确定。容易验证 ∫R T?(χ)Tj(χ)dχ=2ckcjπ∫??1T?(ξ)Tj(ξ)w(ξ)dξ=δkj，
∫R T?′(χ)Tj′(χ)dχ=1ck((4ck?1?ck?2)(k?1)2/16 + (4ck+1?ck+...

数值实验
在本节中，我们通过数值求解NLSW方程来报告方案I和方案II的一些数值结果。所有模拟都是使用Matlab R2018a软件在配备Intel Core i7和16 GB RAM的计算机上进行的。为了说明所提方案的数值精度和离散能量守恒律，我们始终取 rmax=2, C?=2，并通过 |EN??EN?|/EN? 计算相对修改能量误差。L?范数和最大范数误差是通过选择 M=2N 来计算的。

结论
在本文中，我们提出了一种适用于二维无界域中NLSW方程的自适应时间步进、线性隐式和能量守恒方案。所提出的方案被证明是唯一可解的，并且在时间上具有二阶精度。大量的数值比较证实了理论分析的正确性，并表明采用自适应时间步进策略的方案在计算效率上优于均匀时间步进策略。

作者贡献声明
姜慧玲：撰写 – 原始草案、方法论、研究。
胡东东：撰写 – 审稿与编辑、监督、资金获取。

伦理批准
本文不包含任何由作者进行的涉及人类参与者或动物的研究。

利益冲突声明
作者声明他们没有已知的竞争性财务利益或个人关系可能影响本文报道的工作。