Conservative finite difference schemes for the chiral nonlinear Schrödinger equation

Abstract

In this paper, we derive three finite difference schemes for the chiral nonlinear Schrödinger equation (CNLS). The CNLS equation has two kinds of progressive wave solutions: bright and dark soliton. The proposed methods are implicit, unconditionally stable and of second order in space and time directions. The exact solutions and the conserved quantities are used to assess the efficiency of these methods. Numerical simulations of single bright and dark solitons are given. The interactions of two bright solitons are also displayed.

1 Introduction

The chiral nonlinear Schrödinger (CNLS) equation [1, 2] we are going to study is given by

$$ i\frac{\partial\psi}{\partial t}+\frac{1}{2}\frac{\partial^{2}\psi}{ \partial x^{2}}+i\lambda \biggl( \psi^{\ast}\frac{\partial\psi }{\partial x }-\psi\frac{\partial\psi}{\partial x}^{\ast} \biggr) \psi=0, \quad -\infty< x< \infty, $$

(1)

where $\psi(x,t)$ is a complex-valued function, $\psi^{\ast}(x,t)$ denotes the complex conjugate of $\psi(x,t)$, and λ is a nonlinear coupling constant appearing through derivative coupling. This kind of nonlinearity is also known as the current density, unlike the case of cubic nonlinearity which is also known as the Kerr nonlinearity. The chiral nonlinear Schrödinger is a non-integrable equation by the classical method of inverse scattering. The single bright soliton solution of Eq. (1) is given by Biswas [1, 3] as

$$ \psi(x,t)=A\operatorname{sech}\bigl(\beta(x-ct)\bigr)\exp \bigl(i(cx+\omega t)\bigr), $$

(2)

where

$$ \beta=A\sqrt{2\lambda c}\quad \text{and}\quad \omega= \frac{c(2\lambda A^{2}-c)}{2}, $$

(3)

A is the amplitude of the soliton, β is the inverse width of the soliton, c is the soliton velocity, and ω is the wave number. The bright soliton solution (2) exists for $\lambda c>0$. The dark soliton solution exists for $\lambda c<0$ and has the form

$$ \psi(x,t)=A\tanh\bigl(\beta(x-ct)\bigr)\exp\bigl(i(cx+\omega t) \bigr), $$

(4)

where

$$ \beta=A\sqrt{-2\lambda c}\quad \text{and}\quad \omega= \frac{c(4\lambda A^{2}-c)}{2}. $$

(5)

Thus Eq. (1) has bright or dark soliton solutions that are given by (2) or (4), respectively, depending on the sign of λc. This phenomenon makes the solitons chiral.

Eq. (1) has at least four integrals of motion of conserved quantities, namely [1, 2, 4]

$$\begin{aligned}& I_{1}(t)=\int_{-\infty}^{\infty}\bigl\vert \psi(x,t)\bigr\vert ^{2}\, dx, \end{aligned}$$

(6)

$$\begin{aligned}& I_{2}(t)=i\int_{-\infty}^{\infty}\bigl(\bar{ \psi}(x,t)\psi_{x}(x,t)-\psi (x,t)\bar{\psi}_{x}(x,t)\bigr) \, dx, \end{aligned}$$

(7)

$$\begin{aligned}& I_{3}(t)=\int_{-\infty}^{\infty}\bigl[\lambda \bigl\vert \psi(x,t)\bigr\vert ^{4}+2i\bigl(\bar{\psi }(x,t) \psi_{x}(x,t)-\psi(x,t)\bar{\psi}_{x}(x,t)\bigr)\bigr]\, dx \end{aligned}$$

(8)

and

$$ I_{4}(t)=\int_{-\infty}^{\infty} \bigl[\lambda\bigl\vert \psi_{x}(x,t)\bigr\vert ^{2}+2i \bigl(\bar {\psi}(x,t)\psi_{x}(x,t)-\psi(x,t)\bar{\psi}_{x}(x,t) \bigr)+\lambda\bigl\vert \psi (x,t)\bigr\vert ^{6}\bigr]\, dx. $$

(9)

Due to the exponential decay of the bright soliton solution [5] when $|x|\rightarrow\infty$, the conserved quantities (6)-(9) are well defined. By substituting the bright soliton solution (2) into the conserved quantities (6)-(9), we obtain

$$\begin{aligned}& I_{1}=\frac{\beta}{\lambda c},\qquad I_{2}=\frac{2\beta}{\lambda c}, \qquad I_{3}=\frac{4A^{2}}{3\beta}\bigl(\lambda A^{2}+6c\bigr), \\& I_{4}=\frac{4A^{2}}{ 3b}\bigl(15c^{2}+5 \beta^{2}+40\lambda cA^{2}+8\lambda^{2}A^{4} \bigr). \end{aligned}$$

