Explicit High-Order Method to Solve Coupled Nonlinear Schrödinger Equations

Khadijah Alamoudi, Mohmmad Said Hammoudeh

Department of Math, King Abdulaziz University, Jeddah, Saudi Arabia.

DOI: 10.4236/apm.2021.115033                

Abstract

Models of the coupled nonlinear Schrödinger equations submit various critical physical phenomena with a typical equation for optical fibres with linear refraction. In this article, we will presuppose the Compact Finite Difference method with Runge-Kutta of order 4 (explicit) method, which is sixth-order and fourth-order in space and time respectively, to solve coupled nonlinear Schrödinger equations. Many methods used to solve coupled nonlinear Schrödinger equations are second order in time and need to use extra-technique to rise up to fourth-order as Richardson Extrapolation technique. The scheme obtained is immediately fourth-order in one step. This approach is a conditionally stable method. The conserved quantities and the exact single soliton solution indicate the competence and accuracy of the article’s suggestion schemes. Furthermore, the article discusses the two solitons interaction dynamics.

Keywords

Coupled Nonlinear Schrodinger Equations, Sixth Order Method, Interaction of Two Solitons, Compact Finite Difference, Runge-Kutta of Order 4 Method

1. Introduction

An enormous assortment of physical situations can be represented by the coupled nonlinear Schrödinger equations (CNLSEs). These equations have been exhibited in a fibre communication system to capture vibration diffusion over perpendicular polarization axes in nonlinear optical fibres and wavelength division of multiplexed systems. CNLSEs also model radiation inside photorefractive matter or crystals by modelling the radiation as water wave collision. Singular waves in these equations are overwhelmingly called vector solitons in the writing as singular waves, in general, consist of two components. In the previous physical situations, clashing vector solitons is a significant case. In the last decade, the system has been examined seriously. Studies have found that, in addition to crossing through the clash, vector solitons can also rebound each other or hinder one another. This article will present (CNLSEs) as the follows [1] [2]:

$i{\Phi }_{1t}+\frac{1}{2}{\Phi }_{1xx}+\left({\left|{\Phi }_1\right|}^2+c{\left|{\Phi }_2\right|}^2\right){\Phi }_1=0,-\infty <x<\infty$ (1)

$i{\Phi }_{2t}+\frac{1}{2}{\Phi }_{2xx}+\left(c{\left|{\Phi }_1\right|}^2+{\left|{\Phi }_2\right|}^2\right){\Phi }_2=0,-\infty <x<\infty$ (2)

where ${\Phi }_1$ and ${\Phi }_2$ are the wave amplitudes in pair polarizations, and c is the cross-phase modulation coefficient. The initial conditions are:

${\Phi }_1\left(x,0\right)={\zeta }_1\left(x\right),{\Phi }_2\left(x,0\right)={\zeta }_2\left(x\right)$ (3)

and the boundary conditions are

${\Phi }_1\left(x,t\right)={\Phi }_2\left(x,t\right)=\mathrm{0,}\text{as}\left|x\right|\to \infty$ (4)

For linearly birefringent fibres $c=\frac{2}{3}$, and for elliptically birefringent fibres, c

is a positive number. When $c=0$, this system enhances two (decoupled) nonlinear Schrödinger equations (NLS), and when $c=1$, the system is a Manakov equations. In both situations, we get an integrable system. In particular, the solitons of one polarization should pass through the pulses of the opposite interaction’s polarization in the absence of forming shadows, creating an elastic collision. We have a nonintegrable system for other values of c. The clash may be nontrivial, and several complex phenomena like reflection, transmission, trapping, and creation of other solitary waves can happen. Solitons in optical fibres can be described as nonlinear beats that propagate nearly distortion-free, extending into the distance and undergoing (elastic) collision. The sixth-order compact finite difference method is perfect to derive the analytical solution of the CNLS system. We determined the exact solution of the CNLSEs as follows:

${\Phi }_1\left(x,t\right)=\sqrt{\frac{2\alpha }{1+c}}\sec h\left(x-v_{1}t\right)\exp i\left\{v_{1}x-\left[\frac{v_1^2}{2}-\alpha \right]t\right\},$

${\Phi }_2\left(x,t\right)=\pm \sqrt{\frac{2\alpha }{1+c}}\sec h\left(x-v_{1}t\right)\exp i\left\{v_{1}x-\left[\frac{v_1^2}{2}-\alpha \right]t\right\},$

