ABSTRACT

In this paper, we have introduced fourth-order compact finite difference, the time splitting method and the Crank-Nicolson method to solve the nonlinear Schrödinger-Poisson equations. Based on fast Sine transform, we construct a fast solver for the fully discretized system. The presented numerical algorithm has been used to solve one-dimensional, two-dimensional and three-dimensional nonlinear Schrödinger-Poisson equations. We provide specific numerical examples. Through the MATLAB software, we write matlab programs based on the presented numerical algorithm, calculate approximated error and draw the approximated numerical solution. The numerical results prove that the presented algorithm has spectral accuracy in space direction. They also confirm its efficiency and stability.

Keywords:Nonlinear Schrodinger-Poisson Equations, High Order Compact Difference Scheme, Sine Transform, Time Splitting Method

求解薛定谔–泊松方程组的时间分裂紧致 差分格式

姜珊，刘荣华，马秀，王汉权

云南财经大学统计与数学学院，云南 昆明

收稿日期：2018年12月11日；录用日期：2019年1月2日；发布日期：2019年1月9日

摘 要

本文通过应用高阶紧致差分格式，时间分裂法与Crank-Nicolson法等方法求解非线性薛定谔–泊松方程组。基于快速Sine变换，我们设计了一种求解离散系统的快速算法。我们用该算法分别求解一维、二维、三维的非线性薛定谔–泊松方程组。我们列举具体的数值例子，使用MATLAB软件编写出算法程序，并且通过程序计算近似误差与画出近似解的图像。数值计算结果证实了该算法在空间方向具有谱精度，也证实了该算法的高效性与稳定性。

关键词 :非线性薛定谔–泊松方程组，四阶紧致差分，Sine变换，时间分裂法

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

带电粒子的运动规律通常用耦合的动力学方程(例如弗拉索夫方程或波尔兹曼方程)与泊松方程进行描述 [1] [2] [3] [4] 。在考虑粒子的量子隧道效应情形下，带电粒子的运动规律也可用非线性薛定谔泊松方程组来描述，这主要是因为耦合的动力学-泊松方程组可近似化简为非线性薛定谔泊松方程组 [3] [4] 。非线性薛定谔泊松方程组最早由Hatree得到。他在研究原子物理学中的电子对原子核的库伦势的自洽作用(self-consistenteffect)时，发现带电粒子的运动可由非线性薛定谔泊松方程组 [5] 来描述。后来此方程组被广泛应用到量子等离子体 [3] 、量子输运 [6] 、量子半导体 [7] 、分子链 [8] 等物理现象中。本文主要想通过非线性薛定谔泊松方程组的数值解来模拟相关量子物理现象(例如纳米级带电粒子的运动规律)。

据我们所知，前人提出了多种数值方法来求解非线性薛定谔泊松方程组。为了模拟半导体纳米结构，Chen等人提出了一种新奇的快速谱元法 [9] ，他们在求解该方程组的过程中，未知函数的空间变量由高阶勒让德多项式展开；为了计算球对称的薛定谔泊松方程组的基态和动力学现象，Dong研究了一种简化伪谱方法 [10] ；为研究一种非磁化作用的稠密热等离子的物理现象，Shaikh等人针对薛定谔泊松方程组提出了一种傅里叶谱方法，并进行了三维非线性流体在纳米级尺度下湍流的模拟 [11] ；Mauser等人为薛定谔泊松方程组提出了一个向后欧拉法和半隐式/蛙跳正弦伪谱方法，并用该方法计算相关粒子的基态和动力学 [12] ；Zhang等人为该方程组提出了一种紧致差分格式，并研究了该数值方法的最优误差估计；我们还注意到前人研究了球对称形式的薛定谔泊松方程组的稳态解 [13] 。

关于非线性薛定谔泊松方程组的数学理论研究，前人已经开展了许多工作：Ben Abdallahd等人研究了薛定谔泊松方程组的稳态解与动态解的弱解存在性 [2] ；Brezzi等人基于非线性薛定谔泊松方程组的解证明了量子Vlasov-Poisson方程组的全局经典解的存在性 [14] ；Castella等人研究了薛定谔泊松方程组在L²范数意义下解的存在性、唯一性 [15] ；Lange等人显示了Wigner-Poisson方程组与薛定谔泊松方程组的解在不同边界条件(周期性边界条件、混合边界条件)下的存在性、唯一性与正则性 [6] 。

本文将为非线性薛定谔泊松方程组的求解提出一种时间分裂紧致差分格式，重点讨论如何基于Sine变换实现该方法的快速求解过程，并简单讨论该方法的数值特点。本文的结构安排如下：在第1.1节至第1.3节中，对薛定谔泊松方程组的物理模型做了简单介绍，同时也对紧致差分格式的构造、离散Sine变换做了相关阐述。在第2节中，详细推导了用时间分裂紧致差分格式来求解一维、二维、三维的薛定谔泊松方程组的数值计算过程，重点讨论如何基于Sine变换实现离散系统的快速求解过程。在第3节中，通过几个具体的数值例子的计算结果，证实了在前两节中介绍的方法可实现一维、二维、三维薛定谔泊松方程组方程组的数值计算，并验证了数值方法的守恒性。

1.1. 1.1. 问题背景

多粒子量子系统是由多个粒子通过相互作用形成的微观系统，它广泛存在于半导体物理、等离子物理、凝聚态物理和分子动力学等领域，对它的科学研究具有重要意义。本文研究的薛定谔–泊松方程组是对多粒子非相对论量子系统的一种单粒子波函数近似。薛定谔–泊松方程组是对含时Hartree-Fock方程组的一种单粒子波函数局部近似，它的无量纲形式为 [16] [17] [18] [19] [20] [21] [22]

$i{\partial }_t\psi \left(x,t\right)=-\frac{1}{2}\Delta \psi +V\left(x\right)\psi +\varphi \left(x,t\right)\psi +\beta {\left|\psi \right|}^{\frac{4}{d}}\psi ,x\in R^d,t>0,$ (1.1)

