摘要:在提出一种基于时间分裂格式的纯无网格有限点集(split-step finite pointset method, SS-FPM)法的基础上, 数值模拟了含孤立波的二维非线性薛定谔 (nonlinear Schr?dinger, NLS) / (Gross-Pitaevskii, GP) 方程. SS-FPM的构造过程为: 1) 基于时间分裂的思想将非线性薛定谔方程分成线性导数项和非线性项; 2) 采用基于Taylor展开和加权最小二乘法的有限点集法, 借助Wendland权函数, 对线性导数项进行数值离散. 随后, 模拟了带有Dirichlet和周期性边界条件的NLS方程, 将所得结果与解析解做对比. 数值结果表明: 给出的SS-FPM粒子法的优点是在粒子分布非均匀情况下仍具有近似二阶精度, 且较网格类有限差分算法实施容易, 较已有改进的光滑粒子动力学方法计算误差小. 最后, 运用SS-FPM对无解析解的二维周期性边界NLS方程和Dirichlet边界玻色-爱因斯坦凝聚二分量GP方程进行了数值预测, 并与其他数值结果进行对比, 准确展现了非线性孤立波奇异性现象和量子化涡旋过程.
关键词: 非线性孤立波/
有限点集法/
分裂格式/
Gross-Pitaevskii方程

English Abstract

全文HTML

--> --> -->

含孤立波解的非线性问题常见于非线性光学和晶体的热脉冲物理现象^[1—3], 以及玻色-爱因斯坦凝聚态 (BEC)动力学特性^[3—6]中. 该非线性物理现象的研究通常会涉及非线性薛定谔(nonlinear Schr?dinger, NLS)方程或(Gross-Pitaevskii, GP)方程的求解^[5—7], 然而NLS/GP方程中含有的非线性项或旋转角动量算子使得多分量或高维情况下的孤立波形解难以用解析手段精确地得到^[7,8]. 目前已提出了多种数值方法对NLSE/GPE进行求解或对复杂孤立波传播过程进行预测, 如有限元法^[9]、有限差分法^[10,11]、蒙特卡罗法^[12,13]、修正的欧拉算法^[14]、时间分裂谱方法^[11,15]和基于背景网格无网格方法^[16—18]等. 然而上述基于网格的方法在多维复杂区域或非均匀节点分布情况下编程模拟实现都较复杂. 近些年来, 纯无网格方法(粒子方法)以其完全不依赖于网格的优势在实数域偏微分方程的模拟中得到了许多应用, 如光滑粒子动力学方法(smoothed particle hydrodynamics, SPH)方法^[19—23]和有限点集法(finite pointset method, FPM)^[24—26]等. 但上述纯无网格方法对复数域上含孤立波解非线性问题的模拟研究还处于起步阶段, 特别是FPM方法对非线性薛定谔方程的模拟在国际上还鲜有研究.
不依赖于背景网格的FPM方法具有粒子方法的特点, 其在非线性薛定谔方程的模拟应用上还处于试探性研究阶段, 这与其在数值稳定和数值精度方面还待进一步研究相关. FPM方法精度和稳定性的提高, 以及非线性薛定谔方程的数值模拟研究均是国际上的两个研究热点. FPM方法模拟复杂非线性薛定谔方程较网格类方法的主要优点在于: 可以任意布点, 不受区域复杂性限制; 计算空间导数可采用二阶精度显格式且不依赖于网格; 模拟程序易于实现且便于并行实施.
本文针对非线性薛定谔方程的特点, 为提高直接推广FPM方法模拟带非线性项问题的精度和稳定性, 首先将非线性薛定谔方程进行时间分裂分为非线性项和线性导数项; 其次, 引入Wendland权函数^[26], 采用基于Taylor展开和移动最小二乘思想的显式FPM格式对线性导数部分进行离散; 从而得到一种能够准确、高效地求解非线性薛定谔方程的基于时间分裂有限点集法 (split-step finite pointset method, SS-FPM). 所得数值结果表明, 提出的SS-FPM能够有效地求解带有Dirichlet和周期性边界条件的二维非线性薛定谔方程且具有二阶精度, 并能够准确预测出周期边界条件下孤立波变化奇异现象和Dirichlet边界条件下带外旋转BEC二分量的孤立波值随时间变化的量子化涡旋过程.

Gross-Pitaevskii方程常被用来描述量子力学中BEC动力学特性^[1,11], 本文考虑如下带角动量旋转项的无量纲化GP方程:

