本文考虑如下的Schrödinger方程初值问题[1],

$ \begin{equation} \begin{aligned} {i}\varepsilon{u}_t^\varepsilon+\frac{\varepsilon^2}{2}u_{xx}^\varepsilon-V{u}^\varepsilon -f(|{u}^\varepsilon|^2){u}^\varepsilon-\varepsilon\beta\arg({u}^\varepsilon){u}^\varepsilon=0, t>0, x\in R. \end{aligned} \end{equation} $

(1.1)

$ \begin{equation} \begin{aligned} {u}^\varepsilon(x, t=0)={u}_0^\varepsilon(x), x\in R. \end{aligned} \end{equation} $

(1.2)

其中$ V=V(x)>0 $为已知的实函数, $ \varepsilon $为普朗克常数, $ f $为足够光滑的实函数, $ \beta\geq0 $为常数, $ {u}^\varepsilon={u}^\varepsilon(x, t) $为复波函数, 并有

$ \begin{align} \varepsilon\arg({u}^\varepsilon(x, t))={S}^\varepsilon(x, t). \end{align} $

(1.3)

方程(1.1) 包含很多的Schrödinger方程[1].例如, 当$ f\equiv0 $且$ \beta=0 $时, (1.1) 变为线性Schrödinger方程; 当$ V\equiv0 $, $ f(\rho)\equiv\gamma_\varepsilon\rho $且$ \beta=0 $时, 变为广义Schrödinger方程; 当$ f(\rho)=\alpha(\rho-1) $且$ \beta=0 $时, 变为潜水波方程. 方程$ (1.1)\sim(1.2) $遵守以下守恒律,

$ \begin{align} \int_{-\infty}^{+\infty}|{u}^\varepsilon(x, t)|^2dx=const, \end{align} $

(1.4)

如果$ \beta=0 $还遵守下面的守恒律,

$ \begin{align} \int_{-\infty}^{+\infty}E^\varepsilon(x, t)dx=const, \end{align} $

(1.5)

其中

$ \begin{align} E^\varepsilon(x, t)=\frac{\varepsilon^2}{2}|{u}_x^\varepsilon(x, t)|^2+F(|{u}^\varepsilon(x, t)|^2)+V(x)|{u}^\varepsilon(x, t)|^2, \end{align} $

(1.6)

并且

$ F(s)=\int_0^sf(z)dz. $

在文[1]中作者对方程(1.1)–(1.2)的周期初边值问题研究了两个分裂的谱格式, 文[2]对一类半经典长短波方程构造了分裂龙格-库塔格式, 文[3]将分裂的龙格-库塔格式用于耦合Schrödinger方程, 都取得了很好的数值结果. 本文将分裂的三点、四点龙格-库塔格式和一个分裂谱格式应用于半经典Schrödinger方程, 数值结果表明本文的算法是有效的和可靠的. 周期边值问题如下:

$ {u}^\varepsilon(x_L+x, t)={u}^\varepsilon(x_R+x, t), \ {u}_x^\varepsilon(x_L+x, t)={u}_x^\varepsilon(x_R+x, t), $

方程(1.1)–(1.2)变为

$ \begin{align} {i}\varepsilon{u}_t^\varepsilon+\frac{\varepsilon^2}{2}{u}_{xx}^\varepsilon-V(x){u}^\varepsilon -f(|{u}^\varepsilon|^2){u}^\varepsilon-\varepsilon\beta\arg({u}^\varepsilon){u}^\varepsilon=0, \ x_L<x<x_R, \ t>0. \end{align} $

(1.7)

$ \begin{align} {u}^\varepsilon(x, t=0)={u}_0^\varepsilon(x), \ x_L<x<x_R, \ t>0. \end{align} $

(1.8)

$ \begin{align} {u}^\varepsilon(x_L+x, t)={u}^\varepsilon(x_R+x, t), \ {u}_x^\varepsilon(x_L+x, t)={u}_x^\varepsilon(x_R+x, t), \ t>0. \end{align} $