$\Delta \varphi \left(x,t\right)=-{\left|\psi \right|}^2,x\in R^d,$ (1.2)

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

其中d (=1, 2, 3)为空间维数，β为一已知常数，波函数 $\psi \left(x,t\right)$ 是复值函数，且满足无穷远条件

$\underset{\left|x\right|\to \infty }{lim}\left|\psi \left(x,t\right)\right|=\mathrm{0,}$ (1.3)

$V=V\left(x\right)$ 为外势函数，形式已知。ϕ为泊松位势，通常写成卷积形式。

薛定谔–泊松方程组(1.1)~(1.2)本身有许多守恒量，例如：粒子数N(t)，粒子冲量P(t)，粒子角动量L(t)与粒子总能量E(t)均不随时间变化而变化，其中，

粒子数

$N\left(t\right)={\displaystyle {\int }_{{ℝ}^d}}{\left|\psi \left(x\mathrm{,}t\right)\right|}^2\text{d}x\mathrm{,}$ (1.4)

粒子冲量

$P\left(t\right)=-i{\displaystyle {\int }_{{ℝ}^d}}{\psi }^{\mathrm{<mo>\; *\; </mo>}}\left(x\mathrm{,}t\right)\nabla \psi \left(x\mathrm{,}t\right)\text{d}x\mathrm{,}$ (1.5)

粒子旋转量

$L\left(t\right)=-i{\displaystyle {\int }_{{ℝ}^d}}{\psi }^{*}\left(x\mathrm{,}t\right)x\times \nabla \psi \left(x\mathrm{,}t\right)\text{d}x,$ (1.6)

与粒子总能量

$E\left(t\right)={\displaystyle {\int }_{{ℝ}^d}}\left[\frac{1}{2}{\left|\nabla \psi \right|}^2+V{\left|\psi \right|}^2+\frac{{\left|\nabla \varphi \right|}^2}{2}+\frac{d}{2+d}{\left|\psi \right|}^{2+4/d}\right]\text{d}x.$ (1.7)

1.2. 紧致差分格式设计

有限差分法是计算机数值模拟最早采用的方法，至今仍被广泛使用，其主要思想就是在空间和时间两个方向上将问题离散化。传统的差分格式是显式地通过各点函数值的线性组合来表示函数的导函数，而紧致差分格式是隐式地利用网格点的函数值的线性组合来表示导函数的线性组合，需要通过求解方程组来得到每个网格点的导数值。这两种方法相比较，当达到相同精度时，紧致差分格式需要的网格点少，计算精度高。

紧致差分格式的构造最早由 紧致差分格式的构造最早由Lele提出 [23] ，在文 [24] 中也有介绍。在本节中，我们简要地叙述紧致差分格式的构成。通常把未知函数u(x)的定义区间[a, b]作M等分， $h=\left(b-a\right)/M$ 称为步长，

$x_i=a+ih\left(0<i<M\right)$ 称为节点，记 $u_i:=u\left(x_i\right)$ ， ${u^{\prime }}_i:=\left(\frac{\text{d}u}{\text{d}x}\right)\left(x_i\right)$ ， ${{u^{\prime }}^{\prime }}_i:=\left(\frac{{\text{d}}^{2}u}{\text{d}{x}^2}\right)\left(x_i\right)$ ， ${{{u^{\prime }}^{\prime }}^{\prime }}_i:=\left(\frac{{\text{d}}^{3}u}{\text{d}{x}^3}\right)\left(x_i\right)$ ， ${{{{u^{\prime }}^{\prime }}^{\prime }}^{\prime }}_i:=\left(\frac{{\text{d}}^{4}u}{\text{d}{x}^4}\right)\left(x_i\right)$ ，其中i与M是正整数，a与b是给定的区间端点。

紧致差分格式想从下式中隐式地得到函数u(x)在节点x_i处的近似值 ${{u^{\prime }}^{\prime }}_i$ ( $i=0,1,\cdots ,M$ )：

$\begin{array}{l}\beta {{u^{\prime }}^{\prime }}_{i-2}+\alpha {{u^{\prime }}^{\prime }}_{i-1}+{u^{\prime }}_i+\alpha {{u^{\prime }}^{\prime }}_{i+1}+\beta {{u^{\prime }}^{\prime }}_{i+2}\\ =c\frac{u_{i+3}-2u_i+u_{i-3}}{9h^2}+b\frac{u_{i+2}-2u_i+u_{i-2}}{4h^2}+a\frac{u_{i+1}-2u_i+u_{i-1}}{h^2},\end{array}$ (1.8)

其中α，β，a，b，c是一些待定系数。由泰勒展开式可知由下式得到的近似值 ${{u^{\prime }}^{\prime }}_i$ ( $i=0,1,\cdots ,M$ )具有四阶精度

$\frac{1}{10}{{u^{\prime }}^{\prime }}_{i-1}+{{u^{\prime }}^{\prime }}_i+\frac{1}{10}{{u^{\prime }}^{\prime }}_{i+1}=\frac{6}{5}\cdot \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}.$ (1.9)

上式可称为一种四阶紧致差分格式。

另外由泰勒展开式可知由下式得到的近似值 ${{u^{\prime }}^{\prime }}_i$ ( $i=0,1,\cdots ,M$ )具有六阶精度

$2{{u^{\prime }}^{\prime }}_{i-1}+11{{u^{\prime }}^{\prime }}_i+2{{u^{\prime }}^{\prime }}_{i+1}=\frac{3}{4}\cdot \frac{u_{i+2}-2u_i+u_{i-2}}{h^2}+12\cdot \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$ (1.10)

上式可称为一种六阶紧致差分格式。

1.3. 一维的离散Sine变换与逆变换

为了实现高效的数值算法，首先看一维的离散Sine变换：

$u_i={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\hat{u}}_k\sin \left(\frac{ik\text{π}}{M}\right),i=1,2,\cdots ,M-1,$ (1.11)

一维的离散Sine逆变换：

