Hao Ji¹, Yinghong Xu^,1,∗, Chaoqing Dai², Lipu Zhang³¹Department of Mathematics, Zhejiang Sci-Tech University, Hangzhou, 310018, China
²Department of Physics, Zhejiang A&F University, Hangzhou, 311300, China
³College of Media Engineering, Communication University of Zhejiang, Hangzhou, 310018, China

First author contact: ∗Author to whom any correspondence should be addressed.
Received:2021-07-18Revised:2021-08-31Accepted:2021-09-24Online:2021-10-26


Abstract
We consider the (2+1)-dimensional nonlinear Schrödinger equation with power-law nonlinearity under the parity-time-symmetry potential by using the Crank–Nicolson alternating direction implicit difference scheme, which can also be used to solve general boundary problems under the premise of ensuring accuracy. We use linear Fourier analysis to verify the unconditional stability of the scheme. To demonstrate the effectiveness of the scheme, we compare the numerical results with the exact soliton solutions. Moreover, by using the scheme, we test the stability of the solitons under the small environmental disturbances.
Keywords： nonlinear Schrödinger equation;localized spatial solitons;PT-symmetric potential;ADI difference scheme;stability


PDF (1889KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite
Cite this article
Hao Ji, Yinghong Xu, Chaoqing Dai, Lipu Zhang. Propagation of local spatial solitons in power-law nonlinear PT-symmetric potentials based on finite difference. Communications in Theoretical Physics, 2021, 73(12): 125002- doi:10.1088/1572-9494/ac29b6

1. Introduction

An optical spatial soliton has important applications in the connection of integrated optical components and information transmission, electro-optical information processing, optical computing, optical storage and optical display [1–13]. The propagation of solitary waves in strongly nonlocal nonlinear media satisfies the strongly nonlocal nonlinear Schrödinger equation (NLSE). Therefore, the study of the NLSE are of great significance for studying the propagation of optical spatial soliton.

Recently, parity-time (PT)-symmetric potentials were introduced into the field of optics [14, 15], various nonlinear-localized structures in PT-symmetric potentials have been extensively studied, especially in the structure of the NLSE [16, 17]. Dai et al [18] and Xu and Dai [19] studied the dynamical behaviors of nonautonomous solitons in PT-symmetric potentials. Wang et al [20] studied the stable localized spatial solitons in PT-symmetric potentials with power-law nonlinearity. Dai et al [21] studied the dynamics of light bullets in inhomogeneous cubic-quintic-septimal nonlinear media with PT-symmetric potentials. Dai et al [22] studied the reconstruction of stability for Gaussian spatial solitons in quintic-septimal nonlinear materials under PT-symmetric potentials. Zhou et al [23] also studied the solving forward and inverse problems of the logarithmic NLSE with PT-symmetric harmonic potential via deep learning.

Considering that the analytical solutions of the NLSE are complex and even unsolvable, numerical method [24–32] offers another choice for the NLSE. Fourier spectral method is generally used to solve NLSE in optics. It ensures high accuracy and low computation, but it also has the disadvantage that it is only suitable for differential equations under periodic boundary conditions. In contrast, the generality of the finite difference method is much better. In this paper, we use the Crank–Nicolson alternating direction implicit (ADI) scheme to investigate the (2+1)-dimensional NLSE with power-law nonlinearity under the PT-symmetry potentials. To reduce the influence of the initial value on the solution, we also discuss the stability of the scheme. Furthermore, to study the effectiveness of the scheme in terms simulating the evolution of optical soliton, we compare the numerical results of the ADI scheme with the exact soliton solutions. Finally, we test the stability of the solitons under the small environmental disturbances by using the ADI scheme.

This paper is organized as follows. In section 2, we introduce the (2+1)-dimensional NLSE with power-law nonlinearity under the PT-symmetry potentials and its exact soliton solutions. In section 3, we give the derivation process and stability analysis of the second-order ADI scheme. Numerical experiment results are presented in section 4. We study the influence of small environmental disturbances on the evolution of solitons in section 5. We give some conclusions in section 6.

2. The NLSE and its exact spatial soliton solutions

The propagation of an optical beam in a nonlinear medium of non-Kerr index, which is perturbed along the transverse x-and y-directions by a complex profile with n_p = n₀[1 + δn_R(x, y) + iδn_I(x, y)], is governed by the following (2+1)-dimensional NLSE:(1)$\begin{eqnarray}\begin{array}{l}{\rm{i}}\displaystyle \frac{\partial u}{\partial z}+\beta \left(\displaystyle \frac{{\partial }^{2}u}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}u}{\partial {y}^{2}}\right)+{\gamma }_{m}| u{| }^{2m}u\\ +\,[V(x,y)+\ {\rm{i}}W(x,y)]u=0,\\ (x,y,z)\in {\rm{\Omega }}\times (0,Z],\end{array}\end{eqnarray}$where ω is a rectangular domain in ${{\mathbb{R}}}^{2}$, β is the coefficient of the diffraction, γ_m for m = 1, ⋯ , n stand for the nonlinearities of orders up to 2n + 1, u(x, y, z) is the complex envelope of the electrical field, z is the longitudinal coordinate and x, y are the transverse coordinates. Functions V(x, y) and W(x, y) are the real and imaginary components of the complex PT-symmetric potential and correspond to the index guiding and the gain or loss distribution of the optical potential, respectively. V and W are even and odd functions which regard to {x, y}.