$\begin{split}{\rm{i}}\frac{{\partial \psi \left( {{{x}},t} \right)}}{{\partial t}} = & - \frac{1}{2}\Delta \psi + {V_d}\left( {{{x}},t} \right)\psi + {\beta _d}{\left| \psi \right|^2}\psi \\& -\theta {L_z}\psi,\;\;{{x}} \in {\mathbb{R}^d},\;t \geqslant 0,\end{split}$

初值条件

$\psi \left( {{{x}},0} \right) = \varphi \left( {{x}} \right),\;{{x}} \in \mathbb{R}{^d},$

边值条件

$\mathop {\lim }\nolimits_{\left| x \right| \to \infty } \psi \left( {{{x}},t} \right) = 0,\;t \geqslant 0,$

或周期边值条件

$\begin{split}\psi \left( {x,y,t} \right) = \psi \left( {x + 2{\text{π}}y,t} \right),\\\psi \left( {x,y,t} \right) = \psi \left( {x,y + 2{\text{π}}t} \right),\end{split}$

其中空间变量${{x}} = \left( {x,y} \right)$, 或$\left( {x,y,z} \right)$; $\psi \left( {{{x}},t} \right)$是d维空间${\mathbb{R}^d}\left( {d = 2,3} \right)$中的复值波函数; $\Delta $是d维空间的拉普拉斯算子, ${\rm{i}} = \sqrt { - 1} $是虚数单位; $\varphi \left( {{x}} \right)$是一个给定初值的复值函数. 无量纲实数${\beta _d}$通常用来描述三次非线性项的相互强度${\left| \psi \right|^2}\psi $, ${L_z} =$$ - {\rm{i}}\left( {x{\partial _y} - y{\partial _x}} \right)$是无量纲旋转角速度θ下的角动量算子的z分量, 势${V_d}\left( {x,t} \right)$是一个实值函数, 其形式为

