# An Efficient Parallel Laplace Transform Spectral Element Method for the Two-Dimensional Schrödinger Equation

Authors: Wenting Shao
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, a parallel Laplace transform spectral element method for the linear Schrödinger equation and the cubic nonlinear Schrödinger equation is introduced. The main advantage of this new scheme is that it has the high accuracy of spectral methods and the parallel efficiency of the Laplace transform method. The key idea is twofold. First, Laplace transform is utilized to eliminate time dependency and the subsequent boundary value problems are discretized by using spectral element method. Second, a parallel Laplace transform method is developed to improve the computational efficiency. For the cubic nonlinear Schrödinger equation, the increment linearization technique is employed to deal with the nonlinear terms. Numerical experiments are addressed to demonstrate the accuracy and effectiveness of the proposed method. Compared with the Crank–Nicolson scheme, the CPU time is greatly saved by the parallel Laplace transform for both linear and nonlinear problems.


---

## Body

## 1. Introduction

In physics, especially in quantum mechanics, the Schrödinger equation is used to describe how the quantum state of a physical system changes in time. Currently, this equation has been widely applied in many areas, for example, optics, underwater acoustics, material science, plasma physics, electromagnetic wave propagation, design of optoelectronic devices, and Bose–Einstein condensation [1, 2].Since in general, finding the analytical solutions of the two-dimensional Schrödinger equation can be a hard task, one has to make use of numerical procedures and extensive numerical methods haven been studied in recent years [3]. Concerning the time discretization, different directions can be investigated, for example, the finite difference scheme [4–6], the Crank–Nicolson (CN) scheme [7–9], the time-splitting scheme [10, 11], the Runge–Kutta (RK) scheme [1, 2, 12–14], and the integrating factor method [15]. Among these schemes, the CN scheme is implicit and unconditionally stable. The main disadvantage of the CN scheme is that it needs to solve a nonlinear system at each time step for the nonlinear problem. Although the nonlinearity can be dealt with using some fixed-point iteration method, it might be rather taxing and time consuming. A possible alternative way of handling the nonlinearity is the time-splitting scheme. Its basic idea is to split NLSE into a linear part and a nonlinear part. The linear subproblem can be solved by the CN method in time. The nonlinear subproblem in some situations where the potential and nonlinearity can be integrated analytically can be solved exactly. Thus, the scheme is explicit and unconditionally stable. These temporal schemes can be coupled with high-order spatial discretizations like spectral or pseudo-spectral methods. Nonetheless, most of the schemes are low-order accuracy in time which limits the use of large time steps in numerical implementation and the global accuracy. It is therefore necessary to develop high-order temporal discretizations to improve the accuracy and efficiency.The Laplace transform (LT) represents a very effective tool for solving time-dependent problems in science and engineering [16]. It removes the time dependency and converts the original problem into a series of boundary value problems (BVPs). Chen and Wu et al. combined LT with the rational spectral collocation method to solve the advection-diffusion equation [17] and one-dimensional Burgers’ equation (18). A coupled method of LT and Legendre wavelets was presented in [19] for one-dimensional Klein–Gordon equations. In [20], some different types of unsteady flows of the third grade fluid were studied by the use of LT in time. However, LT never really becomes popular in computational work, especially in conjunction with spectral element method. There are two reasons: one is the numerical inversion of LT that is an ill-posed problem when the transform is known only as a real-valued function, and the other one is the nonlinear term that makes the direct use of LT challenging. For the first issue, the idea proposed by Talbot [21] is adopted. The technique combines the trapezoid rule or the midpoint rule with the contour integration by using the complex arithmetic. The convergence rate for the optimized quadrature formula could reach a geometric rate Oe−cM [22], where M is the number of sample points. For the second issue, several techniques have been developed to overcome the limitation, such as the homotopy perturbation method [23] and Adomian decomposition method [24, 25]. The Adomian decomposition method is done by expanding the nonlinear term as an infinite series of the Adomian polynomials, which can be calculated by a recursive relationship. However, it causes a lot of extra cost for solving time-dependent problems. The increment linearization method was first proposed by Wu et al. [26] and then integrated with LT to solve Burgers’ equation (18). Its key idea is to define an increment function at each time step: Δux,t=ux,t−ux,tn, tn≤t≤tn+1, linearize the equation by omitting the nonlinear terms of the order Oτγ,γ≥2, and finally get the approximation by ux,tn+1≈ux,tn+Δux,tn+1. In this work, we will employ the increment linearization method to deal with the nonlinear terms.Concerning the spatial discretization, many techniques have been developed including the finite difference method [11], finite element method [8, 27], Legendre spectral method [13], Fourier/Chebyshev pseudo-spectral methods [1, 7, 10, 12, 15, 28], Chebyshev wavelet method [29], radial basis meshless collocation method [4], radial basis function-finite difference method [2, 14], and Sinc collocation and Sinc Galerkin methods [5]. It is well known that the spectral method is a very useful tool for the solution of PDEs because of its spectral accuracy, when the geometry of the problem is smooth and regular. In [13, 28], Galerkin–Legendre spectral method and Chebyshev pseudo-spectral method were combined with the RK method to solve the two-dimensional linear Schrödinger equation. However, there have been few references about integrating spectral element method (SEM) with LT to solve time-dependent PDEs. This motivates our study to introduce SEM for solving two-dimensional Schrödinger equations. SEM was first proposed by Patera [30] for the numerical solution of the incompressible Navier–Stokes equations, and it combines the high accuracy of spectral methods and the flexibility of finite element method. The approximate result of SEM provides global spectral convergence, not at the expense of “uncontrollable” resolution. Compared with conventional spectral or pseudo-spectral methods, SEM requires weaker smoothness in operations and the fine grid can be obtained by increasing either the interpolation degree (p-convergence) in each element or the amount of elements (h-convergence). In this paper, the choice of SEM makes high-order accuracy in space and it is compatible with the time-discrete method LT.In this paper, we consider the two-dimensional linear Schrödinger equation (LSE) [28].(1)i∂Φx,t∂t+∇⋅K∇Φx,t+wxΦx,t+fx,t=0,and the cubic nonlinear Schrödinger equation (NLSE) [4].(2)i∂Φx,t∂t+∇⋅K∇Φx,t+βΦx,t2Φx,t+wxΦx,t=0,in the rectangle domainΩ and 0<t≤T with the initial condition(3)Φx,0=Φ0x,inΩ,and the homogenous Dirichlet boundary conditions(4)Φx,t=0,on∂Ω×0,T.Here,i2=−1, x=x,y, Φx,t is the complex-valued function, wx is a given real-valued potential function, and K and β are given real constants. The cubic NLSE admits the mass conservation law [3](5)Mt=∫ΩΦx,t2dσ=∫ΩΦ0x2dσ=M0.and the energy conservation law(6)Et=∫ΩK∇Φx,t2−wxΦx,t2−β2Φx,t4dσ=∫ΩK∇Φ0x2−wxΦ0x2−β2Φ0x4dσ=E0.for t>0.A coupled method of Laplace transform and spectral element method, named LTSEM, is presented for solving the LSE and the cubic NLSE. The application of LT converts the original time-dependent problem into a series of BVPs with different LT parameters. SEM is applied to solve the subsequent equations. For the cubic NLSE, the increment linearization technique is introduced to deal with the nonlinear terms. It removes the stumbling block of using LT to solve the nonlinear problem. Unlike the conventional time-marching method, the LT method has an advantage of excellent parallel efficiency, but it is not addressed in detail by numerical experiments in literature. Inspired by this, we develop a parallel LTSEM (pLTSEM) to improve the computational efficiency of LTSEM and it can be easily implemented with the MATLAB parallel computing toolbox. This is the main contribution of this work.The paper is organized as follows. In Section2, we consider LSE and convert it into a series of BVPs by using LT. In Section 3, Talbot’s method is briefly reviewed for the numerical inversion of LT. In Section 4, SEM is utilized to solve the BVPs generated by using LT and a detailed description of pLTSEM is given for LSE. In Section 5, the cubic NLSE is considered. The increment linearization method is introduced to deal with the nonlinearity. The description of pLTSEM for this case is given. In Section 6, numerical examples are tested to show the performance of our proposed method. Finally, in Section 7, a brief conclusion is drawn and some comments are included.


