Vol. 44 No. 6Nov. 2020
第44卷第6期2020年11月江西师范大学学报(自然科学版)
Journal of Jiangxi Normal University ( Natural Science)
文章编号：1000-5862 (2020) 06-0599-05
2维Gross-Pitaevskii 方程的分裂高阶紧致差分格式
贺增甲1，孔令华"，符芳芳2
(1.江西师范大学数学与统计学院，江西南昌330022；2，南昌工学院基础部，江西南昌330108)
摘要：该文为带有旋转角动量的Gross-Pitaevskii 方程构造了分裂高阶紧致差分格式.首先通过时间分裂 将其分为线性方程和非线性方程，非线性方程可以通过质量守恒定律进行精确求解，线性方程通过高阶
紧致格式和局部1维方法进行离散，最终得到的格式时间方向2阶收敛和空间方向4阶收敛，并保持质 量守恒.最后用数值算例验证了格式的收敛阶以及质量守恒性.
关键词:Gross-Pitaevskii 方程；旋转效应；分裂方法；高阶紧致格式；质量守恒中图分类号：0 241.8 文献标志码:A DOI ：10.16357/j. cnki. issnlOOO-5862.2020.06.09
o 引言
20世纪20年代，萨特延德拉•纳特•玻色和
阿尔伯特•爱因斯坦预言了一种新的物态即Bo-
Einstein condensate ( BEC )的存在，直到 1995 年，人
们才在碱金属原子稀薄气体中真正获得了
BEC [1-2] ,Bo-Einstein 碱原子与氢的缩合反应已在
实验室中进行了广泛研究，使人们对宏观量子世界 有了新的认识•在温度T 远小于临界冷凝温度T c =
-273. 15七的情况下，根据平均场理论，旋转框架 下的BEC 可以用带有旋转效应的Gross-Pitaevskii (GP)方程来模拟，具体形式如下：
询= ( - 7]2v 2/(2m) + v,%) +
惊慌
2
NU 0 | ijj | - QL z )tfj{x ,t) ,x U R 2 ,t > 0,
其中m 为原子的质量，"为普朗克常数,N 为凝聚态 中的原子数，匕(％)是与外部势阱相对应的一个实
值函数是旋转激光束的角速度，乩表现为旋转
BEC 中粒子之间的相互作用,Q 是角动量在z 轴方
向上的分量，其定义为L* = - i(xdy - yto).
为了方便起见，将上述GP 方程做无量纲化处 理，得到如下无量纲形式2维GP 方程
i ◎妙二(- V 2/2 + V(x) +/3| if/1 - QLz)屮，策 e德国游记
〃 U RS > 0,
(1)
其中0为无量纲实数，兀二(％』)，*%)二(亡/ + 厂紀)/2几迅为频率常数.
在本文中，考虑齐次Dirichlet 边界条件
&(%,/) = 0,% g dU,t M 0, (2)同时给出初值
0(%,0) = G U,
(3)
通过计算，初边值问题(1)〜(3)保持质量守恒
Q (0(・』))=f I 0(兀丿)「血=0(0(・,0))=
Ju
Q o ,t^o
和能量守恒
E (申(• ,t))=[(」V (/r I 2/2 + V(x) | (/r | 2 +
旅游图片0 | & | /2 - Q 屮)ck 三 E Q ,t M 0,
其中/表示函数/的共辄.
在过去十几年中，一些学者对无旋GP 方程(12 =
0)进行了研究，构造了一些数值方法.如Bao Weizhu 等⑶构造了 4阶时间分裂正弦或傅里叶的
拟谱方法•符芳芳等⑷构造了其辛格式，因为是全
隐格式从而计算效率不高•对于计算无旋BEC 问题 的动力学行为，他们能够较为准确地模拟且非常有
效.对于带旋的GP 方程，由于旋转角动量显式含有 空间变量，导致在数值求解时存在一些困难.Wang
Lan 等⑸讨论了其保持能量结构特征的保能量式,
但是这种数值方法计算效率不高，难以用于求解大
区域细网格上较为准确的数值解•分裂步方法的基 本思想是把原问题分解成若干个更容易求解的子问
题，通过适当的组合方式，得到原问题的近似问题. 通过近似问题来得到原问题的近似解•同时，为
了提高空间方向的离散精度，又不增加太大的计算
收稿日期：2020-08-11
基金项目：国家自然科学基金(11961036)和江西省教育厅基金(GJJ200310)资助项目.
通信作者：孔令华(1977-)，男，江西石城人，教授，博士，博士生导师，主要从事微分方程数值方法的研究.E-mail :konglh@
mail. ustc. edu. cn
600江西师范大学学报（自然科学版）2020年
复杂度，充分利用节点信息，高阶紧致方法比通常有
限差分法能够更精确、更高效地模拟微分方程，对于
具有一定光滑性的问题，经常采用高阶紧致
方法金in.
为此，本文将对带有旋转效应的GP方程构造
出一种时间分裂高阶紧致差分格式，在时间上运用
分裂方法卩问把复杂问题分解成由线性方程和非
线性方程组成的更为简单的问题•对于非线性方程
可以精确求解，不存在离散误差，而对于线性问题使
用局部1维法以及高阶紧致格式进行求解•该
数值格式在时间方向上达到2阶精度，在空间方向
上达到4阶精度，并保持质量守恒且无条件稳定.
1高阶紧致格式
首先回顾1阶导数/；'=扌/（%）高阶紧致方
d%x=x
法的构造.对于1阶导数，假设有离散形式⑺创：
+/'+昭=a（加-几）/（2心，其中a、a为待
定系数.将以上各点的导数值或者函数值在勺处泰
勒展开，并根据误差阶合并同类项，得到系数满足的
线性代数方程组a=1+2a,a=6a，解得a=1/4,
a=3/2.则1阶导数的4阶格式能够整理为
（命+耀'+兀）/6=（Z+1-4J/C2A）,
因此1阶导数在齐次边界条件下的逼近能用矩阵表
示为Af'=昭其中
(4 11
41
14
1
1
4> 1
(0
-11 B=2h
-10
-11 0丿
同理,2阶导数力〃=
d2能够离散⑴]
X=Xj
为隊；+旷+戲+；二6（力+]-窈+久）/亿其中0、厶为待定系数.通过泰勒展开舍去高阶项得到系数关系满足的方程组6=1+20,6=120，解得0=1/10, b=6/5.因此,2阶导数的4阶格式能够写成
（启+1助"+昴）/12=（加-窈+几）/心
则2阶导数在齐次边界条件下的逼近能用矩阵表示为Cf"=巧，其中
(101A
1101
°•°•°・,
1101
I110丿
/-21\
1-21
<1-2丿
巧的反义词引理1⑸设A、c、D是同阶循环矩阵，则他们满足以下性质：
（i）CD=DC且A"、C"和CD也为循环矩阵;
（ii）对于A、C、D存在唯一对称正定矩阵乩、乩、尺3使得A=酣,C=疋,Q=Rj.
2 分裂高阶紧致格式
本节将采用时间分裂技巧、局部1维方法、高阶紧致格式等数值方法为2维GP方程构造出一个时间方向2阶精度、空间方向4阶精度且无条件稳定的数值格式•为此，首先对时空区域九，力]x[九，%]x[0，口进行网格剖分，并给出本文常用的定义及符号，用网格[亏，亏+J x[y4,n+1]x将区域进行剖分，将％方向分成J等份,y方向分成K等份，时间方向分成N等份，“=（X R-x L）/J,h y=（%-yJ/K,T=77N，贝呵=x L+jh x,y k=y L+kh y, t n=riT,j=0,1,2,…,J,k=0,1,2,■■■,K,t=0,1, 2,…,N；用坯=&（亏,:nJ”）和恃〜gjgtj 分别表示方程的精确解和数值解•为了表述方便，本文采用以下差分符号：
=Wj+l,k1Wj,k）/h*,6旳=（tpj.k1屮脱=（畅+i-0；”）/心松=（畅-
>n\/7,n a a,n<>2>n a<n
0/，盒-1"心玄0/，盒-§花北花松屮j,k=歹咖，盒，砥0/上-（*Aj+l,k1WU/Qh*）,§2沖;,*=（0；*+110；*一1）/（2心）.
为构造2维GP方程的分裂步高阶紧致格式，首先将方程（1）分裂为如下形式：
砂，=-（〃”+码）/2,（4）
砂，=02（xdy-（5）
2
妙，=+0|护|护.（6）
对于分裂后的子问题（6）满足点点的质量守恒律，即V%,t,
22
I0（%丿）I=I I-
事实上，在子问题（6）的两边同乘以広得
连接不可用—24
咖妙（光丿）=V（x）|（//（%,?）|+0|&（%』）I .
第6期贺增甲，等：2维Gross-Pitaevskii方程的分裂高阶紧致差分格式601
2
对上式两端取虚部得d|g,t)|/At=0.
将其代入(6)式中有
2
id妙(％,t)=(V(x)+0|0(%,t”)|)&(%,/).(7)再对(7)式从到i…+1上积分可得
=严⑴小旳心”汗“(吋”).
即得出了非线性子问题(6)的精确解.
对分裂后的GP方程(4)~(6)通过Strang组合为时间2阶精度的格式，再运用局部1维法和高阶紧致法离散为
(8)
(9)
(10)
(11)
(12)
因此对于分裂高阶紧致差分格式(8)~(12)，在每个时间步上仅需求解一些线性方程组，易知此格式在空间方向上具有4阶收敛精度，在时间方向上具有2阶收敛精度.在实际计算中比其他进行非线性迭代的方法成本更低，下面给出离散形式下的质量守恒.
引理2同对于任意的正整数p(N1)以及序列有
p P_1
T,a k8^b k=a p b p-a x b0-8x a k b k,
k=l k=l
特别地当％=a p ,b0=b p时有丫a i8i b k=-》8x a k b k.
k=1k=0
/八1©I-定理1定义W"||=丨如，
y j=i k=i
则分裂高阶紧致差分格式(8)~(12)满足离散的质量守恒，即110"*||=||||=…=||.
证将(8)式中的第1式与0⑴+做内积
丄〈〃⑴-0"妙⑴+〃"〉=-*〈C；&(〃⑴+
折纸飞机教案T O
0")妙⑴+F〉，(13)根据引理1和引理2得
〈C：&帥⑴+&"),〃“+0"〉=-〈尽戈辆⑴+ 0"),«2&(〃⑴+〃")〉，
此式为实数从而对(13)式两边取虚部有W⑴||=
II II-
将(9)式中的第1式与〃⑶+〃⑵作内积
丄〈0⑶-&⑵妙⑶+0⑵〉=¥〈曲；他(屮⑶+
T4
〃⑵)妙⑶+〃⑵〉，(14)根据引理1和引理2，得
〈A获2®⑶+〃⑵),〃⑶+〃⑵〉=
-〈艮(0⑶+0⑵)，乩％(0⑶+0⑵)〉，
则对(14)式两边取实部有||0⑶||=||||.
对于(10)式，两边取范数显然||0⑸||=
|助⑷II.同理可证，忖⑵||=W⑴II,W⑷II=
||沪II,n(6)II=w⑸II,II〃⑺II=II严,
II〃⑻II=II〃⑺II,II严II=II〃⑻II•
因此有||<A°+1||=||0"II•定理1得证.
3 数值实验
本节主要通过数值实验来检验所构造出的时间分裂高阶紧致格式的有效性，并且给出一些数值结果，通过表格及生成的图像直观地展示格式的收敛阶及离散的质量守恒等.
取[/=[_8,8]x[_8,8]，在方程(1)中取势
函数V(%)=(X2+y2)/2,/3=10，取初值0o(%)= 2(%+iy)e_"W/馅?.
由于方程的复杂性，因此不能给出其精确解的解析表达式•为便于分析比较，取非常精细的网格和较小的时间步长，如九=h y=1/64,r=0.001，本文构造的分裂高阶紧致格式(8)~(12)计算得到的数值解作为“精确解”•令e(h,T)表示在网格长度人和时间步长丁下误差的无穷范数||e"||“格式的收敛阶用如下公式进行计算：
_ln(e(/ii,T)/e(/i2,T))_ln(eG,m)/e(h,T2)) 01ln(hi/fi2)S ln(T/T2)'下面考察分裂高阶紧致格式(8)~(12)时间和空间方向的误差和收敛阶.为了计算空间方向的收敛阶，选取不同的空间步长％，并固定足够小的时间步长丁=0.001来得出相对于空间离散产生的误差，由时间离散产生的误差可以忽略不计.同理，为了计算时间方向的收敛阶，选取不同的时间步长T，并选取足够小的空间步长亿=h y=1/64来得出相对于时间离散产生的误差，由空间离散产生的误差可以忽略不计.以上计算都取时间/=0.5.由表1~表2可以看出，分裂高阶紧致格式(8)~(12)在时间方向上具有2阶收敛速度，在空间方向上具有4阶收敛速度
.
602江西师范大学学报（自然科学版）2020年
7=0.001时，空间方向的收敛速度和误差
表1当
2范数无穷范数
1/lb
II e"H201II e"||.01
1/4 5.6139x10-3 1.4175x10_3
n1/8 3.3307x10-4 4.0758.4808x IO-5 4.063
U
1/16 2.0257x10-5 4.039 5.1171x10" 4.051
1/32 1.1847x10-6 4.096 2.9932x IO-7 4.096
1/4 5.5721x10-3 1.4086x IO-3
1/8 3.3090x10-4 4.074&4222x IO-5 4.064
0.1
1/16 2.0129x10-5 4.039 5.0862x10" 4.050
1/32 1.1772x10" 4.096 2.9755x IO-7 4.095
表2当h x=h y=1/64时，时间方向的收敛速度和误差
2范数无穷范数
J
II e"||2°2II e"||.°2
1/32 1.6588x10-3 4.1264x10_4
n1/64 4.1379x10-4 2.003 1.0064x IO-4 2.040
1/128 1.0212x10-4 2.019 2.4699x IO-5 2.027
1/256 2.4251x10-5 2.074 5.8570x10" 2.076
1/32 1.6598x10-3 4.1432x IO-4
1/64 4.4194x10-4 2.000 1.0135x IO-4 2.030
0.1
1/128 1.0319x10-4 2.008 2.5039x IO-5 2.017
1/256 2.5066x10-5 2.041 6.0072x10" 2.059
下面讨论2维GP方程（1）所描述的Gaussian 脉冲随时间的演化关系，同时考虑格式对质量的保持情况•用分裂高阶紧致格式⑻〜（⑵在九1/32，丁二0.01下模拟此问题•数值解在f二0、5、10、15、20、25时的波形图如图1所示.图2描述了质量的误差随时间的演化关系.
=h y
0.07
0.06
M0.05
W0.04
M0.03
一0.02
0.01
图1
心0.15
全国211大学排名名单0.20-t=25数值解I就I
在不同时刻下的波形
第6期贺增甲，等：2维Gross-Pitaevskii 方程的分裂高阶紧致差分格式603
.A
〒O
I X
、«5S *
胆
6
120
5
10
15
20 25
图2离散质量的误差随时间的演化关系
由图1可以观察到Gaussian 脉冲的波形和振幅 随时间呈周期性变化，而由图2可以看出格式保持 系统的质量不变.
4结论
对于带有旋转角动量作用的GP 方程，为克服 高维、非线性在数值模拟中带来的困难和不足，采用
分裂技术，把原问题分解成若干个易于计算的局部
1维的子问题和可以精确求解的子问题•为提高计
算效率采用高精度且计算复杂度较低的高阶紧致方 法对空间导数进行离散.通过理论推导可知这样构
造的数值格式能够保持原问题的质量守恒，使得数 值计算更加稳定高效•数值实验表明所构造的数值 格式具有高精度高效率的优点，同时能够在较长时
间内模拟原问题所描述的物理现象•本文所构造的 数值方法可以直接扩展到其他类型微分方程的数值 求解，如3维GP 方程、多维Maxwell 方程等.
5参考文献
[1 ] Anderson M H ,Ensher J R ,Matthewa M R,et al. Obr
vation of Bo-Einstein condensation in a dilute atomic va
por [J], Science ,1995,269(5221) ： 198-201.
[2] Davis K B ,Mewes M 0,Andrews M R,et al. Bo-Einstein
condensation in a gas of sodium atoms [ J ]. Phys Rev Lett, 1995,75(22) ：3969-3973.
[3] Bao Weizhu,Jaksch D ,Markowich P A. Numerical solution
of the Gross-Pitaevskii equation for Bo-Einstein conden sation [J]. J Comput Phys,2003,187(1) ：318-342.
[4] 符芳芳，周媛兰.2维Gross-Pitae 询di 方程的辛格式[J].
江西师范大学学报：自然科学版,2016,40(6) ：599-602.
[5] Wang Lan,Cai Wenjun, Wang Yushun. An energy-prer
ving scheme for the coupled Gross-Pitaevskii equations [J], Adv Appl Math Meeh,2021,13(1) ：203-231.
[6] Kong Linghua , Hong Yuqi , Tian Nana, et al. Stable and
efficient numerical schemes for two-dimensional Maxwell equations in lossy medium [ J ]. J Comput Phys, 2019, 397：108703.
[7] Wang Hanquan. A time-splitting spectral method for cou
pled Gross-Pitaevskii equations with applications to rota ting Bo-Einstein condensates [ J]. J Comput Appl Math , 2007,205(1) ：88-104.
⑻ 陈萌，孔令华，王兰.Burge 第方程的跳点紧致格式[J].江
西师范大学学报：自然科学版,2017,41(5) ：526-530.
[9]翟步祥，聂涛，薛翔.5次非线性Schrodinger 方程的一
个线性化4层紧致差分格式[J].江西师范大学学报: 自然科学版,2019,43(1):35-38.
[10] Li Jichun ,Chen Yitung ,Liu Guoqing. High order compact
ADI methods for the parabolic equations [ J ]. Comput Math Appl,2006,52(8/9) ： 1343-1356.
[11 ] Lele S K. Compact finite difference schemes with spectral-
like solution [ J]. J Comput Phys ,1992,103(1) : 1642.[12] 孔令华，田娜娜，张鹏.2维Maxwell 方程的局部1维高
阶紧致格式[J]-江西师范大学学报：自然科学版, 2019,43(1) ：31-34,
[13] Bao Weizhu,Cai Yongyong. Optimal error estimates of fi
nite difference methods for the Gross-Pitaevskii equation with angular momentum rotation [ J]. Math Comp,2013, 82：99-128.
The Splitting High-Order Compact Difference Scheme for
Two-Dimensional Gross-Pitaevskii Equation
HE Zengjia 1, KONG Linghua 1 * , FU Fangfang 2
(1. School of Mathematics and Statistics , Jiangxi Normal University ,Nanchang Jiangxi 330022,China ；
2. Department of Fundamental Education , Nanchang Institute of Science and Technology , Nanchang Jiangxi 330108 , China)
Abstract : The splitting high-order compact difference scheme for the Gross-Pitaevskii equation with angular momen tum rotation term is constructed. Firstly , the equation is divided into linear equations and nonlinear equations by time splitting method. Secondly ,the nonlinear equations can be accurately solved by the conrvation law of mass, and the linear equation is discretized by a high-order compact scheme and a local one -dimensional
method. The re sulting scheme converges cond-order in time and fourth-order in space while maintaining mass conrvation. Final ly ,numerical experiments verify the convergence orders and mass conrvation of the scheme.
Key words ： Gross-Pitaevskii equation ； rotating effect ； splitting method ； high-order compact scheme ； mass conrva tion law (责任编辑：曾剑锋
)野景