$\begin{split}& {V_d}\left( {{{x}},t} \right) \\=& \left\{ \begin{array}{*{20}{l}}&\!\!\!\!\!\!\!\!\dfrac{\gamma _x^2 x^2}{2}, & \;d = 1, \\&\!\!\!\!\!\!\!\!\dfrac{\gamma _x^2 x^2 + \gamma _y^2 y^2}{2},& d = 2, \\&\!\!\!\!\!\!\!\!\dfrac{\gamma _x^2{x^2} + \gamma _y^2{y^2} + \gamma _z^2{z^2}}{2},& d = 3,\end{array} \right.\end{split}$

其中${\gamma _x},{\gamma _y},{\gamma _z}$为实常数; ${V_d}$还可以是周期性的函数.

本文针对超冷原子BEC动力学特性问题的粒子法模拟, 将时间分裂格式与FPM进行耦合, 得到能够准确模拟GP方程的SS-FPM. 本文提出的SS-FPM方法较直接拓展应用FPM方法^[24,25]模拟GP方程具有较高精度和更好的长时间模拟稳定性, 这是由于长时间模拟NLS/GP方程时, 方程中的非线性项对算法精度和稳定性的要求较高. 本文提出的SS-FPM方法的思想是: 首先, 应用时间分裂法将方程分解成带非线性和线性导数项的两个方程; 其次, 采用具有较好稳定性的Wendland权函数^[23]拓展应用FPM法对带线性导数方程进行二阶显格式离散; 最后, 采用二阶时间分裂伪谱法对两个方程进行交替求解.
2

-->

3.1.时间分裂格式

引入文献[7,10,15]中时间分裂 (time-split, TS)法, GP方程(1)可被写为

${\rm{i}}\frac{{\partial \psi \left( {{{x}},t} \right)}}{{\partial t}} = \left( {A + B} \right)\psi \left( {{{x}},t} \right),\;{{x}} \in {\mathbb{R}^d},$

其中$A = \!- (1/2)\Delta \!-\! \theta {L_z}$是线性微分算子, $B ={V_d}\left( {{{x}},t} \right) +$$ {\beta _d}{\left| {\psi \left( {{{x}},t} \right)} \right|^2}$是非线性算子. 于是上述方程可以被分解为如下的线性方程

${\rm{i}}\frac{{\partial \psi }}{{\partial t}} = A\psi $

和非线性方程

${\rm{i}}\frac{{\partial \psi }}{{\partial t}} = B\psi .$

本文采用的二阶分裂步法^[7]: 首先求解(8)式, 其次以(8)式的解作为初始条件求解(7)式, 最后以(7)式的解作为初始条件求解(8)式, 从而得到下一时间层的解.
2 -->

3.2.FPM离散格式

对方程(7)的求解采用FPM法, 时间上采用二阶龙格-库塔离散格式, 函数导数项近似采用显式FPM离散格式, FPM格式求解函数一阶导数和二阶导数的基本思想^[24,25]是: 第一步, 将求解域任意离散成有限点并赋初始值; 第二步, 每个点x在权函数支持域内有相邻粒子x_j ($j = 1,2, \cdots,n,$ n为支持域内相邻粒子数), 每个x_j在x上进行Taylor展开并保留到二阶导数项; 第三步, 对Taylor级数余项进行整理, 引入最小二乘法得到一个以x处一阶/二阶导数为未知函数的线性方程组; 第四步, 求解涉及局部矩阵的线性方程组得到每个点x处函数导数近似值; 第五步, 采用二阶龙格-库塔时间离散格式得到下一个时间层函数值.
本文以二维空间上均匀布点为例, 采用具有较好稳定性Wendland权函数^[23,27]. Wendland权函数^[23,27]形式如下:

$\begin{split}&{\omega _{ij}} = \omega \left( {r,h} \right) =\\ &{w_0}\left\{ {\begin{array}{*{20}{l}}{{{\left[ {2 - \left( {r/h} \right)} \right]}^4}\left[ {\left( {2r/h} \right) + 1} \right],}&{0 \leqslant \dfrac{r}{h} < 2,}\\{0,}&{r/h \geqslant 2,}\end{array}} \right. \end{split}$

其中$r = \left| {{{{x}}_j} - {{x}}} \right|$为支持域半径; ${\omega _0}$是一个正常数, ${\omega _0}$在二维空间中为$7/\left( {64{\text{π}}{h^2}} \right)$, h为光滑长度, 此处取$h \approx 1.1{\lambda _0}$($\lambda _0$为分布粒子初始间距).
对任意点x处函数$\psi \left( {{{x}},t} \right)$, 其支持域内相邻粒子为x_j, 函数$\psi \left( {{{{x}}_j},t} \right)$在x点的Taylor级数展开

$\begin{split}& \psi \left( {{{{x}}_j},t} \right) \\=& \psi \left( {{{x}},t} \right) + {\psi _x}\left( {{x_j} \!-\! x} \right) \!+\! {\psi _y}\left( {{y_j} - y} \right)\\ & + \frac{1}{2}{\psi _{xx}}{\left( {{x_j} - x} \right)^2} + {\psi _{xy}}\frac{1}{2}{\psi _{xx}}\left( {{x_j} - x} \right)\left( {{y_j} - y} \right)\\& +\frac{1}{2}{\psi _{yy}}{\left( {{y_j} - y} \right)^2} + {e_j},\end{split}$

其中e_j为Taylor级数展开时的余项误差. 将(10)式右端第一项放到等式左边, 则(10)式可以写成

${{e}} = {{Ma}} - {{b}},$

其中
${{M}} \!\!=\!\! \left(\!\!\!\!\! {\begin{array}{*{20}{c}}{{\rm{d}}{x_1}}&{{\rm{d}}{y_1}}&{\begin{array}{*{20}{c}}\!{\dfrac{1}{2}{\rm{d}}x_1^2}&\;\;\;\;{{\rm{d}}{x_1}{\rm{d}}{y_1}}&\;\;\;\;{\dfrac{1}{2}{\rm{d}}y_1^2}\end{array}}\\{\begin{array}{*{20}{c}}{{\rm{d}}{x_2}}\\ \vdots \\{{\rm{d}}{x_n}}\end{array}}&{\begin{array}{*{20}{c}}{{\rm{d}}{y_2}}\\ \vdots \\{{\rm{d}}{y_n}}\end{array}}&{\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}{\dfrac{1}{2}{\rm{d}}x_2^2}\\ \vdots \\{\dfrac{1}{2}{\rm{d}}x_n^2}\end{array}}&{\begin{array}{*{20}{c}}{{\rm{d}}{x_2}{\rm{d}}{y_2}}\\ \vdots \\{{\rm{d}}{x_n}{\rm{d}}{y_n}}\end{array}}&{\begin{array}{*{20}{c}}{\dfrac{1}{2}{\rm{d}}y_2^2}\\ \vdots \\{\dfrac{1}{2}{\rm{d}}y_n^2}\end{array}}\end{array}}\end{array}} \!\!\!\!\!\!\!\!\!\right),$