## 2. Laplace Transform Method

In this section, we will introduce the LT method to convert the LSE into a series of BVPs.In practice, we decompose the complex-valued functionΦx,t into its real and imaginary parts by writing Φx,t=ux,t+ivx,t, and set fx,t=fux,t+ifvx,t, where ux,t, vx,t, fux,t and fvx,t are real-valued functions. Problem (1) is transformed into the following coupled real-valued differential equations(7)utx,t+∇⋅K∇vx,t+wxvx,t+fvx,t=0,vtx,t−∇⋅K∇ux,t−wxux,t−fux,t=0,inΩ.with the initial condition(8)ux,0=u0x,vx,0=v0x,inΩ.and the homogenous boundary condition(9)ux,t=0,vx,t=0,on∂Ω×0,T.Letux,t be a piecewise continuous function and s be a complex variable. The Laplace transform u^x,s of ux,t is defined by(10)u^x,s=ℒux,t=∫0∞e−stux,tdt.and it satisfies(11)ℒutx,t=su^x,s−u0x.Taking the Laplace transform of (7) gives(12)su^x,s−u0x+∇⋅K∇v^x,s+wxv^x,s+f^vx,s=0,sv^x,s−v0x−∇⋅K∇u^x,s−wxu^x,s−f^ux,s=0,inΩ.with the boundary conditions(13)u^x,s=0,v^x,s=0,on∂Ω,where f^ux,s and f^vx,s are the Laplace transforms of fux,t and fvx,t, respectively.Next, we set the decompositionss=sR+isI, u^x,s=u^Rx+iu^Ix, v^x,s=v^Rx+iv^Ix, f^ux,s=f^uRx+if^uIx, and f^vx,s=f^vRx+if^vIx, and then system (12) can be recast to the following linear system satisfied by u^R, u^I, v^R, and v^I(14)sRu^R−sIu^I+∇⋅K∇v^R+wv^R=−f^vR+u0,sIu^R+sRu^I+∇⋅K∇v^I+wv^I=−f^vI,−∇⋅K∇u^R−wu^R+sRv^R−sIv^I=f^uR+v0,−∇⋅K∇u^I−wu^I+sIv^R+sRv^I=f^uI,inΩ.with the boundary conditions(15)u^R=u^I=v^R=v^I=0,on∂Ω.Usually, one solves equations (14) and (15) with different LT parameters sm and gets the approximations u^x,sm and v^x,sm. Then, by doing a numerical inversion of LT, the numerical solutions in the original space could be obtained.


## 3. Numerical Inversion of Laplace Transform

In the past decade, several results [21, 22, 31] have been obtained for numerical inversion of LT, among which Talbot’s method [21] has received a lot of attention and constitutes the mainstream thanks to the high accuracy and relatively low computational complexity. In this section, we employ Talbot’s method to realize the numerical inversion of LT.The inversion Laplace transform ofu^x,s is defined by(16)ux,t=12πi∫Bestu^x,sds,which is called the Bromwich integral, with B as the Bromwich line Res=σ>σ0. σ0 is the maximum value of the real part of all singularities of u^x,s.In Talbot’s method, the Bromwich line is deformed into a curveΓ that begins and ends in the left half-plane, such that Res⟶−∞ on the contour. Owing to the exponential factor est, the integrand decays rapidly on it. It has been proved [21] that Talbot’s method is applicable to both a wide range of Laplace transform functions and values of t where we need to compute ux,t.Talbot’s contour is parameterized by(17)Γ:sθ=σ+μθcotθ+νiθ,−π≤θ≤π,where σ, μ, and ν are real parameters that determine the geometry of the curve. Both μ and ν are positive. Assume that Γ lies in the region of analyticity of u^x,s and use the Cauchy integral theorem; the Bromwich integral (16) with Talbot’s contour becomes(18)ux,t=12πi∫−ππesθtu^x,sθs′θdθ.which can be approximated by the midpoint formula. Define integral points as(19)θm=2m+1π2M,m=−M,…,M−1.and allow sm=sθm,sm′=s′θm; the approximation to (18) is given by(20)ux,t≃12Mi∑m=−MM−1esmtsm′u^x,sm=−1MIm∑m=0M−1esmtsm′u^x,sm.where the second equation is deduced by symmetry.Remark 1.
It was analyzed in [21] that Talbot’s method with fixed parameters converges at a subgeometric rate of Oe−cM. Further, Weideman [22] used the saddle point method to give the error estimate, that is, EMt=Oeσ+μt−2πμtM,M⟶∞. This estimate suggests the strategy of choosing σ, and μ∝M/t could lead to an improved convergence rate Oe−cM. Besides, Weideman found optimal parameters for some parabolic problems: when the following Talbot’s contour(21)Γ:sθ=Mt−0.4814+0.4774θcot0.7409θ+0.2698iθ.is used as the Bromwich line, the convergence rate is given by Oe−1.90M.
Generally, we give a uniform partition of the time interval and denote discrete timetn=nτ, where τ is the time step. At each time step, LT with its numerical inversion is used for the temporal discretization. Accordingly, t will be replaced by τ in (21).
Next, we plan to discretize equations (14) and (15) by using SEM in space and get the approximations of u^x,sm and v^x,sm.


## 4. Spatial Discretization and Parallel Computation

In this section, the Legendre SEM is briefly recalled for the spatial discretization and a detailed description of pLTSEM is given for LSE.