Wang et al [20] considered the following PT-symmetric potentials(2)$\begin{eqnarray}\begin{array}{rcl}V(x,y)&=&{V}_{0}[{{\rm{sech}} }^{2}(x)+{{\rm{sech}} }^{2}(y)]\\ & & +{V}_{1}{{\rm{sech}} }^{2}(x){{\rm{sech}} }^{2}(y)\\ & & +{V}_{2}[{{\rm{sech}} }^{2k}(x)+{{\rm{sech}} }^{2k}(y)],\end{array}\end{eqnarray}$(3)$\begin{eqnarray}\begin{array}{rcl}W(x,y)&=&{W}_{0}[{{\rm{sech}} }^{k}(x)\tanh (x)\\ & & +{{\rm{sech}} }^{k}(y)\tanh (y)],\end{array}\end{eqnarray}$where V₀, V₁, V₂ and W₀ are real, and k is the competing parameter. Setting V₀ = (m + 1)β/m², ${V}_{2}={m}^{2}{W}_{0}^{2}/[{\left({km}+2\right)}^{2}\beta ],$ m = 1, 2 and the exact soliton solutions of equation (1) can be obtained as

Case 1. 2D Hyperbolic Scarf II potential (k = 1)(4)$\begin{eqnarray}\begin{array}{l}u(x,y,z)={\left[\sqrt{-\displaystyle \frac{{V}_{1}}{\gamma }}{\rm{sech}} (x){\rm{sech}} (y)\right]}^{\tfrac{1}{m}}\\ \times \,\exp \left({\rm{i}}\displaystyle \frac{2\beta }{{m}^{2}}z+{\rm{i}}\displaystyle \frac{{{mW}}_{0}}{(m+2)\beta }\right.\\ \left.\times \ \{\arctan [\sinh (x)]+\arctan [\sinh (y)]\}\right).\end{array}\end{eqnarray}$

Case 2. 2D PT-symmetric potential with k = 2(5)$\begin{eqnarray}\begin{array}{l}u(x,y,z)={\left[\sqrt{-\displaystyle \frac{{V}_{1}}{\gamma }}{\rm{sech}} (x){\rm{sech}} (y)\right]}^{\tfrac{1}{m}}\\ \times \,\exp \left\{{\rm{i}}\displaystyle \frac{2\beta }{{m}^{2}}z+{\rm{i}}\displaystyle \frac{{{mW}}_{0}}{(m+2)\beta }\times [\tanh (x)\right.\\ \left.+\,\tanh (y)]\right\}.\end{array}\end{eqnarray}$

3. Second-order ADI scheme for the NLSE and its stability analysis

In the following, we will use the second-order ADI scheme to solve the equation (1), and verify the effectiveness of the second-order ADI scheme by comparing the numerical results with the exact soliton solutions (4)–(5).

To apply for a wider range of nonlinear scenarios, we replace the nonlinear term γ_m∣u∣^2m in (1) with f(∣u∣²) and obtain(6)$\begin{eqnarray}\begin{array}{l}{\rm{i}}\displaystyle \frac{\partial u}{\partial z}+\beta \left(\displaystyle \frac{{\partial }^{2}u}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}u}{\partial {y}^{2}}\right)\\ +\ f(| u{| }^{2})u+[V(x,y)+{\rm{i}}W(x,y)]u=0,\\ (x,y,z)\in {\rm{\Omega }}\times (0,Z],\end{array}\end{eqnarray}$the initial and boundary conditions are(7)$\begin{eqnarray}u{| }_{z=0}=\varphi (x,y),\end{eqnarray}$(8)$\begin{eqnarray}u{| }_{{\rm{\Gamma }}}=g(x,y,z),\end{eqnarray}$where f, φ and g are smooth functions. We define v = V(x, y) + iW(x, y) and suppose ω = (a₁, b₁) × (a₂, b₂). Let Γ denotes the boundary of ω.

3.1. Discretization and notations

The following is preliminaries before deriving the difference Crank–Nicolson implicit finite difference scheme.

