非线性Schrödinger方程的算子分裂配点格式

姚梦丽, 田朝薇, 翁智峰

(华侨大学 数学科学学院, 福建 泉州 362021)

Schrödinger方程是量子力学领域的基本方程,对物理领域的研究具有深远的意义.非线性Schrödinger方程(NLS方程)在量子力学、非线性光学、电磁学、等离子理论、固体物理等领域中有着广泛的应用.考虑以下NLS方程的初边值问题,即

(1)

式(1)中：i为虚数单位;a,b,ρ,β均为常数;α(x),v(x)均为已知函数;u为未知函数.

众所周知,NLS方程满足质量守恒和能量守恒,有

(2)

(3)

式(2),(3)中：M(t)为质量函数;E(t)为能量函数.

Schrödinger方程作为一类重要的量子力学方程,在实际复杂的系统中,由于包含复数且为耦合的问题,其精确解往往不容易求得.因此,对其进行数值求解具有重要意义.许多研究者设计了有效的数值算法求解非线性Schrödinger方程.孙传志等[1]研究NLS方程的几种差分格式.王廷春等[2]用紧致有限差分格式求解非线性耗散Schrödinger方程.Gong等[3]采用Fourier拟谱方法构造二维NLS方程的能量和质量守恒格式.Cui等[4]用标量辅助变量(SAV)方法构造NLS方程的任意高阶保结构指数Runge-Kutta方法.Feng 等[5]研究NLS方程的高阶守恒SAV-Gauss配置有限元格式.张法勇等[6]研究带五次项的NLS方程的守恒差分格式.Wang等[7]用两层网格有限元方法对NLS方程进行超收敛性分析.Wang等[8]采用伽辽金有限元法求解阻尼NLS方程.Hu等[9]采用牛顿迭代Crank-Nicolson有限元方法求解NLS方程.Chen等[10]采用两层网格有限体积元法求解NLS方程.以上方法都是基于网格剖分思想求解NLS方程数值解问题.最近,一种新型的无网格方法重心插值配点法广泛应用于求解微分方程初边值问题.重心插值配点法具有计算格式简单、精度高、程序实施方便、节点适应性好等特点,特别使用Chebyshev节点的重心Lagrange插值公式可以有效克服Runge现象.目前,重心插值配点法已被推广到求解高维Fredholm积分方程[11]、耦合Burgers方程[12]、粘弹性波方程[13]、Allen-Cahn方程[14-16]、Cahn-Hilliard方程[17]、Black-Scholes方程[18]、分数阶电报方程[19]等.

算子分裂法[20-22]将复杂的问题分解成几个简单的子问题,并根据子问题的性质采用不同的求解方法,用算子分裂格式将子问题结合起来得到原问题的解,可降低原问题的求解规模,此方法广泛应用于复杂非线性模型的数值求解中.基于此,本文提出一种非线性Schrödinger方程的算子分裂配点格式.

1 重心Lagrange插值配点格式

1.1 重心Lagrange插值

设有插值节点xj和一组对应的实数gj(j=0,1,…,m),采用多项式插值,则在次数不超过m的多项式空间可以找到唯一的插值多项式p(x),满足p(xj)=gj.将p(x)写成Lagrange插值公式,即

(4)

式(4)中：Lj(x)为Lagrange插值基函数,满足

(5)

以及

(6)

令

l(x)=(x-x0)(x-x1)…(x-xm),

(7)

定义重心插值权wj为

(8)

则Lagrange插值基函数Lj(x)可表示为

(9)

将式(9)代入式(4),可得

(10)

当p(x)=1时,有恒等式

(11)

用式(10)除以式(11),可得重心Lagrange插值公式为

(12)

由式(8)可知：重心插值权仅依赖插值节点的分布,因此,对于一些特殊的节点分布可以简化重心插值权.特别地,当选择Chebyshev点族时,可保证重心Lagrange插值具有良好的数值稳定性.

利用重心插值逼近未知函数时均采用第二类Chebyshev节点,即

(13)

其重心Lagrange插值的重心插值权为

(14)

1.2 重心Lagrange插值微分矩阵

对区间[a,b]上的节点a=x1

(15)

由式(15),未知函数g(x)在节点x0,x1,…,xm处的p阶导数表示为

(16)