Several improved numerical methods have been presented for the nonlinear Schrödinger equations [3] [4] [5] [6] [7]. Among the numerical techniques for CNLSEs, the majority of the improved methods apply finite difference methods [1] [2] [8] [9] [10] [11] [12]. Multi-symplectic methods have also been applied for resolving the CNLS systems recently [13] [14] [15] [16]. Good outcomes have boon obtained in terms of the solutions and the conserved quantities. Wang [17] presents a numerical solution of the single and coupled nonlinear Schrödinger equation using split-step finite difference method. Galerkin method is used to solve CNLSE in [1] and [18], and the outcomes are superior in comparison with those of previous numerical works. Ismail using Collocation method, produced very accurate results. Concerning the finite element method, many articles have been published on nonlinear Schrödinger equations [19] [20]. Many articles as [21] and [22] used Runge-Kutta of order 4 to solve CNLSEs, and they obtained lower accuracy in space direction (second and fourth order) without using the Compact Finite Difference method which gives sixth order.

The conserved quantities of Equations (1) and (2) are follows [1]:

1) Mass conservations:

$I_1={\displaystyle {\int }_{-\infty }^{\infty }}{\left|{\Phi }_1\right|}^2\text{d}x,$ (5)

$I_2={\displaystyle {\int }_{-\infty }^{\infty }}{\left|{\Phi }_2\right|}^2\text{d}x$ (6)

2) Energy conservations:

$I_3={\displaystyle {\int }_{-\infty }^{\infty }}\left\{{\displaystyle \underset{j=1}{\overset{2}{\sum }}}\frac{1}{2}{\left|\frac{\partial {\Phi }_j}{\partial x}\right|}^2-\frac{1}{2}{\displaystyle \underset{j=1}{\overset{2}{\sum }}}{\left|{\Phi }_j\right|}^4-c{\left|{\Phi }_1\right|}^2{\left|{\Phi }_2\right|}^2\right\}\text{d}x.$ (7)

2. Numerical Method

We suppose that, outside the interval $x_0\le x\le x_M$, the solution of the system (1) and (2) is negligible. Therefore, to study this system numerically, we prefer to consider the system

$i\frac{\partial {\Phi }_1}{\partial t}+\frac{1}{2}\frac{{\partial }^2{\Phi }_1}{\partial x^2}+\left({\left|{\Phi }_1\right|}^2+c{\left|{\Phi }_2\right|}^2\right){\Phi }_1=0,x_0<x<x_M$ (8)

$i\frac{\partial {\Phi }_2}{\partial t}+\frac{1}{2}\frac{{\partial }^2{\Phi }_2}{\partial x^2}+\left(c{\left|{\Phi }_1\right|}^2+{\left|{\Phi }_2\right|}^2\right){\Phi }_2=0,x_0<x<x_M$ (9)

where ${\Phi }_1$ and ${\Phi }_2$ are the wave gage in two polarizations and the cross-phase modulation coefficient is c, with

${\Phi }_1\left(x,0\right)={\zeta }_1\left(x\right),{\Phi }_2\left(x,0\right)={\zeta }_2(\; x\; )$

as initial conditions and

${\Phi }_1\left(x,t\right)={\Phi }_2\left(x,t\right)=0,\text{at}x=x_0\text{and}x=x_M$

as boundary conditions. To conclude the numerical method to solve 8 and 9, the rectangular mesh point with coordinates

$x_m=x_0+m\left(h\right),m=0,1,\cdots ,M,$

$t_n=nk,n=0,1,2,\cdots$

which cover the region $R=\left[x_0<x<x_M\right]\times \left[t>0\right]$. Numerically, we prefer to disjoint the complex functions ${\Phi }_1$ and ${\Phi }_2$ into their real and imaginary portions by

${\Phi }_1\left(x,t\right)=u_1\left(x,t\right)+i{u}_2\left(x,t\right),{\Phi }_2\left(x,t\right)=u_3\left(x,t\right)+i{u}_4\left(x,t\right)$ (10)

where $u_i\left(x,t\right),i=1,2,3,4$ are real functions, and we write,

$\begin{array}{l}{u}_1\left(x,0\right)={\zeta }_{1R}\left(x\right),u_2\left(x,0\right)={\zeta }_{1I}\left(x\right)\\ u_3\left(x,0\right)={\zeta }_{2R}\left(x\right),u_4\left(x,0\right)={\zeta }_{2I}\left(x\right)\end{array}$ (11)