We discretize the domain $\{(x,y,z)| (x,y,z)\in \overline{{\rm{\Omega }}}\,\times (0,Z]\}$ into grids which is described by the set {(x_p, y_q, z_j)} of nodes, in which x_p = a₁ + ph_x, y_q = a₂ + qh_y, p = 0, 1, ⋯ , M_x, q = 0, 1, ⋯ , M_y, and z_j = jτ, j = 0, 1, ⋯ , N = Z/τ, where h_x, h_y and τ are the uniform step size in the spatial. Let ${u}_{p,q}^{j}=u({x}_{p},{y}_{q},{z}_{j}),$ ω_h = {(x_p, y_q)∣p = 1, ⋯ , M_x − 1, q = 1, ⋯ , M_y − 1}, Γ_h denote the set of nodes on Γ, and let ${\overline{{\rm{\Omega }}}}_{h}={{\rm{\Omega }}}_{h}\cup {{\rm{\Gamma }}}_{h}$. We use the following notations for difference operators:(9)$\begin{eqnarray}{\delta }_{z}{u}_{{pq}}^{j+\tfrac{1}{2}}=\displaystyle \frac{1}{\tau }({u}_{{pq}}^{j+1}-{u}_{{pq}}^{j}),\end{eqnarray}$(10)$\begin{eqnarray}{\delta }_{x}^{2}{u}_{{pq}}^{j}=\displaystyle \frac{1}{{h}_{x}^{2}}({u}_{p-1,q}^{j}-2{u}_{{pq}}^{j}+{u}_{p+1,q}^{j}),\end{eqnarray}$(11)$\begin{eqnarray}{\delta }_{y}^{2}{u}_{{pq}}^{j}=\displaystyle \frac{1}{{h}_{y}^{2}}({u}_{p,q-1}^{j}-2{u}_{{pq}}^{j}+{u}_{p,q+1}^{j}).\end{eqnarray}$

3.2. The establishment of Crank–Nicolson implicit finite difference scheme

Let $h=\max \{{h}_{x},{h}_{y}\}$, by using the Taylor expansion in (10)–(11), we obtain(12)$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}{u}_{{pq}}^{j}+\displaystyle \frac{{\partial }^{2}}{\partial {y}^{2}}{u}_{{pq}}^{j}\\ =\ {\delta }_{x}^{2}{u}_{{pq}}^{j}+{\delta }_{y}^{2}{u}_{{pq}}^{j}+O({h}^{2}).\end{array}\end{eqnarray}$Substituting (12) into (6) and using the Taylor expansion one more time in (9) yields(13)$\begin{eqnarray}\begin{array}{l}{\rm{i}}{\delta }_{z}{u}_{{pq}}^{j+\tfrac{1}{2}}+\beta {\delta }_{x}^{2}{u}_{{pq}}^{j+\tfrac{1}{2}}+\beta {\delta }_{y}^{2}{u}_{{pq}}^{j+\tfrac{1}{2}}\\ +\ (f(| {u}_{{pq}}^{j+\tfrac{1}{2}}{| }^{2})+{v}_{{pq}}){u}_{{pq}}^{j+\tfrac{1}{2}}\\ +\ O({h}^{2}+{\tau }^{2})=0,\end{array}\end{eqnarray}$where ${u}_{{pq}}^{j+\tfrac{1}{2}}=({u}_{{pq}}^{j+1}+{u}_{{pq}}^{j})/2.$

Since the nonlinear term $f(| {u}_{{pq}}^{j+\tfrac{1}{2}}{| }^{2})$ in (13) will add the computational cost, we construct a linearized difference scheme, i.e. applying extrapolation technology to the real coefficient $| {u}_{{pq}}^{j+\tfrac{1}{2}}{| }^{2}$ of the nonlinear term in (13). Setting(14)$\begin{eqnarray}| {u}_{{pq}}^{j+\tfrac{1}{2}}{| }^{2}=\displaystyle \frac{3}{2}| {u}_{{pq}}^{j}{| }^{2}-\displaystyle \frac{1}{2}| {u}_{{pq}}^{j-1}{| }^{2}+O({\tau }^{2}),\end{eqnarray}$we rewrite (13) as the Crank–Nicolson implicit finite difference scheme with accuracy O(h² + τ²)(15)$\begin{eqnarray}\begin{array}{l}{\rm{i}}{\delta }_{z}{u}_{{pq}}^{j+\tfrac{1}{2}}+\beta {\delta }_{x}^{2}{u}_{{pq}}^{j+\tfrac{1}{2}}+\beta {\delta }_{y}^{2}{u}_{{pq}}^{j+\tfrac{1}{2}}\\ +\ (f\left(\displaystyle \frac{3}{2}| {u}_{{pq}}^{j}{| }^{2}-\displaystyle \frac{1}{2}| {u}_{{pq}}^{j-1}{| }^{2}\right)+{v}_{{pq}}){u}_{{pq}}^{j+\tfrac{1}{2}}=0.\end{array}\end{eqnarray}$