$ \begin{aligned} & {{a}} = {\left( {{\psi _x},{\psi _y},{\psi _{xx}},{\psi _{xy}},{\psi _{yy}}} \right)^{\rm{T}}},\\ & {{b}} = {\left( {{\psi _1} - \psi,{\psi _2} - \psi,{\psi _3} - \psi, \cdots,{\psi _n} - \psi } \right)^{\rm T}},\\ & {{e}} = {\left( {{e_1},{e_2},{e_3}, \cdots,{e_n}} \right)^{\rm{T}}}, \\ & {\rm d}x_j=x_j-x,\; {\rm d}y_j=y_j-y,\;(j=1,2,\cdots,n). \end{aligned} $

引入加权最小二乘法^[24]可以得到

$J = \mathop \sum \limits_{j = 1}^n {\omega _j}e_j^2.$

(12)式可以写为

${{J}} = {\left( {{{Ma}} - {{b}}} \right)^{\rm{T}}}{{W}}\left( {{{Ma}} - {{b}}} \right),$

其中

${{W}} = \left( {\begin{array}{*{20}{c}}{{\omega _1}}&0& \cdots &0\\0&{{\omega _2}}& \cdots &0\\ \vdots & \vdots & \ddots & \vdots \\0&0& \cdots &{{\omega _n}}\end{array}} \right),$

根据极值原理, J取最小可得

$\left( {{{{M}}^{\rm{T}}}{{WM}}} \right){{a}} = \left( {{{{M}}^{\rm{T}}}W} \right){{b}}.$

通过(13)式包含5 × 5局部系数矩阵的线性方程组求解可得x处一、二阶导数值.
2 -->

3.3.SS-FPM离散格式

结合3.1节二阶时间分裂格式与3.2节FPM离散格式, 对GP方程(1)进行离散, 可得如下SS-FPM离散格式.
对非线性方程(8)通过求解常微分方程的方式进行求解^[11],

$\psi _j^* = {{\rm{e}}^{ - {\rm{i}}\left( {{V_j} + \beta {{\left| {\psi _j^m} \right|}^2}} \right){\rm{d}}t/2}}\psi _j^m;$

将(14)式求解得到的结果代入有限点集法中进行计算, 再将求得的结果代入线性方程中, 得到

${\rm{i}}\left( {\frac{{\partial \psi }}{{\partial t}}} \right)_j^{**} = - \frac{1}{2}\Delta {\psi ^*} - \theta {L_z}{\psi ^*};$

最后将方程(15)得到的解代入非线性方程(8)中进行求解, 得到

$\psi _j^{m + 1} = {{\rm{e}}^{ - {\rm{i}}\left( {{V_j} + \beta {{\left| {\psi _j^{**}} \right|}^2}} \right){\rm{d}}t/2}}\psi _j^{**};$

其中$j = 0,1,2, \cdots,N$, N是区域上离散化的粒子数; $m = 0,1,2, \cdots $为时间层; (15)式通过3.2节FPM离散格式进行求解得到$\psi _j^{**}$.
时间步长为dt, 为保证数值模拟稳定性, 时间步长的选择通常需要满足限制性条件 (见文献[19,24,25]), 本文选取${\rm{d}}t \leqslant 0.1\lambda _0^2$ ($\lambda _0$为粒子初始间距).
初边值条件施加: 在 (14) 式和(15) 式计算中, 初边值条件的施加也是至关重要的, 初始条件(2)式可以在计算 (14) 式前进行准确施加$\psi _j^0 = \varphi \left( {x} \right)$. 边值条件(3)式采用文献[1,7,10]处理方式, 将其近似为齐次Dirichlet边值条件施加${\left. {{{\left. \psi \right|}^{m + 1}}} \right|_{\partial {{\varOmega }}}} = 0$($\partial {{\varOmega }}$为区域边界). 周期边值条件(4)式的施加中为保证边界上粒子的不足, 采用SPH粒子方法的施加方式 (见文献[19—21]), 在物理量值更新前需要施加周期性条件的边界上取大于等于支持域尺寸2h的粒子数及物理量赋到相对的边界上.