${\hat{u}}_k=\frac{2}{M}{\displaystyle \underset{i=1}{\overset{M-1}{\sum }}}{u}_i\sin \left(\frac{ik\text{π}}{M}\right),k=1,2,\cdots ,M-1,$ (1.12)

其中，记 $u\left(x_i\right)=u_i$

通过方程(1.11)，我们可写出 $u_{i+1}$ ， $u_{i-1}$ ， ${{u^{\prime }}^{\prime }}_i$ ， ${{u^{\prime }}^{\prime }}_{i+1}$ ， ${{u^{\prime }}^{\prime }}_{i-1}$ 的Sine变换：

$u_{i+1}={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\hat{u}}_k\sin \left(\frac{\left(i+1\right)k\text{π}}{M}\right),u_{i-1}={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\hat{u}}_k\sin \left(\frac{\left(i-1\right)k\text{π}}{M}\right),$

${{u^{\prime }}^{\prime }}_i={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{{{\hat{u}}^{\prime }}^{\prime }}_k\sin \left(\frac{ik\text{π}}{M}\right),$ (1.13)

${{u^{\prime }}^{\prime }}_{i+1}={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{{{\hat{u}}^{\prime }}^{\prime }}_k\sin \left(\frac{\left(i+1\right)k\text{π}}{M}\right),{{u^{\prime }}^{\prime }}_{i-1}={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{{{\hat{u}}^{\prime }}^{\prime }}_k\sin \left(\frac{\left(i-1\right)k\text{π}}{M}\right).$

把上述等式带入与四阶紧致差分格式相关的方程(1.9)，得

$\begin{array}{l}{\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{{{\hat{u}}^{\prime }}^{\prime }}_k\left\{\frac{1}{10}\sin \left(\frac{\left(i-1\right)k\text{π}}{M}\right)+\sin \left(\frac{ik\text{π}}{M}\right)+\frac{1}{10}\sin \left(\frac{\left(i+1\right)k\text{π}}{M}\right)\right\}\\ =\frac{6}{5h^2}{\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\hat{u}}_k\left\{\sin \left(\frac{\left(i-1\right)k\text{π}}{M}\right)+\sin \left(\frac{\left(i+1\right)k\text{π}}{M}\right)-2\sin \left(\frac{ik\text{π}}{M}\right)\right\}.\end{array}$

利用三角公式化简为：

${\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{{{\hat{u}}^{\prime }}^{\prime }}_k\left(\frac{1}{5}\cos \left(\frac{k\text{π}}{M}\right)+1\right)\sin \left(\frac{ik\text{π}}{M}\right)=\frac{6}{5h^2}{\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\hat{u}}_k\left(2\cos \left(\frac{k\text{π}}{M}\right)-2\right)\sin \left(\frac{ik\text{π}}{M}\right).$

根据对应项分别相等, 得

${\hat{u}}_k=-{{{\hat{u}}^{\prime }}^{\prime }}_k{\left(\frac{24si{n}^2\left(\frac{k\text{π}}{2M}\right)}{h^2}\right)}^{-1}\left(cos\left(\frac{k\text{π}}{M}\right)+5\right)\mathrm{,}k=\mathrm{1,2,}\cdots \mathrm{,}M-1.$ (1.14)

方程(1.14)显示了 ${\hat{u}}_k$ 与 ${{{\hat{u}}^{\prime }}^{\prime }}_k$ 之间的联系。而式子(1.11)与(1.13)建立了函数u(x)在节点x_i处的值与二阶导数 ${u^{\prime }}^{\prime }\left(x\right)$ 在节点x_i处的值之间的联系。

重要的是，上述离散Sine变换(1.11)与离散Sine逆变换(1.12)可借助快速的离散傅里叶变换实现，快速的离散傅里叶变换将离散Sine变换(1.11)与逆变换(1.12)的计算量从 $O\left(M^2\right)$ 减少至 $O\left(M\log \left(M\right)\right)$ 。

同理，如果把相关等式代入与六阶紧致差分格式相关的方程(1.10)，那么也可得与(1.14)类似的等式：

${\hat{u}}_k=-{{{\hat{u}}^{\prime }}^{\prime }}_k\frac{\frac{24}{h^2}\left(\cos \left(k\text{π}/M\right)-1\right)+\frac{6}{4h^2}\left(\cos \left(2k\text{π}/M\right)-1\right)}{4\cos \left(j\text{π}/M\right)+11},k=1,2,\cdots ,M-1.$ (1.15)

方程(1.15)也显示了 ${\hat{u}}_k$ 与 ${{{\hat{u}}^{\prime }}^{\prime }}_k$ 之间的联系。

2. 快速高效的数值方法

在本节中，我们先分别介绍如何利用时间分裂紧致差分格式法离散一维非线性薛定谔泊松方程组、二维非线性薛定谔泊松方程组、三维非线性薛定谔泊松方程组，然后介绍如何基于一维的快速Sine变换、二维的快速Sine变换、三维的快速Sine变换来求解上述方法离散后的方程组，最后写出详细的离散求解过程。

2.1. 一维非线性薛定谔泊松方程组的求解方法

考虑如下一维非线性薛定谔泊松方程组初边值问题：

$i{\psi }_t=-\frac{1}{2}{\psi }_{xx}+V\left(x\right)\psi +\varphi \psi +\beta {\left|\psi \right|}^2\psi ,a<x<b,0<t<T,$ (2.16)

$-{\varphi }_{xx}={\left|\psi \right|}^2,$

$\psi \left(a,t\right)=0,\psi \left(b,t\right)=0,$ (2.17)

$\psi \left(x,0\right)={\varphi }_0\left(x\right),$

其中 $\Delta ={\partial }_{xx}$ ， $a<x<b$ 。 $a\mathrm{,}b\mathrm{,}\beta \mathrm{,}T$ 为常数， $V\left(x\right)\mathrm{,}{\varphi }_0\left(x\right)$ 为已知函数， $\psi \left(x\mathrm{,}t\right)\mathrm{,}\varphi \left(x\mathrm{,}t\right)$ 为未知函数。

先将[a,b]等分成M份，则空间步长为 $h_x=\frac{b-a}{M}$ ，节点为 $x_j=a+j{h}_x,j=0,1,\cdots ,M$ ；再将[0,T]等分成N份，则时间步长为 $\tau =\frac{T}{L}$ ，节点为 $t_n=n\tau ,\left(n=0,1,\cdots ,L\right)$ ，记 ${\psi }_j^n\mathrm{,}{v}_j\mathrm{,}{\varphi }_j^n$ 分别为 $\psi \left(x_j\mathrm{,}{t}_n\right)\mathrm{,}v\left(x_j\right)\mathrm{,}\varphi \left(x_j\mathrm{,}{t}_n\right)$ 的

近似值。

通过观察会发现：方程(2.16)类似于非线性薛定谔方程，(2.17)是一个泊松方程。对于(2.16)方程，在一个时间步长 $\left[t_n,t_n{}_{+1}\right]$ 上，我们利用时间分裂法将该方程分两步求解：

先求解非线性部分 $i{\psi }_t=V\left(x\right)\psi +\beta {\left|\psi \right|}^2\psi +\varphi \psi$ 。可巧妙地利用共轭与积分性质：等式两边同时乘ψ的共轭 $\bar{\psi }$ ，得： $i{\psi }_t\bar{\psi }=V\left(x\right){\left|\psi \left(x,t\right)\right|}^2+\beta {\left|\psi \left(x,t\right)\right|}^4+\varphi {\left|\psi \left(x,t\right)\right|}^2$ ，再对它两边同时取共轭，得：

$-i{\bar{\psi }}_t\psi =V\left(x\right){\left|\psi \left(x,t\right)\right|}^2+\beta {\left|\psi \left(x,t\right)\right|}^4+\varphi {\left|\psi \left(x,t\right)\right|}^2,$ (2.18)

注意到ϕ是一个实值函数，将以上两式相减得： $i\left(\psi {\bar{\psi }}_t+\bar{\psi }{\psi }_t\right)=0$ ，也即 ${\left({\left|\psi \right|}^2\right)}_t=0$ 。我们可以得到 ${\left|\psi \left(x,t\right)\right|}^2={\left|\psi \left(x,t_n\right)\right|}^2\left(t_n\le t\le t_{n+1}\right)$ 。再从泊松方程 $-\Delta \varphi ={\left|\psi \right|}^2$ 可知 $\varphi \left(x,t\right)=\varphi \left(x,t_n\right)$ $\left(t_n\le t\le t_{n+1}\right)$ 。因此，非

线性部分方程的求解变为：

$i{\psi }_t\left(x,t\right)=V\left(x\right)\psi \left(x,t\right)+\beta {\left|\psi \left(x,t_n\right)\right|}^2\psi \left(x,t\right)+\varphi \left(x,t_n\right)\psi \left(x,t\right),$ (2.19)

$i\frac{{\psi }_t}{\psi }=V\left(x\right)+\beta {\left|\psi \left(x,t_n\right)\right|}^2+\varphi \left(x,t_n\right),$ (2.20)

$\psi \left(x,t_{n+1}\right)=\psi \left(x,t_n\right){\text{e}}^{-i\tau \left(V\left(x\right)+\beta {\left|\psi \left(x,t_n\right)\right|}^2+\varphi \left(x,t_n\right)\right)},$ (2.21)

记 $\psi \left(x,t_{*}\right)=\psi \left(x,t_{n+1}\right)$ ， $\psi \left(x,t_{*}\right)$ 作为时间分裂法的中间变量，但在方程(2.21)中，仍有未知量 $\varphi \left(x\mathrm{,}{t}_n\right)$ 。由已知条件 $-\Delta \varphi ={\left|\psi \right|}^2$ ，即 $-{\varphi }_{xx}\left(x_j,t_n\right)=f_j$ ，其中记 $f_j={\left|\psi \left(x_j\mathrm{,}{t}_n\right)\right|}^2$ 。可先利用四阶紧致差分格式离散一维泊松方程 $-\Delta \varphi ={\left|\psi \right|}^2$ ，然后对离散后的方程组做Sine变换得 $-{\hat{\varphi }}_{xx}={\hat{f}}_j$ (详细过程见 [25] )。该式子右边 ${\hat{f}}_j$

是由f_j做Sine变换得到；同时由基于Sine变换的四阶紧致差分格式，得到关系式：

${\hat{\varphi }}_j^{xx}=-{\hat{\varphi }}_j\left(\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2}\right){\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}^{-1}\mathrm{,}$

$-{\hat{\varphi }}_j^{xx}={\hat{f}}_j\mathrm{,}$

${\hat{\varphi }}_j={\hat{f}}_j{\left(\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2}\right)}^{-1}\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)\mathrm{,}$ (2.22)

其中 $1\le j\le M-1$ ，这样就可以通过式子(2.22)求出 ${\hat{\varphi }}_j$ 。再利用一次Sine变换这样就求出 ${\varphi }_j^n$ 的值，将其带入到方程(2.21)，就可以求出中间变量 $\psi \left(x\mathrm{,}{t}_{\mathrm{<mo>\; *\; </mo>}}\right)$ 。

再考虑线性部分： $i{\psi }_t=-\frac{1}{2}\Delta \psi$ 的求解。此时以 ${\psi }_j^{\mathrm{<mo>\; *\; </mo>}}$ 做为初始条件，对于该方程时间方向上采用

Crank-Nicolson公式离散：

$i\frac{{\psi }_j^{n+1}-{\psi }_j^{\mathrm{<mo>\; *\; </mo>}}}{\tau }=-\frac{1}{4}\left({\left({\psi }_j^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}+{\left({\psi }_j^{n+1}\right)}_{xx}\right)\mathrm{,}$

$\left(1-\frac{i\tau }{4}\Delta \right){\psi }_j^{n+1}=\left(1+\frac{i\tau }{4}\Delta \right){\psi }_j^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}$

${\psi }_j^{n+1}-\frac{i\tau }{4}{\left({\psi }_j^{n+1}\right)}_{xx}={\psi }_j^{*}+\frac{i\tau }{4}{\left({\psi }_j^{*}\right)}_{xx},$ (2.23)

对上述方程，空间上采用四阶紧致差分格式来离散。由于等式右边 ${\psi }_j^{\mathrm{<mo>\; *\; </mo>}}$ 已知，对上述离散后的新系统做Sine逆变换得 ${\hat{\psi }}_j^{\mathrm{<mo>\; *\; </mo>}}$ 。再结合关系式：

${\left({\hat{\psi }}_j^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}=-{\hat{\psi }}_j^{\mathrm{<mo>\; *\; </mo>}}\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}\mathrm{,}$ (2.24)

得 ${\left({\hat{\psi }}_j^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}$ 。再做Sine变换就得 ${\left({\psi }_j^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}$ 。这样方程(2.23)的右边函数值全部已知，记 $g_j={\psi }_j^{*}+\frac{i\tau }{4}{\left({\psi }_j^{*}\right)}_{xx}$ ，对(2.23)两边同时做Sine变换 ${\hat{\psi }}_j^{n+1}-\frac{i\tau }{4}{\left({\hat{\psi }}_j^{n+1}\right)}_{xx}={\hat{g}}_j$ ，其中 ${\hat{g}}_j$ 是 $g_j$ 经过一维的Sine逆变换得到的。对方程

(2.23)左边，我们利用基于Sine变换的四阶紧致差分格式的关系式：

${\left({\hat{\psi }}_j^{n+1}\right)}_{xx}=-{\hat{\psi }}_j^{n+1}\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}$ (2.25)

得到

${\hat{\psi }}_j^{n+1}=\frac{{\hat{g}}_j}{1+\frac{i\tau }{4}\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}}\mathrm{,}$ (2.26)

这样 ${\hat{\psi }}_j^{n+1}$ 可以被求出，利用一次Sine变换，就可求出 ${\psi }_j^{n+1}$ 。在时间方向上就形成了一次完整的迭代过程。体算法可以总结如下：

1) 输入a，b，T，M，L，β，以及 ${\varphi }_0\left(x\right)$ ；

2) 计算 $x_j=a+j{h}_x$ ， $j=0,1,\cdots ,M$ ， $t_n=n\tau$ ， ${\psi }_0^n={\psi }_M^n=0$ ， $n=0,1,\cdots ,L$ 。从 $n=0$ 出发；

