# Parallel Methods and Higher Dimensional NLS Equations

Authors: M. S. Ismail; T. R. Taha
Journal: Abstract and Applied Analysis
 (2013)

---

## 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.


---

## Body

## 1. Introduction

In this paper, we consider the coupled nonlinear Schrödinger equation(1)iψt+δ(ψxx+ψyy)+(|ψ|2+α|ϕ|2)ψ=0,iψt+δ(ψxx+ψyy)+(x,y,t)∈Ω×(0,T],iϕt+δ(ϕxx+ϕyy)+(α|ψ|2+|ϕ|2)ϕ=0
with initial conditions
(2)ψ(x,y,0)=f1(x,y),ϕ(x,y,0)=f2(x,y),ψ(x,y,0)=f1(x,y),1ϕ(x,y,0)(x,y)∈Ω
and Dirichlet boundary conditions
(3)ψ(x,y,t)=g1(x,y,t),ϕ(x,y,t)=g2(x,y,t),ψ(x,y,t)=g1(x,y,t),(x,y,t)∈∂Ω×(0,T],
where Ω is a rectangular domain in ℝ2. We assume that Ω=[a,b]×[c,d],∂Ω is the boundary of Ω,(0,T] is the time interval, and f1(x,y), f2(x,y), g1(x,y,t), and g2(x,y,t) are given sufficiently smooth functions. The two functions ψ(x,y,t),ϕ(x,y,t) 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
(4)ψ(x,y,t)=Asech(β1x+β2y-vt)exp[i(-k1x-k2y+ωt+θ)],(5)ϕ(x,y,t)=Bsech(β1x+β2y-vt)exp[i(-k1x-k2y+ωt+θ)],
where
(6)ω=[δ(β12+β22)-(k12+k22)],v=-2δ(k1β1+k2β2)A=B=[2δ(β12+β22)(1+α)]1/2,
and α, β1, β2, k1, and k2 are arbitrary constants (see the appendix). The system has the conserved quantities [15]
(7)I1=∬-∞∞|ψ|2dxdy=constant,(8)I2=∬-∞∞|ϕ|2dxdy=constant.To avoid the complex computation, we assume(9)ψ=u1+iu2,ϕ=u3+iu4.The system (1) can be written as
(10)∂u1∂t+δ(∂2u2∂x2+∂2u2∂y2)+F(u)u2=0,∂u2∂t-δ(∂2u1∂x2+∂2u1∂y2)-F(u)u1=0,∂u3∂t+δ(∂2u4∂x2+∂2u4∂y2)+G(u)u4=0,∂u4∂t-δ(∂2u3∂x2+∂2u3∂y2)-G(u)u3=0,
where
(11)F(u)=u12+u22+α(u32+u42),G(u)=α(u12+u22)+u32+u42.System (10) can be written in matrix vector form as
(12)∂u∂t+δA(∂2u∂x2+∂2u∂y2)+B(u)u=0,
where
(13)A=[0100-1000000100-10],B=[0f100-f1000000f200-f20],f1=(u12+u22)+α(u32+u42),f2=α(u12+u22)+(u32+u42).The paper is organized as follows. In Section2 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 Ω=[a,b]×[a,b], such that (x,y)∈Ω. The domain is divided by a uniform mesh in each direction such that
(14)xl=a+lh,l=0,1,2,…,L,ym=a+mh,m=0,1,2,…,M,tn=nk,n=0,1,2,…,N,
where h and k are the space and time step sizes, respectively. We denote ul,mn and Ul,mn to be the exact and the numerical solutions at the point (xl,ym,tn), respectively.


### 2.1. Second-Order ADI Method