为验证提出的SS-FPM法模拟NLS方程的准确性和预测无解析解NLS/GP问题孤立波随时间演化特性的可靠性, 本节首先通过对带两种边值条件的NLS方程的求解, 与解析解做比较, 对SS-FPM法数值精度和收敛速度进行分析; 其次运用提出的粒子方法对周期边界条件下孤立波奇异特性和BEC中孤立波量子化涡旋过程进行数值预测, 并与其他数值结果(SS-FDM (split-step finite difference method)法^[4,7]和SS-ICPSPH (split-step implicit corrected parallel smoothed particle hydrodynamics)法)^[20]进行对比. 为分析所提方法的精度和收敛性, 本文定义如下误差 (精确解与数值解最大误差范数)和收敛阶为:

${e_r} = {\left\| {{\psi _{{\rm{exact}}}} - {\psi _{{\rm{num}}}}} \right\|_{\infty}},$

$Co = \frac{{{\rm{log}}\left[ {{e_r}\left( {{\lambda _{02}}} \right)/{e_r}\left( {{\lambda _{01}}} \right)} \right]}}{{{\rm{log}}\left( {{\lambda _{02}}/{\lambda _{01}}} \right)}},$

这里的${\lambda _{01}}$和${\lambda _{02}}$分别表示不同的粒子初始间距.
2

-->

4.1.SS-FPM方法精度和收敛性分析

运用SS-FPM对两个不同边值条件下NLS方程进行求解, 分析了所提方法模拟NLS方程的精度和收敛速度, 讨论了粒子分布非均匀情况下的数值误差.
3 -->

4.1.1.Dirichlet边值NLS方程

考虑正方形区域$\varOmega:\left[ {0,2{\text{π}}} \right] \times \left[ {0,2{\text{π}}} \right]$的Dirichlet边值条件的NLS方程^[7,23], 其对应的方程和初边值条件分别为

${\rm{i}}\frac{{\partial \psi }}{{\partial t}} = - \frac{1}{2}\left( {{\psi _{xx}} + {\psi _{yy}}} \right) + V\left( {x,y} \right)\psi + {\left| \psi \right|^2}\psi,$

其中$V\left( {x,y} \right) = 1 - {\sin ^2}x{\rm{si}}{{\rm{n}}^2}y$.
初值条件为

${\psi _0}\left( {x,y} \right) = \sin x\sin y,$

边值条件为

$\psi \left( {0,y} \right) = \psi \left( {2{\text{π}}y} \right) = \psi \left( {x,0} \right) = \psi \left( {x,2{\text{π}}} \right) = 0.$

对应该问题的解析解为

$\psi \left( {x,y,t} \right) = \sin x\sin y\exp \left( { - 2{\rm{i}}t} \right).$

该算例模拟中采用粒子初始间距为${{\text{π}}}/64$, 时间步长为dt = 10^–4. 图1给出了几个时刻两种粒子分布情况下SS-FPM法得到沿$y = 0.5{\text{π}} $处波函数变化曲线, 并与解析解进行对比. 图2展示了本算例模拟中采用的均匀和非均匀两种粒子分布情况, 其中图2(a)的粒子是均匀分布方式, 分别沿x, y方向每隔${\text{π}} /64$的距离分布一个粒子; 图2(b)中的粒子是非均匀分布方式, 以$\left({{\text{π}},{\text{π}}} \right)$为圆心, 靠近圆心第一层布6个粒子, 第二层布12个粒子, 逐步成等差数列方式向外扩展沿圆形分布粒子, 相邻两个圆形层之间距离为相等均为${\text{π}}/64$, 但边界上仍采用均匀布点方式. 由图1可知, 随时间演化提出的粒子方法得到的波函数与解析解吻合, 即使粒子分布非均匀下得到SS-FPM结果仍与解析解一致.
图 1 几个不同时刻处沿$y = 0.5{\text{π}}$的$\left| \psi \right|$变化曲线　(a) 粒子均匀分布; (b) 粒子非均匀分布
Figure1. The change curve of $\left| \psi \right|$ along $y = 0.5{\text{π}}$ at different time: (a) Uniform mode; (b) non-uniform mode.

图 2 两种不同的粒子分布　(a) 均匀粒子分布; (b) 非均匀粒子分布
Figure2. Two kinds of particle distribution: (a) Uniform mode; (b) non-uniform mode