3) 计算 $f_j={\left|{\psi }_j^n\right|}^2$ ，通过Sine逆变换得 ${\hat{f}}_j$ ，利用 ${\hat{\varphi }}_j^n={\hat{f}}_j{\left(\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2}\right)}^{-1}\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)$ ，再对 ${\hat{\varphi }}_j^n$ 做Sine 变换得 ${\varphi }_j^n$ ；

4) 由式子 ${\psi }_j^{\mathrm{<mo>\; *\; </mo>}}={\psi }_j^n{\text{e}}^{-i\tau \left(V\left(x_j\right)+\beta {\left|\psi \left(x_j\mathrm{,}{t}_n\right)\right|}^2+{\varphi }_j^n\right)}$ ，计算出 ${\psi }_j^{\mathrm{<mo>\; *\; </mo>}}$ ( $j=0,1,\cdots ,M$ )；

5) 通过Sine逆变换求出 ${\hat{\psi }}_j^{\mathrm{<mo>\; *\; </mo>}}$ ，再利用关系式(2.24)，求出 ${\left({\hat{\psi }}_j^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}$ 。利用Sine变换求出 ${\psi }_j^{\mathrm{<mo>\; *\; </mo>}}{}_{xx}$ ，再从 $g_j$ 的Sine逆变换求出 ${\hat{g}}_j$ ；

6) 由 ${\hat{\psi }}_j^{n+1}=\frac{{\hat{g}}_j}{1+\frac{i\tau }{4}\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}}$ 得到 ${\hat{\psi }}_j^{n+1}$ ，再对其做一次Sine变换就可得到 ${\psi }_j^{n+1}$ 。再令 $n=n+\text{1}$ ，重复3~6直到 $n=L$ 为止。