The conserved quantities (6)-(9) using dark soliton solution (4) are not well defined [5] due the nonzero boundary condition as $|x|\rightarrow \infty$. To overcome this difficulty, we present a Robin-type boundary condition [5], which will lead us to a modified form of the conserved $I_{1}(t)$ which is conserved exactly as we will see in the Numerical results section.

By assuming $\psi(x,t)=u(x,t)+iv(x,t)$, where $u(x,t)$, $v(x,t)$ are real functions, CNLS equation (1) can be written as the nonlinear coupled system [6–8] as follows:

$$\begin{aligned}& \frac{\partial u}{\partial t}+\frac{1}{2}\frac{\partial^{2}v}{\partial x^{2}}+2\lambda\biggl(u \frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x}\biggr)v=0, \end{aligned}$$

(10)

$$\begin{aligned}& \frac{\partial v}{\partial t}-\frac{1}{2}\frac{\partial^{2}u}{\partial x^{2}}-2\lambda\biggl(u \frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x}\biggr)u=0. \end{aligned}$$

(11)

The resulting system (10)-(11) can be displayed in a matrix vector form as

$$ \frac{\partial\mathbf{w}}{\partial t}+\frac{1}{2}A\frac{\partial^{2}\mathbf{w}}{\partial x^{2}}+G( \mathbf{w})A\mathbf{w}=\mathbf{0}, $$

(12)

where

$$ \mathbf{w}=\left [ \textstyle\begin{array}{@{}c@{}} u \\ v\end{array}\displaystyle \right ] , \qquad A=\left ( \textstyle\begin{array}{@{}c@{\quad}c@{}} 0 & 1 \\ -1 & 0\end{array}\displaystyle \right ) , \qquad G( \mathbf{w})=2\lambda \biggl( u\frac{\partial v}{\partial x}-v \frac{\partial u}{\partial x} \biggr) . $$

(13)

There are many theoretical and numerical studies in the literature about the nonlinear Schrödinger equations (NLS). Most of these works are motivated to study single NLS and coupled NLS (see [6–9] and references therein). However, to the authors’ knowledge, there are few numerical studies for the CNLS equation. In this paper we have derived three conservative finite difference schemes for the CNLS equation. The paper is organized as follows. In Section 2, three conservative schemes are proposed for the numerical solution of the chiral NLS. In Section 3, theoretical and numerical conservation properties are proved. Accuracy of the proposed schemes is studied in Section 4. Stability analysis is given in Section 5. Numerical results are presented in Section 6. Finally, some conclusions are drawn in Section 7.

2 Numerical methods

We will consider the numerical solution of the nonlinear system (10)-(11) in a finite interval $[x_{L},x_{R}]$. We assume $x_{m}=x_{L}+mh$, where $m=1,2,\ldots,M-1$, and h is called the space grid size, also we assume $t_{n}=nk$, k is the time step size. We denote the exact and numerical solutions at the grid point $(x_{m},t_{n})$ by $\mathbf{w}_{m}^{n}$ and $\mathbf{W}_{m}^{n}$, respectively. In this work we will present the following numerical schemes for solving (12).

2.1 Scheme 1 (nonlinear implicit scheme)

In this scheme, we will use Crank Nicholson like approach [6, 7, 10]. The scheme we propose can be given as