### 4.1. Notations and Discrete Schemes

We denote by⋅,⋅ and ⋅ the inner product and the norm of the space L2Ω. Denote HαΩα≥0 to be the Sobolev space with the norm ⋅HαΩ and introduce the space: H01Ω=ϕx∈H1Ω:ϕx|∂Ω=0.The weak form of equations (14) and (15) can be written as follows: for any given sm, find u^R,u^I,v^R,v^I∈H01Ω4, such that(22)smRu^R,ϕ−smIu^I,ϕ−K∇v^R,∇ϕ+wv^R,ϕ=−f^vR+u0,ϕ,smIu^R,ϕ+smRu^I,ϕ−K∇v^I,∇ϕ+wv^I,ϕ=−f^vI,ϕ,K∇u^R,∇ϕ−wu^R,ϕ+smRv^R,ϕ−smIv^I,ϕ=f^uR+v0,ϕ,K∇u^I,∇ϕ−wu^I,ϕ+smIv^R,ϕ+smRv^I,ϕ=f^uI,ϕ,∀ϕ∈H01Ω.We takeTh be a uniform regular partition of quadrilateral element of domain Ω: Ω¯=∪i=1neΩ¯i, ∩i=1neΩi=∅. Let I=−1,1 be the reference interval and denote Ωst=−1,12. Define ℙNI to be the space of polynomials with degree at most N on the interval I, and denote QNΩst=ℙNI⊗ℙNI. The modal bases of QNΩst are defined by(23)φ˜j,nξ,η=φjξφnη,0≤j,n≤N,with(24)φ0ξ=1−ξ2,φNξ=1+ξ2,φjξ=Lj−1ξ−Lj+1ξ,1≤j≤N−1,where Ljξ is j th Legendre polynomial and φnη has the same definition with φjξ.For any elementΩi=a,b×c,d, we can use the following mapping functions to transfer Ωi to Ωst and its inverse(25)ξx=2x−ab−a−1,ηy=2y−cd−c−1,x,y∈Ωi.The modal bases ofQNΩi can be easily defined by using (23) and (25).Define the approximation space(26)WN0Ω=ϕNx∈H01Ω:ϕNx|Ωi∈QNΩi,∀Ωi∈Th.The discrete weak form of (22) is as follows: for any given sm, find u^NR,u^NI,v^NR,v^NI∈WN0Ω4, such that(27)smRu^NR,ϕN−smIu^NI,ϕN−K∇v^NR,∇ϕN+wv^NR,ϕN=−f^vR+u0,ϕN,smIu^NR,ϕN+smRu^NI,ϕN−K∇v^NI,∇ϕN+wv^NI,ϕN=−f^vI,ϕN,K∇u^NR,∇ϕN−wu^NR,ϕN+smRv^NR,ϕN−smIv^NI,ϕN=f^uR+v0,ϕN,K∇u^NI,∇ϕN−wu^NI,ϕN+smIv^NR,ϕN+smRv^NI,ϕN=f^uI,ϕN,∀ϕN∈WN0Ω.On any elementΩi, we assume the approximations(28)u^NRx=∑j=0N∑n=0Nu˜j,nRφ˜j,nx,u^NIx=∑j=0N∑n=0Nu˜j,nIφ˜j,nx,v^NRx=∑j=0N∑n=0Nv˜j,nRφ˜j,nx,v^NIx=∑j=0N∑n=0Nv˜j,nIφ˜j,nx.with the expansion coefficientsu˜j,nRj,n=0N, u˜j,nIj,n=0N, v˜j,nRj,n=0N, and v˜j,nIj,n=0N. Here, we still use the notation φ˜ to represent the modal bases.The local linear system arising from (27) can be written in a partitioned matrix form as(29)smRMi−smIMi−KAi+CiOismIMismRMiOi−KAi+CiKAi−CiOismRMi−smIMiOiKAi−CismIMismRMiU˜i=b˜i.where the element vector U˜i is composed of expansion coefficients(30)U˜i=Vecu˜j,nRj,n=0NT,Vecu˜j,nIj,n=0NT,Vecv˜j,nRj,n=0NT,Vecv˜j,nIj,n=0NTT.and “Vec” is the matrix vectorization operator. The element matrices Mi, Ai, and Ci are of dimension N+12×N+12 and could be calculated as follows:(31)Mi=φ˜j′,n′,φ˜j,nΩij,n,j′,n′,Ai=∇φ˜j′,n′,∇φ˜j,nΩij,n,j′,n′,Ci=wφ˜j′,n′,φ˜j,nΩij,n,j′,n′.and Oi=0N+12×N+12. The element vectors(32)b˜i=b˜i1T,b˜i2T,b˜i3T,b˜i4TT.with the components(33)b˜i1=Vec−f^vR+u0,φ˜j′,n′Ωij′,n′,b˜i2=Vec−f^vI,φ˜j′,n′Ωij′,n′,b˜i3=Vecf^uR+v0,φ˜j′,n′Ωij′,n′,b˜i4=Vecf^uI,φ˜j′,n′Ωij′,n′.Here, the Gauss integral formula is adopted in the quadratures (31) and (33).Assembling the contributions from all elements leads to the global matricesM, A, C and vector b˜, and finally, arrive at the following matrix form(34)smRM−smIM−KA+COsmIMsmRMO−KA+CKA−COsmRM−smIMOKA−CsmIMsmRMU˜=b˜,where the vector U˜ contains the expansion coefficients of all elements. The right-hand side(35)b˜=b˜1T,b˜2T,b˜3T,b˜4TT.where b˜1,2,3,4 could be calculated by (33) defined on the whole domain.In the final of this subsection, we present the following theorem that will be used in the analysis of the existence and uniqueness of the solution of (34).Theorem 1 (see [32]).
The matrixL∈ℝn×n is positive definite if and only if its symmetric part W=L+LT/2 is positive definite. Besides, W is positive definite if and only if its eigenvalues are positive.
Working back from the matrix formulation to the functional form of the integral operator, it is easy to show that matricesA and C are symmetric, and the mass matrix M is symmetric and positive definite [33]. Denote the coefficient matrix of (34) by L; then, we have(36)L+LT2=smRMOOOOsmRMOOOOsmRMOOOOsmRM,
and it is a block diagonal matrix.
It is a practical choice that avoids sampling the integrand atsmR=0. If smR>0, since the mass matrix M is positive definite, we obtain that its eigenvalues are all positive, then the block diagonal matrix L+LT/2 is also positive definite. From Theorem 1, it is direct to get that the coefficient matrix L is positive definite; thus, the existence and uniqueness of the solution of (34) are guaranteed. If smR<0, we consider the coefficient matrix L˜=−L and obtain the same conclusion.
Moreover, since the modal bases (23) are constructed by using the orthogonal Legendre polynomials, the matrices A and M are typically sparse although they may have a full bandwidth [33]. Thus, (34) is large but very sparse. In this paper, we adopt the built-in algorithm in MATLAB to solve it, which is already optimized for solving the system of linear equations.