(1.9)

本文作如下的安排：第2节提出了三个新的格式, 第3节研究了格式的守恒律, 第4节研究了$ \beta=0 $时的平面波解, 第5节给出了稳定性条件, 第6节给出了模型问题, 第7节给出了数值实验.

下面对方程(1.7)–(1.9)构造几个差分格式, 首先引入如下的符号:

$ x_j=x_L+jh, \ t^n=n\tau, \ 0\leq\ j \leq J=[\frac{x_R-x_L}{h}], \ n=0, 1, 2, ..., [\frac{T}{\tau}], $

$ u^\varepsilon(j, n)\equiv u^\varepsilon(x_j, t^n), \quad u_j^{\varepsilon, n}\sim u^\varepsilon(j, n), \quad V_j=V(x_j), $

$ (w_j^n)_x=\frac{w_{j+1}^n-w_j^n}{h}, \quad\quad(w_j^n)\bar{x}=\frac{w_j^n-w_{j-1}^n}{h}, $

$ (w_j^n)_t=\frac{w_j^{n+1}-w_j^n}{\tau}, \quad\quad r\equiv \frac {\tau}{h^2}, $

$ \|w^n\|_2^2=h\sum\limits_{j=1}^J|w_j^n|^2, \quad\quad\|w^n\|_\infty=\sup\limits_{1\leq j\leq J}|w_j^n|. $

其中C表示一个常数. 从$ t^n $到$ t^{n+1} $, 方程(1.7) 可分为两步求解^[1], 首先求解

$ \begin{align} {i}\varepsilon{u}_t^\varepsilon-V{u}^\varepsilon-f(|{u}^\varepsilon|^2){u}^\varepsilon-\varepsilon\beta\arg({u}^\varepsilon){u}^\varepsilon=0, \end{align} $

(2.1)

接下来求解

$ \begin{align} {i}\varepsilon{u}_t^\varepsilon+\frac{\varepsilon^2}{2}{u}_{xx}^\varepsilon=0. \end{align} $

(2.2)

龙格-库塔方法是一种具有很高精度的算法, 我们用三点格式离散(2.2)中的导数项$ u_{xx}^\varepsilon $, 然后将离散后的方程用龙格-库塔方法求解, 格式如下