$$\begin{aligned}& \frac{\mathbf{W}_{m}^{n+1}-\mathbf{W}_{m}^{n}}{k}+\frac{1}{2}A\delta _{x}^{2} \biggl[ \frac{\mathbf{W}_{m}^{n+1}+\mathbf{W}_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}+ \biggl\{ \frac{G(\mathbf{W}_{m}^{n+1})+G(\mathbf{W}_{m}^{n})}{2} \biggr\} A \biggl[ \frac{\mathbf{W}_{m}^{n+1}+\mathbf{W}_{m}^{n}}{2} \biggr] =0, \quad m=1,2,\ldots,M-1. \end{aligned}$$

(14)

The scheme in (14) is a nonlinear implicit difference scheme. This will lead us to a block nonlinear tridiagonal system. This system can be solved by using any iterative method such as Newton’s method or fixed point method. In this work we adopt the latter. The fixed point iterative method we used to solve (14) can be displayed as

$$\begin{aligned}& \frac{\mathbf{W}_{m}^{n+1,(s+1)}-\mathbf{W}_{m}^{n}}{k}+\frac {1}{2}A\delta _{x}^{2} \biggl[ \frac{\mathbf{W}_{m}^{n+1,(s+1)}+\mathbf {W}_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}+ \biggl\{ \frac{G(\mathbf{W}_{m}^{n+1,(s)})+G(\mathbf {W}_{m}^{n})}{2} \biggr\} A \biggl[ \frac{\mathbf{W}_{m}^{n+1,(s+1)}+\mathbf{W}_{m}^{n}}{2} \biggr] =0,\quad m=1,2,\ldots,M-1 \end{aligned}$$

(15)

for $s=0,1,\ldots$ , where the initial approximation is taken as $\mathbf {W}^{n+1,(0)}=\mathbf{W}^{n}$. In each iteration, a block tridiagonal system is solved by Crout’s method. The iteration continues until the condition

$$ \bigl\Vert \mathbf{W}_{m}^{n+1,(s+1)}-\mathbf{W}_{m}^{n,(s)} \bigr\Vert \leq10^{-6} $$

is satisfied, and the value $W_{m}^{n+1,(s+1)}$ is used as $W_{m}^{n+1}$. The iteration procedure is repeated at each time level. The scheme conserves the discrete analog of the conserved quantity (6). The scheme is of second order accuracy in time and space, and it is unconditionally stable.

2.2 Scheme 2 (linearly implicit scheme)

The second scheme we present in this work is the three time level scheme [6–8, 10]

$$\begin{aligned}& \frac{\mathbf{W}_{m}^{n+1}-\mathbf{W}_{m}^{n-1}}{2k}+\frac{1}{2}A\delta _{x}^{2} \biggl[ \frac{\mathbf{W}_{m}^{n+1}+\mathbf{W}_{m}^{n-1}}{2h^{2}} \biggr] \\& \quad {}+G\bigl(\mathbf{W}_{m}^{n}\bigr)A \biggl[ \frac{\mathbf{W}_{m}^{n+1}+\mathbf {W}_{m}^{n-1}}{2} \biggr] =0, \end{aligned}$$

(16)

where

$$\begin{aligned}& G\bigl(\mathbf{W}_{m}^{n}\bigr) = 2\lambda \biggl[ U_{m}^{n}\frac{(V_{m+1}^{n}-V_{m-1}^{n})}{2h}-V_{m}^{n} \frac {(U_{m+1}^{n}-U_{m-1}^{n})}{2h} \biggr] , \end{aligned}$$

(17)

$$\begin{aligned}& \delta_{x}^{2}\mathbf{W}_{m}^{n} = \mathbf{W}_{m+1}^{n}-2\mathbf{W}_{m}^{n}+ \mathbf{W}_{m-1}^{n}. \end{aligned}$$

(18)

On expansion of the central difference operator, one can end with the following equation:

$$\begin{aligned}& pA\mathbf{W}_{m-1}^{n+1}+\bigl(I-2pA+kG\bigl( \mathbf{W}_{m}^{n}\bigr)A\bigr)\mathbf{W}_{m}^{n+1}+pA \mathbf{W}_{m+1}^{n+1} \\& \quad =-pA\mathbf{W}_{m-1}^{n-1}+\bigl(I+2pA-kG\bigl( \mathbf{W}_{m}^{n}\bigr)A\bigr)\mathbf{W}_{m}^{n-1}-pA \mathbf{W}_{m+1}^{n-1}, \end{aligned}$$

(19)

$p=\frac{k}{2h^{2}}$. The proposed scheme (19) forms a block linear tridiagonal system in the unknown vector $\mathbf{W}^{n+1}$ and can be solved directly using Crout’s method. The scheme conserves the discrete analog of the conserved quantity (6). The scheme is of second order accuracy in time and space, it is unconditionally stable according to von Neumann stability analysis. In order to start the iteration in this scheme, we need the solution at $t=0$ and $t=k$, this can be easily obtained from the initial condition for $t=0$ and any two-level scheme, like Scheme 1 for $t=k$.

2.3 Scheme 3 (linearly implicit scheme 3)

In order to overcome the difficulty of solving the nonlinear block tridiagonal system obtained in Scheme 1, we present the linearized implicit scheme [8]

$$\begin{aligned}& \frac{\mathbf{W}_{m}^{n+1}-\mathbf{W}_{m}^{n}}{k}+\frac{1}{2}A\delta _{x}^{2} \biggl[ \frac{\mathbf{W}_{m}^{n+1}+\mathbf{W}_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}+ \biggl\{ \frac{3G(\mathbf{W}_{m}^{n})-G(\mathbf{W}_{m}^{n-1})}{2} \biggr\} A \biggl[ \frac{\mathbf{W}_{m}^{n+1}+\mathbf{W}_{m}^{n}}{2} \biggr] =0, \end{aligned}$$

(20)

where

$$ G\bigl(\mathbf{W}_{m}^{n}\bigr)=2\lambda \biggl[ U_{m}^{n}\frac{(V_{m+1}^{n}-V_{m-1}^{n})}{2h}-V_{m}^{n} \frac {(U_{m+1}^{n}-U_{m-1}^{n})}{2h} \biggr] , $$

(21)

where $G(\mathbf{w})$ is approximated by the extrapolation formula $\frac{3G(\mathbf{W}_{m}^{n})-G(\mathbf{W}_{m}^{n-1})}{2}$.

The resulting system in (20) and (21) is a linear block tridiagonal system for the unknown numerical solution $\mathbf{W}^{n+1}$, which can be easily solved by Crout’s method. The scheme conserves the discrete analog of the conserved quantity (6). The accuracy of Scheme 3 is of second order in time and space, it is unconditionally stable according to the von Neumann stability analysis. The scheme is a three-level scheme and the solutions at $t=0$ and $t=k$ are required in order to get the solution at $t=nk$, $n=2,3,\ldots$ , the same procedure in Scheme 2 can be adopted.

3 Conserved quantity

To prove that the decomposed system

$$\begin{aligned}& \frac{\partial u}{\partial t}+\frac{1}{2}\frac{\partial^{2}v}{\partial x^{2}}+2\lambda \biggl(u\frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x}\biggr)v=0, \end{aligned}$$

(22)

$$\begin{aligned}& \frac{\partial v}{\partial t}-\frac{1}{2}\frac{\partial^{2}u}{\partial x^{2}}-2\lambda \biggl(u\frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x}\biggr)u=0 \end{aligned}$$

(23)

satisfies the conserved quantity (6), we multiply Eq. (22) and Eq. (23) by 2u and 2v, respectively, this will lead us to the following system:

$$\begin{aligned}& 2u\frac{\partial u}{\partial t}+u\frac{\partial^{2}v}{\partial x^{2}}+4\lambda\biggl(u \frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x} \biggr)uv=0, \end{aligned}$$

