Parallel Methods and Higher Dimensional NLS Equations

Abstract

Alternating direction implicit (ADI) schemes are proposed for the solution of the two-dimensional coupled nonlinear Schrödinger equation. These schemes are of second- and fourth-order accuracy in space and second order in time. The resulting schemes in each ADI computation step correspond to a block tridiagonal system which can be solved by using one-dimensional block tridiagonal algorithm with a considerable saving in computational time. These schemes are very well suited for parallel implementation on a high performance system with many processors due to the nature of the computation that involves solving the same block tridiagonal systems with many right hand sides. Numerical experiments on one processor system are conducted to demonstrate the efficiency and accuracy of these schemes by comparing them with the analytic solutions. The results show that the proposed schemes give highly accurate results.

1. Introduction

In this paper, we consider the coupled nonlinear Schrödinger equation 𝑖𝜓_𝑡+𝛿𝜓_𝑥𝑥+𝜓_𝑦𝑦+𝜓²+𝛼𝜙²(],𝜓=0,𝑥,𝑦,𝑡∈Ω×0,𝑇𝑖𝜙_𝑡+𝛿𝜙_𝑥𝑥+𝜙_𝑦𝑦+𝛼𝜓²+𝜙²𝜙=0 (1) with initial conditions 𝜓𝑥,𝑦,0=𝑓₁𝑥,𝑦,𝜙𝑥,𝑦,0=𝑓₂𝑥,𝑦,𝑥,𝑦∈Ω (2) and Dirichlet boundary conditions 𝜓𝑥,𝑦,𝑡=𝑔₁𝑥,𝑦,𝑡,𝜙𝑥,𝑦,𝑡=𝑔₂,(],𝑥,𝑦,𝑡𝑥,𝑦,𝑡∈𝜕Ω×0,𝑇 (3) where Ω is a rectangular domain in ℝ². We assume that Ω=[𝑎,𝑏]×[𝑐,𝑑],𝜕Ω is the boundary of Ω,(0,𝑇] is the time interval, and 𝑓₁(𝑥,𝑦), 𝑓₂(𝑥,𝑦), 𝑔₁(𝑥,𝑦,𝑡), and 𝑔₂(𝑥,𝑦,𝑡) are given sufficiently smooth functions. The two functions 𝜓(𝑥,𝑦,𝑡),𝜙(𝑥,𝑦,𝑡) are representing the amplitudes of the two circularly polarized waves. The values of 𝛼 vary over a wide range; that is, 𝛼≥2/3 and 𝛼≤7 correspond to kerr type electronic nonlinearity [1]. The physical significance of system (1) can be seen in the transverse effects in nonlinear optics. Since solitons interact like particles, the studies for the interactions of solitons have been of both experimental and theoretical interest. Many numerical methods have been developed for solving the coupled nonlinear Schrödinger [2–6]. Many published works for solving the two-dimensional nonlinear Schrödinger equation are given in [7–11]. In this work, we are going to derive an ADI method for solving the two-dimensional coupled nonlinear Schrödinger equation.

In this paper, we derive two ADI schemes for solving the two-dimensional coupled nonlinear Schrödinger system, one is of second order in space and time directions, and the other one is of fourth order in space and second order in time. Both methods are producing a block nonlinear tridiagonal system; a fixed point method has been developed to solve this system. The proposed schemes are unconditionally stable using Fourier stability analysis.

The ADI method [8–15], which replaces the solution of multidimensional problems by sequences of one-dimensional cases, only needs to solve tridiagonal linear system or block tridiagonal systems, and the resulting schemes are unconditionally stable and received much attention in recent years.

