A combination of multiscale time integrator and two-scale formulation for the nonlinear Schrödinger equation with wave operator

Accepted Manuscript A combination of multiscale time integrator and two-scale formulation for the nonlinear Schrödinger equation with wave operator Xiaofei Zhao

PII: DOI: Reference:

S0377-0427(17)30306-0 http://dx.doi.org/10.1016/j.cam.2017.06.006 CAM 11184

To appear in:

Journal of Computational and Applied Mathematics

Received date : 21 December 2016 Revised date : 27 April 2017 Please cite this article as: X. Zhao, A combination of multiscale time integrator and two-scale formulation for the nonlinear Schrödinger equation with wave operator, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.06.006 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Manuscript Click here to view linked References

A combination of multiscale time integrator and two-scale formulation for the nonlinear Schr¨odinger equation with wave operator Xiaofei Zhaoa,b a

INRIA-Rennes Bretagne Atlantique, IPSO Project, France b IRMAR, Universit´ e de Rennes 1, Rennes, France

Abstract In this paper, we consider the nonlinear Schr¨ odinger equation with wave operator (NLSW), which contains a dimensionless parameter 0 < ε ≤ 1. As 0 < ε  1, the solution of the NLSW propagates fast waves in time with wavelength O(ε2 ) and the problem becomes highly oscillatory in time. The oscillations come from two parts. One part is from the equation and another part is from the initial data. For the ill-prepared initial data case as described in [2] which brings inconsistency in the limit regime, standard numerical methods have strong convergence order reduction in time when ε becomes small. We review two existing methods to solve the NLSW: an exponential integrator and a two-scale method. We comment on their order reduction issues. Then we derive a multiscale decomposition two-scale method for solving the NLSW by first performing a multiscale decomposition on the NLSW which decomposes it into a well-behaved part and an energy-unbounded part, and then applying an exponential integrator for the well-behaved part and a twoscale approach for the energy-unbounded part. Numerical experiments are conducted to test the proposed method which shows uniform second order accuracy without significant order reduction for all 0 < ε ≤ 1. Comparisons are made with the existing methods. Keywords: nonlinear Schr¨odinger equation with wave operator, highly oscillatory, ill-prepared initial data, multiscale decomposition, two-scale formulation, uniformly accurate, order reduction