To derive the first scheme, we approximate the space derivative using the central difference formulae(15)∂2u∂x2=1h2δx2ul,m+O(h2)=1h2(ul+1,m-2ul,m+ul-1,m)+O(h2),∂2u∂y2=1h2δy2ul,m+O(h2)=1h2(ul,m+1-2ul,m+ul,m-1)+O(h2).We apply (15) to system (12); this will lead us to the following first order differential system in time
(16)ut+1h2A[δx2+δy2]ul,m+B(ul,m)ul,m=0.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 accuracyO(k2+h2):
(17)1k[Ul,mn+1-Ul,mn]+12h2A[δx2+δy2][Ul,mn+1+Ul,mn]+B(U-)U-=0,
where
(18)U-=12[Ul,mn+1+Ul,mn].This scheme can be written as(19)[I+12rAδx2+12rAδy2]Ul,mn+1=[I-12rAδx2-12rAδy2]Ul,mn-B(U-)U-
which can be approximated by using the factored form
(20)[I+12rAδx2][I+12rAδy2]Ul,mn+1=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(U-)U-+r24{Aδx2}{Aδy2}(Ul,mn+1-Ul,mn).By Taylor’s expansion, the last term in (20) can be written as
(21)r24{Aδx2}{Aδy2}(Ul,mn+1-Ul,mn)=-k24δx2δx2h4(k(∂u∂t)+O(k3))=-k34[∂5u∂2x∂2y∂t]+O(k3h2)
which is of the same order of accuracy of the truncation error, and we ignore it to obtain
(22)[I+12rAδx2][I+12rAδy2]Ul,mn+1=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(U-)U-.By introducing a new intermediate vectorUl,m*, we propose a D’Yakonov [12, 17] ADI-like scheme for the coupled system
(23)[I+12rAδx2]Ul,m*=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(U-)U-(24)[I+12rAδy2]Ul,mn+1=Ul,m*,
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
(25)[I+12rAδx2]Ul,m*=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(Ul,mn+1,s+Ul,mn2)×[Ul,mn+1,s+Ul,mn]2,[I+12rAδy2]Ul,mn+1,s+1=Ul,m*,
where the superscript s denotes the sth 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 Ul,mn+1,0 is chosen as
(26)Ul,mn+1,0=Ul,mn.The iteration continues until the condition(27)∥Ul,mn,s+1-Ul,mn,s∥∞≤10-6
is satisfied.The boundary value of the intermediate variableU˙l,m can be extracted from (24) and is given by the following formulas:
(28)U0,m*=[I+r2Aδy2]U0,mn+1,UL+1,m*=[I+r2Aδy2]UL,mn+1.


### 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(29)∂2u∂x2=1h2δx21+(1/12)δx2ul,m+O(h4),∂2u∂y2=1h2δy21+(1/12)δy2ul,m+O(h4).By using these approximations together with Crank-Nicolson for the time direction, we get the numerical scheme(30)1k[Ul,mn+1-Ul,mn]+12h2A[δx21+(1/12)δx2+δy21+(1/12)δy2]×[Ul,mn+1+Ul,mn]+B(U-)U-=0
which can be written as
(31)[I+r2A(δx21+(1/12)δx2+δy21+(1/12)δy2)]Ul,mn+1=[I-r2A(δx21+(1/12)δx2+δy21+(1/12)δy2)]Ul,mn-kB(U-)U-,
and this can be written in the factored form as
(32)[I+r2A(δx21+(1/12)δx2)]=×[I+r2A(δy21+(1/12)δy2)]Ul,mn+1=[I-r2A(δx21+(1/12)δx2)]=×[I-r2A(δy21+(1/12)δy2)]Ul,mn-kB(U-)U-.The difference between (31) and (32) is
(33)-r24[δx21+(1/12)δx2][δy21+(1/12)δy2](Ul,mn+1-Ul,mn)=-k24h4[δx21+(1/12)δx2][δy21+(1/12)δy2]=×(k(∂u∂t)+O(k3))=-k34[∂5u∂2x∂2y∂t]+O(k3h4)
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)δx2)(1+(1/12)δy2)(34)[(1+112δx2)I+r2Aδx2][(1+112δy2)I+r2Aδy2]Ul,mn+1=[(1+112δx2)I-r2Aδx2]=×[(1+112δy2)I-r2Aδy2]Ul,mn-k(1+112δx2)(1+112δy2)B(U-)U-.The fourth-order D’Yakonov ADI-like scheme in this case can be displayed as(35)[(1+112δx2)I+r2Aδx2]U˙=[(1+112δx2)I-r2Aδx2][(1+112δy2)I-r2Aδy2]Ul,mn-k(1+112δx2)(1+112δy2)B(U-)U-(36)[(1+112δy2)I+r2Aδy2]Ul,mn+1=U˙.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 variableU˙l,m in (35) can be given by the following formulas:
(37)U˙0,m=[(1+112δy2)I+r2Aδy2]U0,mn+1,U˙L+1,m=[(1+112δy2)I+r2Aδy2]UL+1,mn+1.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
(38)[(1+σδx2)I+r2Aδx2]U˙=[(1+σδx2)I-r2Aδx2][(1+σδy2)I-r2Aδy2]Ul,mn-k(1+σδx2)(1+σδy2)B(U-)U-,[(1+σδy2)I+r2Aδy2]Ul,mn+1=U˙
for arbitrary value of σ, and for σ=0, σ=1/12 we recover the second- and fourth-order schemes, respectively.