Following Biswas [16], we derive the soliton solution of the system (1) which can be written as 𝜓𝛽𝑥,𝑦,𝑡=𝐴sech₁𝑥+𝛽₂𝑦−𝑣𝑡exp𝑖−𝑘₁𝑥−𝑘₂𝑦+𝜔𝑡+𝜃,𝜙𝛽𝑥,𝑦,𝑡=𝐵sech₁𝑥+𝛽₂𝑦−𝑣𝑡exp𝑖−𝑘₁𝑥−𝑘₂𝑦+𝜔𝑡+𝜃, (4)(5) where 𝜔=𝛿𝛽²₁+𝛽²₂−𝑘²₁+𝑘²₂,𝑣=−2𝛿𝑘₁𝛽₁+𝑘₂𝛽₂𝛽𝐴=𝐵=2𝛿²₁+𝛽²₂(1+𝛼)^1/2, (6) and 𝛼, 𝛽₁, 𝛽₂, 𝑘₁, and 𝑘₂ are arbitrary constants (see the appendix). The system has the conserved quantities [15] 𝐼₁=^∞_−∞𝜓²𝑑𝑥𝑑𝑦=constant,𝐼₂=^∞_−∞𝜙²𝑑𝑥𝑑𝑦=constant. (7)(8)

To avoid the complex computation, we assume 𝜓=𝑢₁+𝑖𝑢₂,𝜙=𝑢₃+𝑖𝑢₄. (9)

The system (1) can be written as 𝜕𝑢₁𝜕𝜕𝑡+𝛿²𝑢₂𝜕𝑥²+𝜕²𝑢₂𝜕𝑦²+𝐹(𝐮)𝑢₂=0,𝜕𝑢₂𝜕𝜕𝑡−𝛿²𝑢₁𝜕𝑥²+𝜕²𝑢₁𝜕𝑦²−𝐹(𝐮)𝑢₁=0,𝜕𝑢₃𝜕𝜕𝑡+𝛿²𝑢₄𝜕𝑥²+𝜕²𝑢₄𝜕𝑦²+𝐺(𝐮)𝑢₄=0,𝜕𝑢₄𝜕𝜕𝑡−𝛿²𝑢₃𝜕𝑥²+𝜕²𝑢₃𝜕𝑦²−𝐺(𝐮)𝑢₃=0, (10) where 𝐹(𝐮)=𝑢²₁+𝑢²₂+𝛼𝑢²₃+𝑢²₄,𝐺(𝐮)=𝛼𝑢²₁+𝑢²₂+𝑢²₃+𝑢²₄. (11)

System (10) can be written in matrix vector form as 𝜕𝐮𝜕𝜕𝑡+𝛿𝐴²𝐮𝜕𝑥²+𝜕²𝐮𝜕𝑦²+𝐵(𝐮)𝐮=𝟎, (12) where 𝐴=0100−1000000100−10,𝐵=0𝑓₁00−𝑓₁000000𝑓₂00−𝑓₂0,𝑓₁=𝑢²₁+𝑢²₂+𝛼𝑢²₃+𝑢²₄,𝑓₂=𝛼𝑢²₁+𝑢²₂+𝑢²₃+𝑢²₄. (13)

The paper is organized as follows. In Section 2 we give two ADI schemes for solving the two-dimensional coupled nonlinear Schrödinger equation, and in Section 3, we present the von Neumann stability analysis for the proposed schemes. Numerical experiments for several problems are presented in Section 4. A parallel algorithm for the proposed ADI schemes is given in Section 5. Conclusions are given in Section 6.

2. Numerical Method

To derive the numerical schemes for solving system (1), we consider the domain of interest Ω=[𝑎,𝑏]×[𝑎,𝑏], such that (𝑥,𝑦)∈Ω. The domain is divided by a uniform mesh in each direction such that 𝑥_𝑙𝑦=𝑎+𝑙ℎ,𝑙=0,1,2,…,𝐿,_𝑚𝑡=𝑎+𝑚ℎ,𝑚=0,1,2,…,𝑀,_𝑛=𝑛𝑘,𝑛=0,1,2,…,𝑁, (14) where ℎ and 𝑘 are the space and time step sizes, respectively. We denote 𝑢^𝑛_𝑙,𝑚 and 𝑈^𝑛_𝑙,𝑚 to be the exact and the numerical solutions at the point (𝑥_𝑙,𝑦_𝑚,𝑡_𝑛), respectively.

