含阻尼效应的非线性薛定谔方程的共形分裂高阶紧致差分格式

罗奕杨,王 兰*,万 隆,孔令华

(1.江西师范大学数学与统计学院,江西 南昌 330022;2.豫章师范学院小学教育学院,江西 南昌 330103)

0 引言

非线性薛定谔方程是在数学物理中非常重要的模型,它在材料科学、等离子体物理学、量子物理学、非线性光学、双分子动力学等领域中得到了广泛应用.含有阻尼效应的非线性薛定谔方程描述了在非线性介质中的几种谐振现象,特别是在垂直振荡水槽中的非线性法拉第谐振[1]和相敏放大器对光纤中孤子的影响[2].本文考虑如下带有阻尼效应的非线性薛定谔(DNLS)方程:

iψt+ψxx+β|ψ|2ψ+iαψ=0,x∈Ω,0

(1)

的周期初边值问题,取初始条件为

ψ(x,0)=ψ0(x),x∈Ω,

(2)

其中ψ(x,t)是定义在Ω×[0,T]上的未知复函数,i是虚数单位,β是不为0的常数(β>0和β<0分别表示方程是聚焦和散焦的情况),α(>0)表示系统耗散率.若α=0,则方程变为一般的薛定谔方程,其质量与能量以及动量是守恒的.

根据R. McLachlan等[3]的研究,初边值问题(1)～(2)保持共形质量守恒律

以及共形动量守恒律

近年来,科研工作者对DNLS方程进行了一些研究,构造了一些数值格式,如Bao Weizhu等[4]提出了2阶时间分裂谱方法,Jiang Chaolong等[5-6]利用傅里叶伪谱法构造了非线性隐式与线性隐式2种方法,Cai Jiaxiang等[7]提出了紧致隐式积分算子方法,Hu Weipeng等[8]提出了DNLS方程的保辛格式等.上述工作虽然有很多优点,但是部分方法依然存在精度不高、计算效率较低或共形质量守恒律不保持等问题.

基于此,本文提出了一种高精度、高效率且保持共形质量守恒律的有限差分方法,该方法结合了时间分裂方法[9-10]与高阶紧致方法[11](这2种方法已经得到广泛运用[12-15]).其中时间分裂方法可以将复杂的原问题分裂成几个易于求解的子问题,便于构造简便的数值格式;且高阶紧致方法是一种利用了更少的模板节点达到了更高精度的方法,它适合于构造高精度、高效率的数值格式.

1 共形分裂高阶紧致格式的构造

本部分先对Ω×[0,T]进行网格剖分.设Ω=[a,b],将区间[a,b]作J等分,将区间[0,T]作N等分,记空间步长与时间步长分别为h、τ,则h=(b-a)/J,τ=T/N,可以得到时空区域Ω×[0,T]的一致网格剖分

Ωhτ={(xj,tn)|xj=a+jh,tn=nτ,j=0,1,…,J,n=0,1,…,N}.

为构造DNLS方程(1)的高阶精度格式,首先回顾一下高阶紧致方法.对于具有一定光滑性的函数u(x),其2阶导数u″(x)可以由u(xj+1)、u(xj)、u(xj-1) 3点近似[16],即

(u(xj+1)-2u(xj)+u(xj-1))/h2=u″(xj)+h2u(4)(xj)/12+O(h4)=u″(xj)+O(h2),

若略去高阶项O(h2),则得到的格式具有2阶精度.若再将上式中的4阶导数d4u(x)/dx4写成d2(d2u(x)/

dx2)/dx2,对前一2阶微分采用2阶中心差商离散,后一2阶微分保持不变,则得到如下离散格式:

(u″(xj+1)+10u″(xj)+u″(xj-1))/12=(u(xj+1)-2u(xj)+u(xj-1))/h2+O(h4).

忽略高阶项O(h4),得到截断误差为4阶的紧致格式:

(3)

此格式只用了3个节点xj-1、xj、xj+1就得到了4阶精度的差分格式,若采用一般的差分离散,则这是不可能达到的.在周期边界条件下,将其写成矩阵形式u″=P-1Qu,其中

为构造DNLS方程(1)高效的分裂步高阶紧致格式,先把它分裂成如下形式:

iψt+ψxx=0,

(4)

iψt+β|ψ|2ψ=0,

(5)

iψt+iαψ=0,

(6)

由用Crank-Nicolson格式对式(4)进行时间离散和对高阶紧致格式(3)进行空间离散可得

(7)

这说明各个点在各个时刻处的质量|ψ(x,t)|2保持不变,即

|ψ(x,tn+1)|2=|ψ(x,tn)|2.

这说明方程(5)满足逐点的质量守恒律.这使得非线性方程(5)可以精确求解.将其代入式(5),则可得方程(5)的精确解

ψ(xj,tn+1)=eiβ|ψ(xj,tn)| 2τψ(xj,tn).

(8)

对于式(6),由观察可知,此方程为线性常系数常微分方程,可以用分离变量法得其精确解

ψ(xj,tn+1)=e-ατψ(xj,tn).

(9)

最后,将方程(4)～(6)的解(7)～(9)利用Strang组合,得到在时间上2阶精度的格式:

(10)