where we have assumed

${\zeta }_1\left(x\right)={\zeta }_{1R}+i{\zeta }_{1I}\text{and}{\zeta }_2\left(x\right)={\zeta }_{2R}+i{\zeta }_{2I}$

and ${\zeta }_{1R}\mathrm{,}{\zeta }_{1I}\mathrm{,}{\zeta }_{1R}$ and ${\zeta }_{1I}$ are all real functions.

By substituting (10) into (8) and (9), the next system is obtained:

$\frac{1}{2}\frac{{\partial }^2{u}_2}{\partial x^2}=-\frac{\partial u_1}{\partial t}-G_1{u}_2$ (12)

$\frac{1}{2}\frac{{\partial }^2{u}_1}{\partial x^2}=\frac{\partial u_2}{\partial t}-G_1{u}_1$ (13)

$\frac{1}{2}\frac{{\partial }^2{u}_4}{\partial x^2}=-\frac{\partial u_3}{\partial t}-G_2{u}_4$ (14)

$\frac{1}{2}\frac{{\partial }^2{u}_3}{\partial x^2}=\frac{\partial u_4}{\partial t}-G_2{u}_3$ (15)

The previous system in a matrix-vector is

$\frac{\partial u}{\partial t}+\frac{1}{2}A\frac{{\partial }^{2}u}{\partial x^2}+F\left(u\right)u=0$ (16)

where

$u=\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]\mathrm{,}A=\left[\begin{array}{cccc}0& 1& 0& 0\\ -1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& -1& 0\end{array}\right]\mathrm{,}F\left(u\right)=\left[\begin{array}{cccc}0& G& 0& 0\\ -G& 0& 0& 0\\ 0& 0& 0& G\\ 0& 0& -G& 0\end{array}\right]\mathrm{.}$

At the grid point $\left(x_m\mathrm{,}{t}_n\right)$, the exact solution is denoted by $u_{i\mathrm{,}m}^n$ while the numerical one is denoted with $U_{i\mathrm{,}m}^n$. The Compact Finite Difference (CFD) formula is derived by [23]:

$\alpha {\left(\frac{{\partial }^2{u}_i}{\partial x^2}\right)}_{m-1}+{\left(\frac{{\partial }^2{u}_i}{\partial x^2}\right)}_m+\alpha {\left(\frac{{\partial }^2{u}_i}{\partial x^2}\right)}_{m+1}=\frac{b}{4h^2}{\delta }_{\hat{x}}^2{U}_{i,m}+\frac{a}{h^2}{\delta }_x^2{U}_{i,m}$ (17)

where

${\delta }_x^2{U}_{i,m}=U_{i,m+1}-2U_{i,m}+U_{i,m-1},$

${\delta }_{\hat{x}}^2{U}_{i,m}=U_{i,m+2}-2U_{i,m}+U_{i,m-2},i=1,2,3,4.$

Now, by Taylor Expansion, we can have the truncation error as the following [23]:

$\begin{array}{l}R\equiv \left(2\alpha +1-b-a\right){\left(\frac{{\partial }^2{u}_i}{\partial x^2}\right)}_m+\left(\alpha -\frac{1}{3}b-\frac{1}{12}a\right)h^2{\left(\frac{{\partial }^4{u}_i}{\partial x^4}\right)}_m\\ +\left(\frac{1}{12}\alpha -\frac{2}{45}b-\frac{1}{360}a\right)h^4{\left(\frac{{\partial }^6{u}_i}{\partial x^6}\right)}_m\end{array}$

if we solve $\left(2\alpha +1-b-a\right)=0$ and $\left(\alpha -\frac{1}{3}b-\frac{1}{12}a\right)=0$, we get

$a=\frac{4}{3}\left(1-\alpha \right)\text{and}b=\frac{1}{3}\left(10\alpha -1\right)$

so the truncation error becomes

$R\equiv -\frac{4}{6!}\left(11\alpha -2\right)h^4{\left(\frac{{\partial }^6{u}_i}{\partial x^6}\right)}_m+O(\; h\; 6\; )$

If $\alpha =0$ then $a=\frac{4}{3}$ and $b=-\frac{1}{3}$, which gives the explicit scheme in order four for the second derivative. Furthermore, when $\alpha =\frac{2}{11}$, the scheme becomes sixth-order accurate, in this case $a=\frac{12}{11}$ and $b=\frac{3}{11}$.