2.1. Second-Order ADI Method

To derive the first scheme, we approximate the space derivative using the central difference formulae 𝜕²𝑢𝜕𝑥²=1ℎ²𝛿²_𝑥𝑢_𝑙,𝑚+𝑂ℎ²=1ℎ²𝑢_𝑙+1,𝑚−2𝑢_𝑙,𝑚+𝑢_{𝑙−1,𝑚}+𝑂ℎ²,𝜕²𝑢𝜕𝑦²=1ℎ²𝛿²_𝑦𝑢_𝑙,𝑚+𝑂ℎ²=1ℎ²𝑢_𝑙,𝑚+1−2𝑢_𝑙,𝑚+𝑢_{𝑙,𝑚−1}+𝑂ℎ². (15)

We apply (15) to system (12); this will lead us to the following first order differential system in time 𝐮_𝑡+1ℎ²𝐴𝛿²_𝑥+𝛿²_𝑦𝐮_𝑙,𝑚+𝐵𝐮_𝑙,𝑚𝐮_𝑙,𝑚=𝟎. (16)

Equation (16) is of second order in space.

By applying the Crank-Nicolson method for the temporal discretization, we get the following difference scheme with accuracy 𝑂(𝑘²+ℎ²): 1𝑘𝐔^𝑛+1_𝑙,𝑚−𝐔^𝑛_𝑙,𝑚+12ℎ²𝐴𝛿²_𝑥+𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚+𝐔^𝑛_𝑙,𝑚+𝐵𝐔𝐔=𝟎, (17) where 1𝐔=2𝐔^𝑛+1_𝑙,𝑚+𝐔^𝑛_𝑙,𝑚. (18)

This scheme can be written as 1𝐼+2𝑟𝐴𝛿²_𝑥+12𝑟𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=1𝐼−2𝑟𝐴𝛿²_𝑥−12𝑟𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝐵𝐔𝐔 (19) which can be approximated by using the factored form 1𝐼+2𝑟𝐴𝛿²_𝑥1𝐼+2𝑟𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=1𝐼−2𝑟𝐴𝛿²_𝑥1𝐼−2𝑟𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘𝐵𝐔𝑟𝐔+²4𝐴𝛿²_𝑥𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚−𝐔^𝑛_𝑙,𝑚. (20)

By Taylor’s expansion, the last term in (20) can be written as 𝑟²4𝐴𝛿²_𝑥𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚−𝐔^𝑛_𝑙,𝑚𝑘=−²4𝛿²_𝑥𝛿²_𝑥ℎ⁴𝑘𝜕𝑢𝑘𝜕𝑡+𝑂³𝑘=−³4𝜕⁵𝑢𝜕²𝑥𝜕²𝑘𝑦𝜕𝑡+𝑂³ℎ² (21) which is of the same order of accuracy of the truncation error, and we ignore it to obtain 1𝐼+2𝑟𝐴𝛿²_𝑥1𝐼+2𝑟𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=1𝐼−2𝑟𝐴𝛿²_𝑥1𝐼−2𝑟𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘𝐵𝐔𝐔. (22)