3.3. ADI scheme for Crank–Nicolson implicit finite difference scheme

In the following, we will rewrite (15) as its ADI scheme.

We define(16)$\begin{eqnarray}{\omega }_{{pq}}^{j}=f\left(\displaystyle \frac{3}{2}| {u}_{{pq}}^{j}{| }^{2}-\displaystyle \frac{1}{2}| {u}_{{pq}}^{j-1}{| }^{2}\right)+{v}_{{pq}}.\end{eqnarray}$Substituting (16) into (15), we transform (15) into(17)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j+1}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}\\ +\ \displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}({u}_{{pq}}^{j+1}+{u}_{{pq}}^{j}).\end{array}\end{eqnarray}$Since(18)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j+1}\\ =\ \left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j+1}+O({\tau }^{2}),\end{array}\end{eqnarray}$and(19)$\begin{eqnarray}\begin{array}{l}\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}+O({\tau }^{2}),\end{array}\end{eqnarray}$we ignore the truncation error O(τ²) in (18) and (19), the equation (17) can be expressed as(20)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j+1}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}\\ +\ \displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}({u}_{{pq}}^{j+1}+{u}_{{pq}}^{j}).\end{array}\end{eqnarray}$

To replaces the solution of two-dimensional problems by sequences of one-dimensional cases, we introduce a new intermediate variable u^(*) into (20), and get the ADI scheme for the NLSE (6), as(21)$\begin{eqnarray}\begin{array}{rcl}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){u}_{{pq}}^{j(* )}&=&\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}+\displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}{u}_{{pq}}^{j},\end{array}\end{eqnarray}$(22)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j+1}-\displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}{u}_{{pq}}^{j+1}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){u}_{{pq}}^{j(* )},\end{array}\end{eqnarray}$where j ≥ 2.

The boundary of intermediate variable u^(*) in (21)–(22) is$\begin{eqnarray*}\begin{array}{l}{u}_{0q}^{j(* )}=\displaystyle \frac{1}{2}({g}_{0q}^{j}+{g}_{0q}^{j+1})-\displaystyle \frac{{\rm{i}}\tau \beta }{4{h}_{y}^{2}}{\delta }_{y}^{2}({g}_{0q}^{j+1}-{g}_{0q}^{j})\\ -\ \displaystyle \frac{{\rm{i}}\tau }{4}\left[f\left(\displaystyle \frac{3}{2}| {g}_{0q}^{j}{| }^{2}-\displaystyle \frac{1}{2}| {g}_{0q}^{j-1}{| }^{2}\right)+{v}_{0q}\right]({g}_{0q}^{j+1}-{g}_{0q}^{j}),\\ {u}_{{M}_{x}q}^{j(* )}=\displaystyle \frac{1}{2}({g}_{{M}_{x}q}^{j}+{g}_{{M}_{x}q}^{j+1})-\displaystyle \frac{{\rm{i}}\tau \beta }{4{h}_{y}^{2}}{\delta }_{y}^{2}({g}_{{M}_{x}q}^{j+1}-{g}_{{M}_{x}q}^{j})\\ -\ \displaystyle \frac{{\rm{i}}\tau }{4}\left[f\left(\displaystyle \frac{3}{2}| {g}_{{M}_{x}q}^{j}{| }^{2}-\displaystyle \frac{1}{2}| {g}_{{M}_{x}q}^{j-1}{| }^{2}\right)+{v}_{{M}_{x}q}\right]({g}_{{M}_{x}q}^{j+1}-{g}_{{M}_{x}q}^{j}).\end{array}\end{eqnarray*}$

3.4. Stability analysis for Crank–Nicolson implicit finite difference scheme

We set ${\omega }_{{jk}}^{n}=c$, where c is a constant, the scheme (21)–(22) can be rewritten as(23)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left[\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right)-\displaystyle \frac{{\rm{i}}\tau }{2}c\right]{u}_{{pq}}^{j+1}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left[\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right)+\displaystyle \frac{{\rm{i}}\tau }{2}c\right]{u}_{{pq}}^{j}.\end{array}\end{eqnarray}$