### 4.2. Parallel LTSEM for the Linear Problem

In Algorithm 1, a description of the pLTSEM for solving LSE is given. The underlying idea of the pLTSEM is straightforward, and the calculation of (34) with different LT parameters sm is independent. This step can be handled more efficiently by the parallel procedure that means the computing procedure can be implemented simultaneously in different CPU cores. Algorithm 1 is divided into three stages of initialization, two-level iteration followed by a solution update step. In the initialization stage, since the matrices M, A, and C depend only on the spatial grid size and the order of polynomial expansion, they can be precomputed before entering into the time loop. The iterative procedure is composed of two loops: the outer loop and the inner loop. The inner loop is a parallel computing process. Its main purpose is to solve the discretized problems on the current time level. The parallel computation is carried out by using the MATLAB parallel computing toolbox, where the parfor loop is a tool for parallel computing single loop. The degree of parallelism is related to the number of CPU cores. The computational time can be saved significantly under the condition of the balance between the computation cost and the data communication. The outer loop is used to forward the solution process to the next time step until the final time is reached.Step 0. for time level n=0 (i.e., t=0),generate the initial mesh, establish matrices M, A, C.Step 1. for time level n≥1,%parallel computing:parfor (m=0; m<M; m++){for each given sm,establish vectors b˜1,2,3,4,solve the global system equation (30) to get coefficient U˜,determineu^Nx,sm, v^Nx,sm.}make numerical inversion of LT to getuNx,tn+1, vNx,tn+1.Step 2. set n=n+1,update the right-hand side (33),goto Step 1end for


## 4.1. Notations and Discrete Schemes

We denote by⋅,⋅ and ⋅ the inner product and the norm of the space L2Ω. Denote HαΩα≥0 to be the Sobolev space with the norm ⋅HαΩ and introduce the space: H01Ω=ϕx∈H1Ω:ϕx|∂Ω=0.The weak form of equations (14) and (15) can be written as follows: for any given sm, find u^R,u^I,v^R,v^I∈H01Ω4, such that(22)smRu^R,ϕ−smIu^I,ϕ−K∇v^R,∇ϕ+wv^R,ϕ=−f^vR+u0,ϕ,smIu^R,ϕ+smRu^I,ϕ−K∇v^I,∇ϕ+wv^I,ϕ=−f^vI,ϕ,K∇u^R,∇ϕ−wu^R,ϕ+smRv^R,ϕ−smIv^I,ϕ=f^uR+v0,ϕ,K∇u^I,∇ϕ−wu^I,ϕ+smIv^R,ϕ+smRv^I,ϕ=f^uI,ϕ,∀ϕ∈H01Ω.We takeTh be a uniform regular partition of quadrilateral element of domain Ω: Ω¯=∪i=1neΩ¯i, ∩i=1neΩi=∅. Let I=−1,1 be the reference interval and denote Ωst=−1,12. Define ℙNI to be the space of polynomials with degree at most N on the interval I, and denote QNΩst=ℙNI⊗ℙNI. The modal bases of QNΩst are defined by(23)φ˜j,nξ,η=φjξφnη,0≤j,n≤N,with(24)φ0ξ=1−ξ2,φNξ=1+ξ2,φjξ=Lj−1ξ−Lj+1ξ,1≤j≤N−1,where Ljξ is j th Legendre polynomial and φnη has the same definition with φjξ.For any elementΩi=a,b×c,d, we can use the following mapping functions to transfer Ωi to Ωst and its inverse(25)ξx=2x−ab−a−1,ηy=2y−cd−c−1,x,y∈Ωi.The modal bases ofQNΩi can be easily defined by using (23) and (25).Define the approximation space(26)WN0Ω=ϕNx∈H01Ω:ϕNx|Ωi∈QNΩi,∀Ωi∈Th.The discrete weak form of (22) is as follows: for any given sm, find u^NR,u^NI,v^NR,v^NI∈WN0Ω4, such that(27)smRu^NR,ϕN−smIu^NI,ϕN−K∇v^NR,∇ϕN+wv^NR,ϕN=−f^vR+u0,ϕN,smIu^NR,ϕN+smRu^NI,ϕN−K∇v^NI,∇ϕN+wv^NI,ϕN=−f^vI,ϕN,K∇u^NR,∇ϕN−wu^NR,ϕN+smRv^NR,ϕN−smIv^NI,ϕN=f^uR+v0,ϕN,K∇u^NI,∇ϕN−wu^NI,ϕN+smIv^NR,ϕN+smRv^NI,ϕN=f^uI,ϕN,∀ϕN∈WN0Ω.On any elementΩi, we assume the approximations(28)u^NRx=∑j=0N∑n=0Nu˜j,nRφ˜j,nx,u^NIx=∑j=0N∑n=0Nu˜j,nIφ˜j,nx,v^NRx=∑j=0N∑n=0Nv˜j,nRφ˜j,nx,v^NIx=∑j=0N∑n=0Nv˜j,nIφ˜j,nx.with the expansion coefficientsu˜j,nRj,n=0N, u˜j,nIj,n=0N, v˜j,nRj,n=0N, and v˜j,nIj,n=0N. Here, we still use the notation φ˜ to represent the modal bases.The local linear system arising from (27) can be written in a partitioned matrix form as(29)smRMi−smIMi−KAi+CiOismIMismRMiOi−KAi+CiKAi−CiOismRMi−smIMiOiKAi−CismIMismRMiU˜i=b˜i.where the element vector U˜i is composed of expansion coefficients(30)U˜i=Vecu˜j,nRj,n=0NT,Vecu˜j,nIj,n=0NT,Vecv˜j,nRj,n=0NT,Vecv˜j,nIj,n=0NTT.and “Vec” is the matrix vectorization operator. The element matrices Mi, Ai, and Ci are of dimension N+12×N+12 and could be calculated as follows:(31)Mi=φ˜j′,n′,φ˜j,nΩij,n,j′,n′,Ai=∇φ˜j′,n′,∇φ˜j,nΩij,n,j′,n′,Ci=wφ˜j′,n′,φ˜j,nΩij,n,j′,n′.and Oi=0N+12×N+12. The element vectors(32)b˜i=b˜i1T,b˜i2T,b˜i3T,b˜i4TT.with the components(33)b˜i1=Vec−f^vR+u0,φ˜j′,n′Ωij′,n′,b˜i2=Vec−f^vI,φ˜j′,n′Ωij′,n′,b˜i3=Vecf^uR+v0,φ˜j′,n′Ωij′,n′,b˜i4=Vecf^uI,φ˜j′,n′Ωij′,n′.Here, the Gauss integral formula is adopted in the quadratures (31) and (33).Assembling the contributions from all elements leads to the global matricesM, A, C and vector b˜, and finally, arrive at the following matrix form(34)smRM−smIM−KA+COsmIMsmRMO−KA+CKA−COsmRM−smIMOKA−CsmIMsmRMU˜=b˜,where the vector U˜ contains the expansion coefficients of all elements. The right-hand side(35)b˜=b˜1T,b˜2T,b˜3T,b˜4TT.where b˜1,2,3,4 could be calculated by (33) defined on the whole domain.In the final of this subsection, we present the following theorem that will be used in the analysis of the existence and uniqueness of the solution of (34).Theorem 1 (see [32]).
The matrixL∈ℝn×n is positive definite if and only if its symmetric part W=L+LT/2 is positive definite. Besides, W is positive definite if and only if its eigenvalues are positive.
Working back from the matrix formulation to the functional form of the integral operator, it is easy to show that matricesA and C are symmetric, and the mass matrix M is symmetric and positive definite [33]. Denote the coefficient matrix of (34) by L; then, we have(36)L+LT2=smRMOOOOsmRMOOOOsmRMOOOOsmRM,
and it is a block diagonal matrix.
It is a practical choice that avoids sampling the integrand atsmR=0. If smR>0, since the mass matrix M is positive definite, we obtain that its eigenvalues are all positive, then the block diagonal matrix L+LT/2 is also positive definite. From Theorem 1, it is direct to get that the coefficient matrix L is positive definite; thus, the existence and uniqueness of the solution of (34) are guaranteed. If smR<0, we consider the coefficient matrix L˜=−L and obtain the same conclusion.
Moreover, since the modal bases (23) are constructed by using the orthogonal Legendre polynomials, the matrices A and M are typically sparse although they may have a full bandwidth [33]. Thus, (34) is large but very sparse. In this paper, we adopt the built-in algorithm in MATLAB to solve it, which is already optimized for solving the system of linear equations.