将式(16)写成矩阵形式,即

g(p)=D(p)g.

(17)

由文献[23]可知,一阶和二阶微分矩阵元素的计算公式为

(18)

通过数学归纳法,可得p阶微分矩阵的递推公式为

(19)

2 基于重心Lagrange插值的NLS方程算子分裂格式

2.1 Strang算子分裂法

将式(1)改写为

iut=L(u)+N(u).

(20)

式(20)中：L(u)=-ρΔu;N(u)=-v(x)u-β|u|2u.

算子分裂法的本质是将方程(20)分为两个子问题,即

线性部分SA：iut=L(u)=-ρΔu,

(21)

非线性部分SB：iut=N(u)=-v(x)u-β|u|2u.

(22)

式(21),(22)中：SA和SB分别为方程(21)和方程(22)的精确解,可得Strang算子分裂格式为

(23)

2.2 SA和SB的数值逼近

设空间节点xj(j=0,1,…,m)为Chebyshev节点,时间节点tk=kτ(k=0,1,2,…,n).

首先,在空间方向上采用重心Lagrange插值配点法,可得半离散格式为

(24)

将式(24)写成矩阵形式,即

i(uh(t))t=-ρD(2)uh(t),

(25)

式(25)中：uh(t)=[u0(t),u1(t),…,um(t)]′.

(26)

因此,有

(27)

方程(22)采用解析求解,即

(28)

则算子分裂格式为

(29)

3 相容性分析

考虑线性子问题的空间半离散格式(24)的相容性分析.设p(x)为u(x)的Lagrange插值函数,满足p(xj)=u(xj),j=0,1,…,m.定义误差函数r(x)为

r(x)=u(x)-p(x).

(30)

引理1[19]假设u∈Cm+1(a,b),则有

式中：lx表示区间[a,b]的一半.类似地,可得

为了给出线性子问题SA的空间半离散格式(24)的相容性分析,首先,定义算子Γ为

(31)

定理1如果u∈C1([0,T])∪Cm+1([a,b]),有

证明：由式(21)和式(24),可得

Γu(x,t)-Γu(xj,t)=iut(x,t)+ρuxx(x,t)-iut(xj,t)-ρuxx(xj,t)=
iut(x,t)-iut(xj,t)+ρuxx(x,t)-ρuxx(xj,t)=R1+R2.

(32)

式(32)中：

R1=iut(x,t)-iut(xj,t),R2=ρuxx(x,t)-ρuxx(xj,t).

(33)

对于R1,有

R1=iut(x,t)-iut(xj,t)=i(ut(x,t)-ut(xj,t))=irt(xj,t).

(34)

由引理1可得

(35)

类似地,可得

(36)

将式(35),(36)代入式(32)即可,证毕.

4 数值算例

定义L∞误差和L2误差分别为

式中：ui为点xi处的精确解;Ui为点xi处的数值解.

4.1 算例1

比较重心Lagrange插值配点格式和有限差分格式的数值结果,以验证插值配点格式的有效性及高精度.ρ=1,β=1,选取带右端项的真解u(x,t)=costsin πx,0≤x≤1,右端项为f(x)=-isintsin πx-π2costsin πx+|costsin πx|2costsin πx,初始条件为u(x,0)=sin πx,0≤x≤1, 边界条件为u(0,t)=u(1,t)=0,0≤t≤1,势函数v(x)=0.

固定时间步长(τ=10-4),改变空间节点数(M),重心Lagrange插值配点格式和有限差分格式的误差,如表1所示.由表1可知：算例1的重心Lagrange插值配点格式选取6个节点差值即可达到10-5量级,而有限差分格式选取256个节点差值才能达到10-5量级,这表明重心Lagrange插值配点格式在空间上用较少的点可达到更高的精度.

表1 重心Lagrange插值配点格式和有限差分格式的误差(算例1)Tab.1 Errors of barycentric Lagrange interpolation collocation scheme and finite difference scheme (example 1)

固定空间节点数(M=8),改变时间步长,重心Lagrange插值配点格式的误差及时间收敛阶,如表2所示.表2中：Rt,1,Rt,2分别为L∞误差和L2误差的时间收敛阶.由表2可知：基于重心Lagrange插值配点格式的NLS方程在时间上是二阶精度.