$ \begin{equation} u_j^{\varepsilon, *}=\left\{ \begin{array}{cc} \begin{aligned} e^{-\frac{i\tau}{\varepsilon}(V_j+f(|u_j^{\varepsilon, n}|^2))}u_j^{\varepsilon, n}, \beta=0 \\ \\ e^{-\frac{i}{\varepsilon\beta}(V_j+f(|u_j^{\varepsilon, n}|^2))(1-e^{-\tau\beta})+iS_j^{\varepsilon}(u_j^{\varepsilon, n}, \tau)}|u_j^{\varepsilon, n}|, \beta\neq0 \\ \end{aligned} \end{array} \right. \end{equation} $

(2.3)

$ \begin{align} u_j^{\varepsilon, n+1}=u_j^{\varepsilon, *}+(K1_j+2K2_j+2K3_j+K4_j)/6, \end{align} $

(2.4)

$ K1_j=G(u_j^{\varepsilon, *}), K2_j=G(u_j^{\varepsilon, *}+0.5K1_j), K3_j=G(u_j^{\varepsilon, *}+0.5K2_j), K4_j=G(u_j^{\varepsilon, *}+K3_j), $

$ \begin{align} G(\eta_j^n)=\frac{i\varepsilon\tau}{2h^2}(\eta_{j+1}^n-2\eta_j^n+\eta_{j-1}^n). \end{align} $

(2.5)

$ S^\varepsilon(u^\varepsilon, t) $如下

$ \begin{align} S^\varepsilon(u^\varepsilon, t)=e^{-\beta(t-t_n)}Im(\int_a^x\frac{u_v^\varepsilon(v, t_n)}{u^\varepsilon(v, t_n)}dv), \end{align} $

(2.6)

离散后为

$ \begin{align} S_j^\varepsilon(u^\varepsilon, t)=e^{-\beta t}\sum\limits_{l=1}^j\frac{h}{2}Im(\frac{D_x^su|_x=x_l-1}{u_l-1}+\frac{D_x^su|_x=x_l}{u_l}), \end{align} $

(2.7)

$ \begin{align} j=1, 2..., J, \ s_0(u, \tau)=0. \end{align} $

(2.8)

$ D_x^s $为导数项$ \partial_x $的谱格式:

$ \begin{align} D_x^su|_{x=x_j}=\frac{1}{J}\sum\limits_{l=-J/2}^{J/2-1}i\mu_l\hat u_le^{i\mu_l(x_j-a)}, \end{align} $

(2.9)

其中

$ \begin{align} \hat u_l=\sum\limits_{j=0}^{J-1}u_je^{-i\mu_l(x_j-a)}, l=-\frac{J}{2}, ..., \frac{J}{2}-1, \end{align} $

(2.10)

格式(2.3)–(2.5)的截断误差为$ O(\tau^2+h^2) $.

类似(2.3)–(2.5)我们考虑如下的五点龙格-库塔格式,

$ \begin{equation} u_j^{\varepsilon, *}=\left\{ \begin{array}{cc} \begin{aligned} e^{-\frac{i\tau}{\varepsilon}(V_j+f(|u_j^{\varepsilon, n}|^2))}u_j^{\varepsilon, n}, \beta=0 \\ \\ e^{-\frac{i}{\varepsilon\beta}(V_j+f(|u_j^{\varepsilon, n}|^2))(1-e^{-\tau\beta})+iS_j^{\varepsilon}(u_j^{\varepsilon, n}, \tau)}|u_j^{\varepsilon, n}|, \beta\neq0 \\ \end{aligned} \end{array} \right. \end{equation} $

(2.11)

$ \begin{align} u_j^{\varepsilon, n+1}=u_j^{\varepsilon, *}+(K1_j+2K2_j+2K3_j+K4_j)/6, \end{align} $

(2.12)

$ K1_j=G(u_j^{\varepsilon, *}), K2_j=G(u_j^{\varepsilon, *}+0.5K1_j), K3_j=G(u_j^{\varepsilon, *}+0.5K2_j), K4_j=G(u_j^{\varepsilon, *}+K3_j), $

$ \begin{align} G(\eta_j^n)=\frac{i\varepsilon\tau}{2h^2}(-\eta_{j+2}^n/12+4\eta_{j+1}^n/3-5\eta_{j}^n/2+4\eta_{j-1}^n/3-\eta_{j-2}^n/12), \end{align} $

(2.13)

格式(2.11)–(2.13) 的精度为$ O(\tau^2+h^4) $.

方程(2.2)是一个线性方程可以用傅里叶方法精确求解, 我们构造下面的格式,

$ \begin{equation} u_j^{\varepsilon, *}=\left\{ \begin{array}{cc} \begin{aligned} e^{-\frac{i\tau}{\varepsilon}(V_j+f(|u_j^{\varepsilon, n}|^2))}u_j^{\varepsilon, n}, \beta=0 \\ \\ e^{-\frac{i}{\varepsilon\beta}(V_j+f(|u_j^{\varepsilon, n}|^2))(1-e^{-\tau\beta})+iS_j^{\varepsilon}(u_j^{\varepsilon, n}, \tau)}|u_j^{\varepsilon, n}|, \beta\neq0 \\ \end{aligned} \end{array} \right. \end{equation} $

(2.14)

$ \begin{align} u_j^{\varepsilon, n+1}=\frac{1}{J}\sum\limits_{l=-J/2}^{J/2-1}e^{-i\varepsilon\mu_l^2\tau/2}\hat u_l^{\varepsilon, *}e^{i\mu_l(x_j-x_L)}, \ \ j=0, 1, 2, ..., J-1, \end{align} $

(2.15)

$ \begin{align} \mu_l=\frac{2\pi l}{x_R-x_L}, \hat u_l^{\varepsilon, *}=\sum\limits_{j=0}^{J-1}u_j^{\varepsilon, *}e^{-i\mu_l(x_j-x_L)}, l=-\frac{J}{2}, ..., \frac{J}{2}-1, \end{align} $

(2.16)

格式(2.14)–(2.16)的精度为$ O(\tau^2+h^m) $, $ m $为方程(1.7)–(1.9)的光滑度.

作为比较我们给出文[1]中的两个格式.

(1) BSP1

$ \begin{align} u_j^{\varepsilon, *}=\frac{1}{J}\sum\limits_{l=-J/2}^{J/2-1}e^{-i\varepsilon\mu_l^2\tau/2}\hat u_l^{\varepsilon, n}e^{i\mu_l(x_j-x_L)}, j=0, 1, 2, ..., J-1, \end{align} $

(2.17)

$ \begin{equation} u_j^{\varepsilon, n+1}=\left\{ \begin{array}{cc} \begin{aligned} e^{-\frac{i\tau}{\varepsilon}(V_j+f(|u_j^{\varepsilon, *}|^2))}u_j^{\varepsilon, *}, \beta=0 \\ \\ e^{-\frac{i}{\varepsilon\beta}(V_j+f(|u_j^{\varepsilon, *}|^2))(1-e^{-\tau\beta})+iS_j^{\varepsilon}(u_j^{\varepsilon, *}, \tau)}|u_j^{\varepsilon, *}|, \beta\neq0 \\ \end{aligned} \end{array} \right. \end{equation} $

(2.18)

$ \begin{align} \mu_l=\frac{2\pi l}{x_R-x_L}, \hat u_l^{\varepsilon, n}=\sum\limits_{j=0}^{J-1}u_j^{\varepsilon, n}e^{-i\mu_l(x_j-x_L)}, l=-\frac{J}{2}, ..., \frac{J}{2}-1. \end{align} $

(2.19)

(2) BSP2

$ \begin{equation} u_j^{\varepsilon, *}=\left\{ \begin{array}{cc} \begin{aligned} e^{-\frac{i\tau}{2\varepsilon}(V_j+f(|u_j^{\varepsilon, n}|^2))}u_j^{\varepsilon, n}, \beta=0 \\ \\ e^{-\frac{i}{2\varepsilon\beta}(V_j+f(|u_j^{\varepsilon, n}|^2))(1-e^{-\tau\beta})+iS_j^{\varepsilon}(u_j^{\varepsilon, n}, \tau/2)}|u_j^{\varepsilon, n}|, \beta\neq0 \\ \end{aligned} \end{array} \right. \end{equation} $

(2.20)

$ \begin{align} u_j^{\varepsilon, **}=\frac{1}{J}\sum\limits_{l=-J/2}^{J/2-1}e^{-i\varepsilon\mu_l^2\tau/2}\hat u_l^{\varepsilon, *}e^{i\mu_l(x_j-x_L)}, \ \ j=0, 1, 2, ..., J-1, \end{align} $

(2.21)

$ \begin{equation} u_j^{\varepsilon, n+1}=\left\{ \begin{array}{cc} \begin{aligned} e^{-\frac{i\tau}{2\varepsilon}(V_j+f(|u_j^{\varepsilon, **}|^2))}u_j^{\varepsilon, **}, \beta=0 \\ \\ e^{-\frac{i}{2\varepsilon\beta}(V_j+f(|u_j^{\varepsilon, **}|^2))(1-e^{-\tau\beta})+iS_j^{\varepsilon}(u_j^{\varepsilon, **}, \tau/2)}|u_j^{\varepsilon, **}|, \beta\neq0 \\ \end{aligned} \end{array} \right. \end{equation} $

(2.22)

格式(2.17)–(2.19)和(2.20)–(2.22)的精度为$ O(\tau^2+h^m) $, $ m $为方程(1.7)–(1.9)的光滑度.

考虑如下的两个守恒律：

$ \begin{align} H_1^n=\|u^{\varepsilon, n}\|^2, \end{align} $

(3.1)

$ \begin{align} H_2^n=\frac{\varepsilon^2}{2}h\sum\limits_{j=0}^{j=J-1}|\frac{u_{j+1}^{\varepsilon, n}-u_j^{\varepsilon, n}}{h}|^2 +F(\|u^{\varepsilon, n}\|^2)+h\sum\limits_{j=0}^{J-1}V_j|u_j^{\varepsilon, n}|^2. \end{align} $

(3.2)

定理3.1   格式SP0, BSP1和BSP2满足守恒律$ \|u^{\varepsilon, n}\|^2=const. $

证   由(2.3)得$ |u_j^{\varepsilon, *}|=|u_j^{\varepsilon, n}|, $对SP0由巴什瓦定理得

$ \|u^{\varepsilon, n+1}\|^2=(x_R-x_L)\sum\limits_{l=-J/2}^{J/2-1}(\hat u_l^{\varepsilon, *})^2=h\sum\limits_{j=1}^{J}|u_j^{\varepsilon, *}|^2=h\sum\limits_{j=1}^{J}|u_j^{\varepsilon, n}|^2=\|u^{\varepsilon, n}\|^2, $

于是有

$ \|u^{\varepsilon, n+1}\|^2=...=\|u^{\varepsilon, 0}\|^2=\|u_0^{\varepsilon}\|^2=const, $

对于格式BSP1有

$ \|u^{\varepsilon, n+1}\|^2=\|u^{\varepsilon, *}\|^2=(x_R-x_L)\sum\limits_{l=-J/2}^{J/2-1}(\hat u_l^{\varepsilon, n})^2=h\sum\limits_{j=1}^{J}|u_j^{\varepsilon, n}|^2=\|u^{\varepsilon, n}\|^2, $

所以

$ \|u^{\varepsilon, n+1}\|^2=...=\|u^{\varepsilon, 0}\|^2=\|u_0^{\varepsilon}\|^2=const, $

与格式SP0和BSP1类似, 可以证明格式BSP2有

$ \|u^{\varepsilon, n+1}\|^2=\|u^{\varepsilon, **}\|^2=(x_R-x_L)\sum\limits_{l=-J/2}^{J/2-1}(\hat u_l^{\varepsilon, *})^2=h\sum\limits_{j=1}^{J}|u_j^{\varepsilon, *}|^2=h\sum\limits_{j=1}^{J}|u_j^{\varepsilon, n}|^2=\|u^{\varepsilon, n}\|^2, $

所以

$ \|u^{\varepsilon, n+1}\|^2=...=\|u^{\varepsilon, 0}\|^2=\|u_0^{\varepsilon}\|^2=const. $

尽管三点龙格-库塔格式(2.3)–(2.5)和五点龙格-库塔格式(2.11)–(2.15)不能够满足守恒律(3.1) 和(3.2), 在数值实验中我们可以看到, 它们也可以很好地满足两个守恒律.

考虑当$ \beta=0 $时方程(1.7)–(1.9)的平面波解, 令

$ \begin{align} V=d, u^\varepsilon=Ae^{i(k\pi x/p-\omega t)}, \end{align} $

(4.1)

其中d, A, k和$ \omega $为常数. 将(4.1) 代入到(1.7)中, 得如下等式

$ \begin{align} i\varepsilon(-i\omega)u^{\varepsilon}+\frac{\varepsilon^2}{2}(i\cdot\frac{k\pi}{p})^2u^\varepsilon-(f(A^2)+d)u^\varepsilon=0, \end{align} $

(4.2)

得

$ \omega=\frac{\frac{\pi^2\varepsilon^2k^2}{2p^2}+f(A^2)+d}{\varepsilon}. $

在下面的论述中我们考虑如下离散形式的平面波解

$ \begin{align} V=d, u^\varepsilon=Ae^{i(k\pi jh/p-\omega n\tau)}. \end{align} $

(4.3)

将(4.3) 代入到(2.3)–(2.5), 得

$ u_j^{\varepsilon, *}=Ae^{-\frac{i\tau}{\varepsilon}(d+f(A^2))}\cdot e^{i(k\pi jh/p-\omega n\tau)}, u_j^{\varepsilon, n+1}=u_j^{\varepsilon, *}+\frac{1}{6}(K1_j+2K2_j+2K3_j+K4_j), $

整理得,

$ u_j^{\varepsilon, n+1}=(1+q+\frac{q^2}{2}+\frac{q^3}{6}+\frac{q^4}{24})u_j^{\varepsilon, *}, $

其中

$ q=\frac{i\varepsilon\tau}{h^2}(\cos\frac{k\pi h}{p}-1)=\frac{i\varepsilon\tau}{h^2}(-\frac{k^2\pi^2h^2}{2p^2}+O(h^4))=-\frac{i\varepsilon k^2\pi^2\tau}{2p^2}+O(h^2), $

于是有

$ u_j^{\varepsilon, n+1}=u_j^{\varepsilon, *}(1-\frac{i\varepsilon k^2\pi^2\tau}{2p^2}+O(\tau^2+h^2)), $

因为$ u_j^{\varepsilon, n+1}=Ae^{k\pi jh/p-\omega(n+1)\tau}, $所以有

$ \omega=\frac{\frac{\pi^2\varepsilon^2k^2}{2p^2}+f(A^2)+d}{\varepsilon}+O(\tau^2+h^2), $

所以三点龙格-库塔方法的精度为$ O(\tau^2+h^2) $, 同理可得五点龙格-库塔格式(2.11)–(2.13)的精度为$ O(\tau^2+h^4) $.

将(4.3)代入到格式SP0中, 由离散的傅里叶变换,

$ u_j^{\varepsilon, *}=e^{-\frac{i\tau}{\varepsilon}(d+f(A^2))}\cdot Ae^{i(\frac{k\pi jh}{p}-\omega n\tau)}=Ae^{-\frac{i\tau}{\varepsilon}(d+f(A^2))+i\frac{k\pi jh}{p}-i\omega n\tau}, $

$ \begin{eqnarray} \hat u_j^{\varepsilon, *} &=& \sum\limits_{j=0}^{J-1}u_j^{\varepsilon, *}e^{-i\frac{l\pi jh}{p}}=Ae^{[-\frac{i\tau}{\varepsilon}(d+f(A^2))-i\omega n\tau]}\cdot \sum\limits_{j=0}^{J-1}e^{i(\frac{k\pi jh}{p}-\frac{l\pi jh}{p})}\\ &=& \left\{ \begin{array}{cc} \begin{aligned} 0, l\neq k \\ J\cdot Ae^{[-\frac{i\tau}{\varepsilon}(d+f(A^2))-i\omega n\tau]}, l=k\nonumber \end{aligned} \end{array} \right. \end{eqnarray} $

并且

$ u_j^{\varepsilon, n+1}=\frac{1}{J}\sum\limits_{l=-\frac{J}{2}}^{\frac{J}{2}-1}e^{-\frac{i\varepsilon \pi^2l^2\tau}{2p^2}}\cdot \hat u_j^{\varepsilon, *}\cdot e^{i\frac{\pi ljh}{p}}=Ae^{-\frac{i\tau}{\varepsilon}(d+f(A^2))-i\omega n\tau}\cdot e^{-\frac{i\varepsilon \pi^2k^2\tau}{2p^2}}\cdot e^{i\frac{\pi kjh}{p}}, $

又因为

$ u_j^{\varepsilon, n+1}=Ae^{i(\frac{k\pi jh}{p}-\omega(n+1)\tau)}, $

整理得

$ \omega=\frac{\frac{\pi^2\varepsilon^2k^2}{2p^2}+f(A^2)+d}{\varepsilon}, $

所以分裂谱格式可以精确满足平面波解, 同理BSP1和BSP2也精确满足平面波解.

关于格式的稳定性有如下的定理:

定理5.1   三点龙格-库塔格式(2.3)–(2.5)、五点龙格-库塔格式(2.11)–(2.13)的稳定性条件分别为$ r\leq \frac{\sqrt{2}}{\varepsilon} $, $ r\leq\frac{12\sqrt{2}}{7\varepsilon} $, 其中$ r=\tau/h^2 $.

证明   假设$ V(x)=d $, 其中$ d $是个常数. 下面用线性化的方法证明三点龙格-库塔格式和五点龙格-库塔格式的稳定性. 考虑如下的解的形式

$ \begin{align} u_j^{\varepsilon, n}=\xi^{\varepsilon, n}e^{ikjh}, \end{align} $

(5.1)

将(5.1)式代入到(2.3)–(2.5), 我们得

$ u_j^{\varepsilon, *}=e^{-\frac{i\tau}{\varepsilon}(d+f(|\xi^n|^2))}\xi^ne^{ikjh}, $

$ u_j^{\varepsilon, n+1}=u_j^{\varepsilon, *}+\frac{1}{6}(K1_j+2K2_j+2K3_j+K4_j), $

$ \begin{eqnarray} K1_j&=&G(u_j^{\varepsilon, *})=\frac{i\varepsilon\tau}{2h^2}\cdot(e^{ikh}-2+e^{-ikh})u_j^{\varepsilon, *}\\ &=&i\varepsilon r(coskh-1)u_j^{\varepsilon, *}\equiv qu_j^{\varepsilon, *}, \end{eqnarray} $

其中

$ q=i\varepsilon r(coskh-1), $

$ K2_j=G((1+\frac{q}{2})u_j^{\varepsilon, *})=q(1+\frac{q}{2})u_j^{\varepsilon, *}, $

$ K3_j=q(1+\frac{q}{2}(1+\frac{q}{2}))u_j^{\varepsilon, *}, $

$ K4_j=q(1+q(1+\frac{q}{2}(1+\frac{q}{2})))u_j^{\varepsilon, *}, $

因为

$ u_j^{\varepsilon, n+1}=(1+q+\frac{q^2}{2}+\frac{q^3}{6}+\frac{q^4}{24})u_j^{\varepsilon, *}, $

又因为

$ u_j^{\varepsilon, n+1}=\xi^{\varepsilon, n+1}e^{ikjh}, $

所以有

$ \xi^{\varepsilon, n+1}=(1+q+\frac{q^2}{2}+\frac{q^3}{6}+\frac{q^4}{24})e^{-\frac{i\tau}{\varepsilon}(d+f(|\xi^n|^2))}\xi^{\varepsilon, n}. $

令

$ \begin{align} G(\tau, k)=(1+q+\frac{q^2}{2}+\frac{q^3}{6}+\frac{q^4}{24})e^{-\frac{i\tau}{\varepsilon}(d+f(|\xi^n|^2))}, p=\varepsilon r(coskh-1), \end{align} $

(5.2)

则

$ \begin{align} G(\tau, k)=(1+ip-\frac{p^2}{2}-i\frac{p^3}{6}+\frac{p^4}{24})e^{-\frac{i\tau}{\varepsilon}(d+f(|\xi^n|^2))}, \end{align} $

(5.3)

格式(2.3)–(2.5)是稳定的, 则$ G(\tau, k) $满足

$ |G(\tau, k)|^2=(1-\frac{p^2}{2}+\frac{p^4}{24})^2+(p-\frac{p^3}{6})^2\leq 1, $

于是有

$ \varepsilon^2r^2\sin^4{\frac{kh}{2}}-2\leq 0, $

三点龙格-库塔格式(2.3)–(2.5)的稳定性条件为$ r\leq \frac{\sqrt{2}}{\varepsilon}. $同理可证, 五点龙格-库塔格式的稳定性条件为$ r\leq \frac{12\sqrt{2}}{7\varepsilon}. $另外, 由定理3.1可知谱格式SP0, BSP1和BSP2是无条件稳定的^[1].

令

$ u_0^{\varepsilon}(x)\equiv A_0(x)e^{\frac{i}{\varepsilon}s_0(x)}. $

模型Ⅰ   平面波解的初值问题：当$ \beta=0 $时, 我们给出了(1.7)–(1.9)的平面波解, 令初值$ A_0(x)=A, S_0(x)=\frac{\varepsilon k\pi x}{p}, V_0^{\varepsilon}(x)=d. $精确解为

$ u^{\varepsilon}(x, t)=Ae^{i(\frac{k\pi x}{p}-\omega t)}, V_0^{\varepsilon}(x)=d, $

并且

$ \omega=\frac{\frac{\pi^2\varepsilon^2k^2}{2p^2}+f(A^2)+d}{\varepsilon}, $

在计算中取如下参数: $ A=5, k=1, p=10, d=1. $

模型Ⅱ   当$ \beta\neq 0 $时, 我们用第2节中的格式来计算密度函数$ \rho^{\varepsilon}(x, t)=|u^{\varepsilon}(x, t)|^2 $. 初值条件为$ A_0(x)=e^{-x^2} $, $ S_0(x)=\frac{1}{e^{3x}+e^{-3x}} $, $ \beta=1.0 $, $ V(x)=0 $, $ f(\rho)=-\rho $.

本节用第2节的数值格式来计算模型Ⅰ和模型Ⅱ. 在计算模型Ⅰ时, 固定精度($ L_\infty $), 让步长$ h $和$ \tau $自由. 对每个格式, 令步长$ h $和$ \tau $尽量大, 既要达到给定的精度, 又能使格式保持稳定, 比较格式的计算时间. 计算的量如下:

守恒律

$ H_1^n=\|u^{\varepsilon, n}\|^2, $

和

$ H_2^n=\frac{\varepsilon^2}{2}h\sum\limits_{j=0}^{j=J-1}|\frac{u_{j+1}^{\varepsilon, n}-u_j^{\varepsilon, n}}{h}|^2 +F(\|u^{\varepsilon, n}\|^2)+h\sum\limits_{j=0}^{J-1}V_j|u_j^{\varepsilon, n}|^2, $

误差为

$ E_1=\frac{H_1^n-H_1^0}{H_1^0}, E_2=\frac{H_2^n-H_2^0}{H_2^0}, $

精确解的误差

$ L_\infty=\max\limits_{j}|u_j^{\varepsilon, n}-u^{\varepsilon}(x_j, t^n)|, $

其中$ u^{\varepsilon}(x_j, t^n) $为在点$ (x_j, t_n) $处的精确解. 模型Ⅰ的计算结果见表 1.

文[4][5]中对小$ \varepsilon $的线性schrödinger方程的有限差分进行了研究, 结果表明为了保证良好的近似, 时间步长$ \tau $和空间步长$ h $满足如下关系：

$ h=o(\varepsilon), \tau=o(\varepsilon). $

对于模型Ⅱ, 我们取$ \varepsilon=0.1 $, 为了画出密度函数的较精确图像, 令空间步长$ h=0.015625 $和一个较小的时间步长$ \tau=0.002 $. 图 1和图 2给出了用五点龙格-库塔格式和Bsp2计算的结果, 表 2给出了所有格式的计算时间.

由表 1可知, 当$ \varepsilon $=1, 0.1, 0.01, 0.001, 0, 0001时, 所有的格式都能够很好的计算平面波解. 其中, 分裂的五点龙格-库塔格式和三点龙格-库塔格式用的时间较少. 表 2中数据也表明五点龙格-库塔格式和三点龙格-库塔格式用的时间较少, 具有较高的计算效率. 由计算结果可知, 所有的格式都可以很好的计算, 其中分裂的五点龙格-库塔格式精度高, 计算时间少, 是最有效的和可靠的.