By using (12)-(15), Equations (17) can be written as

$\begin{array}{l}-\alpha {\left(\frac{\partial u_1}{\partial t}+Gu\right)}_{m-1}-{\left(\frac{\partial u_1}{\partial t}+Gu\right)}_m-\alpha {\left(\frac{\partial u_1}{\partial t}+Gu\right)}_{m+1}\\ =\frac{1}{2}\left[\frac{b}{4h^2}{\delta }_{\hat{x}}^2{U}_{2,m}+\frac{a}{h^2}{\delta }_x^2{U}_{2,m}\right]\end{array}$ (18)

$\begin{array}{l}\alpha {\left(\frac{\partial u_2}{\partial t}-G{u}_1\right)}_{m-1}+{\left(\frac{\partial u_2}{\partial t}-G{u}_1\right)}_m+\alpha {\left(\frac{\partial u_2}{\partial t}-G{u}_1\right)}_{m+1}\\ =\frac{1}{2}\left[\frac{b}{4h^2}{\delta }_{\hat{x}}^2{U}_{1,m}+\frac{a}{h^2}{\delta }_x^2{U}_{1,m}\right]\end{array}$ (19)

$\begin{array}{l}-\alpha {\left(\frac{\partial u_3}{\partial t}+G{u}_4\right)}_{m-1}-{\left(\frac{\partial u_3}{\partial t}+G{u}_4\right)}_m-\alpha {\left(\frac{\partial u_3}{\partial t}+G{u}_4\right)}_{m+1}\\ =\frac{1}{2}\left[\frac{b}{4h^2}{\delta }_{\hat{x}}^2{U}_{4,m}+\frac{a}{h^2}{\delta }_x^2{U}_{4,m}\right]\end{array}$ (20)

$\begin{array}{l}\alpha {\left(\frac{\partial u_4}{\partial t}-G{u}_3\right)}_{m-1}+{\left(\frac{\partial u_4}{\partial t}-G{u}_3\right)}_m+\alpha {\left(\frac{\partial u_4}{\partial t}-G{u}_3\right)}_{m+1}\\ =\frac{1}{2}\left[\frac{b}{4h^2}{\delta }_{\hat{x}}^2{U}_{3,m}+\frac{a}{h^2}{\delta }_x^2{U}_{3,m}\right]\end{array}$ (21)

The above we can rewrite as

$\begin{array}{l}-\alpha {\stackrel{\dot{}}{U}}_{1,m-1}-{\stackrel{\dot{}}{U}}_{1,m}-\alpha {\stackrel{\dot{}}{U}}_{1,m+1}-\alpha {\left(G{U}_2\right)}_{m-1}-{\left(G{U}_2\right)}_m-\alpha {\left(G{U}_2\right)}_{m+1}\\ =p_1\left(U_{2,m-2}-2U_{2,m}+U_{2,m+2}\right)+p_2\left(U_{2,m-1}-2U_{2,m}+U_{2,m+1}\right)\end{array}$ (22)

$\begin{array}{l}\alpha {\stackrel{\dot{}}{U}}_{\mathrm{2,}m-1}+{\stackrel{\dot{}}{U}}_{\mathrm{2,}m}-\alpha {\stackrel{\dot{}}{U}}_{\mathrm{2,}m+1}-\alpha {\left(G{U}_1\right)}_{m-1}-{\left(G{U}_1\right)}_m-\alpha {\left(G{U}_1\right)}_{m+1}\\ =p_1\left(U_{\mathrm{1,}m-2}-2U_{\mathrm{1,}m}+U_{\mathrm{1,}m+2}\right)+p_2\left(U_{\mathrm{1,}m-1}-2U_{\mathrm{1,}m}+U_{\mathrm{1,}m+1}\right)\end{array}$ (23)

$\begin{array}{l}-\alpha {\stackrel{\dot{}}{U}}_{\mathrm{3,}m-1}-{\stackrel{\dot{}}{U}}_{\mathrm{3,}m}-\alpha {\stackrel{\dot{}}{U}}_{\mathrm{3,}m+1}-\alpha {\left(G{U}_4\right)}_{m-1}-{\left(G{U}_4\right)}_m-\alpha {\left(G{U}_4\right)}_{m+1}\\ =p_1\left(U_{\mathrm{4,}m-2}-2U_{\mathrm{4,}m}+U_{\mathrm{4,}m+2}\right)+p_2\left(U_{\mathrm{4,}m-1}-2U_{\mathrm{4,}m}+U_{\mathrm{4,}m+1}\right)\end{array}$ (24)