表2 重心Lagrange插值配点格式的误差及时间收敛阶(算例1)Tab.2 Errors and time convergence orders of barycentric Lagrange interpolation collocation scheme (example 1)

不同数值格式的空间收敛阶,如图1所示.图1中：Rs为空间收敛阶.由图1可知：有限差分格式的空间收敛阶是二阶精度;重心Lagrange插值配点格式的空间收敛阶满足指数收敛的性质.

图1 不同数值格式的空间收敛阶(算例1)Fig.1 Spatial convergence orders for different numerical schemes (example 1)

4.2 算例2

算例2的重心Lagrange插值配点格式和有限差分格式的误差(τ=10-4),如表3所示.由表3可知：算例2的重心Lagrange插值配点格式选取8个节点差值即可达到10-4量级,而有限差分格式选取128个节点差值也只能达到10-4量级,故重心Lagrange插值配点格式具有较高的精度.

表3 重心Lagrange插值配点格式和有限差分格式的误差(算例2)Tab.3 Errors of barycentric Lagrange interpolation collocation scheme and finite difference scheme (example 2)

重心Lagrange插值配点格式的数值解和精确解(M=64,t=1),如图2所示.不同数值格式的空间收敛阶,如图3所示.重心Lagrange插值配点格式守恒量的保持情况,如图4所示.

(a) 数值解 (b) 精确解图2 重心Lagrange插值配点格式的数值解和精确解Fig.2 Numerical solution and exact solution of barycentric Lagrange interpolation collocation scheme

图3 不同数值格式的空间收敛阶(算例2)Fig.3 Spatial convergence orders for different numerical schemes (example 2)

(a) 总质量 (b) 质量误差

由图2可知：NLS方程的重心Lagrange插值配点格式是稳定的.由图3可知：有限差分格式是二阶精度,重心Lagrange插值配点格式呈指数递减.由图4可知：基于重心Lagrange插值配点格式的NLS方程满足质量守恒和能量守恒.

4.3 算例3

通过算例3进行孤立波的演化实验.选取空间节点M=248,时间步长τ=10-4,对应的计算区间为[-10π,10π]×[0,T],ρ=1,β=2,选取真解u(x,t)=exp(i(2x-3t))sech (x-4t),初始条件为u(x,0)=exp(2ix)sechx,边界条件为u(a,t)=u(b,t)=0,势函数v(x)=0.

算例3的重心Lagrange插值配点格式的误差及时间收敛阶(M=248),如表4所示.由表4可知：算例3验证了时间收敛阶为二阶精度.

表4 重心Lagrange插值配点格式的误差及时间收敛阶(算例3)Tab.4 Errors and time convergence orders of barycentric Lagrange interpolation collocation scheme (example 3)

初始孤立波波形,如图5所示.

图5 初始孤立波波形Fig.5 Initial isolated wave shape

重心Lagrange插值配点格式的孤立波波形,如图6所示.由图6可知：算例3的重心Lagrange插值配点格式的孤立波波形随时间的变化而变化,符合物理意义,而且每个时刻的波形都没有出现震荡,这说明该格式是稳定的.

(a) t=1 (b) t=3 (c) t=5图6 重心Lagrange插值配点格式的孤立波波形Fig.6 Isolated wave shapes of barycentric Lagrange interpolation collocation scheme

重心Lagrange插值配点格式守恒量的保持情况,如图7所示.由图7可知：算例3的重心Lagrange插值配点格式保持质量守恒和能量守恒.

(a) 总质量 (b) 总能量图7 重心Lagrange插值配点格式守恒量的保持情况(算例3)Fig.7 Conservation of conserved quantities of barycentric Lagrange interpolation collocation scheme (example 3)

5 结束语

将算子分裂格式与重心Lagrange插值配点格式相结合,提出一种求解NLS方程的有效数值格式.此外,对方程的非线性部分采用解析求解,降低非线性迭代的计算量,提高计算效率,并给出线性子问题空间半离散格式的相容性分析.最后,通过数值算例验证格式的高精度和有效性.与经典的有限差分格式进行比较可知,文中格式可用较少的点达到较高的精度.今后,可将该方法推广到其他非线性微分方程,为解决同类问题提供一种高精度数值求解方法,如耦合的非线性Schrödinger方程、时间分数阶非线性Schrödinger方程等.