## 4.2. Parallel LTSEM for the Linear Problem

In Algorithm 1, a description of the pLTSEM for solving LSE is given. The underlying idea of the pLTSEM is straightforward, and the calculation of (34) with different LT parameters sm is independent. This step can be handled more efficiently by the parallel procedure that means the computing procedure can be implemented simultaneously in different CPU cores. Algorithm 1 is divided into three stages of initialization, two-level iteration followed by a solution update step. In the initialization stage, since the matrices M, A, and C depend only on the spatial grid size and the order of polynomial expansion, they can be precomputed before entering into the time loop. The iterative procedure is composed of two loops: the outer loop and the inner loop. The inner loop is a parallel computing process. Its main purpose is to solve the discretized problems on the current time level. The parallel computation is carried out by using the MATLAB parallel computing toolbox, where the parfor loop is a tool for parallel computing single loop. The degree of parallelism is related to the number of CPU cores. The computational time can be saved significantly under the condition of the balance between the computation cost and the data communication. The outer loop is used to forward the solution process to the next time step until the final time is reached.Step 0. for time level n=0 (i.e., t=0),generate the initial mesh, establish matrices M, A, C.Step 1. for time level n≥1,%parallel computing:parfor (m=0; m<M; m++){for each given sm,establish vectors b˜1,2,3,4,solve the global system equation (30) to get coefficient U˜,determineu^Nx,sm, v^Nx,sm.}make numerical inversion of LT to getuNx,tn+1, vNx,tn+1.Step 2. set n=n+1,update the right-hand side (33),goto Step 1end for


## 5. The Increment Linearization Method for the Nonlinear Problem