为进一步体现SS-FPM方法模拟Dirichlet边值NLS方程的数值精度和收敛性, 表1和表2分别给出了提出方法得到数值结果的误差和收敛阶. 通过观察表1, 给出的粒子方法在粒子分布均匀和非均匀情况下得到的数值误差差距不大, 从而表明粒子方法模拟方程时无论区域是否规则都可以方便处理且具有较高的数值精度, 较网格类方法具更好的灵活推广应用性. 由表2可知, SS-FPM法具有二阶收敛速度, 与文献[23]中SS-ICPSPH法和文献[7]中网格类SS-FDM法的收敛阶基本一致, 但本文提出的粒子方法得到的数值误差较SS-ICPSPH法和SS-FDM法的误差小. SS-FPM方法相较于直接采用FPM方法, 前者具有较小误差和较快的收敛速度, 这是因SS-FPM方法对非线性薛定谔方程中非线性项采用了较准确的二阶精度分裂格式. 通过粒子法SS-FPM与网格类SS-FDM的构造过程以及粒子分布非均匀情况下图2和表1结果可知, 本文粒子方法不受网格限制, 可以在模拟区域任意布点情况下对问题进行纯无网格方法的模拟实现, 较网格类方法具有更好的灵活应用推广性, 易被推广应用于复杂非规则区域上薛定谔问题的模拟.

t	均匀分布/10–4	非均匀分布/10–4
0.5	2.48	3.22
1.0	4.94	6.12
1.5	7.40	9.26
2.0	9.88	16.50

表1粒子分布均匀/非均匀两种情况下的最大误差e_r
Table1.Maximum error e_r under uniform/non-uniform particles distribution

	粒子间距	误差	收敛阶
SS-ICPSPH	${\lambda _0} = {\text{π}}/32$	8.99×10–2
	${\lambda _0} = {\text{π}}/64$	2.23×10–3	2.007
	${\lambda _0} = {\text{π}}/128$	5.52×10–4	2.017
SS-FDM	${\lambda _0} = {\text{π}}/16$	2.016×10–2
	${\lambda _0} = {\text{π}}/32$	5.045×10–3	1.9986
	${\lambda _0} = {\text{π}}/64$	1.262×10–3	1.9997
FPM	${\lambda _0} = {\text{π}}/32$	4.6×10–3
	${\lambda _0} = {\text{π}}/64$	1.35×10–3	1.768
	${\lambda _0} = {\text{π}}/128$	3.86×10–4	1.806
SS-FPM	${\lambda _0} = {\text{π}}/32$	3.95×10–3
	${\lambda _0} = {\text{π}}/64$	9.88×10–4	2.000
	${\lambda _0} = {\text{π}}/128$	2.46×10–4	2.006

表2四种不同方法在t = 2时的数值收敛阶
Table2.The rate of convergence obtained using four different methods at t = 2

3 -->

4.1.2.周期边值NLS方程

为体现提出的SS-FPM方法对带周期性边值NLS方程模拟的准确性, 考虑正方形区域$\varOmega:$$\left[ {0,2{\text{π}}\left] \times \right[0,2{\text{π}}} \right]$的周期边界条件的非线性薛定谔方程, 其对应的方程和初值条件^[11]为

${\rm{i}}{\psi _t} + {{\Delta }}\psi + \beta {\left| \psi \right|^2}\psi = 0,$

初值条件

$\psi \left( {x,y,0} \right) = A{\rm{exp}}\left[ {{\rm{i}}\left( {{k_1}x + {k_2}y} \right)} \right],$

周期边界条件

$\begin{split}&\psi \left( {x,y,t} \right) = \psi \left( {x + 2{\text{π}}y,t} \right),\\&\psi \left( {x,y,t} \right) = \psi \left( {x,y + 2{\text{π}}t} \right),\\&\left( {x,y} \right) \in \mathbb{R} \times \mathbb{R},\;0 < t \leqslant T,\;T = 1.\end{split}$

对应的解析解为

$\psi \left( {x,y,t} \right) = A{\rm{exp}}\left[ {{\rm{i}}\left( {{k_1}x + {k_2}y - wt} \right)} \right],$

其中$w = k_1^2 + k_2^2 - \beta {\left| A \right|^2}$, $A = 1$, $\beta = - 2$.
模拟中, 选取${k_1} = {k_2} = 1$, 粒子初始间距${\text{π}}/64$, 时间步长${\rm{d}}t = {10^{ - 4}}$. 图3给出了两个不同位置、不同时刻SS-FPM法得到的波函数实部变化曲线, 并与解析解进行比较. 通过图3可知, SS-FPM方法准确地展示了波函数实部随时间演化呈现出的周期性波形, 且数值结果与解析解吻合, 从而表明提出的粒子方法可以准确地模拟带周期边值的NLS方程. 表3给出了周期边界条件下当t = 2时刻三种不同方法的数值收敛阶. 由表3可知, 本文提出的SS-FPM法与其他两种数值方法收敛阶接近, 且较SS-ICPSPH方法误差更小; 给出的SS-FPM法模拟周期边界下非线性薛定谔问题具有二阶收敛速度.
图 3 两个不同位置不同时刻ψ实部变化曲线图　(a) 沿对角线; (b) 沿$y = {\text{π}}$
Figure3. The change curve of real part at two positions with different times: (a) Along the diagonal; (b) along $y = {\text{π}}$

	粒子间距	误差/× 10–4	收敛阶