By introducing a new intermediate vector 𝐔^∗_𝑙,𝑚, we propose a D’Yakonov [12, 17] ADI-like scheme for the coupled system 1𝐼+2𝑟𝐴𝛿²_𝑥𝐔^∗_𝑙,𝑚=1𝐼−2𝑟𝐴𝛿²_𝑥1𝐼−2𝑟𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘𝐵𝐔𝐔1𝐼+2𝑟𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=𝐔^∗_𝑙,𝑚, (23)(24) which is a nonlinear scheme. An iterative algorithm of fixed point nature can be used to solve the system of the nonlinear equations (23)-(24). The fixed point that we propose can be given by 1𝐼+2𝑟𝐴𝛿²_𝑥𝐔^∗_𝑙,𝑚=1𝐼−2𝑟𝐴𝛿²_𝑥1𝐼−2𝑟𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚𝐔−𝑘𝐵^{𝑛+1,𝑠𝑙,𝑚}+𝐔^𝑛_𝑙,𝑚2×𝐔^𝑛+1,𝑠_𝑙,𝑚+𝐔^𝑛_𝑙,𝑚2,1𝐼+2𝑟𝐴𝛿²_𝑦𝐔^{𝑛+1,𝑠+1}_𝑙,𝑚=𝐔^∗_𝑙,𝑚, (25) where the superscript 𝑠 denotes the 𝑠th iterate for solving the nonlinear system of equations for each time step. The block tridiagonal matrix equations of (25) can be solved by Crout’s method. The initial iterate 𝐔^𝑛+1,0_𝑙,𝑚 is chosen as 𝐔^𝑛+1,0_𝑙,𝑚=𝐔^𝑛_𝑙,𝑚. (26)

The iteration continues until the condition 𝐔^𝑛,𝑠+1_𝑙,𝑚−𝐔^𝑛,𝑠_𝑙,𝑚∞≤10⁻⁶ (27) is satisfied.

The boundary value of the intermediate variable 𝐔̇_𝑙,𝑚 can be extracted from (24) and is given by the following formulas: 𝐔^∗_0,𝑚=𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_0,𝑚,𝐔^∗_𝐿+1,𝑚=𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_𝐿,𝑚. (28)

2.2. Fourth-Order ADI Method

Now, we want to derive a highly accurate fourth order ADI scheme; to do this, we approximate the space derivatives by the following formulas 𝜕²𝑢𝜕𝑥²=1ℎ²𝛿²_𝑥1+(1/12)𝛿²_𝑥𝐮_𝑙,𝑚+𝑂ℎ⁴,𝜕²𝑢𝜕𝑦²=1ℎ²𝛿²_𝑦1+(1/12)𝛿²_𝑦𝐮_𝑙,𝑚+𝑂ℎ⁴. (29)

By using these approximations together with Crank-Nicolson for the time direction, we get the numerical scheme 1𝑘𝐔^𝑛+1_𝑙,𝑚−𝐔^𝑛_𝑙,𝑚+12ℎ²𝐴𝛿²_𝑥1+(1/12)𝛿²_𝑥+𝛿²_𝑦1+(1/12)𝛿²_𝑦×𝐔^𝑛+1_𝑙,𝑚+𝐔^𝑛_𝑙,𝑚+𝐵𝐔𝐔=𝟎 (30) which can be written as 𝑟𝐼+2𝐴𝛿²_𝑥1+(1/12)𝛿²_𝑥+𝛿²_𝑦1+(1/12)𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=𝑟𝐼−2𝐴𝛿²_𝑥1+(1/12)𝛿²_𝑥+𝛿²_𝑦1+(1/12)𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘𝐵𝐔𝐔, (31) and this can be written in the factored form as 𝑟𝐼+2𝐴𝛿²_𝑥1+(1/12)𝛿²_𝑥×𝑟𝐼+2𝐴𝛿²_𝑦1+(1/12)𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=𝑟𝐼−2𝐴𝛿²_𝑥1+(1/12)𝛿²_𝑥×𝑟𝐼−2𝐴𝛿²_𝑦1+(1/12)𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘𝐵𝐔𝐔. (32)

The difference between (31) and (32) is −𝑟²4𝛿²_𝑥1+(1/12)𝛿²_𝑥𝛿²_𝑦1+(1/12)𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚−𝐔^𝑛_𝑙,𝑚𝑘=−²4ℎ⁴𝛿²_𝑥1+(1/12)𝛿²_𝑥𝛿²_𝑦1+(1/12)𝛿²_𝑦×𝑘𝜕𝑢𝑘𝜕𝑡+𝑂³𝑘=−³4𝜕⁵𝑢𝜕²𝑥𝜕²𝑘𝑦𝜕𝑡+𝑂³ℎ⁴ (33) which is of the same order of accuracy of the truncation error, we ignore it, and then we operate on both sides of (32) by (1+(1/12)𝛿²_𝑥)(1+(1/12)𝛿²_𝑦)11+𝛿12^2𝑥𝑟𝐼+2𝐴𝛿²_𝑥11+𝛿12^2𝑦𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=11+𝛿12^2𝑥𝑟𝐼−2𝐴𝛿²_𝑥×11+𝛿12^2𝑦𝑟𝐼−2𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚1−𝑘1+𝛿12^2𝑥11+𝛿12^2𝑦𝐵𝐔𝐔. (34)