The LT method proves its virtues by transforming any differential or integral equation into the algebraic equation; however, it has the limitation to deal with the nonlinearity. In this context, we follow the idea of the increment linearization method [26] and incorporate it into pLTSEM to solve the cubic NLSE (2). The numerical results in [18, 26] demonstrate its advantages of easy formulation and relatively little computational cost.First, by the decomposition of real and imaginary parts, (2) can be rewritten as(37)utx,t+∇⋅K∇vx,t+wxvx,t+βvx,tu2x,t+v2x,t=0,vtx,t−∇⋅K∇ux,t−wxux,t−βux,tu2x,t+v2x,t=0.with the initial condition 8 and the homogenous boundary condition (9).At each time step, we decomposeux,t and vx,t into two parts(38)ux,t=ux,tn+px,t,vx,t=vx,tn+qx,t,t∈tn,tn+1,where ux,tn and vx,tn are the solutions at tn=nτ. px,t and qx,t are defined as the increment functions and obviously px,tn=0, qx,tn=0.Substituting (38) into (37) and omitting the nonlinear term, which is Oτγ,γ≥2, we get the following linear system which is satisfied by px,t and qx,t(39)ptx,t+∇⋅K∇qx,t+g11tnxpx,t+g12tnxqx,t+g13tnx=0,qtx,t−∇⋅K∇px,t+g21tnxpx,t+g22tnxqx,t+g23tnx=0.in Ω and tn<t≤tn+1 with the initial condition(40)px,tn=qx,tn=0,inΩ.and the boundary condition(41)px,t=qx,t=0,on∂Ω×tn,tn+1,where gijtnx,i=1,2,j=1,2,3 are given by(42)g11tnx=2βux,tnvx,tn,g12tnx=wx+3βv2x,tn+βu2x,tn,g13tnx=KΔvx,tn+wx+βu2x,tn+βv2x,tnvx,tn,g21tnx=−wx−3βu2x,tn−βv2x,tn,g22tnx=−2βux,tnvx,tn,g23tnx=−KΔux,tn−wx+βu2x,tn+βv2x,tnux,tn.In order not to render this paper unduly long, we avoid presenting in detail the derivation of above equations. The interested readers may consult reference [26]. In particular, we mention that the approximation problem (39) is obtained by neglecting the nonlinear terms of the order Oτγ,γ≥2. Therefore, it is expected that the approximations can reach the second-order accuracy in the temporal direction as reported in [18, 26], which will be verified in the numerical experiments.Second, LT is used for the temporal discretization at each time step. Denote the Laplace transform ofpx,t, qx,t by p^x,s, q^x,s, subject to the initial conditions (40), on the system (39) yields: for any given sm(43)smp^x,sm+∇⋅K∇q^x,sm+g11tnxp^x,sm+g12tnxq^x,sm+g13tnxsm=0,smq^x,sm−∇⋅K∇p^x,sm+g21tnxp^x,sm+g22tnxq^x,sm+g23tnxsm=0.with homogenous Dirichlet boundary conditions.Set the decompositionsm=smR+ismI, p^x,sm=p^Rx+ip^Ix, and q^x,sm=q^Rx+iq^Ix, and (43) becomes(44)smR+g11tnp^R−smIp^I+∇⋅K∇q^R+g12tnq^R=−smRg13tnsm2,smIp^R+smR+g11tnp^I+∇⋅K∇q^I+g12tnq^I=smIg13tnsm2,−∇⋅K∇p^R+g21tnp^R+smR+g22tnq^R−smIq^I=−smRg23tnsm2,−∇⋅K∇p^I+g21tnp^I+smIq^R+smR+g22tnq^I=smIg23tnsm2.with homogenous Dirichlet boundary conditions.The discrete weak form is as follows: for any givensm, find p^NR,p^NI,q^NR,q^NI∈WN0Ω4, such that(45)smR+g11tnp^NR,ϕN−smIp^NI,ϕN−K∇q^NR,∇ϕN+g12tnq^NR,ϕN=−smRsm2g13tn,ϕN,smIp^NR,ϕN+smR+g11tnp^NI,ϕN−K∇q^NI,∇ϕN+g12tnq^NI,ϕN=smIsm2g13tn,ϕN,K∇p^NR,∇ϕN+g21tnp^NR,ϕN+smR+g22tnq^NR,ϕN−smIq^NI,ϕN=−smRsm2g23tn,ϕN,K∇p^NI,∇ϕN+g21tnp^NI,ϕN+smIq^NR,ϕN+smR+g22tnq^NI,ϕN=smIsm2g23tn,ϕN,∀ϕN∈WN0Ω.Assume the approximations on any elementΩi as follows:(46)p^NRx=∑j=0N∑n=0Np˜j,nRφ˜j,nx,p^NIx=∑j=0N∑n=0Np˜j,nIφ˜j,nx,q^NRx=∑j=0N∑n=0Nq˜j,nRφ˜j,nx,q^NIx=∑j=0N∑n=0Nq˜j,nIφ˜j,nx.with the expansion coefficients p˜j,nRj,n=0N, p˜j,nIj,n=0N, q˜j,nRj,n=0N, and q˜j,nIj,n=0N.The local linear system arising from (45) has the following partitioned matrix form(47)smRMi+Dip−smIMi−KAi+EiqOismIMismRMi+DipOi−KAi+EiqKAi+EipOismRMi+Diq−smIMiOiKAi+EipsmIMismRMi+DiqU˜i=b˜i,where the element vector U˜i is composed of the expansion coefficients(48)U˜i=Vecp˜j,nRj,n=0NT,Vecp˜j,nIj,n=0NT,Vecq˜j,nRj,n=0NT,Vecq˜j,nIj,n=0NTT.The local element matricesDip, Eip, Diq, and Eiq are computed by(49)Dip=g11tnφ˜j′,n′,φ˜j,nΩij,n,j′,n′,Eip=g21tnφ˜j′,n′,φ˜j,nΩij,n,j′,n′,Diq=g22tnφ˜j′,n′,φ˜j,nΩij,n,j′,n′,Eiq=g12tnφ˜j′,n′,φ˜j,nΩij,n,j′,n′.and the right-hand side(50)b˜i=−smRsm2b˜i1T,smIsm2b˜i1T,−smRsm2b˜i2T,smIsm2b˜i2TT.with the components(51)b˜i1=Vecg13tn,φ˜j′,n′Ωij′,n′,b˜i2=Vecg23tn,φ˜j′,n′Ωij′,n′.Assemble the contributions from all elements, and we obtain the following global linear system with the partitioned matrix form(52)smRM+Dp−smIM−KA+EqOsmIMsmRM+DpO−KA+EqKA+EpOsmRM+Dq−smIMOKA+EpsmIMsmRM+DqU˜=b˜,where the vector U˜ is composed of the expansion coefficients of all elements. The coefficient matrix is sparse as can be appreciated from the fact that M is sparse and O is the zero matrix. The right-hand side(53)b˜=−smRsm2b˜1T,smIsm2b˜1T,−smRsm2b˜2T,smIsm2b˜2TT,where b˜1 and b˜2 could be calculated by (51) defined on the whole domain.Algorithm 2 describes the process of the pLTSEM for solving the cubic NLSE. The procedure is similar to Algorithm 1, which also has two loops. The outer loop is to forward the time step to the final level, and the inner loop is conducted by the parallel computing for each time step. The difference between two algorithms is that the main purpose of the inner loop in Algorithm 2 is to solve the linear system which is satisfied by the increment functionspx,t, qx,t rather than ux,t, vx,t. After that, the computation on the current time level is complete by the addition operation (38) to obtain the approximations uNx,tn+1 and vNx,tn+1. Notice that matrices Dp, Ep, Dq, Eq and vectors b˜1, b˜2 are independent of sm, and they only need to be computed once in the outer loop.Step 0. for time level n=0 (i.e., t=0),generate the initial mesh, establish matrices M, A.Step 1. for time level n≥1,computeg11tnx, g12tnx, g13tnx, g21tnx, g22tnx and g23tnx,establish matrices Dp, Ep, Dq, Eq and vectors b˜1, b˜2.%parallel computing:parfor (m=0; m<M; m++){for each given sm,solve the global system equation (44) to get coefficient U˜,determinep^Nx,sm, q^Nx,sm.}make numerical inversion of LT to getpNx,tn+1, qNx,tn+1,computeuNx,tn+1 and vNx,tn+1 by using (38).Step 2. set n=n+1, goto Step 1endfor


## 6. Numerical Examples