$\begin{array}{l}\alpha {\stackrel{\dot{}}{U}}_{\mathrm{4,}m-1}+{\stackrel{\dot{}}{U}}_{\mathrm{4,}m}+\alpha {\stackrel{\dot{}}{U}}_{\mathrm{4,}m+1}-\alpha {\left(G{U}_3\right)}_{m-1}-{\left(G{U}_3\right)}_m-\alpha {\left(G{U}_3\right)}_{m+1}\\ =p_1\left(U_{\mathrm{4,}m-2}-2U_{\mathrm{4,}m}+U_{\mathrm{4,}m+2}\right)+p_2\left(U_{\mathrm{4,}m-1}-2U_{\mathrm{4,}m}+U_{\mathrm{4,}m+1}\right)\end{array}$ (25)

where ${\stackrel{\dot{}}{U}}_{i,m}=\frac{\partial u_i}{\partial t},i=1,2,3,4$, $p_1=\frac{b}{8h^2}$ and $p_2=\frac{a}{2h^2}$.

This system in a compact form will be written as

$M{\stackrel{\dot{}}{U}}_i=F_i\left(U_1^n\mathrm{,}{U}_2^n\mathrm{,}{U}_3^n\mathrm{,}{U}_4^n\right)\mathrm{,}$ (26)

where $i=\mathrm{1,2,3,4}$,

$M=\frac{1}{11}\left(\begin{array}{cccccccc}11& 2& 0& 0& \cdots & \cdots & \cdots & 0\\ 2& 11& 2& 0& 0& \cdots & \cdots & 0\\ 0& 2& 11& 2& 0& 0& \cdots & 0\\ 0& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & 0& 2& 11& 2\\ 0& \cdots & \cdots & \cdots & 0& 0& 2& 11\end{array}\right)\mathrm{,}$

${\stackrel{\dot{}}{U}}_i=\left[\begin{array}{c}{U}_{i\mathrm{,1}}\\ U_{i\mathrm{,2}}\\ ⋮\\ ⋮\\ U_{i\mathrm{,}M-2}\\ U_{i\mathrm{,}M-1}\end{array}\right]\text{and}{F}_i=\left[\begin{array}{c}{F}_{i\mathrm{,1}}\\ F_{i\mathrm{,2}}\\ ⋮\\ ⋮\\ F_{i\mathrm{,}M-2}\\ F_{i\mathrm{,}M-1}\end{array}\right]$

and

$\begin{array}{c}{F}_{1,m}=-p_1\left(U_{2,m-2}^n-2U_{2,m}^n+U_{2,m+2}^n\right)-p_2\left(U_{2,m-1}^n-2U_{2,m}^n+U_{2,m+1}^n\right)\\ +\alpha G_{m-1}{U}_{2,m-1}^n+G_m{U}_{2,m}^n+\alpha G_{m+1}{U}_{2,m+1}^n\end{array}$ (27)

$\begin{array}{c}{F}_{\mathrm{2,}m}=p_1\left(U_{\mathrm{1,}m-2}^n-2U_{\mathrm{1,}m}^n+U_{\mathrm{1,}m+2}^n\right)+p_2\left(U_{\mathrm{1,}m-1}^n-2U_{\mathrm{1,}m}^n+U_{\mathrm{1,}m+1}^n\right)\\ +\alpha G_{m-1}{U}_{\mathrm{1,}m-1}^n+G_m{U}_{\mathrm{1,}m}^n+\alpha G_{m+1}{U}_{\mathrm{1,}m+1}^n\end{array}$ (28)

$\begin{array}{c}{F}_{\mathrm{3,}m}=-p_1\left(U_{\mathrm{4,}m-2}^n-2U_{\mathrm{4,}m}^n+U_{\mathrm{4,}m+2}^n\right)-p_2\left(U_{\mathrm{4,}m-1}^n-2U_{\mathrm{4,}m}^n+U_{\mathrm{4,}m+1}^n\right)\\ +\alpha G_{m-1}{U}_{\mathrm{4,}m-1}^n+G_m{U}_{\mathrm{4,}m}^n+\alpha G_{m+1}{U}_{\mathrm{4,}m+1}^n\end{array}$ (29)