1. Introduction In this paper, we consider the following nonlinear Schr¨ odinger equation with wave operator (NLSW) in d-dimensions (d = 1, 2, 3) [1, 2, 29, 19]: (
2i∂t uε (x, t) − ε2 ∂tt uε (x, t) + ∆uε (x, t) + f |uε (x, t)|2 uε (x, t) = 0, x ∈ Rd , t > 0, (1.1) uε (x, 0) = ϕ1 (x), ∂t uε (x, 0) = ϕ2 (x), x ∈ Rd . Here t is time, x is the spatial coordinate, uε := uε (x, t) is the unknown complex-valued scalar field, 0 < ε ≤ 1 is a parameter, ϕ1 and ϕ2 are two given initial functions which are uniformly bounded for all 0 < ε ≤ 1, and f (·) : R → R is a given nonlinearity independent of ε. The above NLSW comes widely from both mathematical analysis and physics applications. On the mathematical side, in the theoretical analysis of the nonrelativistic limit of the nonlinear Klein-Gordon equation, the NLSW (1.1) has been considered as an intermediate limit model which has bounded energy as ε → 0 [24, 25, 23, 28]. On the physics applications side, (1.1) has been derived as the model for the envelope approximation of Langmuir turbulence [5, 13] or the modulated planar pulse approximation of the sine-Gordon equation for light bullets [7, 30]. In these cases, the ε represents a dimensionless parameter Email address: [email protected] (Xiaofei Zhao)

Preprint submitted to Elsevier

April 27, 2017

scaled to be inversely proportional to certain physical quantity such as the speed of light or plasma frequency. The NLSW (1.1) is known to conserve the mass Z   ε (1.2) M (t) := |u (x, t)|2 dx − ε2 Im (uε (x, t)∂t uε (x, t)) dx ≡ M (0), t ≥ 0, Rd

and the energy

E(t) :=

Z

Rd

 2  ε |∂t uε (x, t)|2 + |∇uε (x, t)|2 − F (|uε (x, t)|2 ) dx ≡ E(0), t ≥ 0,

(1.3)

Rρ with F (ρ) = 0 f (ρ0 )dρ0 . In the limit regime of the NLSW (1.1), i.e. when ε → 0, the equation has been shown to converge in some energy space to the nonlinear Schr¨ odinger equation [24, 25, 23, 5]: (
2i∂t u(x, t) + ∆u(x, t) + f |u(x, t)|2 u(x, t) = 0, x ∈ Rd , t > 0, (1.4) u(x, 0) = ϕ1 (x), x ∈ Rd . [31] studied the vortex dynamics in (1.1) with a wide range of ε. As indicated from the results, the NLSW (1.1) bridges a phase transition of dynamics between wave equation and pure dispersive equation. As 0 < ε  1, the equation (1.1) propagates waves with wavelength O(ε2 ) in time and O(1) in space of the solution uε (x, t), respectively. For the two given initial data ϕ1 and ϕ2 in (1.1), one can denote that ϕ2 (x) −

 i ∆ϕ1 (x) + f (|ϕ1 (x)|2 )ϕ1 (x) = εp r(x), 2

p ≥ 0,

(1.5)

for some function r(x) = O(1). When p ≥ 2, the pair of initial data ϕ1 and ϕ2 is called as the ‘well-prepared’ initial data, in which case the leading oscillating term in the solution uε is due to the term ε2 ∂tt uε in the equation (1.1) and one can establish the uniform boundedness of ∂tt uε (x, t) as ε → 0 for some t ≥ 0 [23]. On the contrary, when 0 ≤ p < 2, the initial data ϕ1 and ϕ2 due to the inconsistency in the limit regime will introduce extra oscillations apart from the equation itself to the solution uε (x, t) in time, which makes ∂tt uε (x, t) unbounded as ε → 0. In this sense, this type of initial data is called as ‘ill-prepared’ in [2, 1]. The oscillation in the solution in the ill-prepared case is stronger than that in the well-prepared case and that, as a matter of fact, will essentially affect the temporal accuracy of standard numerical methods. On the numerical discretisation aspect of the NLSW (1.1), the work in the literature [19, 29, 13] mainly considered some conservative finite difference time domain methods that could conserve the mass (1.2) and energy (1.3) in the discrete level. [1] investigated the error bounds of those finite difference discretisation for the NLSW ∗ as 0 < ε  1. The results show that the finite difference integrator satisfies two error bounds O(∆t2 /ε4−p ) where p∗ = min{p, 2} and O(∆t2 ) + O(ε2 ). The combination of the two error bounds gives the so-called ∗ uniform convergence of order O(∆t4/(6−p ) ) for all ε ∈ (0, 1] which implies that the finite integrator can at most reach a linear uniform convergence rate for the NLSW (1.1) even with the well-prepared initial data. Later, [2] applied and analysed an exponential integrator to solve the NLSW (1.1), where the uniform ∗ accuracy order has been raised to O(∆t4/(4−p ) ). Hence, in the well-prepared case the exponential integrator manages to push the uniform convergence order to the optimal second order rate, and this fact has found important application in the design of the multiscale time integrator for the nonlinear Klein-Gordon equation [4] and related models [10, 11]. However, the convergence rate of the exponential integrator could drop down to linear in the worst ill-prepared case, i.e. p = 0. Thus, all the classical numerical methods could not reach the optimal uniform second order rate or in other words that the methods have significant order reduction in the integration of the NLSW (1.1) with the ill-prepared initial data as ε → 0. Recently, a two-scale formulation approach has been considered for a class of highly oscillatory equations [12] where a uniform second order error bound has been established rigorously. While on the application to the Klein-Gordon equation in the nonrelativistic limit, unexpected convergence order reduction in the two-scale method has also been observed when the time step is not small enough in practical computing. This work is devoted to concerning the NLSW (1.1) under the ill-prepared initial data case, in particular the worse case  i ϕ2 (x) − ∆ϕ1 (x) + f (|ϕ1 (x)|2 )ϕ1 (x) = r(x), (1.6) 2 2

i.e. (1.5) with p = 0 and some non-zero function r(x) = O(1). We are going to design a second order uniformly accurate numerical scheme without significant order reduction for all ε ∈ (0, 1]. We shall first review the exponential integrator in [2] and the two-scale method in [12] for solving NLSW (1.1), where we remark on the order reduction issue of the two-scale method and give some explanations. Next, by using the spirit of the multiscale time integrator [4, 8, 10, 11], which is to decompose the equation into different-behaving parts and then integrate them with different suitable integrators, we perform a multiscale decomposition on the NLSW by amplitudes which decomposes (1.1) into a well-behaved part and an energyunbounded part. A combined multiscale integrator two-scale method will then be derived by solving the well-behaved part of the NLSW with an exponential integrator and solving the energy-unbounded part with the revised two-scale method. Numerical experiments are done to show the uniform second order accuracy of the proposed method and comparisons are made with the existing direct exponential integrator or two-scale method on (1.1). The rest part of the paper is organised as follows. In section 2, we give a brief review of the exponential integrator and the two-scale method. In section 3, we introduce the combined multiscale integrator two-scale method. Numerical results and discussions are given in section 4. Some conclusions are drawn in section 5. Throughout this paper, we adopt the notation A . B to represent that there exists a generic constant C > 0, which is independent of ∆t (or n), ∆x, ∆τ and ε, such that |A| ≤ CB. 2. Review of existing results In this section, for readers’ convenience and as benchmark for comparisons later, we shall briefly review two recently proposed numerical methods for solving the NLSW (1.1). One is an exponential integrator approach, and the other is the two-scale formulation approach. 2.1. Exponential integrator We start with a brief review of the exponential integrator (or exponential wave integrator) Fourier pseudospectral method (EIFP) proposed in [2] and its convergence results. To present the scheme, we consider the one space dimension case of the NLSW for simplicity, i.e. d = 1 in (1.1). As widely used in the literature [1, 2, 29], the whole space problem is truncated onto a bounded interval Ω = (−L, L) and imposed with periodic boundary condition due to the fact that the solution uε decays fast at far field:
ε  ε 2 ε ε ε 2   2i∂t u (x, t) − ε ∂tt u (x, t) + ∂xx u (x, t) + f |u (x, t)| u (x, t) = 0, x ∈ Ω, t > 0, (2.1) uε (x, 0) = ϕ1 (x), ∂t uε (x, 0) = ϕ2 (x), x ∈ Ω,   ε ε ε ε u (−L, t) = u (L, t), ∂x u (−L, t) = ∂x u (L, t), t ≥ 0.

To avoid significant aliasing error, in practise L > 0 is choosing large enough. We take ∆t > 0 as the time step and denote tn = n∆t for n = 0, 1, . . .. To present the spectral discretisation in space, we introduce the following notations. Choose N > 0 an even integer and define   πl N N XN := span eiµl (x+L) : x ∈ Ω, µl = , l = − , . . . , −1 . L 2 2 For a periodic function v(x) on Ω, let PN : L2 (Ω) → XN be the standard L2 -projection operator [20, 26], N/2−1

(PN v)(x) =

X

l=−N/2

vbl e

iµl (x−a)

1 vbl = b−a

,

Z

b

v(x)e−iµl (x−a) dx.

a

Applying Fourier transformation on both sides of the equation (2.1) and denoting N/2−1 ε

PN u (x, t) =

X

l=−N/2

N/2−1

u bl (t)e

iµl (x+L)

,

ε 2

ε

PN f (|u | )u (x, t) = 3

X

l=−N/2

fbl (t)eiµl (x+L) ,

(2.2)

one has for l = −N/2, . . . , N/2 − 1,

2ib u0l (t) − ε2 u b00l (t) − µ2l u bl (t) + fbl (t) = 0,

t > 0.

Then the variation-of-constant formula or the Duhamel principle shows for some n ≥ 0, Z s 2 0 u bl (tn + s) = al (s)b ul (tn ) + ε bl (s)b ul (tn ) + bl (s − θ)fbl (tn + θ)dθ, s ≥ 0,

(2.3)

0

for l = −N/2, . . . , N/2 − 1, where −

+

isλl λ+ eisλl − λ− l e al (s) := l , + − λl − λl

+

−

eisλl − eisλl bl (s) := i 2 − , ε (λl − λ+ l )

λl±

=

1±

p

1 + ε2 µ2l . ε2

(2.4)

The EIFP as proposed in [2] then takes s = ±∆t in (2.3) and cancels out the derivative term ψbl0 (tn ) from the equation as: Z ∆t u bl (tn+1 ) =νl u bl (tn−1 ) + (al (∆t) − νl al (∆t))b bl (∆t − θ)fbl (tn + θ)dθ ul (tn ) + − νl

Z

0

∆t

0

bl (∆t − θ)fbl (tn − θ)dθ,

n ≥ 1,

l=−

N N ,..., − 1, 2 2

(2.5)

where νl = bl (∆t)b−1 l (−∆t), and hereafter c denotes the complex conjugate of a complex number c. With approximation: i θ hb fbl (tn ± θ) ≈ fbl (tn ) ± fl (tn ) − fbl (tn−1 ) , n ≥ 1, (2.6) ∆t the two unknown trigonometric integrations in (2.5) can be evaluated explicitly, i.e. one can find exactly Z ∆t Z ∆t cl (∆t) = bl (∆t − θ)dθ, dl (∆t) = θbl (∆t − θ)dθ. 0

0

With a starting value at t1 approximated similarly based on (2.3), then one ends up with an exponential integrator spectral method as presented in [2] that proceeds with information of u on three time-levels: ul (tn ) + (cl (∆t) − νl cl (∆t))fbl (tn ) u bl (tn+1 ) ≈νl u bl (tn−1 ) + (al (∆t) − νl al (∆t))b
fbl (tn ) − fbl (tn−1 ) N N , n ≥ 1, l = − , . . . , − 1, + dl (∆t) + νl dl (∆t) ∆t 2 2 u bl (t1 ) ≈al (∆t)b ul (0) + ε2 bl (∆t)b u0l (0) + cl (∆t)fbl (0) + dl (∆t)fbl0 (0).

(2.7)

Under condition (1.5), by utilising the error analysis techniques from asymptotic preserving schemes [27], i.e. analysing the method on the problem itself and analysing the method on the limit model, the following two finite time error bounds have been established under some Sobolev norm for the above EIFP method as  
∆t2 O(∆xm0 ) + O and O(∆xm0 ) + O(∆t2 ) + O ε2 , (2.8) ∗ 2−p ε

where p∗ = min{p, 2} with p ≥ 0 given in (1.5) and m0 > 0 is constant which dependents on the regularity of the solution uε of the NLSW (1.1). As a consequence, a uniform error bound can be archived by taking the minimum of the two independent error bounds as   ∗ O(∆xm0 ) + O ∆t4/(4−p ) . This error estimate of the exponential integrator scheme (2.7) has been shown to be sharp [2]. Thus, in spatial approximation, the method can have uniform spectral accuracy when the solution is smooth. For the temporal approximation, the method can reach the optimal uniform second order rate in the well-prepared case. However, for the ill-prepared initial data case, where 0 ≤ p < 2, the scheme does not have the uniform second order temporal accuracy, or as people often say that the method has order reduction. For the worst case, i.e. p = 0, the method only has first order temporal accuracy for all ε ∈ (0, 1]. 4

2.2. Two-scale formulation approach Next, we review the two-scale method that has been proposed and used to solve the Klein-Gordon (KG) equation in [12]. It is known that the KG and NLSW are connected by a gauge transformation. So in order to apply the two-scale method in [12] for KG to NLSW, we first introduce the gauge transformation on (1.1) by defining: 2 v ε (x, t) := e−it/ε uε (x, t), x ∈ Rd , t ≥ 0. (2.9) Then the NLSW (1.1) can be transformed into the nonlinear KG equation 
1   ε2 ∂tt v ε (x, t) − ∆v ε (x, t) + 2 v ε (x, t) − f |v ε (x, t)|2 v ε (x, t) = 0, ε  ε ε  v (x, 0) = ϕ1 (x), ∂t v (x, 0) = ϕ2 (x) − i ϕ1 (x), x ∈ Rd . ε2

x ∈ Rd , t > 0,

(2.10)

We remark that the gauge transformation is often performed in the literatures [24, 25, 23] to transfer the KG into NLSW in order to eliminate the unbounded term, i.e. ε12 v ε in (2.10), while here we use it in the reverse way. Although (2.10) is mathematically equivalent to (1.1), one can see that the behavior of solution v ε is much worse than uε in terms of the oscillation by noting as ε → 0,
∂t v ε (x, t) = O ε−2 , ∂t uε (x, t) = O(1), x ∈ Rd , t ≥ 0. As derived in [12], introducing the following series of transformations:

z+ (x, t) = v ε (x, t) − iε2 (1 − ε2 ∆)−1/2 ∂t v ε (x, t), z− (x, t) = v ε (x, t) − iε2 (1 − ε2 ∆)−1/2 ∂t v ε (x, t),
  1   √ f 2 (z+ + z− ) z+ −it/ε2 1−ε2 ∆ e
z(x, t), (2.11) f (z) = , z = , w(x, t) = e z− f 12 (z− + z+ )

the KG (2.10) can be written into a first order system   √  √  ∂t w(x, t) = i(1 − ε2 ∆)−1/2 e−it 1−ε2 ∆/ε2 fe eit 1−ε2 ∆/ε2 w(x, t) ,  w(x, 0) = ϕ (x) − i(1 − ε2 ∆)−1/2 ε2 ϕ (x) − iϕ (x)
=: w (x), 1

2

1

x ∈ Rd , t > 0,

(2.12)

0

√ √ where 1 − ε2 ∆ is interpreted as the pseudo-differential operator. In the sense of operator, 1 − ε2 ∆/ε2 can be decomposed into a constant (Fourier frequency) leading part and a bounded (as ε → 0) part: √ √ 1 − ε2 ∆ 1 − ε2 ∆ − 1 1 1 = + D , 0 ≤ D := ≤ − ∆, ε ε 2 2 2 ε ε ε 2 which indicates (2.12) is suitable for the separation of scales. Then the augmented problem is introduced as   ∂ W (x, t, τ ) + 1 ∂ W (x, t, τ ) = F(t, τ, W ), t > 0, x ∈ Rd , τ ∈ T, (2.13a) t τ ε2  W (x, 0, τ ) = W0 (x, τ ), x ∈ Rd , τ ∈ T, (2.13b)

where T := R/2π is the torus and

F(t, τ, W ) := i(1 − ε2 ∆)−1/2 e−iτ e−itDε fe(eiτ eitDε W ).

In the two-scale formulation, τ is interpreted as an extra space variable and the unknown function W (·, ·, τ ) satisfies periodic boundary condition. The initial data W0 (x, τ ) of (2.13) is only specified as W0 (x, 0) = w0 (x) such that W (x, t, t/ε2 ) = w(t). Hence, one has freedom to choose the initial data. With the help of the following operators defined for some periodic function ϕ(τ ) on T as Z 2π Z τ 1 Lϕ(τ ) = ∂τ ϕ(τ ), Πϕ = ϕ(τ )dτ ; L−1 ϕ = (I − Π) ϕ(θ)dθ, if Πϕ = 0; A := L−1 (I − Π), 2π 0 0 5

the following two pairs of initial data for (2.13) are constructed by the Chapman-Enskog expansion [12]: W01st (x, τ ) := v0 (x) + ν1 (τ, v0 (x)) − ν1 (0, v0 (x)),

W02nd (x, τ ) := v0 (x) + ν1 τ, W01st (x, τ ) − ν1 0, W01st (x, τ ) + ν2 (τ, v0 (x)) − ν2 (0, v0 (x)),

(2.14a) (2.14b)

where

ν1 (τ, W ) := ε2 A (F(0, τ, W )) ,

ν2 (τ, W ) := −ε4 A2 [∂t F(0, τ, W ) + ∂W F(0, τ, W )Π (F(t, τ, W ))] .

As rigorously proved in [12] for fixed t ≥ 0 under some suitable Sobolev norm, if

W0 (x, τ ) = W01st (x, τ ),

then ∂tk W (x, t, τ ) = O(1),

k = 0, 1, 2, ε → 0,

(2.15)

k = 0, 1, 2, 3, ε → 0.

(2.16)

and if W0 (x, τ ) = W02nd (x, τ ),

then

∂tk W (x, t, τ ) = O(1),

W02nd

Thus, with the second order initial data, i.e. in (2.14b), a finite difference semi-discretisation is given in [12] to solve (2.13) by denoting W n (x, τ ) ≈ W (x, tn , τ ) as  τ τ  τ ∈ T, (2.17)  W n+1/2 (x, τ ) = W n (x, τ ) + F (tn , τ, W n (x, τ )) − 2 ∂τ W n+1/2 (x, τ ), 2 2ε 
τ   W n+1 (x, τ ) = W n (x, τ ) + τ F tn+1/2 , τ, W n+1/2 (x, τ ) − ∂τ W n+1 (x, τ ) + W n (x, τ ) . 2ε2

By letting τ = tn /ε2 , one can get w(x, tn ) ≈ W n (x, tn /ε2 ). Then by taking the inverse transformations of (2.11) and (2.9), one can obtain the approximation to uε (x, tn ). Considering 1D for instance, together with the truncation of the problem as (2.1) and then using Fourier 9 pseudospectral discretisation [20, 26] in x and τ , under assumption that W0 (x, τ ) ∈ L∞ τ (H ), the method (2.17) that solves (2.13a) for approximated W (x, t, τ ) could give finite time uniform error bound as
O (∆xm0 ) + O (∆τ m1 ) + O ∆t2 ,

where m1 > 0 is another constant that dependents on the regularity of the solution of (2.13). The above rigorous error bound has been established in [12] when the time step ∆t ≤ C0 with C0 > 0 a small constant but independent of ε. However, order reduction in the temporal accuracy of the above finite difference twoscale approach has been observed in [12, 15] in practical computing. Here we can also see the temporal convergence order reduction from the numerical results given in Table 4) when ε and ∆t are in a resonance regime (∆t and ε are comparable). As can be seen that the uniform second order temporal accuracy of the finite difference integrator (2.17) essentially relies on the boundedness of the ∂t3 U as ε → 0 in (2.16). Although one has the the boundedness in (2.16), there is still a period that ∂t3 U grows as ε decreases. Let us illustrate this issue by considering an 1D example in (2.10) with computation domain Ω = (−8, 8) and  2 i f (|v ε |2 )v ε = −|v ε |2 v ε , ϕ1 (x) = e−x , ϕ2 (x) = ∂xx ϕ1 (x) + f (|ϕ1 (x)|2 )ϕ1 (x) − sech(x2 ). 2

We compute the two quantities ∂t2 W (x, t, τ ) and ∂t3 W (x, t, τ ) of system (2.13) till t = 1. Figure 1 shows k∂t2 W kL∞ := k∂t2 W kL∞ (Ω×[0,1]×T) and k∂t3 W kL∞ := k∂t3 W kL∞ (Ω×[0,1]×T) of (2.13) under first order data and second order data, respectively as ε decreases. From the results in Figure 1, we can see that under the first or second order initial data, ∂t2 W and 3 ∂t W are indeed uniformly bounded as ε → 0 and they behaves increasingly as ε decreases. It is well-known [1, 2, 12] that without the corrections to the initial data in (2.13), e.g. W0 (x, τ ) = w0 (x), the time derivatives ∂t2 W and ∂t3 W will grow to infinity as ε → 0. This unbounded growth in the time derivative is the essential reason that leads to convergence order reduction or non-convergence of traditional numerical integrators. Thus, here even though under the first or second order corrected initial data that helps to bound the time derivative ∂t2 W or ∂t3 W for all ε ∈ (0, 1], within the finite increasing period in ∂t2 W or ∂t3 W (with respect to ε), the proposed numerical integrator for the two-scale formulation (2.13) will mimic the behaviour of non-uniformly accurate schemes to some extent, and this would be a reason to induce order reduction in the numerical scheme (2.17) as observed in [12, 15]. 6

1

∞

10

2

t

||∂2W||

L

||∂3t W||L∞

10

−2

−2

−1

10

−1

10

10

10

ε

ε

Figure 1: k∂t2 W kL∞ and k∂t3 W kL∞ as functions of ε in log-log scale: (2.13) under the first order data (left) or under the second order data (right).

3. Combined multiscale integrator two-scale method In this section, we combine the spirit of the multiscale (time) integrator [4, 8, 10, 11], which is to perform some multiscale decomposition on the equation and then integrate different parts with different suitable integrators, with the two-scale method to solve the NLSW. To do so, we first derive a multiscale decomposition by amplitude of the NLSW (1.1) and then propose the numerical scheme for the decomposed system. 3.1. Multiscale decomposition We assume that the solution uε of the NLSW (1.1) satisfying uε (x, t) = ψ(x, t) + ε2 φ(x, t),

x ∈ Rd , t ≥ 0,

(3.1)

with ψ and φ two undetermined functions. Plugging (3.1) into (1.1), we find 2i∂t ψ + 2iε2 ∂t φ − ε2 ∂tt ψ − ε4 ∂tt φ + ∆ψ + ε2 ∆φ + f (|ψ + ε2 φ|2 )(ψ + ε2 φ) = 0. We decompose the above equation as 2i∂t ψ − ε2 ∂tt ψ + ∆ψ + f (|ψ|2 )ψ = 0,

(3.2)

2i∂t φ − ε2 ∂tt φ + ∆φ + g(ψ, φ) = 0,

(3.3)

and where

 1  f (|ψ + ε2 φ|2 )(ψ + ε2 φ) − f (|ψ|2 )ψ . 2 ε Note that the nonlinearity g is bounded, i.e. g(ψ, φ) = O(1), ε → 0, but it is not gauge invariant any more, i.e. g(ψ, eis φ) 6= eis g(ψ, φ). For example, for the cubic nonlinearity f (ρ) = λρ, λ ∈ R, then we have explicitly

g(ψ, φ) = λ|ψ|2 φ + λ φ|ψ|2 + φψ 2 + ελ (φ)2 ψ + |φ|2 ψ + ελ|φ|2 ψ + λε2 |φ|2 φ. g(ψ, φ) :=

Hence, we remark that (3.3) is not a classical NLSW, and its gauge transformation is not a classical KG either. Concerning the initial data for the decomposed system (3.2) and (3.3), we let t = 0 in (3.1) and its time derivative where we find ϕ1 (x) = ψ(x, 0) + ε2 φ(x, 0),

ϕ2 (x) = ∂t ψ(x, 0) + ε2 ∂t φ(x, 0). 7

Here we have some freedom to choose the four unspecified initial values ψ(x, 0), φ(x, 0), ∂t ψ(x, 0), ∂t φ(x, 0) under the above two constraints. By matching the asymptotic order in ε, we have ψ(x, 0) = ϕ1 (x),

φ(x, 0) = 0.

For the initial velocity of (3.2), we choose the well-prepared initial data, i.e. ∂t ψ(x, 0) = and it leaves

 i ∆ψ(x, 0) + f (|ψ(x, 0)|2 )ψ(x, 0) . 2

1 [ϕ2 (x) − ∂t ψ(x, 0)] . ε2 Thus, we have decomposed the NLSW (1.1) into a well-behaved part: 
 2i∂t ψ(x, t) − ε2 ∂tt ψ(x, t) + ∆ψ(x, t) + f |ψ(x, t)|2 ψ(x, t) = 0,    ψ(x, 0) = ϕ1 (x), ∂t ψ(x, 0) = i ∆ϕ1 (x) + f (|ϕ1 (x)|2 )ϕ1 (x) , 2 ∂t φ(x, 0) =

and an energy-unbounded part:   2i∂t φ(x, t) − ε2 ∂tt φ(x, t) + ∆φ(x, t) + g(ψ(x, t), φ(x, t)) = 0,  φ(x, 0) = 0, ∂t φ(x, 0) = 1 [ϕ2 (x) − ∂t ψ(x, 0)] , x ∈ Rd . ε2

x ∈ Rd , t > 0, x ∈ Rd ,

x ∈ Rd , t > 0,

(3.4)

(3.5)

The above two decomposed systems are still oscillatory problems that propagate fast waves in time. (3.4) describes the leading order term in the asymptotic expansion of uε (x, t) and (3.5) describes the high order terms. Comparing to the original system (1.1), the solution of (3.4) shows milder oscillation due to the well-prepared initial data. We have boundedness of ∂tk ψ under some Sobolev norm for k = 1, 2 as ε → 0. On the contrast, (3.5) has worse behaviour than (1.1). The energy of the system (3.5) Eφ is unbounded under the ill-prepared case (1.6), i.e. Z  2  Eφ (t) := ε |∂t φ(x, t)|2 + |∇φ(x, t)|2 − F (|φ(x, t)|2 ) dx = O(ε−2 ), ε → 0, Rd


since the time derivative of the function is now unbounded, i.e. ∂t φ = O ε−2 . The oscillation in the solution of (3.5) is stronger than (1.1), but thanks to the expansion (3.1), the contribution of this highly oscillatory part is in the high order. We are going to solve the well-behaved part (3.4) by an exponential integrator and solve the highly oscillatory part (3.5) by means of the two-scale formulation approach. In the following to present the numerical discretisations, we turn to the 1D truncated problem (2.1). Generalisations to high dimensions are straightforward. 3.2. Exponential integrator for the well-behaved part Under the truncation (2.1), the well-behaved system (3.4) collapses to
 2i∂t ψ(x, t) − ε2 ∂tt ψ(x, t) + ∂xx ψ(x, t) + f |ψ(x, t)|2 ψ(x, t) = 0,     i ψ(x, 0) = ϕ1 (x), ∂t ψ(x, 0) = ∂xx ϕ1 (x) + f (|ϕ1 (x)|2 )ϕ1 (x) ,  2   ψ(−L, t) = ψ(L, t), ∂x ψ(−L, t) = ∂x ψ(L, t), t ≥ 0.

x ∈ Ω, t > 0, x ∈ Ω,

(3.6)

As reviewed in §2.1, the exponential integrator proposed in [2] if applied only solves (3.6) for the function value of ψ. Here we shall modify the scheme in spirit of general exponential integrator [21] in order to cover an approximation for ∂t ψ since the time derivative contributes the kinetic part to the energy (1.3) of the system and also to the mass (1.2). 8

Here we shall keep the form (2.3) and include the derivative equation. By letting s = ∆t in (2.3) or by taking a derivative with respect to s on both sides of (2.3) and then letting s = ∆t, we get Z ∆t 2 0 b b b ψl (tn+1 ) = al (∆t)ψl (tn ) + ε bl (∆t)ψl (tn ) + bl (∆t − θ)fbl (tn + θ)dθ, 0

ψbl0 (tn+1 ) = a0l (∆t)ψbl (tn ) + ε2 b0l (∆t)ψbl0 (tn ) +

Z

∆t

0

b0l (∆t − θ)fbl (tn + θ)dθ,

l=−

N N ,..., − 1. 2 2

Here al and bl are introduced in (2.4). We approximate the above two unknown integrals by a Gautschi-type quadrature [8, 4], e.g. Z ∆t Z ∆t h i b bl (∆t − θ) fbl (tn ) + θfbl0 (tn ) dθ. bl (∆t − θ)fl (tn + θ)dθ ≈ 0

0

Then we have

ψbl (tn+1 ) ≈ al (∆t)ψbl (tn ) + ε2 bl (∆t)ψbl0 (tn ) + cl (∆t)fbl (tn ) + dl (∆t)fbl0 (tn ), ψb0 (tn+1 ) ≈ a0 (∆t)ψbl (tn ) + ε2 b0 (∆t)ψb0 (tn ) + c0 (∆t)fbl (tn ) + d0 (∆t)fb0 (tn ), l

where

l

cl (s) :=

l

Z

0

l

l

s

bl (s − θ)dθ,

dl (s) :=

Z

l

l

s

θbl (s − θ)dθ,

0

s ≥ 0.

Thus, a detailed exponential integrator for (3.6) reads: Denote ψ n (x) ≈ ψ(x, tn ), ψ˙ n (x) ≈ ∂t ψ(x, tn ) and choose ψ 0 (x) = ψ(x, 0), ψ˙ 0 (x) = ∂t ψ(x, 0). Then for n ≥ 0, the scheme updates as bn bn ψbln+1 = al (∆t)ψbln + ε2 bl (∆t)ψ˙ l + cl (∆t)fbln + dl (∆t)f˙l ,

where

b ψ˙ l

n+1

n

n

b b = a0l (∆t)ψbln + ε2 b0l (∆t)ψ˙ l + c0l (∆t)fbln + d0l (∆t)f˙l , X

ψ (x) =

l=−N/2 N/2−1

ψ˙ n (x) =

(3.7b)

N/2−1

N/2−1 n

ψbln eiµl (x+L) ,

X bn ψ˙ eiµl (x+L) , l

n 2

X

n

f (|ψ | )ψ (x) =

l=−N/2

N/2−1

f˙n (x) =

l=−N/2

and

(3.7a) N N l = − ,..., − 1, 2 2

fbln eiµl (x+L) ,

X bn f˙ eiµl (x+L) , l

l=−N/2

x ∈ Ω,

i h f˙n (x) := f 0 (|ψ n (x)|2 ) (ψ n (x))2 ψ˙ n (x) + |ψ n (x)|2 ψ˙ n (x) + f (|ψ n (x)|2 )ψ˙ n (x).

The scheme (3.7) is a two-level explicit scheme. The memory cost is O(N ) and the computational cost is O(N logN ). Denote the error enψt (x) = ∂t ψ(tn , x) − ψ˙ n (x),

enψ (x) = ψ(tn , x) − ψ n (x),

n = 0, 1, . . . .

Assuming that in (3.6) f ∈ C 2 (R),

ϕ1 ∈ H 2 (Ω),

kψkL∞ ([0,T ];H 2 ) + k∂t ψkL∞ ([0,T ];L2 ) . 1,

then with the standard energy method and mathematical induction (or cut-off), we can establish the error bound for the exponential integrator (3.7) till a finite time T > 0 as ε2 kenψt kL2 + kenψ kH 1 . ∆t2 + ∆xm0 ,

0≤n≤

T . ∆t

The proof can be done by a very similar way as in [1, 4, 33, 32] and we omit the lengthy details here. 9

(3.8)

3.3. Two-scale integrator for the energy-unbounded part Now we turn to the energy-unbounded part (3.5) which under truncation (2.1) becomes  2i∂t φ(x, t) − ε2 ∂tt φ(x, t) + ∂xx φ(x, t) + g(ψ(x, t), φ(x, t)) = 0, x ∈ Ω, t > 0,    1 φ(x, 0) = 0, ∂t φ(x, 0) = 2 [ϕ2 (x) − ∂t ψ(x, 0)] , x ∈ Ω,  ε   φ(−L, t) = φ(L, t), ∂x φ(−L, t) = ∂x φ(L, t), t ≥ 0.

(3.9)

Let us put this equation under a suitable form in order to apply the two-scale techniques. We remark that once again, one can transfers the NLSW type equation (3.9) into a KG by the gauge transformation (2.9) and then perform the two-scale formulation based on the KG similarly as reviewed in §2.2, but here the nonlinearity g is not gauge invariant so one needs to put more efforts in writing the nonlinear terms in the two-scale formulation. Here, we present an alternative approach to derive the two-scale formulation by working on the NLSW (3.9) directly instead of going to KG. Setting   iε2 1 ∂ φ(x, t) , (3.10) φ(x, t) + v1 (x, t) = φ(x, t), v2 (x, t) = √ t 2 1 − 4ε2 ∂xx straightforward calculations lead to the system   p 2i     ∂t v1 (x, t) = 2 v1 (x, t) − 1 − 4ε2 ∂xx v2 (x, t) , x ∈ Ω, t > 0, ε i   (∂xx v1 (x, t) + g(ψ(x, t), v1 (x, t))) ,  ∂t v2 (x, t) = √ 2 1 − 4ε2 ∂xx or equivalently in a compact form

∂t v(x, t) = where we denote v = A=

v1 v2




2i Av(x, t) + g(ψ(x, t), v(x, t)), ε2

(3.11)

and ! √ − 1 − 4ε2 ∂xx , 0

1

2 √ ε ∂xx 4 1−4ε2 ∂xx

g(ψ, v) =

The operator matrix A can be diagonalised as   1 + ε2 λε 0 A=P P −1 , 0 −ε2 λε with P =

x ∈ Ω, t > 0,

1

1

2 ε √ −ε λ 1−4ε2 ∂xx

√ 1+ε λ 1−4ε2 ∂xx

2

ε

!

,

λε = −

P −1 =

s





0

i √ g(ψ, v1 ) 2 1−4ε2 ∂xx

1+

√

1 − 4ε2 ∂xx 1 − ε2 ∂xx

.

2∂xx , 1 − 4ε2 ∂xx 2

ε

√ 1+ε λ 1−4ε2 ∂xx 2 ε √ ε λ 1−4ε2 ∂xx

! −1 1

.

Note that in the sense of operator, all the terms in P and P −1 are uniformly bounded, since 0≤ √

ε2 λε 1 < , 4 1 − 4ε2 ∂xx

0< √

1 + ε2 λε ≤ 1. 1 − 4ε2 ∂xx

(3.12)

Then we filter out the main oscillation in (3.11) by letting 2

u(x, t) := e−2itA/ε v(x, t), and we have

  2 2 ∂t u(x, t) = e−2itA/ε g ψ(x, t); e2itA/ε u(x, t) , 10

x ∈ Ω, t > 0.

(3.13)

From now on, we shall omit the x-variable in all the unknowns for simplicity of notation, i.e. u(x, t) = u(t), ψ(x, t) = ψ(t). Let us now introduce the matrices    ε  1 0 λ 0 −1 B=P P , C=P P −1 , 0 0 0 −λε which commutes with each other and implies (A − B)/ε = C. Then we can rewrite (3.13) as ∂t u(t) = G(t, t/ε2 , u(t)), where

t > 0,

(3.14)


G(t, τ, u) := e−2iτ B e−2itC g ψ(t), e2iτ B e2itC u .

The induced function G(t, τ, u) is 2π-periodic in variable τ and it is smooth enough providing the original nonlinear function f (·) is smooth. In particular, together with noting (3.12), for some integers α, β, γ and α β γ s s > 2(α + β) + d/2,
the functional ∂t ∂τ ∂u G is continuous and locally bounded from R+ × T × H to s γ s−2α−2β L (H ) , H uniformly in ε, which based on results in [12] indicates that (3.14) is suitable to apply the scale separation. Now we perform the two-scale formulation on (3.14). Similarly as reviewed in §2.2, we separate the fast time variable t/ε2 out from the equation (3.14) and consider an augmented equation for U (t, τ ) as   ∂ U (t, τ ) + 1 ∂ U (t, τ ) = G(t, τ, U (t, τ )), τ ∈ T, t > 0, τ t ε2 (3.15)  U (0, τ ) = U0 (τ ), τ ∈ T. We require similarly U0 (0, τ ) = u(0) for

U (t, t/ε2 ) = u(t),

t ≥ 0.

Then the prepared initial data U (0, τ ) for (3.15) can be written down in the form as (2.14). Adopting the same operators introduced in §2.2, here we only call for the first order initial data, i.e. U0 (τ ) = u(0) + h1 (τ, u(0)) − h1 (0, u(0)), with

τ ∈ T,

(3.16)

h1 (τ, U ) := ε2 L−1 (I − Π) (G(0, τ, U )) .

The above initial data (3.16) is sufficient for providing ∂t U, ∂t2 U = O(1) as ε → 0 as we mentioned. Unlike in [12, 15] or as reviewed in §2.2, we would not go for higher order prepared initial data, e.g. the second order data (2.14b). The second order initial data (2.14b) which needs the computation of the derivatives of the nonlinearity, is certainly more involved in practise than the first order data (2.14a). On the other hand, the second order initial data theoretical wise requires higher regularity assumption on the chosen data U0 (τ ) in order to have the boundedness of higher order time derivatives. For instance in 1D, under the second order initial data in order to provide ∂t3 U (τ, t) = O(1) as ε → 0 for some t ≥ 0, one needs to assume that 9 U0 (τ ) = U0 (x, τ ) ∈ L∞ τ (H ) [12], which consequently requires more regularity on the given data ϕ1 and ϕ2 . We shall only work with the first order data and construct a uniformly second order accurate scheme for the two-scale problem (3.15) by means of an exponential integrator rather than the finite difference integrator (2.17). Choose Nτ > 0 an even integer. Applying Fourier transform on (3.15) for τ , we have

with

b 0 (t) + il U bl (t) = Gbl (t), U l ε2

l=−

Nτ Nτ ,..., − 1, 2 2

Nτ /2−1

U (t, τ ) =

X

l=−Nτ /2

bl (t)eilτ , U

G(t, τ, U (t, τ )) = 11

Nτ /2−1

X

l=−Nτ /2

Gbl (t)eilτ .

(3.17)

Integrating (3.17) from tn to tn+1 (n ≥ 0), we get bl (tn+1 ) = e− il∆t bl (tn ) + ε2 U U

Z

∆t

0

il

e− ε2 (∆t−s) Gbl (tn + s)ds.

We approximate the above integration by the following quadrature for n ≥ 1,   Z ∆t il d b − (∆t−s) − il∆t 2 2 b b b e ε Gl (tn ) + s Gl (tn ) ds Ul (tn+1 ) ≈ e ε Ul (tn ) + dt 0   il∆t bl (tn ) + pl Gbl (tn ) + ql 1 Gbl (tn ) − Gbl (tn−1 ) , ≈ e− ε2 U ∆t where

(3.18)

 2  il∆t  iε e− ε2 − 1 , l 6= 0, (∆t−s) pl := e ds = l  0 ∆t, l = 0,  2  ε  2 2 − il∆t  Z ∆t ε2 − il∆t  ε − ε e , l 6= 0, il l2 ql := e− ε2 (∆t−s) sds = 2  0   ∆t , l = 0. 2 For the starting value at t1 , we begin with a prediction as Z

∆t

− εil2

Nτ /2−1

b ∗ (t1 ) = e− il∆t bl (t0 ) + pl Gbl (t0 ), ε2 U U l

U ∗ (t1 , τ ) =

X

l=−Nτ /2

then do a correction as

b ∗ (t1 )eilτ , U l

  bl (t1 ) ≈ e− il∆t bl (t0 ) + pl Gbl (t0 ) + ql 1 Gbl∗ (t1 ) − Gbl (t0 ) , ε2 U U ∆t

with

Nτ /2−1

F(t1 , τ, U ∗ (t1 , τ )) =

X

l=−Nτ /2

n

Gbl∗ (t1 )eilτ .

Thus, denoting U (τ ) ≈ U (tn , τ ) for n ≥ 0, a detailed exponential integrator for the two-scale formulation (3.15) reads: choose U 0 (τ ) = U (0, τ ), then   b n+1 = e− il∆t bln + pl Gbln + ql 1 Gbln − Gbn−1 , ε2 U U l l ∆t   il∆t 1 − 1 0 0 bl = e ε2 U bl + pl Gbl + ql U Gbl∗ − Gbl0 , ∆t

b ∗ = e− il∆t b 0 + pl Fb0 and ε2 U where U l l l

Nτ /2−1

n

U (τ ) =

X

l=−Nτ /2 Nτ /2−1

U ∗ (τ ) =

X

l=−Nτ /2

Here for n ≥ 0,

n ≥ 1,

(3.19b)

Nτ /2−1

bln eilτ , U

G (τ, U (τ )) =

bl∗ eilτ , U

G 1 (τ, U ∗ (τ )) =

n

n

X

l=−Nτ /2 Nτ /2−1

X

l=−Nτ /2

Gbln eilτ ,

Gbl∗ eilτ .


G n (τ, u) := e−2iτ B e−2itn C g ψ n , e2iτ B e2itn C u , 12

(3.19a)

n ≥ 0,

with ψ n the numerical solution from (3.7). After obtaining U n (τ ) for n ≥ 1, we set  n   v1 tn 2itn A/ε2 n , v(x, tn ) ≈ := e U ε2 v2n

(3.20)

which consequently gives approximations to φ(x, tn ) and ∂t φ(x, tn ) from (3.10), i.e. i p 2i h φ(x, tn ) ≈ φn := v1n , ∂t φ(x, tn ) ≈ φ˙ n := 2 v1n − 1 − 4ε2 ∂xx v2n . (3.21) ε The above two-scale method (3.19) is explicit and efficient thanks to fast Fourier transform. The total computational cost including the Fourier pseudospectral discretisation in x direction is O(Nτ N log(Nτ N )) 5 per time step, where N is the number of grids in x direction. Under assumption U0 (x, τ ) ∈ L∞ τ (H ), which consequently is to assume sufficient regularity assumptions on the given functions f, ϕ1 , ϕ2 , for example f ∈ C 2 , ϕ1 ∈ H 6 , ϕ2 ∈ H 4 which indicates that in (3.10) v(x, 0) ∈ H 5 , based on results in [12] in 1D, the solution U (t, τ, x) of the (3.15) satisfies sup k∂tk U (t, ·, ·)kL∞ 5−2k ) . 1, τ (H

k = 0, 1, 2.

ε∈(0,1]

Then the error estimate for the Gautschi-type exponential integrator (3.19) on the transport type equation (3.15) can be established in a quite standard way [21, 22]. Global error bound of (3.19) with full spatial discretisation till a finite time T > 0 is T 2 m0 1 . ∆t + ∆x . kenU kL∞ + ∆τ m1 , enU (τ, x) := U (tn , τ, x) − U n (τ, x), 0 ≤ n ≤ τ (H ) ∆t We remark that this exponential integrator (3.19) has also been used to solve a two-scale formulation in our recent work for a kinetic problem [16]. Noting (3.20), (3.21) and (3.10), and denoting enφt (x) = ∂t φ(tn , x) − φ˙ n (x),

enφ (x) = φ(tn , x) − φn (x),

n = 0, 1, . . . .

we further have the error bound for the numerical solution to (3.9) as T . (3.22) ∆t Now, combining the schemes for the well-behaved part (3.7) and for the energy-unbounded part (3.19)(3.21) which we call it as the combined multiscale integrator with two-scale formulation (MITSF) method, we get the numerical approximations based on the ansatz (3.1) to both the function value and the first order time derivative of solution of the original NLSW (1.1): ε2 kenφt kL2 + kenφ kH 1 . ∆t2 + ∆xm0 + ∆τ m1 ,

un (x) := ψ n (x) + ε2 φn (x),

0≤n≤

unt (x) := ψtn (x) + ε2 φnt (x).

(3.23)

Denoting enut (x) = ∂t u(tn , x) − u˙ n (x),

enu (x) = u(tn , x) − un (x),

n = 0, 1, . . . ,

then the error bound (3.22) together with (3.8) shows an error bound of the MITSF method towards solving (1.1) as T ε2 kenut kL2 + kenu kH 1 . ∆t2 + ∆xm0 + ∆τ m1 , 0 ≤ n ≤ . (3.24) ∆t This error estimate further indicates that for the numerical mass and energy: Z  n  n (3.25a) M := |u (x)|2 dx − ε2 Im (un (x)u˙ n (x)) dx, ZΩ  2 n  E n := ε |u˙ (x, t)|2 + |∂x un (x)|2 − F (|un (x)|2 ) dx, n = 0, 1, . . . , (3.25b) Ω

which are approximations to M (t) in (1.2) and E(t) in (1.3) at t = tn resulting from the numerical solution of MITSF, we have uniform error bounds as |M n − M (tn )| . ∆t2 + ∆xm0 + ∆τ m1 ,

|E n − E(tn )| . ∆t2 + ∆xm0 + ∆τ m1 , 13

0≤n≤

T . ∆t

(3.26)

Table 1: Spatial error in x of MITSF for different ε at time T = 1 under ∆t = 10−6 and Nτ = 32. τ e∆t,∆x,N ε ε0 = 0.5 ε0 /21 ε0 /22 ε0 /23 ε0 /24 ε0 /25 ε0 /26 ε0 /28 ε0 /210

∆x = 1 1.74E-1 2.37E-1 2.59E-1 2.64E-1 2.66E-1 2.67E-1 2.67E-1 2.67E-1 2.67E-1

∆x/2 7.70E-3 6.00E-3 6.00E-3 5.50E-3 5.40E-3 5.40E-3 5.30E-3 5.30E-3 5.30E-3

∆x/4 6.23E-6 3.74E-6 1.65E-6 5.46E-7 3.66E-7 4.02E-7 3.76E-7 3.78E-7 3.78E-7

∆x/8 2.53E-10 3.54E-11 1.77E-11 2.12E-11 1.57E-12 1.41E-12 4.28E-13 1.29E-13 1.32E-13

Table 2: Spatial error in τ of MITSF for different ε at time T = 1 under ∆t = 10−6 and ∆x = 1/16. τ e∆t,∆x,N ε ε0 = 0.5 ε0 /21 ε0 /22 ε0 /23 ε0 /24 ε0 /25 ε0 /26 ε0 /28 ε0 /210

Nτ = 2 2.44E-1 2.98E-1 3.10E-1 3.48E-1 3.50E-1 3.84E-1 2.93E-1 3.27E-1 3.31E-1

Nτ = 4 2.37E-2 6.40E-3 1.80E-3 4.69E-4 1.17E-4 2.94E-5 7.35E-6 4.59E-7 2.87E-8

Nτ = 8 3.25E-4 3.74E-6 6.72E-8 1.11E-9 1.93E-11 5.91E-13 2.15E-13 1.72E-13 1.46E-13

Nτ = 16 1.55E-8 3.58E-12 2.60E-13 1.45E-13 1.69E-13 2.05E-13 2.24E-13 1.84E-13 1.22E-13

Nτ = 32 2.25E-13 3.01E-13 2.00E-13 2.95E-13 1.29E-13 4.15E-13 2.18E-13 1.01E-13 1.18E-13

4. Numerical result In this section, we conduct numerical tests on the proposed multiscale integrator with two-scale formulation (MITSF) method and present numerical results. For comparison purpose, we also present results of the direct exponential integrator Fourier pseudospectral (EIFP) method (2.7), the direct two-scale formulation on KG with finite difference (TSFKG) discretisation (2.17). In the end, we apply the MITSF method investigate the convergence rate of the NLSW to NLS in the limit regime. We take d = 1 in the NLSW (1.1) with the cubic nonlinearity f (|uε |2 )uε = −|uε |2 uε , and choose the ill-prepared initial data: 2

ϕ1 (x) = e−x ,

ϕ2 (x) =

 i ∂xx ϕ1 (x) + f (|ϕ1 (x)|2 )ϕ1 (x) − sech(x2 ), 2

i.e. r(x) = −sech(x2 ) in (1.6). We take the computational domain Ω = (−8, 8), i.e. L = 8 in (2.1) and solve the NLSW till T = 1 where the solution is still localised in the domain and is away from boundary. We measure the error uε (tn , x) − un (x) under H 1 -norm and denote  τ τ τ e∆t,∆x,N (T = tn ) := kuε (tn , ·) − un (·)kH 1 , e∆t,∆x,N (T ) := max e∆t,∆x,N (T ) , ε ∞ ε ε

for the methods TSFKG and MITSF. Since EIFP has nothing to do with τ , we denote in the corresponding ∆t,∆x error functions as eε∆t,∆x (T ) and e∞ (T ). We compute the error at T = 1. The ‘exact’ solution here is obtained numerically by the MITSF method with very small step size, e.g. ∆t = 10−6 , ∆x = 1/16, Nτ = 64. 14

Table 3: Temporal error of EIFP for different ε at time T = 1 under ∆x = 1/16.

e∆t,∆x ε ε0 = 0.5 rate ε0 /21 rate ε0 /22 rate ε0 /23 rate ε0 /24 rate ε0 /25 rate ε0 /26 rate ε0 /28 rate ε0 /210 rate e∆t,∆x ∞ rate

∆t = 0.2 4.60E-2 – 1.86E-1 – 8.18E-2 – 2.71E-2 – 2.66E-2 – 6.46E-2 – 1.26E-1 – 3.71E-2 – 4.92E-2 – 1.86E-1 –

∆t/41 3.30E-3 1.90 1.31E-2 1.91 2.62E-2 0.82 5.20E-3 1.19 1.60E-3 2.03 1.40E-3 2.76 2.00E-3 2.97 2.20E-3 2.03 1.50E-3 2.50 2.62E-2 1.41

∆t/42 2.10E-4 1.99 8.19E-4 2.00 3.60E-3 1.43 3.90E-3 0.21 9.11E-4 0.41 1.31E-4 1.71 1.65E-4 1.80 4.36E-4 1.16 1.08E-4 1.90 3.90E-3 1.37

∆t/43 1.32E-5 2.00 5.15E-5 2.00 2.20E-4 2.01 9.09E-4 1.05 9.51E-4 -0.03 2.08E-4 -0.33 3.29E-5 1.16 5.28E-6 3.18 8.94E-6 1.80 9.51E-4 1.01

∆t/44 8.26E-7 2.00 3.22E-6 2.00 1.38E-5 2.00 5.57E-5 2.01 2.28E-4 1.03 2.33E-4 -0.08 5.27E-5 -0.34 2.75E-6 0.47 1.75E-6 1.17 2.33E-4 1.01

∆t/45 5.12E-8 2.01 2.03E-7 1.99 8.64E-7 2.00 3.49E-6 2.00 1.39E-5 2.01 5.72E-5 1.01 5.85E-5 -0.07 1.83E-6 0.29 2.30E-7 1.46 5.85E-5 0.99

Since there are no CFL-type conditions constricted on each method, so we could test the temporal error and spatial error separately. For the spatial error in x (or in τ ), we use the time step ∆t and mesh size ∆τ (or ∆x) sufficiently small such that discretisation error in t and in τ (or in x) is negligible compared to the spatial error in x (or in τ ), e.g. ∆t = 10−5 , ∆τ = π/32 (or ∆x = 1/16). Table 1 shows the spatial discretisation error eε∆t,∆x,Nτ of the MITSF method in x-direction under different ∆x and ε, and Table 2 shows the spatial error of MITSF in τ -direction under different Nτ (= 2π/∆τ ) and ε. Here we omit the corresponding spatial errors of EIFP and TSFKG, since the results of EIFP are similar to MITSF in xdirection and the results of TSFKG are similar to MITSF in both x and τ . For the temporal error, we use the mesh size sufficiently small such that the spatial discretisation error in both x and τ is negligible compared to the temporal discretisation error, e.g. ∆x = 1/16, ∆τ = π/32. Tables 3-5 show the temporal ∆t,∆x,Nτ τ errors e∆t,∆x,N and e∞ (or eε∆t,∆x and e∆t,∆x ) together with convergence rate of the EIFP, TSF-FD ε ∞ and MITSF, respectively under different ∆t and ε. We compare the uniform accuracy of EIFP, TSFKG and τ MITSF by collecting the temporal error e∆t,∆x,N or e∆t,∆x of the three methods, and they are shown with ∞ ∞ computational time (cputime) in Table 6. The schemes are programmed in Matlab and run on a 2.5GHz CPU laptop. Based on the results in Tables 1-6, we can draw the following conclusions: (i) The EIFP, TSFKG and MITSF have uniform spectral accuracy in x for all ε ∈ (0, 1] (c.f. Table 1). TSFKG and MITSF also have uniform spectral accuracy in τ (c.f. Table 2). The smaller the ε is, the less Nτ is needed to reach machine accuracy. (ii) Under the ill-prepared initial data case (1.6), EIFP method shows second order temporal accuracy when ∆t . ε2 (c.f. the upper diagonal part of Table 3) or when ε . ∆t (c.f. the lower diagonal part of Table 3). However in the resonance regime (c.f. the diagonal part of Table 3), the method has significant order reduction. In terms of uniform accuracy for all ε ∈ (0, 1], EIFP only has first order temporal accuracy (c.f. the last row of Table 3). (iii) TSFKG method also has order reduction for some ε and ∆t, in particular for ε in the intermediate 15

Table 4: Temporal error of TSFKG for different ε at time T = 1 under ∆x = 1/16, Nτ = 64.

τ e∆t,∆x,N ε ε0 = 0.5 rate ε0 /21 rate ε0 /22 rate ε0 /23 rate ε0 /24 rate ε0 /25 rate ε0 /26 rate ε0 /28 rate ε0 /210 rate τ e∆t,∆x,N ∞ rate

∆t = 0.2 2.85E-2 – 9.20E-3 – 6.20E-3 – 5.20E-3 – 5.40E-3 – 5.50E-3 – 5.50E-3 – 5.50E-3 – 5.50E-3 – 2.85E-2 –

∆t/41 1.40E-3 2.16 1.20E-3 1.45 2.82E-4 2.23 2.01E-4 2.34 2.06E-4 2.35 2.10E-4 2.35 2.11E-4 2.35 2.11E-4 2.35 2.11E-4 2.35 1.40E-3 2.16

∆t/42 1.00E-4 1.90 3.92E-4 0.81 5.56E-5 1.17 1.51E-5 1.86 1.22E-5 2.03 1.27E-5 2.02 1.29E-5 2.01 1.29E-5 2.01 1.29E-5 2.01 3.92E-4 0.92

∆t/43 6.33E-6 1.99 2.11E-5 2.10 4.78E-5 0.11 2.84E-6 1.21 9.52E-7 1.84 7.69E-7 2.02 7.95E-7 2.01 8.07E-7 2.00 8.08E-7 2.00 4.78E-5 1.51

∆t/44 3.90E-7 2.01 1.31E-6 2.00 2.94E-6 2.01 2.84E-7 1.66 1.65E-7 1.27 5.79E-8 1.86 4.77E-8 2.03 4.96E-8 2.01 4.97E-8 2.01 2.94E-6 2.01

∆t/45 1.83E-8 2.20 6.13E-8 2.20 1.34E-7 2.28 2.25E-7 0.17 1.27E-8 1.85 7.23E-9 1.50 2.42E-9 2.14 2.31E-9 2.21 2.33E-9 2.20 2.25E-7 1.85

regime (c.f. Table 4), but the resonance regime between ε and ∆t is smaller than EIFP. The temporal uniform accuracy of TSFKG is second order after time step is small enough (c.f. the last row of Table 4). (iv) MITSF method shows second order temporal accuracy either in terms of uniform accuracy or for a fixed ε ∈ (0, 1] (c.f. Table 5), which together with spatial results justify the error bound (3.24). The method does not show any significant order reduction for all pairs of ε and ∆t. With fixed ∆t, the approximation of MITSF is much more accurate than EIFP when ε becomes small. Compared to TSFKG, MITSF is much more accurate and efficient in the intermediate regime of ε and in terms of the uniform accuracy (c.f. Table 6). Next, we measure the fluctuation of the numerical mass M n and energy E n in (3.25) from the exact values during the computation of the MITSF method. The ‘exact’ mass (1.2) and energy (1.3) of the problem are computed from the initial data, i.e. Z Z    2  M (0) = |ϕ1 |2 − ε2 Im (ϕ1 ϕ2 ) dx, E(0) = ε |ϕ2 |2 + |ϕ01 |2 − F (|ϕ1 |2 ) dx, R

R

with a fine mesh size for spatial discretisation, e.g. ∆x = 1/8. The errors |M n − M (0)| and |E n − E(0)| of MITSF with different ∆t and ε (under the fixed mesh size ∆x = 1/8, Nτ = 32) till t = 5 are shown in Figure 2. From the results in Figure 2, we can see that under the time discretisation, the numerical mass and energy (3.25) during the computation of the MITSF method, are just small fluctuations from the exact values. As time step decreases, the numerical mass and energy of the MITSF method converges at uniformly second order rate to the exact values, which verifies the estimate (3.26). Lastly, we apply the proposed method to study the convergence of the NLSW (1.1) as ε → 0 in 1D to

16

Table 5: Temporal error of MITSF for different ε at time T = 1 under ∆x = 1/16, Nτ = 64.

τ e∆t,∆x,N ε ε0 = 0.5 rate ε0 /21 rate ε0 /22 rate ε0 /23 rate ε0 /24 rate ε0 /25 rate ε0 /26 rate ε0 /28 rate ε0 /210 rate τ e∆t,∆x,N ∞ rate

∆t = 0.2 1.95E-2 – 2.25E-2 – 2.62E-2 – 2.83E-2 – 2.86E-2 – 3.13E-2 – 2.84E-2 – 3.05E-2 – 3.09E-2 – 3.13E-2 –

∆t/41 1.20E-3 2.02 1.80E-3 1.83 1.80E-3 1.94 1.70E-3 2.04 1.60E-3 2.08 1.60E-3 2.14 1.50E-3 2.12 1.60E-3 2.11 1.60E-3 2.15 1.80E-3 2.06

∆t/42 7.45E-5 2.00 1.08E-4 2.02 1.31E-4 1.89 1.18E-4 1.92 9.86E-5 2.01 9.64E-5 2.03 9.69E-5 1.98 1.00E-4 2.00 9.57E-5 2.03 1.31E-4 1.89

∆t/43 4.67E-6 2.00 6.75E-6 2.00 8.04E-6 2.01 8.59E-6 1.89 6.54E-6 1.96 6.02E-6 2.00 5.97E-6 2.01 5.96E-6 2.03 5.99E-6 2.00 8.59E-6 1.97

∆t/44 2.88E-7 2.01 4.15E-7 2.01 4.95E-7 2.01 5.25E-7 2.01 5.38E-7 1.80 3.90E-7 1.97 3.68E-7 2.01 3.67E-7 2.01 3.58E-7 2.01 5.38E-7 2.00

∆t/45 1.37E-8 2.20 1.95E-8 2.20 2.32E-8 2.20 2.46E-8 2.20 2.52E-8 2.20 2.55E-8 1.97 2.02E-8 2.09 1.72E-8 2.20 1.72E-8 2.21 2.55E-8 2.19

τ Table 6: Comparisons of uniform temporal accuracy of the three methods at T = 1 with computation time (second): e∆t,∆x,N ∞ ∆t,∆x under ∆x = 1/8, Nτ = 16 for TSFKG and MITSF; e∞ under ∆x = 1/8 for EIFP.

EIFP rate cputime TSFKG rate cputime MITSF rate cputime

∆t = 0.2 1.86E-1 – 4.7E-4 2.85E-2 – 1.1E-2 3.13E-2 – 7.8E-3

∆t/41 2.62E-2 1.41 2E-3 1.40E-3 2.16 4.2E-2 1.80E-3 2.06 3.1E-2

∆t/42 3.90E-3 1.37 6.9E-3 3.92E-4 0.92 1.7E-1 1.31E-4 1.89 9.4E-2

∆t/43 9.51E-4 1.01 2.7E-2 4.78E-5 1.51 6.9E-1 8.59E-6 1.97 5.6E-1

the nonlinear Schr¨odinger (NLS) equation ( 2i∂t u(x, t) + ∆u(x, t) + f (|u(x, t)|2 )u(x, t) = 0, u(x, 0) = ϕ1 (x),

x ∈ R.

∆t/44 2.33E-4 1.01 1.2E-1 2.94E-6 2.01 2.7 5.38E-7 2.00 2.3

x ∈ R, t > 0,

∆t/45 5.85E-5 0.99 4.2E-1 2.25E-7 1.85 10.9 2.55E-8 2.19 8.8

(4.1)

Take the initial data as the smooth pair 2

ϕ1 (x) = e−x ,

ϕ2 (x) = sech(x2 ), 17

(4.2)

ε=0.5

−4

x 10

ε=0.05

−5

x 10

ε=0.005

−5

x 10

1.4 8

8

6

6

∆ t=0.02

0.8

error

error

1

∆ t=0.01 0.6

error

1.2

4

4

∆ t=0.005

0.4 2

2

0.2 0 0

1

2

3

4

0 0

5

1

2

t −4

ε=0.5

x 10

3

4

0 0

5

1

2

t ε=0.05

−4

x 10

3

4

5

4

5

t ε=0.005

−4

x 10

1.2

1.2

1

1

0.8

0.8

5

error

error

∆ t=0.02

4

∆ t=0.01 3 ∆ t=0.005

error

6

0.6

0.6

2

0.4

0.4

1

0.2

0.2

0 0

1

2

3

4

0 0

5

t

1

2

3

4

0 0

5

1

2

t

3

t

Figure 2: Mass error |M n − M 0 | (first row) and energy error |E n − E 0 | (second row) of the MITSF method under different ∆t and ε. 0.04

80 ε=0.02

2

ε=0.005 0.02 0.01 0 0

ε=0.02 ε=0.01 ε=0.005

60

ε=0.01

η(t)/ε

η(t)

0.03

40 20

2

4

6

8

0 0

10

2

4

t 0.2

10

8

η(t)/ε

η(t)

8

10 ε=0.02 ε=0.01 ε=0.005

0.15 0.1

6 ε=0.02

4

ε=0.01

0.05 0 0

6

t

2 2

4

6

8

0 0

10

t

ε=0.005 2

4

6

8

10

t

Figure 3: Convergence of NLSW to NLS: η(t) under the smooth data (4.2) (first row) and the non-smooth data (4.3) (second row).

18

or the non-smooth pair

2

ϕ1 (x) = x|x|e−x ,

ϕ2 (x) = sech(x2 ).

(4.3)

Both initial data (4.2) and (4.3) are the ill-prepared case in (1.6). We solve the NLSW (2.1) by the MITSF method for the solution uε , and solve the NLS (4.1) by the standard Strang time splitting spectral method [3] for the solution u, with very fine mesh on a bounded interval Ω = (−128, 128), e.g. ∆t = 10−4 , ∆x = 1/8, Nτ = 32. We define the error function as η(t) := kuε (·, t) − u(·, t)kH 1 .

(4.4)

Figure 2 shows η(t) with the smooth data (4.2) or the non-smooth data (4.3) under different ε. From Figure 2, we can draw the following conclusions: (i) The solution of the NLSW equation (1.1) converges to that of the NLS (4.1) quadratically in ε (not uniformly in time) provided that the initial data in (1.1) is smooth, i.e. kuε (·, t) − u(·, t)kH 1 ≤ (C0 + C1 T )ε2 ,

0 ≤ t ≤ T,

where the constants C0 , C1 > 0 are independent of ε and T (c.f. the first row of Figure 2). (ii) If the regularity of the initial data is weaker, e.g. ϕ1 and ϕ2 ∈ H 2 (Ω), then the convergence rate collapses to linear rate, i.e. kuε (·, t) − u(·, t)kH 1 ≤ (C2 + C3 T )ε,

0 ≤ t ≤ T,

where C3 and C4 are two positive constants which are independent of ε and T (c.f. the second row of Figure 2). Rigorous mathematical justification for these numerical observations is on-going. (iii) The observed convergence results of the NLSW are consistent with the convergence results reported in [10, 11] for the nonlinear Klein-Gordon equation and related models in the nonrelativistic limit. 5. Conclusions In this paper, we solved numerically the nonlinear Schr¨ odinger equation with wave operator (NLSW) with a small parameter ε which goes to zero in the highly oscillatory regime, by means of a two-scale formulation approach with a pre-decomposition by amplitude on the equation. We first reviewed and remarked on the convergence order reduction issue of the existing numerical methods for the NLSW in the highly-oscillatory regime with particularly the ill-prepared initial data [1]. Then we preformed the multiscale pre-decomposition of NLSW by decomposing the equation into a well-behaved part and an energy-unbounded part. We integrated the well-behaved by an exponential integrator and solved the energy-unbounded part by the two-scale formulation. Extensive numerical results showed that the proposed multiscale integrator with two-scale formulation (MITSF) method has uniform second order accuracy in time and uniform spectral accuracy in space for all ε ∈ (0, 1]. For all pairs of time step ∆t . 1 and ε ∈ (0, 1], there are no notable convergence order reductions. Comparisons were made with direct exponential integrator and direct twoscale method which showed the significant improvement of MITSF in uniform accuracy. Applications were made to study the convergence of the NLSW in the limit regime. Acknowledgements This work was supported by the French ANR project MOONRISE ANR-14-CE23-0007-01 and the IPL FRATRES. The author would like to thank Prof. Weizhu Bao, Philippe Chartier, Mohammed Lemou and Florian M´ehats for the stimulating discussions and helpful suggestions.

19

References [1] W. Bao, Y. Cai, Uniform error estimates of finite difference methods for the nonlinear Schr¨ odinger equation with wave operator, SIAM J. Numer. Anal. 50 (2012) pp. 492-521. [2] W. Bao, Y. Cai, Uniform and optimal error estimates of an exponential wave integrator sine pseudospectral method for the nonlinear Schr¨ odinger equation with wave operator, SIAM J. Numer. Anal. 52 (2014) pp. 1103-1127. [3] W. Bao, Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, Kinet. Relat. Models. 6 (2013) pp. 1-135. [4] W. Bao, Y. Cai, X. Zhao, A uniformly accurate multiscale time integrator pseudospectral method for the Klein-Gordon equation in the nonrelativistic limit regime, SIAM J. Numer. Anal. 52 (2014) pp. 2488-2511. ´, T. Colin, A singular perturbation problem for an enveloppe equation in plasma physics, Phys. D, 84 (1995) pp. [5] L. Bere 437-459. [6] W. Bao, X. Dong, Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime, Numer. Math. 120 (2012) pp. 189-229. [7] W. Bao, X. Dong, J. Xin, Comparisons between sine-Gordon and perturbed nonlinear Schr¨ odinger equations for modeling light bullets beyond critical collapse, Phys. D 239 (2010) pp. 1120-1134. [8] W. Bao, X. Dong, X. Zhao, Uniformly correct multiscale time integrators for highly oscillatory second order differential equations, J. Math. Study 47 (2014) pp. 111-150. [9] W. Bao, X. Dong, X. Zhao, An exponential wave integrator pseudospectral method for the Klein-Gordon-Zakharov system, SIAM J. Sci. Comput. 35 (2013) pp. A2903-A2927. [10] W. Bao, X. Zhao, A uniformly accurate (UA) multiscale time integrator Fourier pseoduspectral method for the KleinGordon-Schr¨ odinger equations in the nonrelativistic limit regime, Numer. Math. 135 (2017) pp. 833-873. [11] W. Bao, X. Zhao, A uniformly accurate multiscale time integrator spectral method for the Klein-Gordon-Zakharov system in the high-plasma-frequency limit regime, J. Comp. Phys. 327 (2016) pp. 270-293. ´hats, Uniformly accurate numerical schemes for highly oscillatory [12] P. Chartier, N. Crouseilles, M. Lemou, F. Me Klein-Gordon and nonlinear Schr¨ odinger equations, Numer. Math. 129 (2015) pp. 211-250. [13] T. Colin, P. Fabrie, Semidiscretization in time for Schr¨ odinger-waves equations, Discrete Contin. Dyn. Syst. 4 (1998) pp. 671-690. ´hats, Asymptotic preserving schemes for the highly oscillatory Vlasov-Poisson equa[14] N. Crouseilles, M. Lemou, F. Me tions, J. Comp. Phys. 248 (2013) pp. 287-308. ´hats, X. Zhao, Uniformly accurate forward semi-Lagrangian methods for highly [15] N. Crouseilles, M. Lemou, F. Me oscillatory Vlasov-Poisson equations, Multiscale Model. Simul. 15 (2017) pp. 723-744. ´hats, X. Zhao, Uniformly accurate particle-in-cell method for the long time two[16] N. Crouseilles, M. Lemou, F. Me dimensional Vlasov-Poisson equation with strong magnetic field, Preprint (2016) hal-01418976. [17] X. Dong, X. Xu, X. Zhao, On time-splitting pseudospectral discretization for nonlinear Klein-Gordon equation in nonrelativistic limit regime, Commun. Comput. Phys. 16 (2014) pp. 440-466. [18] W. Gautschi, Numerical integration of ordinary differential equations based on trigonometric polynomials, Numer. Math. 3 (1961) pp. 381-397. [19] B. Guo, H. Liang, On the problem of numerical calculation for a class of systems of nonlinear Schr¨ odinger equations with wave operator, J. Numer. Methods Comput. Appl. 4 (1983) pp. 176-182. [20] J.S. Hesthaven, S. Gottlieb, D. Gottlieb, Spectral Methods for Time-Dependent Problems, Cambridge University Press, Cambridge, New York, 2007 [21] M. Hochbruck, Ch. Lubich, A Gautschi-type method for oscillatory second-order differential equations, Numer. Math. 83 (1999) pp. 403-426. [22] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numer. 19 (2010) pp. 209-286. [23] N. Masmoudi, K. Nakanishi, From nonlinear Klein-Gordon equation to a system of coupled nonlinear Schr¨ odinger equations, Math. Ann. 324 (2002) pp. 359-389. [24] S. Machihara, The nonrelativistic limit of the nonlinear Klein-Gordon equation, Funkcial. Ekvac. 44 (2001) pp. 243-252. [25] S. Machihara, K. Nakanishi, T. Ozawa, Nonrelativistic limit in the energy space for nonlinear Klein-Gordon equations, Math. Ann. 322 (2002) pp. 603-621. [26] J. Shen, T. Tang, L. Wang, Spectral Methods: Algorithms, Analysis and Applications, Springer-Verlag, Berlin Heidelberg, 2011 [27] J. Shi, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations, SIAM J. Sci. Comput. 21 (1999) pp. 441-454. [28] M. Tsutumi, Nonrelativistic approximation of nonlinear Klein-Gordon equations in two space dimensions, Nonlinear Anal. 8 (1984) pp. 637-643. [29] T. Wang, L. Zhang, Analysis of some new conservative schemes for nonlinear Schr¨ odinger equation with wave operator, Appl. Math. Comput. 182 (2006) pp. 1780-1794. [30] J. Xin, Modeling light bullets with the two-dimensional sine-Gordon equation, Phys. D 135 (2000) pp. 345-368. [31] Y. Yu, Vortex dynamics for nonlinear Klein-Gordon equation, J. Differential Equations 251 (2011) pp. 970-994. [32] X. Zhao, An exponential wave integrator pseudospectral method for the symmetric regularized-long-wave equation, J. Comput. Math. 34 (2016) pp. 49-69. [33] X. Zhao, On error estimates of an exponential wave integrator sine pseudospectral method for the Klein-Gordon-Zakharov system, Numer. Methods Partial Differ. Equ. 32 (2016) pp. 266-291.

20