The fourth-order D’Yakonov ADI-like scheme in this case can be displayed as 11+𝛿12^2𝑥𝑟𝐼+2𝐴𝛿²_𝑥𝐔̇=11+𝛿12^2𝑥𝑟𝐼−2𝐴𝛿²_𝑥11+𝛿12^2𝑦𝑟𝐼−2𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚1−𝑘1+𝛿12^2𝑥11+𝛿12^2𝑦𝐵𝐔𝐔11+𝛿12^2𝑦𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=𝐔̇. (35)(36)

Now the systems in (35) and (36) are nonlinear. By the similar approach used in the previous scheme, we can derive a fixed point iterative formulas.

The boundary value of the intermediate variable 𝐔̇_𝑙,𝑚 in (35) can be given by the following formulas: 𝐔̇_0,𝑚=1+1𝛿12^2𝑦𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_0,𝑚,𝐔̇_𝐿+1,𝑚=1+1𝛿12^2𝑦𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_𝐿+1,𝑚. (37)

We can easily write a generalized version of D’Yakonov ADI like method for solving the two-dimensional coupled nonlinear Schrödinger system (1) as 1+𝜎𝛿²_𝑥𝑟𝐼+2𝐴𝛿²_𝑥𝐔̇=1+𝜎𝛿²_𝑥𝑟𝐼−2𝐴𝛿²_𝑥1+𝜎𝛿²_𝑦𝑟𝐼−2𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘1+𝜎𝛿²_𝑥1+𝜎𝛿²_𝑦𝐵𝐔𝐔,1+𝜎𝛿²_𝑦𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=𝐔̇ (38) for arbitrary value of 𝜎, and for 𝜎=0, 𝜎=1/12 we recover the second- and fourth-order schemes, respectively.

3. Stability Analysis

To study the stability of the proposed scheme, we consider the von Neumann stability analysis which can be only applied for the linear finite difference scheme, so we consider the linear version of the generalized scheme 1+𝜎𝛿²_𝑥𝑟𝐼+2𝐴𝛿²_𝑥1+𝜎𝛿²_𝑦𝑟𝐼+2𝐴𝛿²_𝑦𝐔^𝑛+1_𝑙,𝑚=1+𝜎𝛿²_𝑥𝑟𝐼−2𝐴𝛿²_𝑥1+𝜎𝛿²_𝑦𝑟𝐼−2𝐴𝛿²_𝑦𝐔^𝑛_𝑙,𝑚−𝑘2𝜔1+𝜎𝛿²_𝑥1+𝜎𝛿²_𝑦𝐔, (39) where 𝜔 is constant and 𝜎=0,1/12.

Suppose that the numerical solution can be expressed by Fourier series, whose typical term is 𝐔^𝑛_𝑙,𝑚=𝐺^𝑛𝐖^𝑛exp𝑖𝛽₁𝑙ℎ+𝑖𝛽₂𝑚ℎ, (40) where 𝑖=√−1,𝐺^𝑛 is the amplification matrix at time level 𝑛, and 𝛽₁,𝛽₂ are the wave numbers in 𝑥 and 𝑦 directions. Substituting (40) into (39) will lead to the matrix equation 𝐼−2𝑟𝜇_𝑥𝐴𝐼−2𝑟𝜇_𝑦𝐴𝐺=𝐼+2𝑟𝜇_𝑥𝐴𝐼+2𝑟𝜇_𝑦𝐴−𝑘2𝜔[𝐺+𝐼], (41) where 𝜇_𝑥=sin²𝛽₁𝑙ℎ,𝜇_𝑦=sin²𝛽₂𝑚ℎ. (42)