## 2.1. Second-Order ADI Method

To derive the first scheme, we approximate the space derivative using the central difference formulae(15)∂2u∂x2=1h2δx2ul,m+O(h2)=1h2(ul+1,m-2ul,m+ul-1,m)+O(h2),∂2u∂y2=1h2δy2ul,m+O(h2)=1h2(ul,m+1-2ul,m+ul,m-1)+O(h2).We apply (15) to system (12); this will lead us to the following first order differential system in time
(16)ut+1h2A[δx2+δy2]ul,m+B(ul,m)ul,m=0.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 accuracyO(k2+h2):
(17)1k[Ul,mn+1-Ul,mn]+12h2A[δx2+δy2][Ul,mn+1+Ul,mn]+B(U-)U-=0,
where
(18)U-=12[Ul,mn+1+Ul,mn].This scheme can be written as(19)[I+12rAδx2+12rAδy2]Ul,mn+1=[I-12rAδx2-12rAδy2]Ul,mn-B(U-)U-
which can be approximated by using the factored form
(20)[I+12rAδx2][I+12rAδy2]Ul,mn+1=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(U-)U-+r24{Aδx2}{Aδy2}(Ul,mn+1-Ul,mn).By Taylor’s expansion, the last term in (20) can be written as
(21)r24{Aδx2}{Aδy2}(Ul,mn+1-Ul,mn)=-k24δx2δx2h4(k(∂u∂t)+O(k3))=-k34[∂5u∂2x∂2y∂t]+O(k3h2)
which is of the same order of accuracy of the truncation error, and we ignore it to obtain
(22)[I+12rAδx2][I+12rAδy2]Ul,mn+1=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(U-)U-.By introducing a new intermediate vectorUl,m*, we propose a D’Yakonov [12, 17] ADI-like scheme for the coupled system
(23)[I+12rAδx2]Ul,m*=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(U-)U-(24)[I+12rAδy2]Ul,mn+1=Ul,m*,
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
(25)[I+12rAδx2]Ul,m*=[I-12rAδx2][I-12rAδy2]Ul,mn-kB(Ul,mn+1,s+Ul,mn2)×[Ul,mn+1,s+Ul,mn]2,[I+12rAδy2]Ul,mn+1,s+1=Ul,m*,
where the superscript s denotes the sth 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 Ul,mn+1,0 is chosen as
(26)Ul,mn+1,0=Ul,mn.The iteration continues until the condition(27)∥Ul,mn,s+1-Ul,mn,s∥∞≤10-6
is satisfied.The boundary value of the intermediate variableU˙l,m can be extracted from (24) and is given by the following formulas:
(28)U0,m*=[I+r2Aδy2]U0,mn+1,UL+1,m*=[I+r2Aδy2]UL,mn+1.