(24)

$$\begin{aligned}& 2v\frac{\partial v}{\partial t}-v\frac{\partial^{2}u}{\partial x^{2}}-4\lambda\biggl(u \frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x} \biggr)uv=0. \end{aligned}$$

(25)

Adding (24) and (25) will lead us, after some manipulation, to the following equation:

$$ \frac{\partial}{\partial t}\bigl(u^{2}+v^{2}\bigr)+ \frac{\partial}{\partial x}(uv_{x}-vu_{x})=0. $$

(26)

By integrating (26) with respect to x, we get

$$ \frac{\partial}{\partial t}\int_{-\infty}^{\infty } \bigl(u^{2}+v^{2}\bigr)\, dx+(uv_{x}-vu_{x})|_{-\infty}^{\infty}=0. $$

(27)

By imposing the vanishing boundary conditions, Eq. (27) will be reduced to

$$ \int_{-\infty}^{\infty}\bigl(u^{2}+v^{2} \bigr)\, dx=\mathrm{constant}, $$

which is Eq. (6).

To prove that the proposed schemes preserve the discrete analog of invariant (6), we need the following lemma [9, 10].

Lemma 1

For any two discrete functions $\{ u_{m}|m=0,1,\ldots,M \}$ and $\{ v_{m}|m=0,1,\ldots,M \}$, there is the identity

$$ \sum_{m=1}^{M-1}u_{m}(v_{m})_{xx}=-\sum_{m=1}^{M-1}(u_{m})_{x}(v_{m})_{x}-u_{0}(v_{0})_{x}+u_{M}(v_{M})_{\bar {x}}, $$

where

$$\begin{aligned}& (u_{m})_{xx} = \frac{1}{h^{2}}(u_{m+1}-2u_{m}+u_{m-1}), \\& (u_{m})_{x} = \frac{1}{h}(u_{m+1}-u_{m}), \qquad (u_{m})_{\bar{x}}=\frac{1}{h}(u_{m}-u_{m-1}). \end{aligned}$$

We start with Scheme 1. To prove that Scheme 1 preserves the discrete analog of invariant (6), we rewrite (14) in a component-wise form as