SS-ICPSPH	${\lambda _0} = {\text{π}}/32$	75.530
	${\lambda _0} = {\text{π}}/64$	18.280	2.046
	${\lambda _0} = {\text{π}}/128$	4.316	2.082
SS-FDM	${\lambda _0} = {\text{π}}/32$	24.670
	${\lambda _0} = {\text{π}}/64$	6.634	1.895
	${\lambda _0} = {\text{π}}/128$	1.725	1.943
SS-FPM	${\lambda _0} = {\text{π}}/32$	67.840
	${\lambda _0} = {\text{π}}/64$	16.790	2.015
	${\lambda _0} = {\text{π}}/128$	4.128	2.024

表3三种不同方法在t = 2时的数值收敛阶
Table3.The rate of convergence obtained using three different particle methods at t = 2

2 -->

4.2.数值模拟预测

为进一步体现SS-FPM方法模拟NLS/GP方程的能力, 采用SS-FPM对带周期边界NLS方程中孤立波奇异特性和BEC中GP方程描述的孤立波量子化涡旋过程进行数值预测.
3 -->

4.2.1.具有奇异性周期边值NLS方程

考虑正方形区域$\varOmega:\left[ {0,2{\text{π}}\left] \times \right[0,2{\text{π}}} \right]$的周期边界条件的非线性薛定谔方程, 其对应的方程和初值条件^[11]为

${\rm{i}}{\psi _t} + {{\Delta }}\psi + \beta {\left| \psi \right|^2}\psi = 0,$

初值条件

$\psi \left[ {x,y,0} \right] = \left[{1 + \sin \left( x \right)} \right]\left[ {2 + \sin \left( y \right)} \right],$

其中β = 1.
该算例为无解析解的带周期边界NLS方程, 它将描述随时间演化孤立波出现奇异特性, 常被用来验证一种数值方法预测周期性NLS方程出现奇异值现象的可靠性和稳定性^[11]. 本小节运用SS-FPM粒子法对该算例进行了模拟, 并与文献[7]中SS-FDM法结果进行对比(见图4). 观察图4知, SS-FPM方法得到的带周期边界NLS方程的奇异值与SS-FDM结果吻合, 表明提出的粒子方法模拟预测NLS方程描述的奇异特性是可靠的.
图 4 两个不同时刻波函数$\left| \psi \right|$三维图和等值线图　(a1), (a2) t = 0; (b1), (b2) t = 0.0108
Figure4. The 3D graphs and contour of $\left| \psi \right|$ at two different times: (a1), (a2) t = 0; (b1), (b2) t = 0.0108.

3 -->

4.2.2.BEC中角动量旋转GP方程

为验证SS-FPM方法预测GP方程描述BEC动力学特性的有效性, 考虑如下二分量带旋转项GP方程, 并与其他数值结果进行比较. 正方形区域${{\Omega}}:\left[ { - 8,8\left] \times \right[ - 8,8} \right]$上具有二分量角动量旋转GP方程, 及其对应的初边值条件^[8]为:

$\left\{ \begin{array}{l}{\rm{i}}\dfrac{{\partial {\psi _1}}}{{\partial t}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{\psi _1}}}{{\partial {x^2}}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{\psi _1}}}{{\partial {y^2}}} + \sigma \left( {{{\left| {{\psi _1}} \right|}^2} + \varsigma {{\left| {{\psi _2}} \right|}^2}} \right){\psi _1} + {w_1}\left( {x,y} \right){\psi _1} - \varTheta {L_z}{\psi _1} = 0,\\{\rm{i}}\dfrac{{\partial {\psi _2}}}{{\partial t}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{\psi _2}}}{{\partial {x^2}}} + \dfrac{1}{2}\dfrac{{{\partial ^2}{\psi _2}}}{{\partial {y^2}}} + \sigma \left( {{{\left| {{\psi _2}} \right|}^2} + \varsigma {{\left| {{\psi _1}} \right|}^2}} \right){\psi _2} + {w_2}\left( {x,y} \right){\psi _2} - \varTheta {L_z}{\psi _2} = 0,\end{array} \right.$