In this section, numerical examples are presented to witness the performance of our proposed method in terms of the accuracy and the parallel efficiency. In the first linear numerical example, the convergence rate of LTSEM is investigated. The nonlinear numerical example demonstrates the validity and feasibility of applying the increment linearization method to the cubic NLSE. Numerical results show the pLTSEM provides much more efficiency than LTSEM, and the pLTSEM with 4 cores is the best. It confirms the superiority of parallel computation of the LT method. All the calculations are conducted on the computer of Intel(R) Core(TM)i5-8500T CPU 2.10 GHz, 16.0 GB, Windows 10(x64).Example 1.
In this test case, we study the LSE (1) defined on 0,82, where the parameters adopted are K=1/2, wx=−1, and fx,t which are chosen such that the exact solution is(54)Φx,t=e3t8−x8−ysinxy+i8−xytsinx8−y.
The initial condition could be obtained easily from setting the exact solution att=0.
First, we investigate the order of convergence of LTSEM in the spatial direction and the temporal direction. For the spatial direction, the approximation can converge either algebraically by increasing the number of elements with fixed degree of the polynomials or exponentially by increasing the degree of the polynomials with fixed number of elements [30]. Here, we choose the latter path and employ a uniform quadrilateral element mesh with the size of 2. The discrete L2 norm error etL2 and the discrete H1 norm error etH1 of the modulus Φt at t=1 are computed. In Figure 1(a), we plot the errors with various N (the degree of the polynomial) and set M=80 (the number of integral points) at each time step to ignore the error in the temporal direction. For this linear problem, we choose the time step τ=1. One can see that the convergence rate in the spatial direction is Oe−1.6N for both two norms. Table 1 lists the errors with various M, where we set N=26 to ignore the error in the spatial direction. The convergence rate is shown in Figure 1(b), and it is found that the order is Oe−0.75M for both two norms.
Second, we show the excellent performance of the pLTSEM toward LTSEM. Table1 presents the comparison of CPU time and errors of LTSEM and pLTSEM, where Tsm denotes the CPU time required to get the approximations for different values of sm, TILT denotes the CPU time required for the numerical inversion of LT, and Tall denotes the total CPU time required for all the operations. It is clearly observed that for LTSEM, the numerical inversion of LT part accounts a small amount of calculation and the solving system (30) part is much more time consuming. Besides, the errors for LTSEM and pLTSEM are always the same, independently of the number of processors that are used. Hence, the parallel algorithm for the solving of equation (30) with different values of sm does improve the computational efficiency.
For the sake of comparison, CNSEM is also considered, and it is using the CN scheme for the time derivative and SEM for the space derivative. The CN scheme is a kind of time-marching scheme with second-order convergence rate. The numerical results of CNSEM with various time stepsτ are listed in Table 2.
For convenience, we also plot errors versus the total CPU time in Figure2, where x-axis represents the total CPU time Tall and y-axis represents the errors. It is easy to find that pLTSEM works much better than LTSEM and CNSEM; that is, when the comparable magnitude of the error is obtained, pLTSEM requires much less CPU time and CNSEM is the most expensive. More precisely, from the last three columns in Table 1, we could conclude that pLTSEM with 2 cores reduces nearly 38 % of the total CPU time of LTSEM, and pLTSEM with 4 cores reduces nearly 54 %.
In addition, Table3 provides the errors and the total CPU time at different times with M=80 and N=26. The time interval 0,5 is discretized by the time step τ=1 and LT is used at each time step.Figure 1
Numerical results of LTSEM att=1: (a) errors with various N and M=80; (b) errors with various M and N=26, where N is the degree of polynomial and M is the number of integral points at each time step, Example 1.
(a)(b)Table 1
Results of LTSEM and pLTSEM with variousM and N=26 at t=1, Example 1.
LTSEM/pLTSEM (2/4 cores)LTSEMpLTSEM (2 cores)pLTSEM (4 cores)Me1L2e1H1Tsm (s)TILT (s)Tall (s)Tall (s)Tall (s)504.10e−011.66e+0034.141.5937.7323.7117.30521.00e−014.10e−0136.131.5839.7125.0218.25542.43e−029.89e−0238.171.5941.7626.3019.20565.74e−032.33e−0239.991.5843.5727.4119.99581.33e−035.41e−0341.931.5945.4328.5520.89603.03e−041.23e−0344.011.5947.5029.5721.82626.77e−052.75e−0445.951.6049.5431.1622.77641.49e−056.06e−0547.991.5951.4932.4123.66663.23e−061.31e−0550.041.5953.5433.7024.61686.89e−072.80e−0651.991.5855.5734.9225.55701.45e−075.90e−0753.901.5957.4736.1826.43723.00e−081.22e−0755.851.5959.4337.4127.32746.18e−092.52e−0857.991.6061.5738.7728.26761.23e−095.25e−0959.941.5963.5239.9329.18782.38e−104.19e−0962.071.5865.6541.2930.20804.18e−113.63e−0964.131.5967.7042.6431.12Table 2
Results of CNSEM with variousτ and N=26 at t=1, Example 1.
τe1L2ordere1H1orderTall (s)1/243.40e+00—9.48e+00—9.941/258.51e−011.992.37e+001.9919.891/262.12e−011.995.94e−011.9939.701/275.32e−021.991.48e−011.9979.471/281.33e−021.993.71e−022.00156.991/293.32e−032.009.29e−032.00310.471/2108.31e−042.002.32e−032.00625.831/2112.07e−042.005.80e−042.001269.931/2125.19e−052.001.45e−041.992529.301/2131.30e−051.994.74e−051.615078.361/2143.25e−061.995.81e−05—10060.281/2159.03e−071.841.43e−04—20294.981/2166.96e−070.372.59e−04—39789.691/2171.30e−06—5.23e−04—80595.92Figure 2
Errors att=1 versus the total CPU time Tall: (a) LTSEM and pLTSEM with M = 50, 52,…, 80 and N = 26 ; (b) CNSEM with τ=1/24,1/25,…,1/217 and N = 26 , Example 1.
(a)(b)Table 3
Results of LTSEM and pLTSEM at different times withM=80, τ=1, and N = 26, Example 1.
tLTSEM/pLTSEM (2/4 cores)LTSEMpLTSEM (2 cores)pLTSEM (4 cores)etL2etH1Tall (s)Tall (s)Tall (s)14.18e−113.63e−0967.7042.6431.1227.37e−105.68e−08136.6987.0963.0231.77e−081.82e−06206.68131.9595.2642.67e−072.02e−05275.58175.41126.3857.16e−067.22e−04343.49218.89158.25Example 2.
[9] Consider the cubic NLSE (2) defined on Ω=0,π2, where the parameters adopted are K=1, β=1, wx=−sin2xsin2y, and the exact solution is given by(55)Φx,t=e−2itsinxsiny.
The discreteL2 and H1 norm errors at t=1 with various time steps τ are reported in Table 4, where the quadrilateral element mesh with the size of π/8 is employed and M=16, N=6. These results show that the time discretization scheme (43) is second-order accurate for this two-dimensional nonlinear problem. This is consistent with the numerical results reported in [18, 26] and also verifies that the combination of the increment linearization method and LTSEM/pLTSEM is a feasible approach for solving the cubic NLSE.
In order to illustrate the superiority of the parallel LT, we consider the CN scheme in the temporal direction. Using the central difference approximation forut and vt in (37), the temporal discretized form becomes(56)un+1−unτ+12∇⋅K∇vn+1+∇⋅K∇vn+12wvn+1+vn,+β2vn+1un+12+vn+12+vnun2+vn2=0,vn+1−vnτ−12∇⋅K∇un+1+∇⋅K∇un−12wun+1+un,−β2un+1un+12+vn+12+unun2+vn2=0,where un and vn evaluate ux,tn and vx,tn, respectively. Obviously, this scheme is nonlinearly implicit and coupled in the computation, so it requires heavy iterative calculations to solve the resulting implicit algebraic equations at each time step and is not suitable for the parallel computation. Newton’s method of quasilinearization [34] is used to deal with the nonlinearities, and the scheme becomes(57)un+1,k+1−unτ+12∇⋅K∇vn+1,k+1+∇⋅K∇vn+12wvn+1,k+1+vn+β2g1n+1,k+g2n+1,kun+1,k+1−un+1,k+g3n+1,kvn+1,k+1−vn+1,k+β2vnun2+vn2=0,vn+1,k+1−vnτ−12∇⋅K∇un+1,k+1+∇⋅K∇un−12wun+1,k+1+un−β2g4n+1,k+g5n+1,kun+1,k+1−un+1,k+g6n+1,kvn+1,k+1−vn+1,k−β2unun2+vn2=0.with(58)g1n+1,k=vn+1,kun+1,k2+vn+1,k2,g2n+1,k=2un+1,kvn+1,k,g3n+1,k=un+1,k2+3vn+1,k2,g4n+1,k=un+1,kun+1,k2+vn+1,k2,g5n+1,k=3un+1,k2+vn+1,k2,g6n+1,k=2un+1,kvn+1,k,where k denotes the iteration number and un+1,0=un, vn+1,0=vn. At each time level, we iterate (57) until the residual condition is satisfied:(59)maxun+1,k+1−un+1,k∞,vn+1,k+1−vn+1,k∞<10−8.
After this, we can go to the next time level, and the process is iterated until reaching to the desirable time. The numerical results including the errors att=1 and the total CPU time with various τ are shown in Table 5.
For the convenience of comparison, we use the data in Tables4 and 5 to plot Figure 3, which depicts the errors versus the total CPU time of LTSEM, pLTSEM, and CNSEM. From Figure 3, it is definite that for this two-dimensional cubic NLSE, CNSEM is the most expensive. Moreover, from Table 4, we can find that pLTSEM presents an obvious advantage on computational efficiency over LTSEM. Without loss of accuracy, pLTSEM with 2 cores reduces nearly 18 % of the total CPU time of LTSEM, and pLTSEM with 4 cores reduces nearly 36 %. Besides, Table 6 lists the errors, the mass and energy differences, and the total CPU time of LTSEM and pLTSEM at several final times with τ=1/26, M=16, and N=6.Table 4
Results of LTSEM and pLTSEM with variousτ and M=16, N=6 at t=1, Example 2.
τLTSEM/pLTSEM (2/4 cores)LTSEMpLTSEM (2 cores)pLTSEM (4 cores)e1L2Ordere1H1OrderTall (s)Tall (s)Tall (s)1/22.63e−01—4.73e−01—18.6414.9311.291/226.44e−022.021.13e−012.0637.4230.7223.171/231.72e−021.903.01e−021.9174.7761.0847.581/244.51e−031.937.88e−031.93149.80122.9295.341/251.15e−031.972.01e−031.97301.37247.96190.271/262.91e−041.985.09e−041.98603.43495.35379.401/277.21e−052.011.25e−042.021208.171001.66761.611/281.70e−052.082.98e−052.072426.802022.831525.60Table 5
Results of CNSEM with variousτ and N=6 at t = 1 , Example 2.
τe1L2ordere1H1orderTall (s)1/22.74e−01—4.74e−01—105.161/227.58e−021.851.31e−011.85167.571/231.94e−021.963.37e−021.95330.041/244.90e−031.988.49e−031.98500.051/251.22e−032.002.12e−032.00992.901/263.07e−041.995.32e−041.991980.631/277.67e−052.001.33e−042.004030.341/281.92e−051.993.32e−052.007896.50Figure 3
Errors att=1 versus the total CPU time Tall of LTSEM, pLTSEM, and CNSEM with τ=1/2,1/22,…,1/28, Example 2. (a) e1L2. (b) e1H2.
(a)(b)Table 6
Results of LTSEM and pLTSEM at different times withτ=1/26, M=16, and N = 6, Example 2.
tLTSEM/pLTSEM (2/4 cores)LTSEMpLTSEM (2 cores)pLTSEM (4 cores)etL2etH1Mt−M0Et−E0Tall (s)Tall (s)Tall (s)12.91e−045.09e−042.54e−054.70e−05603.43495.35379.4025.75e−041.00e−035.10e−051.01e−041207.73981.90754.5838.49e−041.47e−037.71e−051.39e−041781.711474.361138.2341.11e−031.92e−031.02e−041.56e−042379.171958.601521.5451.36e−032.42e−031.27e−041.81e−042991.632450.681900.3864.51e−037.88e−031.36e−032.07e−033601.522946.132281.71108.55e−031.48e−022.73e−034.30e−036084.164901.333790.74151.58e−022.75e−022.05e−022.90e−029132.547346.745686.56201.07e−021.90e−022.74e−023.88e−0212162.669803.187585.48Example 3.
Consider the cubic NLSE (2) with a rotational term as follows [35](60)i∂Φx,t∂t+12∂2Φx,t∂x2+12∂2Φx,t∂y2−Φx,t2Φx,t+wxΦx,t+iΘy∂Φx,t∂x−x∂Φx,t∂y=0.which describes the dynamic properties of one-component Bose–Einstein condensates under external rotation Θ>0. We use the following initial data(61)Φ0x=1πe−x2+y2/2.and take parameters wx=1.5x2+0.5y2,Θ=0.7 on the bounded domain Ω=−8,82. The problem is solved using pLTSEM with 4 cores. The quadrilateral element mesh with the size of 1 is employed, and the parameters are taken as N=16, τ=0.01, and M=16. Figure 4 shows the numerical solutions ReΦ, ImΦ, and Φ at t=1.Figure 4
Surface plot ofReΦ, ImΦ, and Φ at t=1, Example 3. (a) Re(Φ(x,1)). (b) Im(Φ(x,1)). (c) |Φ(x,1)|.
(a)(b)(c)