$$\begin{aligned}& \frac{U_{m}^{n+1}-U_{m}^{n}}{k}+\frac{1}{2}\delta_{x}^{2} \biggl[ \frac{V_{m}^{n+1}+V_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}+ \biggl[ \frac{G(\mathbf{W}_{m}^{n+1})+G(\mathbf{W}_{m}^{n})}{2} \biggr] \biggl[ \frac{V_{m}^{n+1}+V_{m}^{n}}{2} \biggr] =0 \end{aligned}$$

(28)

and

$$\begin{aligned}& \frac{V_{m}^{n+1}-V_{m}^{n}}{k}-\frac{1}{2}\delta_{x}^{2} \biggl[ \frac{U_{m}^{n+1}+U_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}- \biggl[ \frac{G(\mathbf{W}_{m}^{n+1})+G(\mathbf{W}_{m}^{n})}{2} \biggr] \biggl[ \frac{U_{m}^{n+1}+U_{m}^{n}}{2} \biggr] =0. \end{aligned}$$

(29)

By multiplying (28) by $(U_{m}^{n+1}+U_{m}^{n})$ and (29) by $(V_{m}^{n+1}+V_{m}^{n})$ and by adding the resulting equations, we obtain

$$ \frac{1}{k}\sum_{m=1}^{M-1} \bigl\{ \bigl[ U^{2}+V^{2} \bigr] _{m}^{n+1}- \bigl[ U^{2}+V^{2} \bigr] _{m}^{n} \bigr\} =0 $$

or

$$ \sum_{m=1}^{M-1} \bigl[ U^{2}+V^{2} \bigr] _{m}^{n+1}=\sum _{m=1}^{M-1} \bigl[ U^{2}+V^{2} \bigr] _{m}^{n}, $$

which is the discrete analog of the conserved quantity (6).

To prove that Scheme 2 preserves the discrete analog of invariant (6), we adopt the same procedure as above and write Scheme 2 in a component-wise form as follows:

$$\begin{aligned}& \frac{U_{m}^{n+1}-U_{m}^{n}}{k}+\frac{1}{2}\delta_{x}^{2} \biggl[ \frac{V_{m}^{n+1}+V_{m}^{n-1}}{2h^{2}} \biggr] +G\bigl(\mathbf{W}_{m}^{n}\bigr) \biggl[ \frac{V_{m}^{n+1}+V_{m}^{n-1}}{2} \biggr] =0, \end{aligned}$$

(30)

$$\begin{aligned}& \frac{V_{m}^{n+1}-V_{m}^{n-1}}{2k}-\delta_{x}^{2} \biggl[ \frac{U_{m}^{n+1}+U_{m}^{n-1}}{2h^{2}} \biggr] -G\bigl(\mathbf{W}_{m}^{n}\bigr) \biggl[ \frac{U_{m}^{n+1}+U_{m}^{n-1}}{2} \biggr] =0. \end{aligned}$$

(31)

By multiplying (30) by $(U_{m}^{n+1}+U_{m}^{n-1})$ and (31) by $(V_{m}^{n+1}+V_{m}^{n-1})$ and adding the resulting equations, we get

$$ \frac{1}{k}\sum_{m=1}^{M-1} \bigl\{ \bigl[ U^{2}+V^{2} \bigr] _{m}^{n+1}- \bigl[ U^{2}+V^{2} \bigr] _{m}^{n-1} \bigr\} =0 $$

or

$$ \sum_{m=1}^{M-1} \bigl[ U^{2}+V^{2} \bigr] _{m}^{n+1}=\sum _{m=1}^{M-1} \bigl[ U^{2}+V^{2} \bigr] _{m}^{n-1}, $$

which is the discrete analog of (6).

Finally, to prove the conservation property of Scheme 3, we again write the scheme in a component-wise form as

$$\begin{aligned}& \frac{U_{m}^{n+1}-U_{m}^{n}}{k}+\frac{1}{2}\delta_{x}^{2} \biggl[ \frac{V_{m}^{n+1}+V_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}+ \biggl\{ \frac{3G(\mathbf{W}_{m}^{n})-G(\mathbf{W}_{m}^{n-1})}{2} \biggr\} \biggl[ \frac{V_{m}^{n+1}+V_{m}^{n}}{2} \biggr] =0, \end{aligned}$$

(32)

$$\begin{aligned}& \frac{V_{m}^{n+1}-V_{m}^{n}}{k}-\frac{1}{2}\delta_{x}^{2} \biggl[ \frac{U_{m}^{n+1}+U_{m}^{n}}{2h^{2}} \biggr] \\& \quad {}- \biggl\{ \frac{3G(\mathbf{W}_{m}^{n})-G(\mathbf{W}_{m}^{n-1})}{2} \biggr\} \biggl[ \frac{U_{m}^{n+1}+U_{m}^{n}}{2} \biggr] =0. \end{aligned}$$

(33)

Now by multiplying (32) by $(U_{m}^{n+1}+U_{m}^{n})$ and (33) by $(V_{m}^{n+1}+V_{m}^{n})$ and by addition and taking the summation, we end with

$$ \frac{1}{k}\sum_{m=1}^{M-1} \bigl\{ \bigl[ U^{2}+V^{2} \bigr] _{m}^{n+1}- \bigl[ U^{2}+V^{2} \bigr] _{m}^{n} \bigr\} =0, $$

and this gives

$$ \sum_{m=1}^{M-1} \bigl[ U^{2}+V^{2} \bigr] _{m}^{n+1}=\sum _{m=1}^{M-1} \bigl[ U^{2}+V^{2} \bigr] _{m}^{n}, $$

which is the discrete analog of (6). This is a good indication that all schemes will not blow up for long time integration.

4 Accuracy of Scheme 2

To study the accuracy of the proposed schemes, we study the accuracy of Scheme 2. The same procedure can be adopted for the other schemes. We start by replacing the numerical solution $\mathbf{W}_{m}^{n}$ by the exact solution $\mathbf{w}_{m}^{n}$ in (16) to get

$$ \frac{1}{2k} \bigl( \mathbf{w}_{m}^{n+1}- \mathbf{w}_{m}^{n-1} \bigr) +\frac{1}{2h^{2}}A \bigl[ \mathbf{w}_{m+1}^{\ast}-2\mathbf{w}_{m}^{\ast }+ \mathbf{w}_{m-1}^{\ast} \bigr] +G\bigl( \mathbf{w}_{m}^{n}\bigr)A\mathbf{w}_{m}^{\ast}=0, $$

(34)

where

$$ \mathbf{w}_{m}^{\ast}=\frac{\mathbf{w}_{m}^{n+1}+\mathbf{w}_{m}^{n-1}}{2}. $$

Taylor’s series expansions of all terms in (34) about the grid point $(x_{m},t_{n})$ can be given as follows:

$$\begin{aligned}& \frac{\mathbf{w}_{m}^{n+1}-\mathbf{w}_{m}^{n-1}}{2k}=\frac{\partial \mathbf{w}}{\partial t}+\frac{k^{2}}{6}\frac{\partial^{3}\mathbf{w}}{\partial t^{3}}+ \frac{k^{4}}{120}\frac{\partial^{5}\mathbf{w}}{\partial t^{5}}+\cdots, \end{aligned}$$

(35)

$$\begin{aligned}& \frac{1}{2h^{2}} \bigl[ \mathbf{w}_{m+1}^{\ast}- \mathbf{w}_{m}^{\ast}+\mathbf{w}_{m-1}^{\ast} \bigr] =\frac{\partial^{2}\mathbf{w}}{\partial x^{2}}+\frac{h^{2}}{12}\frac{\partial^{4}\mathbf{w}}{\partial x^{4}}+ \frac{k^{2}}{2}\frac{\partial^{4}\mathbf{w}}{\partial x^{2}\partial t^{2}}+\cdots, \end{aligned}$$

(36)

$$\begin{aligned}& \mathbf{w}^{\ast}=\frac{\mathbf{w}_{m}^{n+1}+\mathbf{w}_{m}^{n-1}}{2}=\mathbf{w}+ \frac{k^{2}}{2}\frac{\partial^{2}\mathbf{w}}{\partial t^{2}} +\cdots, \\& G\bigl(\mathbf{w}_{m}^{n}\bigr)=G\bigl( \mathbf{w}_{m}^{n}\bigr)+\frac{h^{2}}{6} \frac {\partial ^{3}G}{\partial x^{3}}+\cdots. \end{aligned}$$

(37)

Substituting (35)-(37) into Eq. (34) will lead us to the local truncation error (LTE)

$$\begin{aligned} \mathrm{LTE} =& \biggl[ \frac{\partial\mathbf{w}}{\partial t}+\frac{1}{2}A \frac{\partial^{2}\mathbf{w}}{\partial x^{2}}+G(\mathbf{w})A\mathbf{w} \biggr] \\ &{}+\frac{k^{2}}{6}\frac{\partial^{2}\mathbf{w}}{\partial t^{2}}+\frac{h^{2}}{12}A \frac{\partial^{4}\mathbf{w}}{\partial x^{4}}+\frac{k^{2}}{2}G(\mathbf {w})A\frac{\partial^{2}\mathbf{w}}{\partial t^{2}}+\cdots. \end{aligned}$$

(38)

The first quantity in the right-hand side of Eq. (38) is zero by the differential system under consideration, which means that Scheme 2 is of second order in space and time. Similar analysis can be done for the other schemes.

5 Stability of the proposed schemes

Von Neumann stability analysis is used to study the stability of the proposed schemes. This method is only applicable for linear schemes. To apply this method, we assume that

$$ \mathbf{U}_{m}^{n}=\mathbf{U}_{0}e^{i\beta mh} ,\quad i=\sqrt{-1}, $$

(39)

where $\beta\in\mathbb{R}$, $\mathbf{U}_{0}\in\mathbb{R}^{2}$ is substituted into the difference equation for $\mathbf{U}_{m}^{n}$, it is found that $\mathbf{U}_{m}^{n+1}$ is of the same form with $G\mathbf{U}_{0}$ replacing $\mathbf{U}_{0}$. The matrix G is called the amplification matrix. To apply von Neumann stability analysis, we will consider the linearized form of the proposed system

$$ \frac{\partial\mathbf{u}}{\partial t}+\frac{1}{2}A\frac{\partial^{2}\mathbf{u}}{\partial x^{2}}+\alpha A \mathbf{u}=\mathbf{0}, $$

(40)

where α is constant. The linearized version of Scheme 1 for (40) is

$$ \bigl(\mathbf{U}_{m}^{n+1}- \mathbf{U}_{m}^{n}\bigr)+\frac{1}{4}rA \delta_{x}^{2}\bigl(\mathbf{U}_{m}^{n+1}+ \mathbf{U}_{m}^{n}\bigr)+\frac{1}{2}k\alpha A\bigl( \mathbf{U} _{m}^{n+1}+\mathbf{U}_{m}^{n} \bigr)=\mathbf{0}. $$

(41)

By substituting (39) into (41), we get after some manipulation the amplification matrix G, which is given implicitly by

$$ [ G-I]-\biggl(r\mu A-\frac{k}{2}\alpha A\biggr)[G+I]=0, $$

(42)

where

$$ \mu=\sin^{2}\frac{\beta h}{2} . $$

Equation (42) can be written as

$$ [I-\omega A]G-[I+\omega A]=\mathbf{0}, $$

(43)

and this gives us

$$ G=[I-\omega A]^{-1}[I+\omega A]=\mathbf{0}, $$

(44)

where $\omega=\mu-\frac{k}{2}\alpha$.

The von Neumann necessary condition for the stability of a system is

$$ \max_{j} |\lambda_{j}| \leq 1,\quad j=1,2, $$

where $\lambda_{j}$ are the eigenvalues of G. The eigenvalues of the amplification matrix G are

$$\begin{aligned}& \lambda_{1}=\frac{1+i\omega}{1-i\omega}, \\& \lambda_{2}= \frac{1+i\omega}{1-i\omega}, \end{aligned}$$

the modulus of these eigenvalues is equal to 1. This means that Scheme 1 is unconditionally stable, which means that it is stable for all values of k and h, but these values should be small in order to get accurate results. The same analysis can be applied for Schemes 2 and 3.

6 Numerical results

In this section, we will test the efficiency of the numerical schemes presented in this work, by considering different numerical tests. Trapezoidal rule is used to calculate the conserved quantities.

6.1 Bright soliton solution

To study the behavior of a single bright soliton solution, we choose the initial condition

$$ \psi(x,0)=A\operatorname{sech}(\beta x)\exp(ivx) $$

and the homogenous Dirichlet boundary conditions $\psi(x,0)=0$ at $x=x_{L},x_{R}$. The following set of parameters is used:

$$\begin{aligned}& x_{L} = -25,\qquad x_{R}=25.0, \qquad h=0.05,\qquad k=0.01, \\& A = 0.5,\qquad \lambda=0.5, \qquad v=0.5,\qquad t=0,1,2,\ldots,20. \end{aligned}$$

Tables 1-3 display the errors and the conserved quantities for our proposed schemes. It is very easy to see that the results are almost identical. Those numerical data support the theoretical calculations that all schemes preserve the invariant $I_{1}$. Moreover, we see that all schemes preserve the other invariants quite well. In Figure 1, we display the modulus of the numerical solution. We have noticed that the cpu time required for producing the results in Table 1, Table 2 and Table 3 is, respectively, 7.078, 2.343 and 2.313 seconds.

Figure 1

Bright soliton.

Table 1 Conserved quantities (Scheme 1: bright soliton, $\pmb{A=0.5}$ , $\pmb{\lambda=0.5}$ , $\pmb{v=0.5}$ )

Table 2 Conserved quantities (Scheme 2: bright soliton, $\pmb{A=0.5}$ , $\pmb{\lambda=0.5}$ , $\pmb{v=0.5}$ )

Table 3 Conserved quantities (Scheme 3: bright soliton, $\pmb{A=0.5}$ , $\pmb{\lambda=0.5}$ , $\pmb{v=0.5}$ )

6.2 Dark soliton solution

To study the behavior of a dark soliton solution, we choose the initial condition

$$ \psi(x,0)=A\tanh(\beta x)\exp(icx), \quad x_{L}< x< x_{R} $$

(45)

together with the boundary conditions [5]

$$ \psi_{x}(x_{L},t)-ic\psi(x_{L},t)=0, \qquad \psi_{x}(x_{L},t)-ic\psi (x_{L},t)=0,\quad t>0. $$

(46)

By using (47), we can define the modified conserved quantity of $I_{1}(t)$ as

$$ I_{11}(t)=\int_{x_{L}}^{x_{R}} \bigl\vert \psi(x,t)\bigr\vert ^{2}\, dx+c\int_{0}^{t} \bigl[ \bigl\vert \psi (x_{R},t)\bigr\vert ^{2}-\bigl\vert \psi(x_{L},t)\bigr\vert ^{2} \bigr]\, dt,\quad t \geqslant0, $$

(47)

where the proof is given in the Appendix.

The following parameters are used in this test:

$$\begin{aligned}& x_{L} = -25,\qquad x_{R}=25.0, \qquad h=0.05,\qquad k=0.001, \\& A = 0.5,\qquad \lambda=0.5, \qquad c=-0.5,\qquad t=0,1,2,\ldots,20. \end{aligned}$$

The $L_{\infty}$, $L_{2}$ error norms and the modified conserved quantity $I_{11}$ are calculated for Scheme 1, Scheme 2 and Scheme 3 and are displayed in Tables 4, 5, and 6, respectively. It is very clear that the results of the schemes are almost identical and conserved (47) exactly. The cpu time required to produce Tables 3, 4 and 5 is 46.53, 21.78 and 21.84 seconds, respectively; we have noticed that Scheme 1 is expensive and requires double of the time required for Schemes 2 and 3. In Figure 2, we display the numerical solution of the dark soliton solution at $t=0,1,2,\ldots,20$.

Figure 2

Dark soliton (single).

Table 4 Dark soliton (Scheme 1)

Table 5 Dark soliton (Scheme 2)

Table 6 Dark soliton (Scheme 3)

6.3 Interaction of two bright solitons

To study the interaction of two bright solitons, we choose the initial condition as

$$ \psi(x,0)=\psi_{1}(x,0)+\psi_{2}(x,0), $$

where

$$\begin{aligned}& \psi_{1}(x,0) = A_{1}\operatorname{sech}\bigl( \beta_{1}(x+x_{0})\bigr)\exp(iv_{1}x), \\& \psi_{2}(x,0) = A_{2}\operatorname{sech}\bigl( \beta_{2}(x-x_{0})\bigr)\exp(iv_{2}x), \end{aligned}$$

where $x_{0}$ is chosen such that the two bright solitons are initially centered at $x=\pm x_{0 }$ and well separated. In this test we use Scheme 1 and the following parameters are selected:

$$\begin{aligned}& h = 0.1,\qquad k=0.001,\qquad x_{L}=-50,\qquad x_{0}=20, \qquad A_{1}=1.0, \\ & v_{1} = 1.0,\qquad A_{2}=0.5,\qquad v_{2}=0.1, \qquad \lambda=0.5,\qquad t=0,5,10,\ldots,150. \end{aligned}$$

The conserved quantities at different time during the interaction scenario are given in Table 7, it is very clear that the conserved quantity $I_{1}$ is exactly conserved. In Figure 3, we display the interaction scenario of the two bright solitons. We have noticed that the two solitons approach each other, interact and leave the interaction region unchanged in shape.

Figure 3

Interaction of two bright solitons.

Table 7 Interaction of two bright solitons for CNLS

7 Conclusion

In this work we have solved the chiral nonlinear Schrödinger equation numerically by deriving three different finite difference schemes. In Scheme 1, we derived a nonlinear implicit scheme, we have used a fixed point iterative method to solve the nonlinear block tridiagonal system obtained. In Schemes 2 and 3, we have derived two linearly implicit finite difference schemes. Crout’s method is used to solve the resulting linear block tridiagonal system. All numerical schemes we have derived in this work conserve the energy, and this indicates that no blow-up is expected during the simulation, and hence all schemes are stable. Concerning the accuracy, the proposed schemes all are of second order accuracy in both time and space directions. The linearized schemes, Scheme 2 and Scheme 3, are more efficient than Scheme 1 regarding the issue of the cpu execution time required.

Appendix

To prove the modified conserved quantity

$$ I_{11}(t)=\int_{x_{L}}^{x_{R}} \bigl\vert \psi(x,t)\bigr\vert ^{2}\, dx+c\int_{0}^{t} \bigl[ \bigl\vert \psi (x_{R},t)\bigr\vert ^{2}-\bigl\vert \psi(x_{L},t)\bigr\vert ^{2} \bigr] \, dt,\quad t \geqslant0, $$

(48)

we differentiate (48) with respect to t to get

$$\begin{aligned} \frac{\partial}{\partial t}I_{11}(t) =& \int_{x_{L}}^{x_{R}}(2uu_{t}+2vv_{t}) \, dx \\ & {}+c \bigl\{ \bigl[u^{2}(x_{R},t)+v^{2}(x_{R},t) \bigr]-\bigl[u^{2}(x_{L},t)+v^{2}(x_{L},t) \bigr] \bigr\} ,\quad t\geqslant0. \end{aligned}$$

(49)

Now by integration by parts of (49) and making use of the differential system

$$\begin{aligned}& \frac{\partial u}{\partial t}+\frac{1}{2}\frac{\partial^{2}v}{\partial x^{2}}+2\lambda\biggl(u \frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x}\biggr)v=0, \end{aligned}$$

(50)

$$\begin{aligned}& \frac{\partial v}{\partial t}-\frac{1}{2}\frac{\partial^{2}u}{\partial x^{2}}-2\lambda\biggl(u \frac{\partial v}{\partial x}-v\frac{\partial u}{\partial x}\biggr)u=0, \end{aligned}$$

(51)

together with the boundary conditions

$$ \psi_{x}(x_{L},t)-ic\psi(x_{L},t)=0, \qquad \psi_{x}(x_{L},t)-ic\psi (x_{L},t)=0,\quad t>0, $$

we get after some manipulation

$$ \frac{\partial}{\partial t}I_{11}(t)=0,\quad t>0, $$

which implies the conservation of the modified conserved quantity $I_{11}(t)$.