## 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(29)∂2u∂x2=1h2δx21+(1/12)δx2ul,m+O(h4),∂2u∂y2=1h2δy21+(1/12)δy2ul,m+O(h4).By using these approximations together with Crank-Nicolson for the time direction, we get the numerical scheme(30)1k[Ul,mn+1-Ul,mn]+12h2A[δx21+(1/12)δx2+δy21+(1/12)δy2]×[Ul,mn+1+Ul,mn]+B(U-)U-=0
which can be written as
(31)[I+r2A(δx21+(1/12)δx2+δy21+(1/12)δy2)]Ul,mn+1=[I-r2A(δx21+(1/12)δx2+δy21+(1/12)δy2)]Ul,mn-kB(U-)U-,
and this can be written in the factored form as
(32)[I+r2A(δx21+(1/12)δx2)]=×[I+r2A(δy21+(1/12)δy2)]Ul,mn+1=[I-r2A(δx21+(1/12)δx2)]=×[I-r2A(δy21+(1/12)δy2)]Ul,mn-kB(U-)U-.The difference between (31) and (32) is
(33)-r24[δx21+(1/12)δx2][δy21+(1/12)δy2](Ul,mn+1-Ul,mn)=-k24h4[δx21+(1/12)δx2][δy21+(1/12)δy2]=×(k(∂u∂t)+O(k3))=-k34[∂5u∂2x∂2y∂t]+O(k3h4)
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)δx2)(1+(1/12)δy2)(34)[(1+112δx2)I+r2Aδx2][(1+112δy2)I+r2Aδy2]Ul,mn+1=[(1+112δx2)I-r2Aδx2]=×[(1+112δy2)I-r2Aδy2]Ul,mn-k(1+112δx2)(1+112δy2)B(U-)U-.The fourth-order D’Yakonov ADI-like scheme in this case can be displayed as(35)[(1+112δx2)I+r2Aδx2]U˙=[(1+112δx2)I-r2Aδx2][(1+112δy2)I-r2Aδy2]Ul,mn-k(1+112δx2)(1+112δy2)B(U-)U-(36)[(1+112δy2)I+r2Aδy2]Ul,mn+1=U˙.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 variableU˙l,m in (35) can be given by the following formulas:
(37)U˙0,m=[(1+112δy2)I+r2Aδy2]U0,mn+1,U˙L+1,m=[(1+112δy2)I+r2Aδy2]UL+1,mn+1.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
(38)[(1+σδx2)I+r2Aδx2]U˙=[(1+σδx2)I-r2Aδx2][(1+σδy2)I-r2Aδy2]Ul,mn-k(1+σδx2)(1+σδy2)B(U-)U-,[(1+σδy2)I+r2Aδy2]Ul,mn+1=U˙
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(39)[(1+σδx2)I+r2Aδx2][(1+σδy2)I+r2Aδy2]Ul,mn+1=[(1+σδx2)I-r2Aδx2][(1+σδy2)I-r2Aδy2]Ul,mn-k2ω(1+σδx2)(1+σδy2)U-,
where ω is constant and σ=0,1/12.Suppose that the numerical solution can be expressed by Fourier series, whose typical term is(40)Ul,mn=GnWnexp(iβ1lh+iβ2mh),
where i=-1,Gn is the amplification matrix at time level n, and β1,β2 are the wave numbers in x and y directions. Substituting (40) into (39) will lead to the matrix equation
(41)[I-2rμxA][I-2rμyA]G=[I+2rμxA][I+2rμyA]-k2ω[G+I],
where
(42)μx=sin2(β1lh),μy=sin2(β2mh).The eigenvalues of the matrixG are given by
(43)eigenvaluesofG=[1+2irμx][1+2irμy]-(k/2)ω[1-2irμx][1-2irμy]+(k/2)ω.It can be easily shown that the modulus of the maximum eigenvalue of the amplification matrixG 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 theL∞ 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(44)ψ(x,y,t)=Asech(ξ)exp(iη),ϕ(x,y,t)=Bsech(ξ)exp(iη),
where
(45)ξ=β1x+β2y-vt,η=-κ1x-κ2y+ωt+θ,A=B=[2δ(β12+β22)(1+α)]1/2,
at t=0. The following parameters are used
(46)h=0.2,k=0.001,β1=0.5,β2=1.0,k1=0.5,k2=1.0,α=1,δ=12,-10≤x,y≤10.The boundary conditions are extracted from the exact solution (44).Tables1, 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 t=0. Figure 2 shows the numerical solution at t=1.Table 1
Single soliton withσ=0.
Iter
T
ER
I
0
0
0
24.99989
2
0.2
0.00435
25.00086
2
0.4
0.00886
25.00291
2
0.6
0.01365
25.00586
2
0.8
0.01893
25.00973
2
1.0
0.02439
25.01458Table 2
Single soliton withσ=1/6.
Iter
T
ER
I
0
0
0
24.99989
2
0.2
0.00436
24.99882
2
0.4
0.00859
24.99664
2
0.6
0.01277
24.99351
2
0.8
0.01750
24.98938
2
1.0
0.02206
24.98417Table 3
Single soliton withσ=1/12.
Iter
T
ER
I
0
0
0
24.99989
2
0.2
0.00011
24.99989
2
0.4
0.00018
24.99987
2
0.6
0.00025
24.99981
2
0.8
0.00030
24.99971
2
1.0
0.00036
24.99952Figure 1
Single solitonβ1=κ1=0.5, β2=κ2=1, α=1, t=0.Figure 2
Single solitonβ1=κ1=0.5, β2=κ2=1, α=1, t=1.