$\begin{array}{c}{F}_{4,m}=p_1\left(U_{3,m-2}^n-2U_{3,m}^n+U_{3,m+2}^n\right)+p_2\left(U_{3,m-1}^n-2U_{3,m}^n+U_{3,m+1}^n\right)\\ +\alpha G_{m-1}{U}_{3,m-1}^n+G_m{U}_{3,m}^n+\alpha G_{m+1}{U}_{3,m+1}^n\end{array}$ (30)

The solution of (26), using the explicit RK4, can be displayed as follows:

$U_{j\mathrm{,}m}^{n+1}=U_{j\mathrm{,}m}^n+\frac{1}{6}\left[K_{j\mathrm{,1}}+2K_{j\mathrm{,2}}+2K_{j\mathrm{,3}}+K_{j\mathrm{,4}}\right]\mathrm{,}$ (31)

$K_{j\mathrm{,1}}=k{M}^{-1}{F}_{j\mathrm{,}m}\left(U_{\mathrm{1,}m}^n\mathrm{,}{U}_{\mathrm{2,}m}^n\mathrm{,}{U}_{\mathrm{3,}m}^n\mathrm{,}{U}_{\mathrm{4,}m}^n\right)$

$K_{j\mathrm{,2}}=k{M}^{-1}{F}_{j\mathrm{,}m}\left(U_{\mathrm{1,}m}^n+\frac{1}{2}{K}_{\mathrm{1,1}}\mathrm{,}{U}_{\mathrm{2,}m}^n+\frac{1}{2}{K}_{\mathrm{2,1}}\mathrm{,}{U}_{\mathrm{3,}m}^n+\frac{1}{2}{K}_{\mathrm{3,1}}\mathrm{,}{U}_{\mathrm{4,}m}^n+\frac{1}{2}{K}_{\mathrm{4,1}}\right)$

$K_{j\mathrm{,3}}=k{M}^{-1}{F}_{j\mathrm{,}m}\left(U_{\mathrm{1,}m}^n+\frac{1}{2}{K}_{\mathrm{1,2}}\mathrm{,}{U}_{\mathrm{2,}m}^n+\frac{1}{2}{K}_{\mathrm{2,2}}\mathrm{,}{U}_{\mathrm{3,}m}^n+\frac{1}{2}{K}_{\mathrm{3,2}}\mathrm{,}{U}_{\mathrm{4,}m}^n+\frac{1}{2}{K}_{\mathrm{4,2}}\right)$

$K_{j\mathrm{,4}}=k{M}^{-1}{F}_{j\mathrm{,}m}\left(U_{\mathrm{1,}m}^n+K_{\mathrm{1,3}}\mathrm{,}{U}_{\mathrm{2,}m}^n+K_{\mathrm{2,3}}\mathrm{,}{U}_{\mathrm{3,}m}^n+K_{\mathrm{3,3}}\mathrm{,}{U}_{\mathrm{4,}m}^n+K_{\mathrm{4,3}}\right)$

for $j=1,2,3,4$

The proposed scheme is of sixth order in space and fourth order in time. The previous scheme, which applies the Runge-Kutta of order 4 method, is conditionally stable.

3. Numerical Results

This section presents the numerical results for the aimed scheme. This method’s accuracy is tested with by studying its conservation properties and the truncation error using $L_2$ and $L_{\infty }$ norms. Trapezoidal rule is used to calculate the conserved quantities. To reliably estimate the smooth solution of CNLSEs, discrete conservation laws are critical. The accuracy is scaled by using the $L_2$ error norm and $L_{\infty }$ error norm.

${\Vert ER\Vert }_{\infty }=\underset{1<m<N}{\max }\left\{\left|\left(\Vert {\Phi }_1\left(x_m,t_n\right)\Vert -\Vert U_{1,m}^n+i{U}_{2,m}^n\Vert \right)\right|\right\}$

${\Vert ER\Vert }_2={\left[{\displaystyle \underset{m=1}{\overset{N}{\sum }}}{\left(\Vert {\Phi }_1\left(x_m,t_n\right)\Vert -\Vert U_{1,m}^n+i{U}_{2,m}^n\Vert \right)}^2\right]}^{\frac{1}{2}}$

In the following sections, we will analyze the subsequent problems that arose in studying the properties of the obtained scheme.

3.1. Single Soliton

In this experiment, the initial condition is determined by

${\Phi }_1\left(x\mathrm{,0}\right)=\sqrt{\frac{2\alpha }{1+c}}\sec h\left(\sqrt{2\alpha }x\right)\exp i\left\{v_{1}x\right\}\mathrm{,}$

