Nonlinear Schr?dinger equation with a Dirac delta potential: finite difference method
本站小编 Free考研考试/2022-01-02
Bin Cheng(程彬)1, Ya-Ming Chen(陈亚铭),2, Chuan-Fu Xu(徐传福)1, Da-Li Li(李大力)2, Xiao-Gang Deng(邓小刚)21College of Computer, National University of Defense Technology Changsha 410073, China 2College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China
* National Natural Science Foundation of China.11601517 National Natural Science Foundation of China.61772542
Abstract The nonlinear Schrödinger equation with a Dirac delta potential is considered in this paper. It is noted that the equation can be transformed into an equation with a drift-admitting jump. Then following the procedure proposed in Chen and Deng (2018 Phys. Rev. E 98 033302), a new second-order finite difference scheme is developed, which is justified by numerical examples. Keywords：nonlinear Schrödinger equation;delta potential;finite difference method
PDF (956KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Bin Cheng(程彬), Ya-Ming Chen(陈亚铭), Chuan-Fu Xu(徐传福), Da-Li Li(李大力), Xiao-Gang Deng(邓小刚). Nonlinear Schr?dinger equation with a Dirac delta potential: finite difference method*. Communications in Theoretical Physics, 2020, 72(2): 025001- doi:10.1088/1572-9494/ab6185
1. Introduction
We consider in this paper the following nonlinear Schrödinger equation with a Dirac delta potential:$\begin{eqnarray}{\rm{i}}{\partial }_{t}{\rm{\Psi }}(x,t)=\left(-\displaystyle \frac{1}{2}{\partial }_{x}^{2}+\lambda \delta (x)+g| {\rm{\Psi }}(x,t){| }^{2}\right){\rm{\Psi }}(x,t),\end{eqnarray}$where i is the imaginary unit, λ and g are constants, and δ(x) represents the Dirac delta function [1]. While the function $\Psi$(x, t) is continuous at x=0, its derivative is not. Instead, we have the following matching relation$\begin{eqnarray}{\partial }_{x}{\rm{\Phi }}(0+\,,\,t)-{\partial }_{x}{\rm{\Phi }}(0-\,,\,t)=2\lambda {\rm{\Phi }}(0,t),\end{eqnarray}$which can be obtained by firstly integrating equation (1) over the interval [−ϵ, ϵ], and then considering the limit $\varepsilon \to 0$.
Equation (1) provides an idealized model for short-range interactions [2–6]. Due to the jump condition (2), special techniques shall be employed to discretize equation (1) properly. There are some useful numerical methods for solving equation (1). For instance, a finite element method is developed in [7], where the Dirac delta function is incorporated directly into the discretization of the problem by using the inherent integration of the method. In [8], the evolution operator of equation (1) is approximated by means of the Lie evolution operator. Then spectral splitting methods are implemented. Viewing the problem as an interface at the point x=0, a so-called explicit jump immersed interface method is proposed in [9]. A weak multi-symplectic reformulation is proposed in [10, 11] for equation (1), enabling one to get some preserving properties. One can also regularize the Dirac delta function in equation (1) by using a nonsingular function [12], resulting in an equation that can be handled by some common schemes. For instance, replacing the Dirac delta function by a so-called discrete delta function [13], a wavelet collocation method is constructed in [14].
In this paper, we provide a novel method to solve equation (1) numerically. Rather than discretizing equation (1) directly, we first show in section 2 that equation (1) can be transformed into an equation with discontinuous drift, in which case both the solution and flux are continuous. Then we can employ the method proposed in [15] to solve the transformed equation. For the purpose of efficiency, we provide in section 3 a new finite difference scheme following the derivation procedure in [15]. Numerical experiments are conducted in section 4 to validate the proposed method. Finally, we draw conclusions in section 5.
2. Transformation
Consider the following one-dimensional equation:$\begin{eqnarray}a{\partial }_{t}\psi (x,t)=\left[b{\partial }_{x}^{2}+V(x)+g| \psi {| }^{2}\right]\psi (x,t),\end{eqnarray}$where$\begin{eqnarray}V(x)=\lambda \delta (x)+{V}_{0}(x)\end{eqnarray}$with V0(x) representing a regular potential. Here a, b, g and λ are constants. Let$\begin{eqnarray}\psi (x,t)={{\rm{e}}}^{c| x| }\phi (x,t),\end{eqnarray}$where c is a constant that will be determined later to cancel the δ function in equation (3). It is noted that$\begin{eqnarray}{\partial }_{x}\psi ={{\rm{e}}}^{c| x| }[c\mathrm{sgn}(x)\phi +{\partial }_{x}\phi ],\end{eqnarray}$where the ‘$\mathrm{sgn}$’ stands for the sign function. Furthermore, using the relation ${\partial }_{x}\mathrm{sgn}(x)=2\delta (x)$, we obtain immediately from equation (6) that$\begin{eqnarray}{\partial }_{x}^{2}\psi ={{\rm{e}}}^{c| x| }\{{c}^{2}\phi -2c\delta (x)\phi +2c{\partial }_{x}[\mathrm{sgn}(x)\phi ]+{\partial }_{x}^{2}\phi \}.\end{eqnarray}$Thus, it follows from equations (4) and (7) that$\begin{eqnarray}\begin{array}{l}\left[b{\partial }_{x}^{2}+V(x)\right]\psi ={e}^{c| x| }\{[{{bc}}^{2}+{V}_{0}(x)]\phi \\ \quad +\,(\lambda -2{bc})\delta (x)\phi +2{bc}{\partial }_{x}[\mathrm{sgn}(x)\phi ]+b{\partial }_{x}^{2}\phi \}.\end{array}\end{eqnarray}$Therefore, the choice $c=\lambda /(2b)$ cancels the term with the δ function. In this case, we know that the substitution of the transform (5) into equation (3) results in$\begin{eqnarray}\begin{array}{rcl}a{\partial }_{t}\phi & = & [{{bc}}^{2}+{V}_{0}(x)+{{g}{\rm{e}}}^{2c| x| }| \phi {| }^{2}]\phi \\ & & +\lambda {\partial }_{x}[\mathrm{sgn}(x)\phi ]+b{\partial }_{x}^{2}\phi .\end{array}\end{eqnarray}$
To solve equation (9) involving the sign function, two matching conditions should be imposed at x=0, i.e. the continuity of φ and the continuity of the flux, defined by$\begin{eqnarray}f(x,t)=\lambda \mathrm{sgn}(x)\phi +b{\partial }_{x}\phi .\end{eqnarray}$In addition, it is observed that equation (9) is very similar to the Fokker–Planck equation considered in [15], except the additional source term $[{{bc}}^{2}+{V}_{0}(x)+{{g}{\rm{e}}}^{2c| x| }| \phi {| }^{2}]\phi $. Therefore, it is no doubt that the numerical scheme developed in [15] can be applied directly here. However, the numerical scheme is only of second order for the case with drift-admitting jumps although it is designed to be fifth-order for the case with smooth drifts. In order to save some computational resource, we prefer to solve equation (9) by developing in the next section a new scheme of third order for the case with smooth drifts rather than fifth order.
3. Numerical scheme
Following the procedure in [15], we develop here a new scheme for solving equation (9). Using the method of line, we first consider how to discretize properly the derivative of the flux, i.e. ${\partial }_{x}f$ (see equation (10)). The key is to use a so-called staggered grid.
Suppose that the computational domain is truncated to be [xL, xR] with xL<0<xR. Then we divide the domain into two subdomains, denoted by Ω1=[xL, 0] and Ω2=[0, xR], respectively. In each subdomain Ωi, solution points and flux points are placed uniformly, where the solution points are defined by ${x}_{1,j}={x}_{L}+(j-1/2){h}_{1}$ ($1\leqslant j\leqslant {N}_{1}$) for Ω1 and ${x}_{2,j}\,={x}_{d}+(j-1){h}_{2}$ (1 ≤ j ≤ N2) for Ω2. Here hi are the spacial steps for Ωi with ${h}_{1}=-{x}_{L}/\left({N}_{1}-1/2\right)$ and ${h}_{2}={x}_{R}/({N}_{2}\,-1/2)$, where Ni are the numbers of solution points. It is noted that ${x}_{1,{N}_{1}}={x}_{\mathrm{2,1}}=0$, meaning that x=0 is set to be a solution point. We will see later that this setup is very important for ensuring the continuity of the solution at x = 0. In addition, we define the flux points by ${x}_{1,j+1/2}={x}_{L}+{{jh}}_{1}$ ($0\leqslant j\leqslant {N}_{1}-1$) for Ω1 and ${x}_{2,j+1/2}=(j-1/2){h}_{2}$ (1≤ j≤N2) for Ω2, such that the flux points and the solution points are staggered with each other (see also figure 1 in [15]). Then it follows immediately that both the two boundary points are flux points since ${x}_{\mathrm{1,1}/2}={x}_{L}$ and ${x}_{2,{N}_{2}+1/2}={x}_{R}$.
3.1. Spacial discretization
For the introduced staggered grid, we first evaluate the flux at flux points, and then compute the flux derivative at solution points.
3.1.1. Evaluation of the flux
Here we only present the scheme for the subdomain Ω1, while the scheme for the subdomain Ω2 can be obtained immediately by using a mirror-symmetric property of the grid points. Additionally, let us assume b > 0 in equation (9) without loss of generality.
At a certain time level, suppose that the values of φ are known at solution points ${x}_{1,j}$, denoted by ${\phi }_{1,j}$. Then we derive a third-order interpolation scheme to obtain the values at flux points. For stability reason, we approximate the term $\lambda \mathrm{sgn}(x)\phi $ appearing in the flux at flux points ${x}_{1,j+1/2}$ by$\begin{eqnarray}\max \{\lambda ,0\}{\phi }_{1,j+1/2}^{-}+\min \{\lambda ,0\}{\phi }_{1,j+1/2}^{+},\end{eqnarray}$where ${\phi }_{1,j+1/2}^{-}$ represents interpolation values by using left-biased stencils while ${\phi }_{1,j+1/2}^{+}$ are obtained by using right-biased stencils. The interpolation scheme can be described as follows. Let ${{\rm{\Phi }}}_{1}^{\pm }={[{\phi }_{1,1/2}^{\pm },{\phi }_{1,3/2}^{\pm },...,{\phi }_{1,{N}_{1}-1/2}^{\pm }]}^{{\rm{T}}}$ and ${{\rm{\Phi }}}_{1}\,={\left[{\phi }_{\mathrm{1,1}},{\phi }_{\mathrm{1,2}},...,{\phi }_{1,{N}_{1}}\right]}^{{\rm{T}}}$. Then the right-biased scheme can be expressed as ${{\rm{\Phi }}}_{1}^{+}={I}_{1}^{+}{{\rm{\Phi }}}_{1}$, where the N1×N1 interpolation matrix ${I}_{1}^{+}$ reads$\begin{eqnarray}{I}_{1}^{+}=\left(\begin{array}{ccccc}\tfrac{15}{8} & -\tfrac{5}{4} & \tfrac{3}{8} & & \\ \tfrac{3}{8} & \tfrac{3}{4} & -\tfrac{1}{8} & & \\ & \ddots & \ddots & \ddots & \\ & & \tfrac{3}{8} & \tfrac{3}{4} & -\tfrac{1}{8}\\ & & -\tfrac{1}{8} & \tfrac{3}{4} & \tfrac{3}{8}\end{array}\right).\end{eqnarray}$For the left-biased values, the scheme reads ${{\rm{\Phi }}}_{1}^{-}\,={I}_{1}^{-}{[{\phi }_{1,1/2}^{+},{{\rm{\Phi }}}_{1}^{{\rm{T}}}]}^{{\rm{T}}}$, where the ${N}_{1}\times ({N}_{1}+1)$ interpolation matrix is$\begin{eqnarray}{I}_{1}^{-}=\left(\begin{array}{ccccccc}1 & & & & & & \\ -\tfrac{1}{3} & 1 & \tfrac{1}{3} & & & & \\ & -\tfrac{1}{8} & \tfrac{3}{4} & \tfrac{3}{8} & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\tfrac{1}{8} & \tfrac{3}{4} & \tfrac{3}{8} & \\ & & & & -\tfrac{1}{8} & \tfrac{3}{4} & \tfrac{3}{8}\end{array}\right).\end{eqnarray}$
Let ${{\rm{\Phi }}}_{1,x}={\left[{\left({\partial }_{x}\phi \right)}_{\mathrm{1,1}/2},{\left({\partial }_{x}\phi \right)}_{\mathrm{1,3}/2},...,{\left({\partial }_{x}\phi \right)}_{1,{N}_{1}-1/2}\right]}^{{\rm{T}}}$. Then we can get the difference scheme ${{\rm{\Phi }}}_{1,x}={A}_{1}{{\rm{\Phi }}}_{1}/{h}_{1}$, where the N1 × N1 difference matrix A1 is expressed as follows:$\begin{eqnarray}{A}_{1}=\left(\begin{array}{cccccc}-2 & 3 & -1 & & & \\ -1 & 1 & & & & \\ \tfrac{1}{24} & -\tfrac{9}{8} & \tfrac{9}{8} & -\tfrac{1}{24} & & \\ & \ddots & \ddots & \ddots & \ddots & \\ & & \tfrac{1}{24} & -\tfrac{9}{8} & \tfrac{9}{8} & -\tfrac{1}{24}\\ & & & & -1 & 1\end{array}\right).\end{eqnarray}$Finally, we obtain the approximation to the flux as$\begin{eqnarray}\begin{array}{rcl}{f}_{1,j+1/2} & = & \max \{\lambda ,0\}{\phi }_{1,j+1/2}^{-}+\min \{\lambda ,0\}{\phi }_{1,j+1/2}^{+}\\ & & +b{\left({\partial }_{x}\phi \right)}_{1,j+1/2}.\end{array}\end{eqnarray}$As mentioned before, the flux ${f}_{2,j+1/2}$ for the subdomain Ω2 can be obtained by using the symmetry property of the grid. Moreover, we usually apply boundary conditions ${f}_{\mathrm{1,1}/2}\,={f}_{2,{N}_{2}+1/2}=0$ in applications.
3.1.2. Compute the derivative of the flux
Theoretically, the flux is continuous at x = 0. Thus, we can derive a difference scheme for the whole computational domain [xL, xR] directly. To describe the algorithm, let us introduce the flux vector ${\boldsymbol{f}}=\left[{{\boldsymbol{f}}}_{1}^{{\rm{T}}},{{\boldsymbol{f}}}_{2}^{{\rm{T}}}\right]$ with$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{f}}}_{1} & = & {\left[{f}_{\mathrm{1,1}/2},{f}_{\mathrm{1,3}/2},\cdot \cdot \cdot ,{f}_{1,{N}_{1}-1/2}\right]}^{{\rm{T}}},\\ {{\boldsymbol{f}}}_{2} & = & {\left[{f}_{\mathrm{2,3}/2},{f}_{\mathrm{2,5}/2},\cdot \cdot \cdot ,{f}_{2,{N}_{2}+1/2}\right]}^{{\rm{T}}},\end{array}\end{eqnarray}$and the derivative vector ${{\boldsymbol{f}}}_{x}=\left[{{\boldsymbol{f}}}_{1,x}^{{\rm{T}}},{{\boldsymbol{f}}}_{2,x}^{{\rm{T}}}\right]$ with$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{f}}}_{1,x} & = & {\left[{\left({\partial }_{x}f\right)}_{\mathrm{1,1}},{\left({\partial }_{x}f\right)}_{\mathrm{1,2}},\cdot \cdot \cdot ,{\left({\partial }_{x}f\right)}_{1,{N}_{1}}\right]}^{{\rm{T}}},\\ {{\boldsymbol{f}}}_{2,x} & = & {\left[{\left({\partial }_{x}f\right)}_{\mathrm{2,2}},{\left({\partial }_{x}f\right)}_{\mathrm{2,3}},\cdot \cdot \cdot ,{\left({\partial }_{x}f\right)}_{2,{N}_{2}}\right]}^{{\rm{T}}}.\end{array}\end{eqnarray}$Then we express the difference scheme as ${{\boldsymbol{f}}}_{x}=A{\boldsymbol{f}}$, where the entries ${a}_{j,k}$ of the difference matrix A can be obtain by using the following stencils:$\begin{eqnarray}{\left[{{\boldsymbol{f}}}_{x}\right]}_{j}=\left\{\begin{array}{ll}{{\rm{\Sigma }}}_{k=1}^{3}{a}_{j,k}{[{\boldsymbol{f}}]}_{k}, & j=1,\\ {{\rm{\Sigma }}}_{k=1}^{4}{a}_{j,k}{[{\boldsymbol{f}}]}_{j+k-2}, & 2\leqslant j\leqslant {N}_{x}-1,\\ {{\rm{\Sigma }}}_{k=1}^{3}{a}_{j,k}{[{\boldsymbol{f}}]}_{{N}_{x}+k-2}, & j={N}_{x}.\end{array}\right.\end{eqnarray}$Here [f]k denotes the kth entry of the vector f, and ${N}_{x}={N}_{1}+{N}_{2}-1$.
By using the method of Lagrangian interpolation, it is straightforward to determine the entries of A. While the stencils locate in a single domain, we have$\begin{eqnarray}\begin{array}{rcl}\,{\left[{a}_{j,k}\right]}_{1\leqslant k\leqslant 3} & = & \left\{\begin{array}{ll}\tfrac{1}{{h}_{1}}[-1,1,0], & j=1,\\ \tfrac{1}{{h}_{2}}[0,-1,1], & j={N}_{x},\end{array}\right.\\ {\left[{a}_{j,k}\right]}_{1\leqslant k\leqslant 4} & = & \left\{\begin{array}{ll}\tfrac{1}{{h}_{1}}\left[\tfrac{1}{24},-\tfrac{9}{8},\tfrac{9}{8},-\tfrac{1}{24}\right], & 2\leqslant j\leqslant {N}_{1}-2,\\ \tfrac{1}{{h}_{2}}\left[\tfrac{1}{24},-\tfrac{9}{8},\tfrac{9}{8},-\tfrac{1}{24}\right], & {N}_{1}+2\leqslant j\leqslant {N}_{x}-1.\end{array}\right.\end{array}\end{eqnarray}$While the stencils across the interface point x=0, i.e. ${N}_{1}-1\leqslant j\leqslant {N}_{1}+1$, the entries can be determined as ${a}_{j,k}={d}_{j-{N}_{1}+2,k}({h}_{1},{h}_{2})$, where the expressions for ${d}_{j,k}(x,y)$ are presented in table 1.
Table 1. Table 1.Coefficients that determine the difference scheme (18) for the cases ${N}_{1}-1\leqslant j\leqslant {N}_{1}+1$.
After discretizing the derivative of the flux, we obtain from equation (3) an ordinary differential system, which can be solved by a time-marching scheme. In this paper, the third-order Runge–Kutta scheme as used in [15] is employed (see equations (30) and (31) therein). If not stated otherwise, the time step is chosen to be $0.01\min \{{h}_{1}^{2},{h}_{2}^{2}\}$.
4. Numerical examples
In this section, we solve several problems numerically to show the validity of the aforementioned method. To measure the convergence rate, we introduce the L2-norm error, defined by$\begin{eqnarray}{L}^{2}{\rm{error}}={\left[\displaystyle \sum _{j=1}^{{N}_{1}-1}{e}_{1,j}^{2}{h}_{1}+\displaystyle \frac{1}{2}{e}_{1,{N}_{1}}^{2}\left({h}_{1}+{h}_{2}\right)+\displaystyle \sum _{j=2}^{{N}_{2}}{e}_{2,j}^{2}{h}_{2}\right]}^{1/2},\end{eqnarray}$and the ${L}^{\infty }$-norm error$\begin{eqnarray}{L}^{\infty }{\rm{error}}=\mathop{\max }\limits_{i,j}\left|{e}_{i,j}\right|,\end{eqnarray}$where ei,j are the errors between numerical results and analytic solutions. Then the numerical convergence rate is measured by$\begin{eqnarray}{\rm{rate}}=-\mathrm{ln}\left({E}_{M}/{E}_{N}\right)/\mathrm{ln}(M/N),\end{eqnarray}$where EM and EN are the errors corresponding to the cases with M and N solution points, respectively.
4.1. Stationary nonlinear Schrödinger equation
At first, we consider the stationary nonlinear Schrödinger equation$\begin{eqnarray}\left(-\displaystyle \frac{1}{2}{\partial }_{x}^{2}+\lambda \delta (x)+g| \psi {| }^{2}\right)\psi =\mu \psi ,\end{eqnarray}$where λ, g and μ are constants. Rather than solving equation (23) directly, we can solve the following time-dependent nonlinear Schrödinger equation$\begin{eqnarray}{{\rm{\Psi }}}_{t}=\left(\displaystyle \frac{1}{2}{\partial }_{x}^{2}-\lambda \delta (x)+\mu -g| {\rm{\Psi }}{| }^{2}\right){\rm{\Psi }}\end{eqnarray}$until $\Psi$t reaches the machine zero. Then we obtain equivalently the solution of equation (23).
It is obvious that equation (24) is in the form of the general equation considered in equation (3). Thus the algorithm presented in the previous sections can be applied directly. Here we just specify the initial condition for equation (24) to be Gaussian, i.e.$\begin{eqnarray}{\rm{\Psi }}(x,0)=\displaystyle \frac{1}{\sqrt{2\pi {\tau }_{0}}}{{\rm{e}}}^{-{x}^{2}/\left(2{\tau }_{0}\right)}\end{eqnarray}$with τ0 chosen to be 0.01. In the following two cases, we set the computational domains to be large enough, enabling us to impose the zero flux boundary condition$\begin{eqnarray}f\left({x}_{L},t\right)=f\left({x}_{R},t\right)=0.\end{eqnarray}$In addition, the computation is stopped until the residue, defined by$\begin{eqnarray}{{\rm{Res}}}^{n}=\mathop{\max }\limits_{j}| {{\rm{\Psi }}}_{j}^{n}-{{\rm{\Psi }}}_{j}^{n-1}| ,\end{eqnarray}$reaches the machine zero.
4.1.1. Linear case
When g = 0, equation (23) degenerates to a linear equation, which admits a bound state solution [16]$\begin{eqnarray}\psi (x)={{\rm{e}}}^{\lambda | x| }\end{eqnarray}$with μ=−λ2/2. Here we choose λ=−1/2 and truncate the computational domain to be [−20, 20]. Then we solve equation (24) and obtain the numerical solution as shown in figure 1(a). The convergence process of the numerical solution is shown in figure 1(b). The test of accuracy in table 2 demonstrates that second-order accuracy is achieved, as expected.
Figure 1.
New window|Download| PPT slide Figure 1.Numerical results for equation (23) with g=0 and λ=−1/2, where N1=N2=200 is used to compute the numerical solution that matches well with the analytic solution. (a) Numerical solution; (b) residue.
Table 2. Table 2.Accuracy test for equation (23) with g=0 and λ=−1/2. Here, N1=N2=N.
For $g\ne 0$, the nonlinear term in equation (24) comes into play. Here we confine the study to the case with g=1, which admits the analytic solution [17]$\begin{eqnarray}\psi (x)=\left\{\begin{array}{ll}-k\,\mathrm{cosech}\left(k\left(x+{x}_{0}\right)\right), & x\lt 0,\\ k\,\mathrm{cosech}\left(k\left(x-{x}_{0}\right)\right), & x\gt 0\end{array}\right.\end{eqnarray}$if $\mu =-{(\lambda -1/2)}^{2}/2$. Here x0 is determined by the equality$\begin{eqnarray}\tanh \left({{kx}}_{0}\right)=k/\lambda ,\end{eqnarray}$which is a consequence of the jump condition of equation (23) at x=0; see equation (2). In addition, k is set as $k=-\lambda -1/2$ such that ${\int }_{-\infty }^{\infty }| \psi (x){| }^{2}{\rm{d}}{x}=1$.
Here we set λ=−1 and truncate the computational domain to be [−15, 15]. As we can see from figure 2(a), the numerical results match well with the analytic solution. The convergence process of the numerical solution is shown in figure 2(b). The expected second-order convergence rate is also achieved as shown in table 3.
Figure 2.
New window|Download| PPT slide Figure 2.Numerical results for equation (23) with g=1 and $\lambda =-1$. Here, N1=N2=200. (a) Numerical solution; (b) residue.
Table 3. Table 3.Accuracy test for equation (23) with g=1 and $\lambda =-1$. Here, N1=N2=N.
Now we are trying to solve the time-dependent nonlinear Schödinger equation (1) [11]. Since $\Psi$ is a complex number for fixed x and t, we shall assume ${\rm{\Psi }}=p(x,t)+{\rm{i}}q(x,t)$ and plug it into equation (1) to obtain a couple of equations$\begin{eqnarray}{{\rm{\partial }}}_{t}p=\left(-\displaystyle \frac{1}{2}{{\rm{\partial }}}_{x}^{2}+\lambda \delta (x)+g({p}^{2}+{q}^{2})\right)q,\end{eqnarray}$$\begin{eqnarray}-{{\rm{\partial }}}_{t}q=\left(-\displaystyle \frac{1}{2}{{\rm{\partial }}}_{x}^{2}+\lambda \delta (x)+g({p}^{2}+{q}^{2})\right)p.\end{eqnarray}$Then we can use the proposed numerical method to solve this system.
When g=1, equation (1) possesses the analytic solution$\begin{eqnarray}{\rm{\Psi }}(x,t)=\psi (x)\exp \left({\rm{i}}{\left(2\lambda -1\right)}^{2}t/8\right),\end{eqnarray}$where $\Psi$(x) is defined by equation (29). We choose λ=−0.7 and set the initial condition as $\Psi$(x, 0) according to equation (33). The computational domain is truncated to be [−20, 20]. Figure 3 shows that the numerical results of the real and imaginary parts match well with the analytic solutions at time t=1. In this case, second-order accuracy is still observed, as shown in table 4.
Figure 3.
New window|Download| PPT slide Figure 3.Numerical solution of equation (1) with g=1 and λ=−0.7 at time t=1. Here, N1=N2=200. (a) Solution of the real part p; (b) solution of the imaginary part q.
Table 4. Table 4.Accuracy test for equation (1) with g=1 and $\lambda =-0.7$ at time t=1. Here, N1=N2=N.
To conclude, we have proposed in this paper a transform for the one-dimensional nonlinear Schrödinger equation with a Dirac delta potential, such that the annoying Dirac delta function can be canceled out. For the transformed equation without the delta function, we have provided an efficient second-order finite difference method. In the numerical experiments, we have considered both stationary solutions and time-dependent solutions. All the computed results validate the new method, which is second-order.
In the future, we may consider whether the same idea is applicable for the case with double delta potentials [18–20] or delta comb potentials [21–23]. In addition, we may also consider how to improve the convergence rate of the numerical scheme.
,2, Chuan-Fu Xu(徐传福)1, Da-Li Li(李大力)2, Xiao-Gang Deng(邓小刚)21College of Computer, 
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