### 4.2. Example 2

In this example we will choose the initial condition [1]
(47)ψ(x,y,0)=α12exp(-η2)sech×[Re(κ1x+l1y)+η2]exp(iIm(κ1x+l1y))ϕ(x,y,0)=β12exp(-η2)sech×[Re(κ1x+l1y)+η2]exp(iIm(κ1x+l2y)),
where
(48)exp(η)=12|α1|2+|β1|2(κ1+κ1*)2+(l1+l1*)2,
which represent two solitons of different amplitudes. The parameters α1, β1, κ1, and l1 are complex parameters. In this test, we choose the parameters
(49)h=0.2,k=0.01,xl=-10,α1=1,β1=0.5,κ1=l1=1+i.In Figures3 and 4, we display the numerical solution of |ψ|2 at t=0 and t=1, while in Figures 5 and 6 we display the numerical solution of |ϕ|2 at t=0 and t=1.Figure 3
Soliton solution|ψ|2 at t=0.Figure 4
Soliton solution|ψ|2 at t=1.Figure 5
Soliton solution|ϕ|2 at t=0.Figure 6
Soliton solution|ϕ|2 at t=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(50)ψ(x,y,0)=Aexp(-(x-d)2-y2)exp(-iln(cosh(x2+y2))),ϕ(x,y,0)=Bexp(-(x+d)2-y2)exp(-iln(cosh(x2+y2))),
which represent two solitons moving in the opposite direction and centered at d and -d, respectively. In this test we choose the parameters
(51)h=0.2,k=0.01,xl=-10,d=5.0,A=1,B=0.8,t=0,…,1.5.The interaction scenario is given in Figures7 and 8. In Figure 7 we show the two solitons with two different amplitudes at t=0. In Figure 8 we display the interaction where the two solitons interact at t=0.5, and, in Figure 9, we display the two solitons after the interaction at t=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.
T
I
1
I
2
0.0
1.57080
1.00531
0.5
1.57086
1.00533
1.0
1.57065
1.00522
1.5
1.57060
1.00520Figure 7
Interaction of two solitons|ψ|2 and |ϕ|2 at t=0.Figure 8
Interaction of two solitons|ψ|2 and |ϕ|2 at t=0.5.Figure 9
Interaction of two solitons|ψ|2 and |ϕ|2 at t=1.5.


### 4.4. Example 4

The system under consideration generates a progressive plane wave solutions [7]:
(52)ψ(x,y,t)=Aexp(i(κ1x+κ2y+ωt)),ϕ(x,y,t)=Aexp(i(κ1x+κ2y+ωt)),
where
(53)ω=k12+k22-(1+α)|A|2.Our numerical experiments are conducted in the domain[0,2π]×[0,2π] with
(54)A=1,α=1,k1=k2=1,δ=1,k=0.001,h=π50.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 t=0, and Figure 11 displays, the solution at t=1.Table 5
Periodic solutionk1=k2=α=1.
Iter
T
ER
I
0
0.0
0.000000
39.478420
2
0.2
0.000003
39.478420
2
0.4
0.000005
39.478420
2
0.6
0.000007
39.478420
2
0.8
0.000009
39.478420
2
1.0
0.000008
39.478420Figure 10
Initial condition withk1=k2=α=1, at t=0.Figure 11
Numerical solution withk1=k2=α=1, at t=1.Table4 displays the accuracy of the scheme and preserves the conserved quantities.


## 4.1. Example 1