If the numerical solution obtained in the scheme (21)–(22) can be expressed by a Fourier series, whose typical term is(24)$\begin{eqnarray}{u}_{{pq}}^{j}={\xi }^{n}\exp ({\rm{i}}{\alpha }_{1}{{ph}}_{x}+{\rm{i}}{\alpha }_{2}{{qh}}_{y}),\end{eqnarray}$where ξ^j is the amplitude at longitudinal spatial level j, α₁, α₂ are the wave numbers in the x and y directions. Substituting (24) into (23), we obtain(25)$\begin{eqnarray}\xi =\displaystyle \frac{\left(1-\tfrac{2{\rm{i}}\tau \beta }{{h}_{x}^{2}}{\sin }^{2}\tfrac{{\alpha }_{1}{h}_{x}}{2}\right)\left(1-\tfrac{2{\rm{i}}\tau \beta }{{h}_{y}^{2}}{\sin }^{2}\tfrac{{\alpha }_{2}{h}_{y}}{2}-c\tfrac{{\rm{i}}\tau }{2}\right)}{1+\tfrac{2{\rm{i}}\tau \beta }{{h}_{x}^{2}}{\sin }^{2}\tfrac{{\alpha }_{1}{h}_{x}}{2}\left(1+\tfrac{2{\rm{i}}\tau \beta }{{h}_{y}^{2}}{\sin }^{2}\tfrac{{\alpha }_{2}{h}_{y}}{2}+c\tfrac{{\rm{i}}\tau }{2}\right)}.\end{eqnarray}$Since the numerator and denominator of ξ are complex conjugates to each other, we can conclude that ∣ξ∣ = 1, which means the scheme (21)–(22) is unconditionally stable.

4. Numerical experiments

In this section, we use two known solitons (4)–(5) in the case of self-focusing (γ_m > 0) to verify the validity of the ADI scheme to simulate solitons evolution. All the numerical experiments are obtained by Matlab 2017a on a Thunderobot computer with i7-6700HQ CPU and 8 Gbyte of memory.

In PT-symmetric potential (2) and (3), we consider equation (1) with (x, y) ∈ [−2, 2] × [1, 5], β = γ_m = 1, and V1 = −15. When k = 1 and 2, the exact solitons are Cases 1 and 2 in the section 2, respectively. The boundary and initial values are directly taken from (4)–(5).

We set ${h}_{x}={h}_{y}=\tau =\tfrac{1}{80}$. Figures 1 and 2 present the numerical solutions and absolute errors for case 1 at two longitudinal spatial points, z = 1 and 2. Figures 3 and 4 present the numerical solutions and absolute errors for case 2 at two longitudinal spatial points, z = 1 and 2. In figures 1–4, m are 1, 2, 1 and 2, and W₀ are 0.012, 0.12, 0.066 and 0.022, respectively.

Figure 1.

New window|Download| PPT slide
Figure 1.Numerical solutions (left) and the absolute errors (right) of ADI for Case 1 at: (a) z = 1, (b) z = 2.

Figure 2.

New window|Download| PPT slide
Figure 2.Numerical solutions (left) and the absolute errors (right) of ADI for Case 1 at: (a) z = 1, (b) z = 2.

Figure 3.

New window|Download| PPT slide
Figure 3.Numerical solutions (left) and the absolute errors (right) of ADI for Case 2 at: (a) z = 1, (b) z = 2.


The results of figures 1–4 show that, the error levels are all less than or equal to 10⁻⁴, the ADI scheme is no oscillation.

Figure 4.

New window|Download| PPT slide
Figure 4.Numerical solutions (left) and the absolute errors (right) of ADI for case 2 at: (a) z = 1, (b) z = 2.


If the soliton itself is stable, as long as the appropriate value τ is selected in the numerical scheme and propagated to z = 20, 200, etc, the soliton solution is still stable. However, if it is not selected properly, the error will accumulate with the increase of z, and the numerical solution may be unstable.

We set h = τ, and z = 1. Tables 1 and 2 give the L^∞-norm errors, convergence rate and CPU time of the ADI scheme.


Table 1.
Table 1.Maximum errors and convergence orders of the ADI scheme for Case 1 at z = 1.

	Figure 1(a) m = 1, W0 = 0.012			Figure 2(a) m = 2, W0 = 0.120
h, τ	L∞error	Order	CPU time	L∞error	Order	CPU time
$\tfrac{1}{20}$	5.4174e–4	—	4.501	1.3605e–4	—	4.697
$\tfrac{1}{40}$	1.3592e–4	1.995	39.651	3.4048e–5	1.998	42.549
$\tfrac{1}{80}$	3.4318e–5	1.986	310.942	8.5199e–6	1.999	324.669

New window|CSV

The results of tables 1 and 2 show that, the dense griding will lead to better the effect of numerical approximation. The convergence rate is close to 2, which also verifies the accuracy of the ADI scheme is O(h² + τ²).