The eigenvalues of the matrix 𝐺 are given by eigenvaluesof𝐺=1+2𝑖𝑟𝜇_𝑥1+2𝑖𝑟𝜇_𝑦−(𝑘/2)𝜔1−2𝑖𝑟𝜇_𝑥1−2𝑖𝑟𝜇_𝑦+()𝜔.𝑘/2 (43)

It can be easily shown that the modulus of the maximum eigenvalue of the amplification matrix 𝐺 is less than one; hence, the scheme is unconditionally stable according to the von Neumann stability analysis.

4. Numerical Results

In this section, the efficiency and accuracy of the proposed schemes will be tested by comparing with the exact solutions. We will measure the accuracy of the proposed schemes using the 𝐿_∞ norm. We compute the conserved quantity by using the trapezoidal rule.

4.1. Example 1

In this example, we choose the initial conditions from the exact solutions 𝜓(𝜉),𝜙(𝜉),𝑥,𝑦,𝑡=𝐴sechexp𝑖𝜂𝑥,𝑦,𝑡=𝐵sechexp𝑖𝜂 (44) where 𝜉=𝛽₁𝑥+𝛽₂𝑦−𝑣𝑡,𝜂=−𝜅₁𝑥−𝜅₂𝑦+𝜔𝑡+𝜃,𝐴=𝐵=2𝛿(𝛽²₁+𝛽²₂)(1+𝛼)^1/2, (45) at 𝑡=0. The following parameters are used ℎ=0.2,𝑘=0.001,𝛽₁𝛽=0.5,₂=1.0,𝑘₁=0.5,𝑘₂1=1.0,𝛼=1,𝛿=2,−10≤𝑥,𝑦≤10. (46)

The boundary conditions are extracted from the exact solution (44).

Tables 1, 2, and 3 show the errors (ER) and the conserved quantities (I) for 𝜎=0,1/6 and 𝜎=1/12, respectively. We have noticed that the scheme with 𝜎=1/12 produced highly accurate results, and this is due to the fourth order accuracy in space and second-order accuracy in time. The other two methods using 𝜎=0 and 𝜎=1/6 are of second-order accuracy in space and time. Figure 1 shows a single soliton at 𝑡=0. Figure 2 shows the numerical solution at 𝑡=1.

Table 1 

Single soliton with 𝜎=0.

Table 2 

Single soliton with 𝜎=1/6.

Table 3 

Single soliton with 𝜎=1/12.

Figure 1 

Single soliton 𝛽1=𝜅1=0.5, 𝛽2=𝜅2=1, 𝛼=1, 𝑡=0.

Figure 2 

Single soliton 𝛽1=𝜅1=0.5, 𝛽2=𝜅2=1, 𝛼=1, 𝑡=1.

4.2. Example 2

In this example we will choose the initial condition [1] 𝜓=𝛼𝑥,𝑦,0¹2exp−𝜂2×𝜅sechRe₁𝑥+𝑙₁𝑦+𝜂2exp𝑖Im𝜅₁𝑥+𝑙₁𝑦𝜙=𝛽𝑥,𝑦,0¹2exp−𝜂2×𝜅sechRe₁𝑥+𝑙₁𝑦+𝜂2exp𝑖Im𝜅₁𝑥+𝑙₂𝑦, (47) where 𝜂=1exp2𝛼₁²+𝛽₁²𝜅₁+𝜅^∗₁²+𝑙₁+𝑙^∗₁², (48) which represent two solitons of different amplitudes. The parameters 𝛼₁, 𝛽₁, 𝜅₁, and 𝑙₁ are complex parameters. In this test, we choose the parameters ℎ=0.2,𝑘=0.01,𝑥_𝑙=−10,𝛼₁𝛽=1,₁=0.5,𝜅₁=𝑙₁=1+𝑖. (49)