2.2. 二维非线性薛定谔泊松方程组的求解方法

考虑二维零边界的非线性薛定谔泊松方程组的初边值问题：

$i{\psi }_t=-\frac{1}{2}\Delta \psi +V\psi +\beta {\left|\psi \right|}^2\psi +\varphi \left(x,y,t\right)\psi ,a<x<b,c<y<d,t>0,$ (2.27)

$-\Delta \varphi \left(x,y,t\right)={\left|\psi \right|}^2,$ (2.28)

$\psi \left(a,y,t\right)=0;\psi \left(b,y,t\right)=0,$

$\psi \left(x,c,t\right)=0;\psi \left(x,d,t\right)=0,$

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

其中 $\Delta ={\partial }_{xx}+{\partial }_{yy}$ ， $a<x<b$ ， $c<y<d$ ， $a\mathrm{,}b\mathrm{,}c\mathrm{,}d\mathrm{,}\beta$ 均为常数， ${\varphi }_0\left(x\mathrm{,}y\right)$ ， $V=V\left(x,y\right)$ 均是已 知函数，而 $\psi \mathrm{,}\varphi$ 是未知函数。

首先，将区域 $\left[a\mathrm{,}b\right]$ 等分成M份，步长为 $h_x=\frac{b-a}{M}$ ，节点为 $x_j=a+j{h}_x,j=0,1,\cdots ,M$ ；将区域 $\left[c\mathrm{,}d\right]$ 等分成N份，步长为 $h_y=\frac{d-c}{N}$ ，节点为 $y_k=c+k{h}_y,k=0,1,\cdots ,N$ ；将区域 $\left[\mathrm{0,}T\right]$ 等分成L份，步长为 $\tau =\frac{T}{L}$ ，节点为 $t_n=n\tau ,n=0,1,\cdots ,L$ ；我们记 ${\psi }_{jk}^n,{\varphi }_{jk}^n,V_{jk}$ 分别为 $\psi \left(x_j\mathrm{,}{y}_k\mathrm{,}{t}_n\right)\mathrm{,}\varphi \left(x_j\mathrm{,}{y}_k\mathrm{,}{t}_n\right)\mathrm{,}V\left(x_j\mathrm{,}{y}_k\right)$ 的近似值。