Table 2.
Table 2.Maximum errors and convergence orders of the ADI scheme for Case 2 at z = 1.

	Figure 3(a) m = 1, W0 = 0.066			Figure 4(a) m = 2, W0 = 0.022
h, τ	L∞error	Order	CPU time	L∞error	Order	CPU time
$\tfrac{1}{20}$	5.3701e–4	—	4.521	8.7541e–5	—	4.754
$\tfrac{1}{40}$	1.3471e–4	1.995	42.528	2.1691e–5	2.013	45.050
$\tfrac{1}{80}$	3.4013e–5	1.986	324.258	5.3926e–6	2.008	345.450

New window|CSV

5. A small environmental disturbance to the evolution of solitons

In real world practical problems, the environment will affect the evolution of soliton, under this situation we introduce white noise to simulate the limited disturbances. In the following, we will use the ADI scheme to test the stability of the soliton solutions (4)–(5) of the equation (1) under a small perturbations.

5.1. Adding initial perturbations

We use the ADI scheme to rerun the soliton solutions (4)–(5) with a white noise under 5% initial fluctuation. Figures 5(a) and 6(a) reveal the stable evolution of spatial solitons until z = 2 in the self-focusing cubic nonlinear case (m = 1). Figures 5(b) and 6(b) reveal the unstable evolution of spatial solitons (4) And (5) in the self-focusing quintic nonlinear case (m = 2). The result of figures 5 and 6 show that in the case of self-focusing quintic nonlinearity, the solitons are seriously affected by the initial disturbance.

Figure 5.

New window|Download| PPT slide
Figure 5.Numerical solutions (left) and the absolute errors (right) of ADI at z = 2. A white noise with 5% fluctuation is added to initial solution. The parameters of (a) and (b) are the same as those of figures 1 and 2 respectively.

Figure 6.

New window|Download| PPT slide
Figure 6.Numerical solutions (left) and the absolute errors (right) of ADI at z = 2. A white noise with 5% fluctuation is added to initial solution. The parameters of (a) and (b) are the same as those of figures 3 and 4 respectively.

5.2. Adding boundary perturbations

We use the ADI scheme to rerun soliton the solutions (4)–(5) with a white noise under 1% boundary fluctuation. Figures 7 and 8 show that when the boundary perturbations are added, the solitons will evolve unsteadily at m = 1 and 2 until z = 2. Moreover, in the case of cubic nonlinearity, the soliton evolution is less affected by the boundary disturbance.

Figure 7.

New window|Download| PPT slide
Figure 7.Numerical solutions (left) and the absolute errors (right) of ADI at z = 2. A white noise with 1% fluctuation is added to boundary solution. The parameters of (a) and (b) are the same as those of figures 1 and 2 respectively.

Figure 8.

New window|Download| PPT slide
Figure 8.Numerical solutions (left) and the absolute errors (right) of ADI at z = 2. A white noise with 1% fluctuation is added to boundary solution. The parameters of (a) and (b) are the same as those of figures 1 and 2 respectively.

6. Conclusion

We use the ADI scheme to solve the NLSE with power-law nonlinearity under the PT-symmetry potential, and we also use linear Fourier analysis to verify the unconditional stability of the scheme. Numerical experiments show that the Crank–Nicolson ADI scheme is effective in terms of accuracy. Moreover, we find that in the case of self-focusing cubic nonlinearity, a small initial perturbation keeps the spatial soliton stable, but in the case of self-focusing quintic nonlinearity, a small initial perturbation leads to the spatial soliton unstable, in the case of self-focusing cubic or quintic nonlinearity, a small boundary perturbation leads to the spatial soliton unstable.

Appendix

Remark 1 is a derivation of (20) to (21)–(22). When j = 1, u¹ cannot be calculated by (21)–(22), hence, we use remark 2 to calculate u¹.

Remark 1. By eliminating the intermediate variables in (21) and (22), we have$\begin{eqnarray*}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j+1}\\ =\ \left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left[\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){u}_{{pq}}^{j(* )}\right.\\ \left.+\ \displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}{u}_{{pq}}^{j+1}\right]\\ =\ \left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){u}_{{pq}}^{j(* )}\\ +\ \displaystyle \frac{{\rm{i}}\tau }{2}(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}){\omega }_{{pq}}^{j}{u}_{{pq}}^{j\,+\,1}\\ =\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left[\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}+\displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}{u}_{{pq}}^{j}\right]\\ +\ \displaystyle \frac{{\rm{i}}\tau }{2}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){\omega }_{{pq}}^{j}{u}_{{pq}}^{j\,+\,1}\\ =\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right)\left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{j}\\ +\ \displaystyle \frac{{\rm{i}}\tau }{2}{\omega }_{{pq}}^{j}({u}_{{pq}}^{j}+{u}_{{pq}}^{j+1})+O({\tau }^{2}),\end{array}\end{eqnarray*}$which obtain (20).