其中j=1,2,…,J,n=0,1,2,…,N-1.

格式(10)在时间上具有2阶精度、在空间上具有4阶精度;且在计算过程中,格式(10)在每次时间的推进时只需要求解线性方程组以及简单的代数运算,不需要额外的非线性迭代,计算成本较低,计算效率非常高,这克服了已有数值格式的各种不足.不仅如此,格式(10)还保持DNLS方程固有的一些物理性态,如共形质量守恒.

2 数值格式的共形守恒律以及稳定性

为分析格式(10)的共形守恒律,需要运用到循环矩阵的一些性质.

引理1[15]设X、Y为同阶循环矩阵,则其满足如下性质:

(i)XY=YX,X-1、Y-1、XY也都是循环矩阵;

(ii)对于循环对称正定矩阵X,存在唯一的矩阵Z,使得X=ZTZ.

定理1DNLS方程(1)的Strang分裂格式(10)保持如下的离散共形质量守恒律,即

证将式(10)中的第1式与第2式两边取模并平方,关于空间指标j求和,可得

‖ψ(1)‖2=e-ατ‖ψn‖2,‖ψ(2)‖2=‖ψ(1)‖2,

同理可得‖ψ(4)‖2=‖ψ(3)‖2,‖ψn+1‖2=e-ατ·‖ψ(4)‖2.

用ψ(3)+ψ(4)与式(10)中第3式作内积,取虚部并利用引理1可得

‖ψ(3)‖2=‖ψ(2)‖2.

综上所述有‖ψn+1‖2=e-2ατ‖ψn‖2,即格式(10)满足离散共形质量守恒律Mn+1=e-2ατMn.同时这表明了格式(10)是无条件稳定的,原因是

‖ψn+1‖2=e-2ατ‖ψn‖2=e-2αtn+1‖ψ0‖2≤‖ψ0‖2.

3 数值实验

本部分主要通过数值实验来检验所构造的格式(10)的有效性,并验证该格式的收敛精度和对共形质量守恒律的保持情况.此外,考虑到DNLS方程(1)具有共形动量守恒律,定义方程(1)对应的离散动量表达式为

其中Dx为1阶导数ψx的离散算子δ2x所对应的反对称矩阵.通过上述公式来模拟格式的动量变化情况.

取空间为Ω=[-25,25],DNLS方程(1)中的参数β=2、α=0.005,取初值为ψ0(x)=sech(x)·e2ix.由于无法给出其精确解,所以为便于分析,将利用格式(10)在相对小的时间步长与空间步长下得到的数值解作为“精确解”.这里取h=1/32,τ=1×10-4.为了计算收敛阶,使用公式

oτ=ln(‖e(h,τ1)‖γ/‖e(h,τ2)‖γ)/ln(τ1/τ2),

oh=ln(‖e(h1,τ)‖γ/‖e(h2,τ)‖γ)/ln(h1/h2)

来分别计算时间方向与空间方向的收敛阶数,其中‖e(h,τ)‖γ表示在空间步长h以及时间步长τ下误差的γ范数,γ=2或∞.

在计算空间方向的收敛阶时,可以选取不同的空间步长,并固定一个相对小的时间步长来数值模拟,以此来减少时间方向对误差的影响;在计算时间方向的收敛阶时同理.选择τ=1×10-4,h=1/2、1/4、1/8、1/16来计算空间方向的收敛阶,选择h=1/32,τ=1/25、1/50、1/100、1/200来计算时间方向的收敛阶,问题计算到t=1.从表1和表2容易看出:格式(10)在空间方向上有4阶收敛率,在时间方向上有2阶收敛率.这与理论分析结果相符.

表1 空间收敛阶及误差(τ=1×10-4)

表2 时间收敛阶及误差(h=1/32)

(a)α=0 (b)α=0.005图1 在不同时间上的波形

图2与图3分别展示了当α=0和α=0.005时质量与动量随时间的演化情况.图4描绘了质量和动量的衰减率随时间的演化关系,其纵坐标的数据由(ln(Qn/Q0))/2±0.1(其中Qn=Mn或In)计算得到,中间的线条是参照线y=-αt=-0.005t.从图2可以看出:当α=0时,质量与初始值相比是基本保持不变的,符合薛定谔方程的质量守恒律,而动量与初始值相比具有一定的误差.由图3和图4可以看到:质量随时间的衰减率完全与理论分析的结果一样,且动量也随着时间而衰减.

(a)质量 (b)动量图2 当α=0时质量与动量随时间的演化

(a)质量 (b)动量图3 当α=0.005时质量与动量随时间的演化

图4 质量与动量随时间的衰减率

4 结论

本文通过采用分裂步方法和高阶紧致方法,构造了含有阻尼效应的非线性薛定谔方程的共形分裂高阶紧致方法.通过采用分裂步技术把原方程分解成更简单的2个微分方程,其中非线性问题能够精确求解,这可以避免求解非线性代数方程组.对于空间导数的离散采用高效的高阶紧致格式,这样得到的数值格式不仅能够保持原方程的共形结构,而且计算效率还非常高.这种数值格式的构造思路可以推广到其他微分方程保结构格式的构造中,在未来的工作中将深入讨论此类格式.