在下面的计算过程中，我们需要定义关于 ${\hat{u}}_{kl}$ 的离散二维Sine变换以及关于 $u_{ij}$ 的离散二维Sine逆变换，它们的形式分别如下

$u_{ij}={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\displaystyle \underset{l=1}{\overset{N-1}{\sum }}}{\hat{u}}_{kl}\sin \left(\frac{ik\text{π}}{M}\right)\sin \left(\frac{jl\text{π}}{N}\right),1\le i\le M-1,1\le j\le N-1,$ (2.29)

${\hat{u}}_{kl}=\frac{2}{M}\frac{2}{N}{\displaystyle \underset{i=1}{\overset{M-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}{u}_{ij}\sin \left(\frac{ik\text{π}}{M}\right)\sin \left(\frac{jl\text{π}}{N}\right),1\le k\le M-1,1\le l\le N-1.$ (2.30)

类似于一维的求解方法，在一个时间步长 $\left[t_n,t_n{}_{+\text{1}}\right]$ 上，假设 $\psi \left(x_j\mathrm{,}{y}_k\mathrm{,}{t}_n\right)$ 已知。采用时间分裂法，先求解方程(2.27)的非线性部分： $i{\psi }_t=V+\beta {\left|\psi \right|}^2\psi +\varphi \left(x,y,t\right)\psi$ ，利用共轭变换与定积分性质得：

$\psi \left(x,y,t_{*}\right)=\psi \left(x,y,t_n\right){\text{e}}^{-i\tau \left(V\left(x,y\right)+\beta {\left|\psi \left(x,y,t_n\right)\right|}^2+\varphi \left(x,y,t_n\right)\right)},$ (2.31)

$\psi \left(x\mathrm{,}y\mathrm{,}{t}_{\mathrm{<mo>\; *\; </mo>}}\right)$ 作为时间分裂法的中间变量，即为第二部分的初始值。但是，等式右边仍然有未知量 $\varphi \left(x\mathrm{,}y\mathrm{,}{t}_n\right)$ 。为了求解它，我们要利用方程 $-{\left({\varphi }_{jk}^n\right)}_{xx}-{\left({\varphi }_{jk}^n\right)}_{yy}={\left|{\psi }_{jk}^n\right|}^2$ ，记 $f_{jk}={\left|{\psi }_{jk}^n\right|}^2$ 。对于上式做二维的Sine变换得： $-{\left({\hat{\psi }}_{jk}^n\right)}_{xx}-{\left({\hat{\psi }}_{jk}^n\right)}_{yy}={\hat{f}}_{jk}$ ，观察发现：等式右边可由已知函数 $f_{jk}$ 做二维的Sine逆变换求得；等式

左边出现了二阶导，利用关系式：

${\left({\hat{\varphi }}_{jk}^n\right)}_{xx}=-{\hat{\varphi }}_{jk}^n\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}\mathrm{,1}\le j\le M-\mathrm{1,}$ (2.32)

${\left({\hat{\varphi }}_{jk}^n\right)}_{yy}=-{\hat{\varphi }}_{jk}^n\frac{24si{n}^2\left(\frac{k\text{π}}{2N}\right)}{h_y^2\left(cos\left(\frac{k\text{π}}{N}\right)+5\right)}\mathrm{,1}\le k\le N-\mathrm{1,}$ (2.33)

将上述两个关系式带入到方程 $-{\left({\hat{\varphi }}_{jk}^n\right)}_{xx}-{\left({\hat{\varphi }}_{jk}^n\right)}_{yy}={\hat{f}}_{jk}$ ，得

${\hat{\varphi }}_{jk}^n=\frac{{\hat{f}}_{jk}}{\frac{24si{n}^2\left(\frac{j\text{π}}{2M}\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}+\frac{24si{n}^2\left(\frac{k\text{π}}{2N}\right)}{h_y^2\left(cos\left(\frac{k\text{π}}{N}\right)+5\right)}}\mathrm{.}$ (2.34)

这样 ${\hat{\varphi }}_{jk}^n$ 就可求出，再利用一次二维的Sine变换就得到 ${\varphi }_{jk}^n$ 。然后将其带入到方程(2.31)中，等式左边的中间变量 ${\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}$ 就可以表示出来，进而可以作为第二部分的初始值。

考虑线性部分： $i{\psi }_t=-\frac{1}{2}\Delta \psi$ ，对于此方程，我们在时间方向上采用crank-nicolson法，空间上采用四

阶紧致差分格式来离散，具体如下：

$i\frac{{\psi }_{jk}^{n+1}-{\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}}{\tau }=-\frac{1}{4}\left(\Delta \left({\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)+\Delta \left({\psi }_{jk}^{n+1}\right)\right)$

$\left(1-\frac{i\tau }{4}\Delta \right){\psi }_{jk}^{n+1}=\left(1+\frac{i\tau }{4}\Delta \right){\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}$

${\psi }_{jk}^{n+1}-\frac{i\tau }{4}\Delta \left({\psi }_{jk}^{n+1}\right)={\psi }_{jk}^{*}+\frac{i\tau }{4}\Delta \left({\psi }_{jk}^{*}\right)$

${\psi }_{jk}^{n+1}-\frac{i\tau }{4}\left({\left({\psi }_{jk}^{n+1}\right)}_{xx}+{\left({\psi }_{jk}^{n+1}\right)}_{yy}\right)={\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}+\frac{i\tau }{4}\left({\left({\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}+{\left({\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{yy}\right)\mathrm{,}$ (2.35)

对于(2.35)，我们同时对等式两边做Sine变换得到：

${\hat{\psi }}_{jk}^{n+1}-\frac{i\tau }{4}\left({\left({\hat{\psi }}_{jk}^{n+1}\right)}_{xx}+{\left({\hat{\psi }}_{jk}^{n+1}\right)}_{yy}\right)={\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}+\frac{i\tau }{4}\left({\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}+{\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{yy}\right)\mathrm{,}$ (2.36)

等式右边的 ${\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}$ ，可以由 ${\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}$ 做Sine逆变换得到。由于有关系式：

${\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}={\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\frac{12}{{\left(h_x\right)}^2}\frac{cos\left(\frac{j\text{π}}{M}\right)-1}{cos\left(\frac{j\text{π}}{M}\right)+5}\mathrm{,}$ (2.37)

${\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{yy}={\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\frac{12}{{\left(h_y\right)}^2}\frac{cos\left(\frac{k\text{π}}{N}\right)-1}{cos\left(\frac{k\text{π}}{N}\right)+5}\mathrm{,}$ (2.38)

然后分别对 ${\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}$ ， ${\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{yy}$ 进行Sine变换，这样就可以得到 ${\left({\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}$ ， ${\left({\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{yy}$ 。所以等式(2.36)右边就可以当做已知，等式左边也是利用Sine逆变换及二维的Sine变换关系式，将其换算成只含有 ${\hat{\psi }}_{jk}^{n+1}$ 。具体如下：

${\left({\hat{\psi }}_{jk}^{n+1}\right)}_{xx}={\hat{\psi }}_{jk}^{n+1}\frac{12}{{\left(h_x\right)}^2}\frac{cos\left(\frac{j\text{π}}{M}\right)-1}{cos\left(\frac{j\text{π}}{M}\right)+5}\mathrm{,}$ (2.39)

${\left({\hat{\psi }}_{jk}^{n+1}\right)}_{yy}={\hat{\psi }}_{jk}^{n+1}\frac{12}{{\left(h_y\right)}^2}\frac{cos\left(\frac{k\text{π}}{N}\right)-1}{cos\left(\frac{k\text{π}}{N}\right)+5}\mathrm{.}$ (2.40)

将上两式带入到等式(2.36)的左边，得到的就是一个只含有 ${\hat{\psi }}_{jk}^{n+1}$ 的表达式，等式(2.36)的右边记为 ${\hat{g}}_{jk}$ 。 综合等式(2.36)、等式(2.39)、等式(2.40)可得关系式：

${\hat{\psi }}_{jk}^{n+1}=\frac{{\hat{g}}_{jk}}{1-\frac{i\tau }{4}\left(\frac{12\left(cos\left(\frac{j\text{π}}{M}\right)-1\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}+\frac{12\left(cos\left(\frac{k\text{π}}{N}\right)-1\right)}{h_y^2\left(cos\left(\frac{k\text{π}}{N}\right)+5\right)}\right)}\mathrm{.}$ (2.41)

由于等式右边均可知，所以就可求得 ${\hat{\psi }}_{jk}^{n+1}$ ，在利用一次Sine变换得 ${\psi }_{jk}^{n+1}$ 。这样时间方向上 就形成了完整迭代，然后利用MATLAB程序进行编写，称为计算法则2：

1) 输入a，b，c，d，T，M，N，L，β，以及 ${\varphi }_0\left(x\mathrm{,}y\right)$ ， $h_x$ ， $h_y$ ， $\tau$ ；

2) 计算 $x_j=a+j{h}_x$ ， $j=0,1,\cdots ,M$ ， $y_k=c+k{h}_y$ ， $k=0,1,\cdots ,N$ ， $t_n=n\tau$ ， ${\psi }_{0k}^n={\psi }_{Mk}^n=0$ ， ${\psi }_{j0}^n={\psi }_{jN}^n=0$ ， $n=0,1,\cdots ,L$ 。 从 $n=0$ 出发；

3) 计算 $f_{jk}={\left|{\psi }_{jk}^n\right|}^2$ ，通过二维的Sine逆变换得 ${\hat{f}}_{jk}$ ，利用关系式：

${\hat{\varphi }}_{jk}^n=\frac{{\hat{f}}_{jk}}{-\frac{12\left(cos\left(\frac{j\text{π}}{M}\right)-1\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}-\frac{12\left(cos\left(\frac{k\text{π}}{N}\right)-1\right)}{h_y^2\left(cos\left(\frac{k\text{π}}{N}\right)+5\right)}}\mathrm{,}$ (2.42)

再对 ${\hat{\varphi }}_{jk}^n$ 作二维的Sine变换得 ${\varphi }_{jk}^n$

4) 由 ${\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}={\psi }_{jk}^n{\text{e}}^{-i\tau \left(V\left(x_j\mathrm{,}{y}_k\right)+\beta {\left|\psi \left(x_j\mathrm{,}{y}_k\mathrm{,}{t}_n\right)\right|}^2+{\varphi }_{jk}^n\right)}$ ，计算出 ${\psi }_{jk}^{\mathrm{<mo>\; *\; </mo>}}$ ；

5) 通过二维的Sine逆变换求出 ${\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}$ ，再利用关系式(2.39)，(2.40)，求出 ${\left({\hat{\psi }}_{jk}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{xx}$ ， ${\left({\hat{\psi }}_{jk}\right)}_{yy}$ ；

6) 最后由关系式(2.41)，求出 ${\hat{\psi }}_{jk}^{n+1}$ ，对其作二维的Sine变换得 ${\psi }_{jk}^{n+1}$ 。再令 $n=n+\text{1}$ ，重复3~6直到 $n=L$ 为止。

2.3. 三维非线性薛定谔泊松方程组的求解方法

考虑三维薛定谔泊松方程组：

$-i{\psi }_t=-\frac{1}{2}\Delta \psi +V\psi +\beta {\left|\psi \right|}^2+\varphi \psi ,$ (2.43)

$-\Delta \varphi \left(x,y,z,t\right)={\left|\psi \right|}^2,$ (2.44)

$\psi \left(a,y,z,t\right)=\psi \left(b,y,z,t\right)=0,$

$\psi \left(x,c,z,t\right)=\psi \left(x,d,z,t\right)=0,$

$\psi \left(x,y,e,t\right)=\psi \left(x,y,f,t\right)=0,$

$\psi \left(x,y,z,0\right)={\varphi }_0\left(x,y,z\right),$

其中 $\Delta ={\partial }_{xx}+{\partial }_{yy}+{\partial }_{zz}$ ， $a<x<b$ ， $c<y<d$ ， $e<z<f$ ， $t>0$ ，a，b，c，d，e，f，T，β均是给定常数， ${\varphi }_0,v=v\left(x,y,z\right)$ 已知函数，ϕ，ψ为未知函数。

首先，我们将区域[a, b]等分成M份，步长为 $h_x=\frac{b-a}{M}$ ；将区域[c, d]等分成N份，步长为 $h_y=\frac{d-c}{N}$ ；将区域[e, f]等分成P份，步长为 $h_z=\frac{f-e}{P}$ ；将区域[0, T]等分成L份，步长为 $\tau =\frac{T}{N}$ ；记 ${\psi }_{jkl}^n$ ， ${\varphi }_{jkl}^n$ 为 $\psi \left(x_j\mathrm{,}{y}_k\mathrm{,}{z}_l\mathrm{,}{t}_n\right)$ ， $\varphi \left(x_j\mathrm{,}{y}_k\mathrm{,}{z}_l\mathrm{,}{t}_n\right)$ 的近似值，其中 $x_j=a+j{h}_x$ ， $y_k=c+k{h}_y$ ， $z_l=e+l{h}_l$ ， $t_n=n\tau$ ， $0\le j\le M$ ， $0\le k\le N$ ， $0\le l\le P$ ， $0\le n\le L$ 。

在下面的计算过程中，我们需要定义关于 ${\hat{u}}_{klq}$ 的离散三维Sine变换以及关于 $u_{ijp}$ 的离散三维Sine逆变换，它们的形式分别如下

$u_{ijp}={\displaystyle \underset{k=1}{\overset{M-1}{\sum }}}{\displaystyle \underset{l=1}{\overset{N-1}{\sum }}}{\displaystyle \underset{q=1}{\overset{P-1}{\sum }}}{\hat{u}}_{klq}sin\left(\frac{ik\text{π}}{M}\right)sin\left(\frac{jl\text{π}}{N}\right)sin\left(\frac{qp\text{π}}{P}\right)\mathrm{,}$ (2.45)

${\hat{u}}_{klq}=\frac{2}{M}\frac{2}{N}\frac{2}{P}{\displaystyle \underset{i=1}{\overset{M-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{N-1}{\sum }}}{\displaystyle \underset{p=1}{\overset{P-1}{\sum }}}{u}_{ijp}\sin \left(\frac{ik\text{π}}{M}\right)\sin \left(\frac{jl\text{π}}{N}\right)\sin \left(\frac{qp\text{π}}{P}\right),$ (2.46)

下面讨论具体的求解过程。在时间方向上采用时间分裂法，在空间方向上采用基于快速Sine变换的四阶紧致差分格式。首先考虑在一个时间步长 $\left[t_n,t_n{}_{+\text{1}}\right]$ 上，假设 ${\psi }_{jkl}^n$ 已知，然后将方程(2.43)分为两步求解：

1) 考虑非线性部分： $i{\psi }_t=V\psi +\beta {\left|\psi \right|}^2\psi +\varphi \psi$ ，利用共轭变换与积分性质得：

$\psi \left(x,y,z,t_{*}\right)=\psi \left(x,y,z,t_n\right){\text{e}}^{-i\tau \left(V\left(x,y,z\right)+\beta {\left|\psi \left(x,y,z,t_n\right)\right|}^2+\varphi \left(x,y,z,t_n\right)\right)}.$ (2.47)

由于 $\varphi \left(x\mathrm{,}y\mathrm{,}z\mathrm{,}{t}_n\right)$ 是未知量，所以等式右边无法直接求出。考虑方程(2.44)类似于泊松方程， 我们采用基于Sine变换的四阶紧致差分格式来求解：

$-{\varphi }_{xx}-{\varphi }_{yy}-{\varphi }_{zz}={\left|{\psi }_{jkl}^n\right|}^2\mathrm{,}$

$-{\hat{\varphi }}_{xx}-{\hat{\varphi }}_{yy}-{\hat{\varphi }}_{zz}={\left|{\hat{\psi }}_{jkl}^n\right|}^2\mathrm{,}$ (2.48)

记 ${\hat{f}}_{jlk}^n={\left|{\varphi }_{jkl}^n\right|}^2$ ，由基于Sine变换的四阶紧差分格式有如下关系式：

${\left({\hat{\psi }}_{jkl}^n\right)}_{xx}={\hat{\psi }}_{jlk}^n\frac{12\left(cos\left(\frac{j\text{π}}{M}\right)-1\right)}{h_x^2\left(cos\left(\frac{j\text{π}}{M}\right)+5\right)}\mathrm{,}$

${\left({\hat{\psi }}_{jkl}^n\right)}_{yy}={\hat{\psi }}_{jlk}^n\frac{12\left(cos\left(\frac{k\text{π}}{N}\right)-1\right)}{h_y^2\left(cos\left(\frac{k\text{π}}{N}\right)+5\right)}\mathrm{,}$

${\left({\hat{\psi }}_{jkl}^n\right)}_{zz}={\hat{\psi }}_{jlk}^n\frac{12\left(cos\left(\frac{l\text{π}}{P}\right)-1\right)}{h_z^2\left(cos\left(\frac{l\text{π}}{P}\right)+5\right)}\mathrm{.}$ (2.49)