Remark 2. We use the following ADI linear iterative algorithm to calculate u¹:(26)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){\left({u}_{{pq}}^{0(* )}\right)}^{s+1}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){u}_{{pq}}^{0}+\displaystyle \frac{{\rm{i}}\tau }{2}{\left({\omega }_{{pq}}^{1}\right)}^{s}{u}_{{pq}}^{0},\end{array}\end{eqnarray}$(27)$\begin{eqnarray}\begin{array}{l}\left(I-\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{y}^{2}\right){\left({u}_{{pq}}^{1}\right)}^{s+1}-\displaystyle \frac{{\rm{i}}\tau }{2}{\left({\omega }_{{pq}}^{1}\right)}^{s}{\left({u}_{{pq}}^{1}\right)}^{s+1}\\ =\ \left(I+\displaystyle \frac{{\rm{i}}\tau \beta }{2}{\delta }_{x}^{2}\right){\left({u}_{{pq}}^{0(* )}\right)}^{s+1}.\end{array}\end{eqnarray}$Letting ${\omega }_{{pq}}^{1}=f\left(\tfrac{1}{2}| {u}_{{pq}}^{1}{| }^{2}+\tfrac{1}{2}| {u}_{{pq}}^{0}{| }^{2}\right)+{v}_{{pq}}$ as the initial value of iteration and substituting it into (26)–(27), we obtain the u¹.

Reference View Option By original order
By published year
By cited within times
By Impact factor

[1]Mollenauer L F Stolen R H Gordon J P 1980 Experimental observation of picosecond pulse narrowing and solitons in optical fibers
Phys. Rev. Lett. 45 1095
DOI:10.1103/PhysRevLett.45.1095 [Cited within: 1]

[2]Stegeman G I Mamyshev P V Torruellas W E Villeneuve A Aitchison J S 1995 Photonic applications of spatial soliton switches
Proc. SPIE 2481 270
DOI:10.1117/12.212720

[3]Man W S Demokan M S Tam H 1996 All-optical routing switch using a fiber nonlinear directional coupler and controlling solitons
J. Lightwave Technol. 14 2793
DOI:10.1109/50.545799

[4]Wu Y D 1989 New all-optical switch based on the spatial soliton repulsion
Opt. Express 14 4005
DOI:10.1364/OE.14.004005

[5]Islam M N 2006 Ultrafast all-optical logic gates based on soliton trapping in fibers
Opt. Lett. 14 1257
DOI:10.1364/OL.14.001257

[6]Islam M N Soccolich C E Miller D A B 1990 Ultrafast all-optical logic gates based on soliton trapping in fibers
Opt. Lett. 15 909
DOI:10.1364/OL.15.000909

[7]Bigo S Leclerc O Desurvire E 1997 All-optical fiber signal processing and regeneration for soliton communications
IEEE J. Sel. Top. Quantum Electron. 3 1208
DOI:10.1109/2944.658596

[8]Boscolo S Turitsyn S K Blow K J 2008 Nonlinear loop mirror-based all-optical signal processing infiber-optic communications
Opt. Fiber Technol. 14 299
DOI:10.1016/j.yofte.2008.01.003

[9]Tchier F Aslan E C Inc M 2017 Nanoscale waveguides in optical metamaterials: Jacobi elliptic funtion solutions
J. Nanoelectron. Optoelectron. 12 526
DOI:10.1166/jno.2017.2041

[10]Aslan E C Tchier F Inc M 2017 On optical solitons of the Schrödinger–Hirota equation with power law nonlinearity in optical fibers
Superlattices Microstruct. 105 48
DOI:10.1016/j.spmi.2017.03.014

[11]Inc M Aliyu A I Yusuf A Baleanu D 2018 Gray optical soliton, linear stability analysis and conservation laws via multipliers to the cubic nonlinear Schrödinger equation
Optik 164 472
DOI:10.1016/j.ijleo.2018.02.080

[12]Inc M Aliyu A I Yusuf A Baleanu D 2020 Optical solitons and modulation instability analysis to the quadratic-cubic nonlinear Schrödinger equation
Nonlinear Anal.-Model. 24 20
DOI:10.15388/NA.2019.1.2

[13]Korprinar Z Inc M Alishonrani A S Baleanu D 2020 The deterministic and stochastic solutions of the Schrödinger equation with time conformable derivative in birefrigent fibers
AIMS Math. 5 2326
DOI:10.3934/math.2020154 [Cited within: 1]