## 7. Conclusion

This paper presents an efficient parallel Laplace transform spectral element method (pLTSEM) that greatly outperforms the Crank–Nicolson scheme in the temporal direction for solving the linear Schro¨dinger equation and the cubic nonlinear Schro¨dinger equation. The basic idea is to use the Laplace transform method to eliminate the time dependence and apply spectral element method to discretize the subsequent BVPs with different Laplace transform parameters. Since the solving procedures of the discretized systems are independent of each other, they can be handled in the parallel computation for efficiency reason. The parallel computing could be directly realized with the MATLAB parallel computing toolbox. For the cubic nonlinear Schro¨dinger equation, pLTSEM is also robust when used in conjunction with the increment linearization method. Finally, the approximation is inverted from the transformed domain to the time domain by the numerical inversion of Laplace transform using Talbot’s method. The numerical results indicate that our proposed method is competitive, and good accuracy and speed-up are obtained for both linear and nonlinear problems.Here, we restrict our attention to a class of Schro¨dinger equations. This is the first step of applying pLTSEM to solve the time-dependent problems. Actually, the numerical frameworks presented in this work are essential for extension to more complicated cases, for example, the coupled systems of nonlinear Schro¨dinger equations. Discrete schemes can be established following an analogous analysis to this paper with slight technical modifications.



---

*Source: 4806846-2022-07-21.xml*