${\Phi }_2\left(x\mathrm{,0}\right)=\sqrt{\frac{2\alpha }{1+c}}\sec h\left(\sqrt{2\alpha }x\right)\exp i\left\{v_{1}x\right\}$

where $\alpha$, c and $v_1$ are constants.

To estimate the numerical solution, we use the following parameters:

$x_0=-50,x_M=50,h=0.1,k=0.001,v_1=0.5,\alpha =1.0,c=1$

In Table 1, we demonstrate the conserved quantities recovered from the proposed scheme. It is obvious that all conserved quantities are nearly conserved. Table 2 shows the accuracy of the method where $L_i\left(u_1\right)=L_i\left(u_3\right)$ and $L_i\left(u_2\right)=L_i\left(u_4\right)$, $i=\infty ,2$. Figure 1 shows the development of a single soliton running to the right with velocity $v_1=0.5$.

3.2. Interaction of Two Solitons

To examine the interaction of two solitons, we exercise the following initial conditions.

Table 1. Single solitons (conserved quantities).

Table 2. The errors of single soliton.

Figure 1. Single soliton.

${\Phi }_1\left(x,0\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}}\sqrt{\frac{2{\alpha }_j}{1+c}}\sec h\left(\sqrt{2{\alpha }_j}{x}_j\right)\exp i\left\{v_{1j}{x}_j\right\},$

${\Phi }_2\left(x,0\right)={\displaystyle \underset{j=1}{\overset{2}{\sum }}}\sqrt{\frac{2{\alpha }_j}{1+c}}\sec h\left(\sqrt{2{\alpha }_j}{x}_j\right)\exp i\left\{v_{1j}{x}_j\right\}$

These initial conditions draw two waves. The complete sum of the masses asymptotically can be shown as $I=\frac{2}{1+c}{\displaystyle \underset{j=1}{\overset{2}{\sum }}}\sqrt{2{\alpha }_j}$.

In this test we choose the parameters

$\begin{array}{l}{x}_0=-50,x_M=50,h=0.1,k=0.001,x_1=-25,x_2=25,\\ v_{11}=0.5,v_{12}=-0.5,{\alpha }_1=0.5,{\alpha }_2=0.5,c=1\end{array}$

This method is very nearly conserves all conserved quantities (see Table 3). In Figure 2, we demonstrate the interaction of the pair waves moving at the same velocities in opposite directions for ${\Phi }_1$ (the same for ${\Phi }_2$ ), which is called head-on interaction. It is recognized that both interacting and moving the interaction zone unaltered in shape and velocities. In Figure 3, we plot the contour of two solitons with the same parameters.

Figure 2. Interaction of two solitons.

Figure 3. Refection of two solitons.

Table 3. Interaction of two solitons (conserved quantities).

4. Conclusion