[14]Zhong W P Belić M R Huang T W 2012 Two-dimensional accessible solitons in PT-symmetric potentials
Nonlinear Dyn. 70 2027
DOI:10.1007/s11071-012-0594-9 [Cited within: 1]

[15]Musslimani Z H Markris K G El-Ganainy R Christodoulides D N 2008 Optical solitons in PT periodic potentials
Phys. Rev. Lett. 100 030402
DOI:10.1103/PhysRevLett.100.030402 [Cited within: 1]

[16]Shi Z W Jiang X J Zhu X Li H G 2011 Bright spatial solitons in defocusing Kerr media with PT-symmetric potentials
Phys. Lett. A 84 053855
DOI:10.1103/PhysRevA.84.053855 [Cited within: 1]

[17]Midya B Roychoudhury R 2014 Nonlinear localized modes in PT-symmetric optical media with competing gain and loss
Ann. Phys. 341 12
DOI:10.1016/j.aop.2013.11.011 [Cited within: 1]

[18]Dai C Q Wang X G zhou G Q 2014 Stable light-bullet solutions in the harmonic and parity-time-symmetric potentials
Phys. Rev. A 89 013834
DOI:10.1103/PhysRevA.89.013834 [Cited within: 1]

[19]Xu Y J Dai C Q 2014 Nonlinear tunnelling of spatial solitons in PT-symmetric potential
Opt. Commun. 318 112
DOI:10.1016/j.optcom.2013.12.087 [Cited within: 1]

[20]Wang Y Y Dai C Q Wang X G 2014 Stable localized spatial solitons in PT-symmetric potentials with power-law nonlinearity
Nonlinear Dyn. 77 1323
DOI:10.1007/s11071-014-1381-6 [Cited within: 2]

[21]Dai C Q Chen R P Wang Y Y Fan Y 2017 Dynamics of light bullets in inhomogeneous cubic-quintic-septimal nonlinear media with PT-symmetric potentials
Nonlinear Dyn. 87 1675
DOI:10.1007/s11071-016-3143-0 [Cited within: 1]

[22]Dai C Q Wang Y Y Fan Y Yu D G 2018 Reconstruction of stability for Gaussian spatial solitons in quinticCseptimal nonlinear materials under PT-symmetric potentials
Nonlinear Dyn. 92 1351
DOI:10.1007/s11071-018-4130-4 [Cited within: 1]

[23]Zhou Z J Yan Z Y 2021 Solving forward and inverse problems of the logarithmic nonlinear Schrödinger equation with PT-symmetric harmonic potential via deep learning
Phys. Lett. A 387 127010
DOI:10.1016/j.physleta.2020.127010 [Cited within: 1]

[24]Akrivis G 1993 Finite difference discretization of the cubic Schrödinger equation
IMA J. Numer. Anal. 1 115
DOI:10.1093/imanum/13.1.115 [Cited within: 1]

[25]Xu Y Q Zhang L M 2012 Alternating direction implicit method for solving two-dimensional cubic nonlinear Schrödinger equation
Comput. Phys. Commun. 183 1082
DOI:10.1016/j.cpc.2012.01.006

[26]Zhao X Sun Z Z Hao Z P Fourth-order A 2014 Compact ADI scheme for Two-Dimensional Nonlinear Space Fractional Schrödinger Equation
SIAM J. Sci. Comput. 36 A2865
DOI:10.1137/140961560

[27]Qin J G 2010 The new alternating direction implicit difference methods for solving three-dimensional parabolic equations
Appl. Math. Model. 34 890
DOI:10.1016/j.apm.2009.07.006

[28]Tian Z F Yu P X 2010 High-order compact ADI (HOC-ADI) method for solving unsteady 2D Schrödinger equation
Comput. Phys. Commun. 181 861
DOI:10.1016/j.cpc.2009.12.022

[29]Zhang Q F Lin X M Pan K J Ren Y Z 2020 Linearized ADI schemes for two-dimensional space-fractional nonlinear Ginzburg-Landau equation
Comput. Math. Appl. 80 1201
DOI:10.1016/j.camwa.2020.05.027

[30]Inc M 2017 New type soliton solutions for the Zhiber-Shabat and related equations
Optik 138 1
DOI:10.1016/j.ijleo.2017.02.103

[31]Gao Z Xie S 2011 Fourth-order alternating direction implicit compact finite difference schemes for two-dimensional Schrödinger equations
Appl. Numer. Math. 61 593
DOI:10.1016/j.apnum.2010.12.004

[32]Inc M Yusuf A Aliyu A I Baleanu D 2018 Dark and singular optical solitons for the conformable space-time nonlinear Schrödinger equation with Kerr and power law nonlinearity
Optik 162 65
DOI:10.1016/j.ijleo.2018.02.085 [Cited within: 1]