In Figures 3 and 4, we display the numerical solution of |𝜓|² at 𝑡=0 and 𝑡=1, while in Figures 5 and 6 we display the numerical solution of |𝜙|² at 𝑡=0 and 𝑡=1.

Figure 3 

Soliton solution |𝜓|2 at 𝑡=0.

Figure 4 

Soliton solution |𝜓|2 at 𝑡=1.

Figure 5 

Soliton solution |𝜙|2 at 𝑡=0.

Figure 6 

Soliton solution |𝜙|2 at 𝑡=1.

4.3. Example 3

To study the interaction of two solitons, many numerical tests have been conducted with different initial conditions, and, among these tests, we select the initial conditions of the form 𝜓−()𝑥,𝑦,0=𝐴exp𝑥−𝑑²−𝑦²exp−𝑖lncosh𝑥²+𝑦²,𝜙−()𝑥,𝑦,0=𝐵exp𝑥+𝑑²−𝑦²exp−𝑖lncosh𝑥²+𝑦², (50) which represent two solitons moving in the opposite direction and centered at 𝑑 and −𝑑, respectively. In this test we choose the parameters ℎ=0.2,𝑘=0.01,𝑥_𝑙=−10,𝑑=5.0,𝐴=1,𝐵=0.8,𝑡=0,…,1.5. (51)

The interaction scenario is given in Figures 7 and 8. In Figure 7 we show the two solitons with two different amplitudes at 𝑡=0. In Figure 8 we display the interaction where the two solitons interact at 𝑡=0.5, and, in Figure 9, we display the two solitons after the interaction at 𝑡=1.5. In Table 4, we display the conserved quantities, and we see that the numerical method we proposed conserves the conserved quantities almost exactly. It is easy to see that the interaction regiem is an inelastic one. See [18].

Table 4 

The conserved quantities during the interaction scenario.

Figure 7 

Interaction of two solitons |𝜓|2 and |𝜙|2 at 𝑡=0.

Figure 8 

Interaction of two solitons |𝜓|2 and |𝜙|2 at 𝑡=0.5.

Figure 9 

Interaction of two solitons |𝜓|2 and |𝜙|2 at 𝑡=1.5.

4.4. Example 4

The system under consideration generates a progressive plane wave solutions [7]: 𝜓𝑥,𝑦,𝑡=𝐴exp𝑖𝜅₁𝑥+𝜅₂,𝜙𝑖𝜅𝑦+𝜔𝑡𝑥,𝑦,𝑡=𝐴exp₁𝑥+𝜅₂𝑦+𝜔𝑡, (52) where 𝜔=𝑘²₁+𝑘²₂−(1+𝛼)|𝐴|². (53)

Our numerical experiments are conducted in the domain [0,2𝜋]×[0,2𝜋] with 𝐴=1,𝛼=1,𝑘₁=𝑘₂𝜋=1,𝛿=1,𝑘=0.001,ℎ=.50 (54)

Initial and boundary conditions are extracted from the exact solution. In this example we choose 𝜎=1/12, the fourth order ADI method. The numerical results are presented in Table 5. Figure 10 displays the solution at 𝑡=0, and Figure 11 displays, the solution at 𝑡=1.

Table 5 

Periodic solution 𝑘₁=𝑘₂=𝛼=1.

Figure 10 

Initial condition with 𝑘1=𝑘2=𝛼=1, at 𝑡=0.

Figure 11 

Numerical solution with 𝑘1=𝑘2=𝛼=1, at 𝑡=1.

Table 4 displays the accuracy of the scheme and preserves the conserved quantities.

5. Parallel Algorithm for the Proposed ADI Schemes

It is to be noted that the implementation of the ADI schemes requires solving the same block-tridiagonal matrix with different right-hand sides. This can be done efficiently by using the fast parallel algorithm given in [19]. This algorithm [19] is a generalization of the parallel dichotomy algorithm for solving tridiagonal liner system of equations [20].