We have used a sixth-order finite difference method and Runge-Kutta of order 4 method to solve CNLSEs. The RK4 scheme is fourth order in time and sixth order in space, and it remains conditionally stable. Obviously, from the results, the scheme obtains fourth-order instantly without the need for an extra procedure to improve the accuracy in the temporal direction.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ismail, M.S. (2008) Numerical Solution of Coupled Nonlinear Schrödinger Equation by Galerkin Method. Mathematics and Computers in Simulation, 78, 532-547. https://doi.org/10.1016/j.matcom.2007.07.003
[2]	Ismail, M.S. (2010) Collocation Method for Numerical Solution of Coupled Nonlinear Schrödinger Equation. In: AIP Conference Proceedings, American Institute of Physics, College Park, 1429-1432. https://doi.org/10.1063/1.3498015
[3]	Ismail, M.S. (1996) Finite Difference Method with Cubic Spline for Solving Nonlinear Schrödinger Equation. International Journal of Computer Mathematics, 62, 101-112. https://doi.org/10.1080/00207169608804528
[4]	Muslu, G.M. and Erbay, H.A. (2005) Higher-Order Split-Step Fourier Schemes for the Generalized Nonlinear Schrödinger Equation. Mathematics and Computers in Simulation, 67, 581-595. https://doi.org/10.1016/j.matcom.2004.08.002
[5]	Sanz-Serna, J.M. and Verwer, J.G. (1986) Conservative and Nonconservative Schrödinger Equation. IMA Journal of Numerical Analysis, 6, 25-42. https://doi.org/10.1093/imanum/6.1.25
[6]	Shamerdan, A.B. (1990) The Numerical Treatment of the Nonlinear Schrödinger Equation. Computers & Mathematics with Applications, 19, 67-73. https://doi.org/10.1016/0898-1221(90)90195-P
[7]	Sheng, Q., Khaliq, A.Q.M. and Al-Said, E.A. (2001) Solving the Generalized Nonlinear Schrödinger Equation via Quartic Spline Approximation. Journal of Computational Physics, 166, 400-417. https://doi.org/10.1006/jcph.2000.6668
[8]	Ismail, M.S. and Taha, T.R. (2007) A Linearly Implicit Conservative Scheme for the Coupled Nonlinear Schrödinger Equation. Mathematics and Computers in Simulation, 74, 302-311. https://doi.org/10.1016/j.matcom.2006.10.020
[9]	Ismail, M.S. and Alamri, S.Z. (2004) Highly Accurate Finite Difference Method for Coupled Nonlinear Schrödinger Equation. International Journal of Computer Mathematics, 81, 333-351. https://doi.org/10.1080/00207160410001661339
[10]	Ismail, M.S. and Taha, T.R. (2001) Numerical Simulation of Coupled Nonlinear Schrödinger Equation. Mathematics and Computers in Simulation, 56, 547-562. https://doi.org/10.1016/S0378-4754(01)00324-X
[11]	Ismail, M.S. and Taha, T.R. (2000) A Finite Element Solution for the Coupled Schrödinger Equation. System, 1, u1.
[12]	Sonnier, W.J. and Christov, C.I. (2005) Strong Coupling of Schrödinger Equations: Conservative Scheme Approach. Mathematics and Computers in Simulation, 69, 514-525. https://doi.org/10.1016/j.matcom.2005.03.016
[13]	Aydin, A. and Karasozen, B. (2009) Multi-Symplectic Integration of Coupled Nonlinear Schrödinger System with Soliton Solutions. International Journal of Computer Mathematics, 86, 864-882.
[14]	Aydin, A. and Karasozen, B. (2007) Symplectic and Multisymplectic Methods for the Coupled Nonlinear Schrödinger Equations with Periodic Solutions. Computer Physics Communications, 177, 566-583. https://doi.org/10.1016/j.cpc.2007.05.010
[15]	Sun, J.Q., Gu, X.Y. and Ma, Z.Q. (2004) Numerical Study of the Soliton Waves of the Coupled Nonlinear Schrödinger System. Physica D, 196, 311-328. https://doi.org/10.1016/j.physd.2004.05.010
[16]	Sun, J.Q. and Qin, M.Z. (2003) Multi-Symplectic Methods for the Coupled 1D Nonlinear Schrödinger System. Computer Physics Communications, 155, 221-235. https://doi.org/10.1016/S0010-4655(03)00285-6
[17]	Wang, H. (2005) Numerical Studies on the Split Step Finite Difference Method for the Nonlinear Schrödinger Equations. Applied Mathematics and Computation, 170, 17-35. https://doi.org/10.1016/j.amc.2004.10.066
[18]	Xu, Y. and Shu, C. (2005) Local Discontinuous Galerkin Method for Nonlinear Schrödinger Equations. Journal of Computational Physics, 205, 72-97. https://doi.org/10.1016/j.jcp.2004.11.001
[19]	de Frutos, J. and Sanz-Serna, J.M. (1991) An Easily Implementable Fourth Order Method for the Time Integration of Waves Problems. Report 1991/2, Universidad De Valladolid, Valladolid.
[20]	Griffiths, D.F., Mitchell, A.R. and Li Morris, J. (1982) A Numerical Study of the Nonlinear Schrödinger Equation. NA/52, University of Dundee, Dundee.
[21]	Xu, Q.-B. and Chang, Q.-S. (2010) New Numerical Methods for the Coupled Nonlinear Schrödinger Equations. Acta Mathematicae Applicatae Sinica, English Series, 26, 205-218. https://doi.org/10.1007/s10255-007-7098-2
[22]	Ismail, M.S. (2008) A Fourth-Order Explicit Schemes for the Coupled Nonlinear Schrödinger Equation. Applied Mathematics and Computation, 196, 273-284. https://doi.org/10.1016/j.amc.2007.05.059
[23]	Li, J. and Chen, Y. (2008) Computational Partial Differential Equations Using MATLAB. Taylor & Francis Group, Abingdon-on-Thames. https://doi.org/10.1201/9781420089059