In this example, we choose the initial conditions from the exact solutions(44)ψ(x,y,t)=Asech(ξ)exp(iη),ϕ(x,y,t)=Bsech(ξ)exp(iη),
where
(45)ξ=β1x+β2y-vt,η=-κ1x-κ2y+ωt+θ,A=B=[2δ(β12+β22)(1+α)]1/2,
at t=0. The following parameters are used
(46)h=0.2,k=0.001,β1=0.5,β2=1.0,k1=0.5,k2=1.0,α=1,δ=12,-10≤x,y≤10.The boundary conditions are extracted from the exact solution (44).Tables1, 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 t=0. Figure 2 shows the numerical solution at t=1.Table 1
Single soliton withσ=0.
Iter
T
ER
I
0
0
0
24.99989
2
0.2
0.00435
25.00086
2
0.4
0.00886
25.00291
2
0.6
0.01365
25.00586
2
0.8
0.01893
25.00973
2
1.0
0.02439
25.01458Table 2
Single soliton withσ=1/6.
Iter
T
ER
I
0
0
0
24.99989
2
0.2
0.00436
24.99882
2
0.4
0.00859
24.99664
2
0.6
0.01277
24.99351
2
0.8
0.01750
24.98938
2
1.0
0.02206
24.98417Table 3
Single soliton withσ=1/12.
Iter
T
ER
I
0
0
0
24.99989
2
0.2
0.00011
24.99989
2
0.4
0.00018
24.99987
2
0.6
0.00025
24.99981
2
0.8
0.00030
24.99971
2
1.0
0.00036
24.99952Figure 1
Single solitonβ1=κ1=0.5, β2=κ2=1, α=1, t=0.Figure 2
Single solitonβ1=κ1=0.5, β2=κ2=1, α=1, t=1.


## 4.2. Example 2

In this example we will choose the initial condition [1]
(47)ψ(x,y,0)=α12exp(-η2)sech×[Re(κ1x+l1y)+η2]exp(iIm(κ1x+l1y))ϕ(x,y,0)=β12exp(-η2)sech×[Re(κ1x+l1y)+η2]exp(iIm(κ1x+l2y)),
where
(48)exp(η)=12|α1|2+|β1|2(κ1+κ1*)2+(l1+l1*)2,
which represent two solitons of different amplitudes. The parameters α1, β1, κ1, and l1 are complex parameters. In this test, we choose the parameters
(49)h=0.2,k=0.01,xl=-10,α1=1,β1=0.5,κ1=l1=1+i.In Figures3 and 4, we display the numerical solution of |ψ|2 at t=0 and t=1, while in Figures 5 and 6 we display the numerical solution of |ϕ|2 at t=0 and t=1.Figure 3
Soliton solution|ψ|2 at t=0.Figure 4
Soliton solution|ψ|2 at t=1.Figure 5
Soliton solution|ϕ|2 at t=0.Figure 6
Soliton solution|ϕ|2 at t=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(50)ψ(x,y,0)=Aexp(-(x-d)2-y2)exp(-iln(cosh(x2+y2))),ϕ(x,y,0)=Bexp(-(x+d)2-y2)exp(-iln(cosh(x2+y2))),
which represent two solitons moving in the opposite direction and centered at d and -d, respectively. In this test we choose the parameters
(51)h=0.2,k=0.01,xl=-10,d=5.0,A=1,B=0.8,t=0,…,1.5.The interaction scenario is given in Figures7 and 8. In Figure 7 we show the two solitons with two different amplitudes at t=0. In Figure 8 we display the interaction where the two solitons interact at t=0.5, and, in Figure 9, we display the two solitons after the interaction at t=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.
T
I
1
I
2
0.0
1.57080
1.00531
0.5
1.57086
1.00533
1.0
1.57065
1.00522
1.5
1.57060
1.00520Figure 7
Interaction of two solitons|ψ|2 and |ϕ|2 at t=0.Figure 8
Interaction of two solitons|ψ|2 and |ϕ|2 at t=0.5.Figure 9
Interaction of two solitons|ψ|2 and |ϕ|2 at t=1.5.


## 4.4. Example 4

The system under consideration generates a progressive plane wave solutions [7]:
(52)ψ(x,y,t)=Aexp(i(κ1x+κ2y+ωt)),ϕ(x,y,t)=Aexp(i(κ1x+κ2y+ωt)),
where
(53)ω=k12+k22-(1+α)|A|2.Our numerical experiments are conducted in the domain[0,2π]×[0,2π] with
(54)A=1,α=1,k1=k2=1,δ=1,k=0.001,h=π50.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 t=0, and Figure 11 displays, the solution at t=1.Table 5
Periodic solutionk1=k2=α=1.
Iter
T
ER
I
0
0.0
0.000000
39.478420
2
0.2
0.000003
39.478420
2
0.4
0.000005
39.478420
2
0.6
0.000007
39.478420
2
0.8
0.000009
39.478420
2
1.0
0.000008
39.478420Figure 10
Initial condition withk1=k2=α=1, at t=0.Figure 11
Numerical solution withk1=k2=α=1, at t=1.Table4 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.



---

*Source: 497439-2013-09-16.xml*