It has been shown that this dichotomy yields almost a linear speedup on a high performance system with many processors [19]. We expect the same speedup for our proposed methods. The parallel implementation of the proposed ADI method will be reported in our future work.

6. Concluding Remarks

In this present work, a generalized alternating direction implicit methods for solving two-dimensional coupled nonlinear Schrödinger equation have been established, the methods are unconditionally stable. The methods produced schemes of second and fourth order accuracy in space and second order in time according to the selected value of 𝜎. The numerical results have shown that the proposed schemes successfully combine accuracy and efficiency for the two-dimensional coupled nonlinear Schrödinger system. The implementation of the ADI schemes requires solving the same block-tridiagonal matrix with different right-hand sides. This can be done by using the fast parallel algorithm given in [19]. This algorithm yields almost a linear speedup on a high performance system with many processors [19]. We expect the same speedup for our proposed methods. The proposed schemes can be easily extended to solve higher dimensional coupled nonlinear Schrödinger system.

Appendix

The exact solution of the nonlinear system can be derived as follows: we assume [19] 𝜓(𝜉),𝜙(𝜉),𝑥,𝑦,𝑡=𝐴sechexp𝑖𝜂𝑥,𝑦,𝑡=𝐵sechexp𝑖𝜂 (A.1) where 𝜉=𝛽₁𝑥+𝛽₂𝑦−𝑣𝑡,𝜂=−𝜅₁𝑥−𝜅₂𝑦+𝜔𝑡+𝜃. (A.2)

Now from (A.1), we deduce 𝜓_𝑡=(𝜉)(𝜉)(𝜉,𝜓𝑣𝐴tanhsech+𝑖𝜔𝐴sech)exp𝑖𝜂_𝑥𝑥=𝐴𝛽²₁sech(𝜉)−2𝐴𝛽²₁sech³(𝜉)−𝐴𝜅²₁sech(𝜉)+2𝑖𝜅₁𝛽₁(𝜉)(𝜉,𝜓𝐴tanhsech)exp𝑖𝜂_𝑦𝑦=𝐴𝛽²₂sech(𝜉)−2𝐴𝛽²₂sech³(𝜉)−𝐴𝜅²₂sech(𝜉)+2𝑖𝜅₂𝛽₂𝐴tanh(𝜉)sech(𝜉)exp𝑖𝜂, (A.3) and similar expressions can be obtained for the function 𝜙. By substituting these expressions into the given system, and equating the real and imaginary parts, this will produce the following relations: 𝑣=−2𝛿𝜅₁𝛽₁+𝜅₂𝛽₂,−𝜔𝐴sech(𝜉)+𝛿𝐴𝛽²₁sech(𝜉)−2𝛿𝐴𝛽²₁sech³(𝜉)−𝛿𝐴𝜅²₁sech(𝜉)+𝛿𝐴𝛽²₂sech(𝜉)−2𝛿𝐴𝛽²₂sech³(𝜉)−𝛿𝐴𝜅²₂(𝜉)+𝐴sech₂+𝛼𝐵²𝐴sech³(𝜉)=0. (A.4)(A.5) By equating the coefficient of sech³(𝜉) in (A.5), we obtain 𝐴²+𝛼𝐵²=2𝛿𝛽²₁+𝛽²₂, (A.6) and by equating the coefficients of sech(𝜉), we will get the relation 𝜔=𝛿𝛽²₁+𝛽²₂−𝜅²₁+𝜅²₂. (A.7) From the second equation of the proposed system (𝜙), we will get the relation 𝛼𝐴²+𝐵²=2𝛿𝛽²₁+𝛽²₂, (A.8) and from (A.6) and (A.8), we deduce that 𝐴=𝐵, and then we get the relation 𝐴=𝐵=2𝛿(𝛽²₁+𝛽²₂)(1+𝛼)^1/2 (A.9) with 𝛼, 𝛽, 𝜅, and 𝑙 as arbitrary complex parameters.

Acknowledgments

This work was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under Grant no. (130-074-D1433). The authors, therefore, acknowledge with thanks the DSR technical and financial support.