初值条件 ${\psi _1}\left( {x,y,0} \right) \!=\! \dfrac{{x \!+\! {\rm{i}}y}}{{\sqrt {\text{π}} }}\exp \!\!\left[\! \!{-\dfrac{{\left( {{x^2} \!+\! {y^2}} \right)}}{2}}\!\! \right], $${\rm{on}}\;\;\varOmega,$

${\psi _2}\left( {x,y,0} \right) = \frac{{x + {\rm{i}}y}}{{\sqrt {\text{π}}}}\exp \left[ { - \frac{{\left( {{x^2} + {y^2}} \right)}}{2}} \right],\;\;\;\;{\rm{on}}\;\varOmega ;$

边值条件: ${\psi _1}\left( {x,y,t} \right) = 0,\;{\rm{on}}\;\partial \varOmega \times \left[ {0,T} \right]$,

${\psi _2}\left( {x,y,t} \right) = 0,\;{\rm{on}}\;\partial \varOmega \times \left[ {0,T} \right].$

本文模拟中选取${w_1} \!=\! {w_2} = 1.5{x^2} + 0.5{y^2}$, $\varTheta = $$\;0.7,\;\sigma = - 100,\;\varsigma = 0.8$.
图5给出了两种数值方法得到的不同时刻沿x轴变化的波函数曲线. 由图5中两个不同时刻波函数变化曲线可知, BEC动力学特性在角动量旋转下的变化是复杂的; SS-FPM方法得到的数值结果与SS-FDM法^[7]得结果吻合, 表明提出的粒子方法模拟预测GP方程描述BEC动力学性质是可靠的. 图6给出了三个物理量${\rm{Re}}\left( \psi \right){\rm{,Im}}\left( \psi \right),\left| \psi \right|$三维数值结果, 可以明显观察到角动量旋转项影响下随时间演化的量子化涡旋变化情况. 值得注意的是图5和图6仅展示了第一分量的数值结果, 这是因为根据选取模拟参数两个分量变化是一致的.
图 5 两个不同时刻$\left| \psi \right|$沿x轴 (y = 0)变化曲线　(a) t = 0.05; (b) t = 0.25
Figure5. The change curve of $\left| \psi \right|$ along x-axis (y = 0) at two different times: (a) t = 0.05; (b) t = 0.25.

图 6 两个不同时刻${\rm{Re}}\left( \psi \right){\rm{,Im}}\left( \psi \right),\left| \psi \right|$的三维数值结果　(a1), (a2), (a3) t = 0; (b1), (b2), (b3) t = 0.25
Figure6. Three-dimensional numerical results of ${\rm{Re}}\left( \psi \right){\rm{,Im}}\left( \psi \right),\left| \psi \right|$ at two different times: (a1), (a2), (a3) t = 0; (b1), (b2), (b3) t = 0.25

为提高直接推广有限点集法模拟二维非线性薛定谔方程或Gross-Pitaevskii方程的数值精度, 本文给出了基于分裂格式的有限点集法(SS-FPM), 该方法兼顾时间分裂格式和传统FPM方法的优点. 数值算例中, 考虑了粒子分布均匀和非均匀情况下带有不同边界条件的二维非线性薛定谔方程, 并与解析解进行对比, 对SS-FPM方法的精度和收敛性进行了分析, 验证了数值预测的准确性; 随后运用SS-FPM法对无解析解NLS/GP方程进行了数值预测, 与其他数值结果进行了比较. 数值结果表明:
1) SS-FPM方法模拟二维非线性薛定谔方程具有二阶收敛速度, 较已有分裂有限差分法具有较好的灵活应用性, 且在粒子分布非均匀情况下仍具有较高数值精度;
2) SS-FPM法能够成功地预测带周期边界条件下孤立波传播过程的奇异现象, 且与其他数值结果相吻合;
3) SS-FPM法准确地预测了带外旋转BEC二分量的孤立波值量子化涡旋随时间的演变过程.
目前未见文献将时间分裂格式与FPM耦合对非线性薛定谔方程进行模拟研究, 本文针对不同边界条件下二维非线性薛定谔方程给出的SS-FPM法较网格类方法具有更好的灵活推广应用性, 为复杂区域上含孤立波非线性问题的数值预测提供了一种准确有效的粒子方法.