Low Regularity Exponential-Type Integrators for Semilinear Schrödinger Equations

Low Regularity Exponential-Type Integrators for Semilinear Schrödinger Equations Found Comput Math (2018) 18:731–755 https://doi.org/10.1007/s10208-017-9352-1 Low Regularity Exponential-Type Integrators for Semilinear Schrödinger Equations 1 2 Alexander Ostermann · Katharina Schratz Received: 24 March 2016 / Revised: 22 March 2017 / Accepted: 3 May 2017 / Published online: 26 May 2017 © The Author(s) 2017. This article is an open access publication Abstract We introduce low regularity exponential-type integrators for nonlinear Schrödinger equations for which first-order convergence only requires the bound- edness of one additional derivative of the solution. More precisely, we will prove r r +1 first-order convergence in H for solutions in H (with r > d/2) of the derived schemes. This allows us lower regularity assumptions on the data than for instance required for classical splitting or exponential integration schemes. For one-dimensional quadratic Schrödinger equations, we can even prove first-order convergence without any loss of regularity. Numerical experiments underline the favorable error behav- ior of the newly introduced exponential-type integrators for low regularity solutions compared to classical splitting and exponential integration schemes. Keywords Nonlinear Schrödinger equations · Exponential-type time integrator · Low regularity · Convergence Mathematics Subject Classification 65M12 · 65M70 · 35Q41 Communicated by Arieh Iserles. Alexander Ostermann alexander.ostermann@uibk.ac.at Katharina Schratz katharina.schratz@kit.edu Department of Mathematics, University of Innsbruck, Technikerstr. 13, 6020 Innsbruck, Austria Fakultät für Mathematik, Karlsruhe Institute of Technology, Englerstr. 2, 76131 Karlsruhe, Germany 123 732 Found Comput Math (2018) 18:731–755 1 Introduction Semilinear Schrödinger equations, in particular those of type 2 p i ∂ u =−u + μ|u| u, p ∈ N (1) with μ ∈ R are nowadays extensively studied numerically. In this context, split- ting methods (where the right-hand side is split into the kinetic and nonlinear part, respectively) as well as exponential integrators (including Lawson-type Runge–Kutta methods) contribute particularly attractive classes of integration schemes. For an extensive overview on splitting and exponential integration methods, we refer to [14,16,18,26], and for their rigorous convergence analysis in the context of semi- linear Schrödinger equations, we refer to [2,4,6,8,9,11,12,25,31] and the references therein. However, within the construction of all these numerical methods the stiff part (i.e., the terms involving the differential operator ) is approximated in one way or another which generally requires the boundedness of two additional spatial derivatives of the exact solution. In particular, convergence of a certain order only holds under sufficient additional regularity assumptions on the solution. In the following, we will illustrate the local error behavior of classical splitting and exponential integration methods and the thereby introduced smoothness requirements. Splitting schemes The Strang splitting for the cubic Schrödinger equation (i.e., p = 1) n+1/2 i  n u = e u , n+1/2 n+1/2 −iτμ u n+1/2 u = e u , + − n+1/2 n+1 i u = e u , where the right-hand side is split into the kinetic T (u) = i u and nonlinear part V (u) =−i μ|u| u is rigorously analyzed in [25]. In particular, its second-order con- r r +4 r vergence in H for solutions in H and its first-order convergence in H for solutions r +2 in H holds for all r ≥ 0, see [25] and [10], respectively. The result follows from the fact that the local error of the splitting scheme can be expressed through the double Lie commutator [T , [T , V ]](u) and the Lie commutator [T , V ](u) for second- and first-order methods, respectively. The latter reads [T , V ](u) = (∇u ·∇u) u + (∇uu) ·∇u + (u∇u) ·∇u + uu u, 2μ see [25, Section 4.2]. Due to the appearance of u in the local error, the boundedness of at least two additional derivatives of the exact solution is required. Based on [25], fractional error estimates for the Strang splitting were established in [10] which require the boundedness of 2 + 2γ additional derivatives for convergence of order 1 + γ , with 0 <γ < 1. Furthermore, in [19] first-order convergence of a filtered Lie splitting method for Schrödinger equations of type (1) with 1 ≤ 2 p < 4/d was shown in 2 d 2 d L (R ) for solutions in H (R ). 123 Found Comput Math (2018) 18:731–755 733 Classical exponential integrators The classical first-order exponential integrator in the cubic case reads e − 1 n+1 iτ n n 2 n u = e u − i μτ ϕ (iτ) |u | u with ϕ (z) = . (2) 1 1 Its construction is based on Duhamel’s formula iτ i (τ −s) 2 p u(t + τ) = e u(t ) − i μ e |u(t + s)| u(t + s)ds n n n n by applying the approximation u(t + s) ≈ u(t ) (3) n n in the integral terms and solving the remaining integral over the free Schrödinger i (τ −s) group e exactly. Due to the approximation (3), the above scheme introduces a local error of the form 2 2 2 4 2 2 ∂ i |u| u = i ∂ uu + 2i |u| ∂ u = μ|u| u − 2|u| u + u u, t t t see [16]. Hence, first-order convergence also requires the boundedness of at least two derivatives of the exact solution due to the appearance of u and its complex conjugate counterpart. The main novelty in this work lies in the development and analysis of a first-order exponential-type integrator for Schrödinger equations of type (1) e − 1 p+1 p n+1 iτ n n u = e u − i μτ u ϕ (−2iτ) u with ϕ (z) = (4) 1 1 which only requires the boundedness of one additional derivative of the exact solution. The construction is based on Duhamel’s formula looking at the twisted variable v(t ) = −it  2 p e u(t ) and treating the dominant term triggered by u in an exact way. This idea of twisting the variable is widely used in the analysis of partial differential equations in low regularity spaces (see for instance [3]) and also well known in the context of numer- ical analysis, see [24] for the introduction of Lawson-type Runge–Kutta methods. However, instead of approximating the appearing integrals with a Runge–Kutta method (see for instance [20]) we integrate the dominant second-order stiff parts exactly. A similar approach has been successfully applied to the one-dimensional KdV equa- tion, see [17]. However, due to the Burgers-type nonlinearity, additional smoothness assumptions on the exact solution are necessary. More precisely, first-order conver- 1 3 gence in H is only guaranteed for solutions in H in the KdV setting. The introduced exponential integrators for nonlinear Schrödinger equations are in contrast first-order r r +1 convergent in H for solutions in H (with r > d/2). Furthermore, we show 2 2 that our approach yields methods for quadratic nonlinearities of type u and |u| in r r one dimension which are convergent in H for solutions in H , i.e., no additional 123 734 Found Comput Math (2018) 18:731–755 smoothness of the solution is required. Such generous convergence results for non- linear Schrödinger equations are up to our knowledge not yet know in the numerical analysis literature. For practical implementation issues, we impose periodic boundary conditions and refer to [3] for local wellposedness (LWP) results of nonlinear Schrödinger equations r d in H (T ). Note that in the case of a cubic nonlinearity we obtain LWP for r > 1/2 if d ≤ 3. In the case of a quadratic nonlinearity of type u (or u ), LWP can be pushed down to r > −1/2 on the one-dimensional torus, see [21]. In the following let r r d r > d/2. We denote by · the standard H = H (T ) Sobolev norm. In particular, we exploit the well-known bilinear estimate fg ≤ c  f  g (5) r r,d r r which holds for some constant c > 0. Furthermore, we denote by c a generic r,d constant which may depend on r, d, p and μ. 2 Low Regularity Exponential Integrators We consider nonlinear Schrödinger equations (NLS) of type 2 p i ∂ u(t, x ) =−u(t, x ) + μ|u(t, x )| u(t, x ), 0 d u(0, x ) = u (x ), (t, x ) ∈ R × T , p ∈ N (6) d 2 with μ ∈ R and where we have set  = ∂ . Furthermore, we employ the k=1 x so-called ϕ function defined as e − 1 ϕ (z) = . In the following, we derive exponential-type integrators for Schrödinger equations of −it type (6) based on iterating Duhamel’s formula in the twisted variable v(t ) = e u(t ). Note that the twisted variable satisfies −it  it  2 p it  0 i ∂ v(t ) = μe |e v(t )| e v(t ) ,v(0) = v (7) with mild solution given by −i (t +s) i (t +s) 2 p i (t +s) n n n v(t + τ) = v(t ) − i μ e |e v(t + s)| e v(t + s) ds, n n n n (8) where t = nτ . In order to derive our numerical scheme, we proceed as follows. As it  r e is a linear isometry on H for all t ∈ R and the bilinear estimate (5) holds we obtain that 123 Found Comput Math (2018) 18:731–755 735 2 p+1 v(t + s) − v(t ) ≤|μ| v(t + ξ) dξ n n r n r 2 p+1 ≤ s|μ| sup v(t + ξ) , r > d/2. (9) n r 0≤ξ ≤s In this sense, we have for |s|≤ τ v(t + s) ≈ v(t ) (10) n n for a small time step τ . Remark 2.1 Note that the twisted variable allows the formal expansion 2 p v(t + s) − v(t ) = O s|v| v n n such that the approximation (10) holds without additional regularity assumptions (see (9) for the rigorous estimate). The classical exponential integrator (2), in contrast, is based on the approximation (3) for the original solution u. This approximation error is small only under additional smoothness assumptions, since is 2 p u(t + s) − u(t ) = e − 1 u(t ) + O s|u| u . n n n Plugging (10)into(8) yields the approximation −i (t +s) i (t +s) 2 p i (t +s) n n n v(t + τ) ≈ v(t ) − i μ e |e v(t )| e v(t ) ds (11) n n n n which is the basis of our numerical scheme. Hence, we are left with deriving a numerical approximation to the integral τ −i (t +s) i (t +s) 2 p i (t +s) n n n I (w, t ) := e |e w| e w ds. To illustrate the idea, we will first consider the cubic case p = 1 in Sect. 2.1 below. In Sect. 2.2 we will deal with general nonlinearities p ≥ 1. The special case of quadratic 2 2 nonlinearities of type u and |u| , respectively, will be treated in Sect. 4. 2.1 Cubic Nonlinearities p = 1 For notational simplicity, we first illustrate the idea in one dimension. In Sect. 2.1.2 we will give the generalization to arbitrary dimensions d ≥ 1. 123 736 Found Comput Math (2018) 18:731–755 2.1.1 Cubic Nonlinearities p = 1 in Dimension d = 1 2 ikx Let f ∈ L (T). Then, we will denote its Fourier expansion by f (x ) = f e . k∈Z −1 Furthermore, we define a regularization of ∂ through its action in Fourier space by −1 (ik) if k = 0 −1 −1 −1 ikx (∂ ) := , i.e., ∂ f (x ) = (ik) f e , k k x x 0if k = 0 k∈Z\{0} and, by continuity, 2 2 2 it ∂ −itk it ∂ x x e − 1 e − 1 e − 1 = , = 1. (12) 2 2 2 it ∂ −itk it ∂ x x k=0 k=0 In the case of a cubic nonlinearity [i.e., p = 1in (6)] in one spatial dimension, the integral in (11) can be expressed in terms of the Fourier expansion as follows 2 2 2 τ −i (t +s)∂ −i (t +s)∂ i (t +s)∂ n n n x x x I (w, t ) = e e w e w ds 2 2 2 2 i (t +s)  +k −k −k i x 1 2 3 = e ds w ˆ w ˆ w ˆ e k k k 1 2 3 ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 (13) 2 2 2 2 it  +k −k −k i x 1 2 3 = e w ˆ w ˆ w ˆ e k k k 1 2 3 ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 is 2k −2k (k +k )+2k k 1 2 3 2 3 e ds. Naturally, the next step would be the exact integration of the appearing integral is 2k −2k (k +k )+2k k 1 2 3 2 3 e ds. (14) However, the obtained convolution cannot be converted into the physical space straightforwardly and therefore only yields a practical scheme in dimension one as the iteration needs to be carried out fully in Fourier space. In order to obtain an efficient and practical low regularity implementation, we therefore only treat the dominant quadratic term 2k exactly. For this purpose, we write the integrand in (14)as iβ 2 e − 1 is 2k −2k (k +k )+2k k 1 2 3 2 3 2isk γ 1 1 (15) e = e 1 +|β| |β| 123 Found Comput Math (2018) 18:731–755 737 with β = s(−2k (k + k ) + 2k k and some 0 ≤ γ ≤ 1. Inserting (15)into(13) 1 2 3 2 3 and integrating the first term 2isk 2 e ds = τϕ (2i τ k ) yields that 2 2 2 2 τ it  +k −k −k 2 i x 1 2 3 I (w, t ) = τ e ϕ (2i τ k )w ˆ w ˆ w ˆ e n 1 k k k 1 1 1 2 3 ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 (16) + R (w, t ) 2 2 2 −it ∂ it ∂ 2 −it ∂ τ n n n x x x = τ e e w ϕ (−2iτ∂ )e w + R (w, t ), 1 n x 1 where 2 2 2 2 τ 2 2r it  +k −k −k 1 2 3 R (w, t ) = (1 +||) e 2k (k + k ) − 2k k n 1 2 3 2 3 1 r ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 ⎛ ⎞ is −2k (k +k )+2k k τ 1 2 3 2 3 2 e − 1 i x 2isk γ 1 ⎝ ⎠ w ˆ w ˆ w ˆ e e s ds . k k k 1 2 3 γ s |2k (k + k ) − 2k k 0 1 2 3 2 3 −γ iβ Since |β| (e − 1) is uniformly bounded for β ∈ R and 0 ≤ γ ≤ 1, the integral in 1+γ the above formula is of the order τ . This shows that ⎛ ⎞ ⎜   ⎟ τ 2 2+2γ 2r ⎜ ⎟ R (w, t ) ≤ cτ (1 +||) |k k |+|k k |+|k k | |ˆ w ||w ˆ ||w ˆ | n 1 2 1 3 2 3 k k k 1 r 1 2 3 ⎝ ⎠ ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 ⎛ ⎞ ⎜ ⎟ γ γ γ 2+2γ 2r ⎜ ⎟ ≤ cτ (1 +||) 1 +|k | 1 +|k | 1 +|k | |ˆ w ||w ˆ ||w ˆ | 1 2 3 k k k 1 2 3 ⎝ ⎠ ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 (17) for some constant c > 0 and 0 ≤ γ ≤ 1. ikx Next we define the auxiliary function g(x ) = g ˆ e through its Fourier k∈Z coefficients g ˆ = (1 +|k|) |ˆ w | (18) k k which allows us to express the bound on the remainder given in (17)asfollows τ 1+γ 3 R (w, t ) ≤ cτ g  . n r r 123 738 Found Comput Math (2018) 18:731–755 With the aid of the bilinear estimate (5) we can thus conclude that for r > 1/2 3/2 τ 1+γ 3 1+γ 2(r +γ) 2 R (w, t ) ≤ cτ g = cτ (1 +|k|) |ˆ w | , n r k 1 r k∈Z where we have used the definition of the Fourier coefficients in (18). Hence, we can conclude the following bound on the remainder τ 1+γ 3 R (w, t ) ≤ cτ w for 0 ≤ γ ≤ 1, r > 1/2 (19) n r r +γ for some constant c > 0. Next, we replace the integral in (11) with v in place of v(t ) by the term (16) τ n without remainder R and w = v . Twisting the solution back again, i.e., setting n it ∂ n u := e v yields the following exponential-type integration scheme for the cubic NLS n+1 iτ∂ n n 2 2 u = e u − iτμ(u ) ϕ (−2iτ∂ )u . (20) Next we generalize the above scheme to arbitrary dimensions d ≥ 1. 2.1.2 Cubic Nonlinearities p = 1 in Dimension d ≥ 1 In the following we use the notation j := ( j ,..., j ) ∈ Z , j · x := j x + ··· + j x . 1 d 1 1 d d In the case of a cubic nonlinearity p = 1 in dimension d, the integral in (11) can be expressed in terms of the Fourier expansion as follows τ −i (t +s) −i (t +s) i (t +s) n n n I (w, t ) = e e w e w ds 2 2 2 2 i (t +s)  +j −k −l n i ·x = e ds w ˆ w ˆ w ˆ e . j k l d d ∈Z j,k,l∈Z =−j+k+l Note that we have 2 2 2 2 2 + j − k − l = 2j − 2j · (k + l) + 2k · l. In order to obtain an efficient and practical low regularity implementation, we again treat the dominant quadratic term 2j exactly. 123 Found Comput Math (2018) 18:731–755 739 This yields that 2 2 2 2 2 τ it  +j −k −l i ·x 2isj τ I (w, t ) = e w ˆ w ˆ w ˆ e e ds + R (w, t ) n j k l n 1 1 d d ∈Z j,k,l∈Z =−j+k+l 2 2 2 2 it  +j −k −l n i ·x 2 τ = e w ˆ w ˆ w ˆ e ϕ (2i τ j ) + R (w, t ), j k l 1 n d d ∈Z j,k,l∈Z =−j+k+l (21) where similarly to above [see (19)] we have that τ 1+γ 3 R (w, t ) ≤ cτ w for 0 ≤ γ ≤ 1, r > d/2 n r 1 r +γ for some constant c > 0. Inserting the approximation (21)into(11) and twisting the solution back again, i.e., setting n it  n u := e v yields the generalization of the exponential-type integrator (20) to dimensions d ≥ 1 n+1 iτ n n 2 u = e u − iτμ(u ) ϕ (−2iτ)u . Next we consider the general case of (6) with p ∈ N. 2.2 Nonlinearities with Integer p In the case of a general nonlinearity p ∈ N in (6), the integral in (11) can be expressed in terms of the Fourier expansion as follows τ −i (t +s) i (t +s) 2 p i (t +s) n n n I (w, t ) = e |e w| e w ds. Similar observations to above (see also [28]) yield the following numerical scheme p+1 p n+1 n −it  it  n −it n n n n v = v − i μτ e e v ϕ (−2iτ) e v (22) which satisfies n+1 n τ n τ n v = v − i μI (v , t ) − i μR (v , t ), (23) n n p p with 2 p+1 τ n 1+γ n R (v , t ) ≤ cτ v  for 0 ≤ γ ≤ 1, r > d/2(24) n r p r +γ for some constant c > 0. 123 740 Found Comput Math (2018) 18:731–755 In the original variable u, the exponential-type integration scheme for the nonlinear Schrödinger equation (6) reads p+1 p n+1 iτ n n u = e u − i μτ u ϕ (−2iτ) u . (25) In the following section, we give an error analysis for the above scheme. 3 Error Analysis In this section, we carry out the error analysis of the exponential-type integrator (22). In the following, we set p+1 p τ −it  it  −it (w) := w − i μτ e e w ϕ (−2iτ) e w k+1 τ k such that in particular v = (v ). Lemma 3.1 (Stability) Let r > d/2 and f, g ∈ H . Then, for all t ∈ R we have τ τ ( f ) − (g) ≤ (1 + τ |μ|L) f − g , r r t t where L depends on  f  and g . r r Proof Note that for all t ∈ R and w ∈ H it holds by (12) that ϕ (2it )w ≤ cw . 1 r it  r r Thus, as e is a linear isometry on H and H is an algebra for r > d/2 we obtain that p+1 p τ τ it  −it ( f ) − (g) ≤ f − g + τ |μ| e f ϕ (−2iτ) e f r r 1 t t p+1 p it  −it − e g ϕ (−2iτ) e g 2 p 2 p−l ≤ f − g + τ c|μ| f − g  f  g . r r r l=0 This proves the assertion. Lemma 3.2 (Local error) Let r > d/2 and 0 ≤ γ ≤ 1. Assume that v(t + t ) = t r +γ t φ (v(t )) ∈ H for 0 ≤ t ≤ τ , where φ denotes the exact flow of (7). Then, τ τ 1+γ φ (v(t )) − (v(t )) ≤ cτ , k k r where c depends on sup φ (v(t )) . k r +γ 0≤t ≤τ 123 Found Comput Math (2018) 18:731–755 741 Proof Note that by (23)wehavethat τ −i (t +s) i (t +s) 2 p i (t +s) k k k (v(t )) = v(t ) − i μ e |e v(t )| e v(t ) ds k k k k + i μR (v(t ), t ). k k Furthermore, (8) implies that τ −i (t +s) i (t +s) 2 p i (t +s) k k k φ (v(t )) = v(t ) − i μ e |e v(t + s)| e v(t + s) ds. k k k k it  r r Thus as e is a linear isometry on H and H is an algebra for r > d/2 we obtain that 2 p τ τ φ (v(t )) − (v(t )) ≤|μ| sup v(t + s) v(t + s) − v(t ) ds k k r k r k k r 0≤s≤τ 0 +|μ|R (v(t ), t ) . k k r Together with the estimate on the difference v(t + s) − v(t ) given in (9) and the k k estimate on the remainder R (v, t ) given in (24) this yields for r > d/2 and 0 ≤ γ ≤ that 4 p+1 2 p+1 τ τ 2 2 1+γ φ (v(t )) − (v(t )) ≤ τ c|μ| sup v(t + s) + τ c|μ|v(t ) k k r k r k t r +γ 0≤s≤τ which concludes the stated local error bound. These two lemmata allow us to prove the following convergence bound. Theorem 3.3 Let r > d/2 and 0 <γ ≤ 1. Assume that the exact solution of (7) r +γ satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (22) 0 n is bounded by n γ v(t ) − v  ≤ cτ , n r where c depends on sup v(t ) . r +γ 0≤t ≤T Proof The triangle inequality yields that k+1 τ τ k v(t ) − v  =φ (v(t )) − (v ) k+1 r k r τ τ τ τ k ≤φ (v(t )) − (v(t )) + (v(t )) − (v ) . k k r k r t t t k k k k r Thus, we obtain with the aid of Lemmata 3.2 and 3.1 that as long as v ∈ H we have that k+1 1+γ k (26) v(t ) − v  ≤ cτ + (1 + τ |μ|L) v(t ) − v  , k+1 r k r 123 742 Found Comput Math (2018) 18:731–755 where c depends on sup v(t ) and L depends on v(t ) as well as on r +γ k r t ≤t ≤t k k+1 k r +γ k r v  . As long as v(t ) ∈ H and v ∈ H for 0 ≤ k ≤ n we therefore obtain by r k iterating the estimate (26) that n+1 1+γ n v(t ) − v  ≤ cτ + (1 + τ |μ|L) v(t ) − v n+1 r n r 1+γ τ |μ|L 1+γ n−1 ≤ cτ + e cτ + (1 + τ |μ|L) v(t ) − v n−1 r 1+γ nτ |μ|L γ t |μ|L ≤ cnτ e ≤ ct τ e . The assertion then follows by induction, respectively, a Lady Windermere’s fan argu- ment (see, for example [11,15,25]). The above theorem immediately yields a convergence result for the exponential- type integration scheme (22) approximating the solution of the nonlinear Schrödinger equation (6). Corollary 3.4 Let r > d/2 and 0 <γ ≤ 1. Assume that the exact solution of (6) r +γ satisfies u(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (25) 0 n is bounded by n γ u(t ) − u  ≤ cτ , n r where c depends on sup u(t ) . r +γ 0≤t ≤T it  r Proof As e is a linear isometry on H for all t ∈ R we have that n it  n n u(t ) − u  =e v(t ) − v  =v(t ) − v  . n r n r n r Thus, the assertion follows from Theorem 3.3. 4 Quadratic Schrödinger Equations in One Space Dimension In this section, we focus on quadratic nonlinear Schrödinger equations of type (see, e.g., [1]) 2 2 0 i ∂ u =−∂ u + μu , u(0, x ) = u (x ), (t, x ) ∈ R × T,μ ∈ R (27) as well as of type (see, e.g, [21]) 2 2 0 i ∂ u =−∂ u + μ|u| , u(0, x ) = u (x ), (t, x ) ∈ R × T,μ ∈ R, (28) which have recently gained a lot of attention in the literature, see for instance [1,3,5, 13,21–23,27,29,32] and the references therein. Here, the key relations 2 2 2 (k + k ) − k − k = 2k k (29) 1 2 1 2 1 2 123 Found Comput Math (2018) 18:731–755 743 and 2 2 2 (k − k ) − k + k =−2k (k − k ) (30) 1 2 2 1 2 1 2 allow us to derive first-order exponential-type integrators for Schrödinger equations of type (27) and (28) without imposing any additional regularity assumptions on the solution. 4.1 Quadratic Schrödinger Equations of Type (27) −it ∂ Set v(t ) = e u(t ), where u is the solution of the Schrödinger equation (27). Then the twisted variable v(t ) satisfies 2 2 −it ∂ it ∂ x x i ∂ v = μe e v(t ) (31) which by Duhamel’s formula yields that 2 2 −i (t +s)∂ i (t +s)∂ n n x x v(t + τ) = v(t ) − i μ e e v(t + s) ds. (32) n n n To construct our numerical scheme, we proceed as above. We start from the following approximation 2 2 −i (t +s)∂ i (t +s)∂ τ n n x x v(t + τ) ≈ v(t ) − i μ e e v(t ) ds =: (v(t )) (33) n n n n and rewrite the above integral in Fourier space as follows 2 2 τ −i (t +s)∂ i (t +s)∂ n n x x I (w, t ) = e e w ds 2 2 2 −i (t +s)∂ i (t +s)∂ ik x i (t +s)∂ ik x n n 1 n 2 x x x = e e w ˆ e e w ˆ e ds k k 1 2 k k 1 2 2 2 2 i (t +s) (k +k ) −k −k n 1 2 i (k +k )x 1 2 1 2 = e w ˆ w ˆ e ds. k k 1 2 k ,k 1 2 The key relation (29) now shows that 2 2 2 i τ (k +k ) −k −k 1 2 1 2 2 2 2 e − 1 it (k +k ) −k −k τ n 1 2 i (k +k )x 1 2 1 2 I (w, t ) = e w ˆ w ˆ e n k k 1 2 2ik k 1 2 k ,k 1 2 k ,k =0 1 2 ik x 2 + 2τ w ˆ w ˆ e + τ w ˆ 0 k 1 0 k =0 123 744 Found Comput Math (2018) 18:731–755 2 2 2 i τ (k +k ) −k −k 1 2 1 2 2 2 2 i e − 1 it (k +k ) −k −k i (k +k )x n 1 2 1 2 1 2 = e w ˆ w ˆ e k k 1 2 2 (ik )(ik ) 1 2 k ,k 1 2 k ,k =0 1 2 ik x 2 + 2τ w ˆ w ˆ e − τ w ˆ 0 k 1 0 k ∈Z 2 2 i 2 2 2 2 −it ∂ −iτ∂ i (t +τ)∂ −1 it ∂ −1 n n n x x x x = e e e ∂ w − e ∂ w x x + 2τ w ˆ w − τ w ˆ . Together with the approximation in (33) this yields the following integration scheme: n+1 τ n v = (v ) n n n 2 = 1 − 2i μτ v ˆ v + i μτ (v ˆ ) 0 0 (34) 2 2 μ 2 2 2 2 −it ∂ −iτ∂ i (t +τ)∂ −1 n it ∂ −1 n n n n x x x x + e e e ∂ v − e ∂ v . x x In order to obtain an approximation to the original solution u(t ) of the quadratic Schrödinger equation (27) at time t = nτ , we twist the variable back again, i.e., we n it ∂ n set u := e v . This yields the following exponential-type integrator 2 2 2 μ 2 μ 2 n+1 n iτ∂ n n 2 iτ∂ −1 n iτ∂ −1 n x x x u = 1 − 2i μτ u ˆ e u + i μτ (u ˆ ) + e ∂ u − e ∂ u . x x 0 0 2 2 (35) 4.2 Quadratic Schrödinger Equations of Type (28) For the quadratic Schrödinger equation (28), the equation in the twisted variable v(t ) = −it ∂ e u(t ) reads 2 2 −it ∂ it ∂ x x i ∂ v = μe e v(t ) . To construct our numerical scheme, we proceed as above. We start from the following approximation 2 2 −i (t +s)∂ i (t +s)∂ n n x x v(t + τ) ≈ v(t ) − i μ e e v(t ) ds (36) n n n and rewrite the above integral in Fourier space as follows 2 2 τ −i (t +s)∂ i (t +s)∂ n n x x I (w, t ) = e e w ds 2 2 2 i (t +s) (k −k ) −k +k n 1 2 i (k −k )x 1 2 1 2 = e w ˆ w ˆ e ds. k k 1 2 k ,k 1 2 123 Found Comput Math (2018) 18:731–755 745 The key relation (30) now shows that 2 2 2 i τ (k −k ) −k +k 1 2 1 2 2 2 2 e − 1 it (k −k ) −k +k τ n 1 2 i (k −k )x 1 2 1 2 I (w, t ) = e w ˆ w ˆ e n k k 1 2 −2ik (k − k ) 2 1 2 k ,k 1 2 k =k ,k =0 1 2 2 ik x 2 +τ w ˆ w ˆ e + τ w 0 k 2 i 2 2 2 −1 −i (t +τ)∂ i (t +τ)∂ −i (t +τ)∂ −1 n n n x x x = ∂ e e w e ∂ w x x 2 2 2 −1 −it ∂ it ∂ −it ∂ −1 2 n n n x x x − ∂ e e w e ∂ w + τ w ˆ w + τ w . 0 2 x x Together with the approximation in (36) this yields (by twisting the variable back it ∂ again) the following integration scheme in the original solution u(t ) = e v(t ) n 2 n+1 iτ∂ n n 2 u = 1 − i μτ u ˆ e u − i μτ u (37) μ 2 2 2 −1 iτ∂ n −iτ∂ −1 n iτ∂ n −1 n x x x + ∂ e u e ∂ u − e u ∂ u . x x x 4.3 Error Analysis In this section, we carry out the error analysis of the exponential-type integrators (35) and (37). In the following let r > 1/2. We commence with the quadratic Schrödinger equation of type (27). Let φ denote the exact flow of (31), i.e., v(t ) = φ (v(0)). The following lemma provides a local error bound of the scheme (34). t r Lemma 4.1 Let r > 1/2 and assume that v(t + t ) = φ (v(t )) ∈ H for 0 ≤ t ≤ τ , k k where φ denotes the exact flow of (27). Then τ τ 2 φ (v(t )) − (v(t )) ≤ cτ , k k r where c depends on sup φ (v(t )) . k r 0≤t ≤τ it ∂ r Proof As e is a linear isometry on H for all t ∈ R we obtain by taking the difference of (32) with the approximation (33) and using the bilinear estimate (5) that 2 2 2 2 τ τ i (t +s)∂ i (t +s)∂ k k x x φ (v(t )) − (v(t )) ≤|μ|  e v(t + s) − e v(t )  ds k k r k k 2 2 i (t +s)∂ i (t +s)∂ k k x x ≤ c|μ|τ sup e (v(t + s) − v(t ))  e (v(t + s) + v(t )) k k r k k r 0≤s≤τ ≤ c|μ|τ sup v(t + s) − v(t ) k k r 0≤s≤τ 2 i (t +ξ)∂ ≤ c|μ| τ sup  e v(t + ξ)  dξ, 0≤s≤τ 0 where c depends on sup φ (v(t )) . This yields the assertion. k r 0≤s≤τ 123 746 Found Comput Math (2018) 18:731–755 Next we state the stability estimate. Lemma 4.2 Let r > 1/2 and f, g ∈ H . Then, for all t ∈ R we have τ τ ( f ) − (g) ≤ (1 + τ |μ|L)  f − g , r r t t where L depends on  f + g . Proof The assertion follows from the representation of the numerical flow given in (33) together with the bilinear estimate (5). The above lemmata allow us to prove the following global convergence result. Theorem 4.3 Let r > 1/2. Assume that the exact solution of (31) satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (34) is bounded by 0 n v(t ) − v  ≤ cτ, n r where c depends on sup v(t ) . 0≤t ≤T Proof The proof follows the line of argumentation to the proof of Theorem 3.3 thanks to the local error bound in Lemma 4.1 and the stability estimate of Lemma 4.2. The above theorem immediately yields first-order convergence of the exponential- type integrator (35)(in H ) without any additional smoothness assumptions on the exact solution of (27). Corollary 4.4 Let r > 1/2. Assume that the exact solution of (27) satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (35) is bounded by 0 n u(t ) − u  ≤ cτ, n r where c depends on sup u(t ) . 0≤t ≤T it ∂ Proof The bound follows directly from Theorem 4.3 as e is a linear isometry on H for all t ∈ R. Similarly to above we obtain the following first-order convergence result for the exponential-type integrator (37) approximating the solution of the quadratic Schrödinger equation (28). Corollary 4.5 Let r > 1/2. Assume that the exact solution of (28) satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (37) is bounded by 0 n u(t ) − u  ≤ cτ, n r where c depends on sup u(t ) . 0≤t ≤T Proof The proof is similar to that of Corollary 4.4 and will therefore be omitted here. 123 Found Comput Math (2018) 18:731–755 747 5 Numerical Experiments In this section, we compare the newly developed exponential-type integrator (25) applied to Schrödinger equations with low regularity solutions in the energy space with the classical Lie splitting, the classical Strang splitting and the classical first- order exponential integrator. In the numerical experiments, we use a standard Fourier pseudospectral method for the space discretization where we choose the largest Fourier mode K = 2 (i.e., the spatial mesh size x = 0.0061). The numerical experiments show that the new integrator (25) is preferable over the classical splitting and expo- nential integration methods for low regularity solutions. In particular, the experiments indicate that the proposed integrator (25) is convergent of order one even for rough solutions whereas the classical methods suffer from order reduction. Although the new integrator (25) was derived for nonlinearities with integer p ≥ 1, it can be applied to problem (6) with non-integer p as well. Preliminary numerical results for small non-integer values of p are presented in Sect. 5.3. 5.1 Numerical Experiments for Schrödinger Equations of Type (6) In this section, we consider nonlinear Schrödinger equations of type (6). The associated Lie splitting (of classical order one) reads n+1/2 iτ n u = e u , (38) 2 p n+1/2 n+1 −iτμ u n+1/2 u = e u , L L and the associated Strang splitting (of classical order two) is given by n+1/2 i  n u = e u , 2 p n+1/2 n+1/2 −iτμ u n+1/2 − (39) u = e u , + − n+1 n+1/2 u = e u . Furthermore, the classical first-order exponential integrator (of classical order one)is defined through n+1 iτ n n 2 p n u = e u − i μτ ϕ (iτ) |u | u . (40) E E E In our numerical experiments, we choose μ = 1 (unless otherwise noted) and set K 2K U := [u(0, x ), u(0, x ), . . . , u(0, x )]= rand(2K , 1) + i rand(2K , 1) ∈ C , 1 2 2K (41) where x = (−1 + )π and rand(2K , 1) return 2K uniformly distributed random numbers between 0 and 1. We compare the above schemes with the newly derived exponential-type integrator (25) for initial values 123 748 Found Comput Math (2018) 18:731–755 0.17 0.17 0.16 0.16 0.15 0.15 Fig. 1 Initial values (42) normalized in L for two different values of ϑ. Left ϑ = . Right ϑ = 5 -2 -2 -3 -3 -4 -4 -5 -2 -1 -2 -1 10 10 10 10 Fig. 2 (cubic NLS) Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39) (magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type integration 3/2 2 scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. The slope of the continuous lines is one half and one, respectively (Color figure online) −ϑ |k| if k = 0, K −ϑ K −ϑ u =|∂ | U , |∂ | := (42) x ,K x ,K 0if k = 0 for different values of ϑ normalized in L . (For typical initial values, see Fig. 1.) All experiments are carried out with constant step sizes τ = for j = 1,..., 512. 5.1.1 Cubic Schrödinger Equations In this section, we consider the classical cubic Schrödinger equation, i.e., p = 1in (6). The error at time T = 1 measured in a discrete H norm of the exponential-type integrator (25), the Lie splitting scheme (38), the Strang splitting scheme (39)aswell as the classical exponential integrator (40) for the initial value (42) normalized in L is illustrated in Fig. 2 for ϑ = 3/2, respectively, ϑ = 2 and in Fig. 3 for ϑ = 3, respectively, ϑ = 5. As a reference solution, we take the schemes themselves with a very small time step size if ϑ = 5. For ϑ = 5 we take the Strang splitting scheme (39) with a very small time step size as the reference solution for all schemes. The experiments show that Strang splitting is second-order convergent for ϑ ≥ 5, and 123 Found Comput Math (2018) 18:731–755 749 -2 -2 10 10 -3 -3 10 10 -4 -4 10 10 -5 -5 10 10 -6 -6 10 10 -7 -7 10 10 -8 -8 10 10 -9 -9 10 10 -2 -1 -2 -1 10 10 10 10 Fig. 3 (cubic NLS) Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39)(magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type inte- 3 5 gration scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. The slope of the continuous lines is one and two, respectively (Color figure online) that Lie splitting is first-order convergent for ϑ ≥ 3. These observations are in line with the convergence proofs from the literature [10,25], see also our discussion in the introduction. For smaller values of ϑ, both Lie and Strang splitting show a zigzag behavior, depending on whether the local errors accumulate or happy error cancelation occurs. The classical exponential integrator of order one is first-order convergent for ϑ ≥ 3 as expected (see the discussion in the introduction). For smaller values of ϑ, the convergence behavior of the exponential integrator gets less regular and order reduction occurs. On the other hand, the new exponential-type integrator introduced in this paper is first-order convergent for all ϑ ≥ 2. This behavior is in line with Theorem 3.3. The integrator even shows first-order convergence and small errors for ϑ = . 5.1.2 Quintic Schrödinger Equations In this section, we consider the quintic Schrödinger equation, i.e., p = 2in (6), which appears for instance as the mean field limit of a Boson gas with three-body interactions, see [7]. The error at time T = 1 measured in a discrete H norm of the exponential- type integrator (25), the Lie splitting scheme (38), the Strang splitting scheme (39)as well as the classical exponential integrator (40) for the initial value (42) normalized in L is illustrated in Fig. 4 for ϑ = 3/2, respectively, ϑ = 2 and in Fig. 5 for ϑ = 3, respectively, ϑ = 5. As a reference solution we take the schemes themselves with a very small time step size if ϑ = 5. For ϑ = 5 we take the Strang splitting scheme (39) with a very small time step size as the reference solution for all schemes. The outcome of these numerical experiments is almost the same as for the cubic nonlinear Schrödinger equation. The only notable difference is the exponential-type integrator for ϑ = , where the convergence behavior is now less regular. Note, however, that the errors are still much smaller compared to the other methods. 123 750 Found Comput Math (2018) 18:731–755 -2 -3 -3 -4 -4 -5 -5 -6 -6 10 10 -2 -1 -2 -1 10 10 10 10 Fig. 4 (quintic NLS) Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39)(magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type inte- 3/2 2 gration scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. Left picture The slope of the continuous lines is one quarter and one, respectively. Right picture The slope of the continuous lines is one half and one, respectively (Color figure online) -4 -4 -5 -5 10 -6 -6 -7 -7 -8 -8 -9 -9 -2 -1 -2 -1 10 10 10 10 Fig. 5 (quintic NLS) in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39) (magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type integra- 3 5 tion scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. The slope of the continuous lines is one and two, respectively (Color figure online) 5.2 Numerical Experiments for Quadratic Schrödinger Equations In this section, we consider quadratic Schrödinger equations of type (27). In order to derive the splitting methods, we split the right-hand side of (27) into the linear and nonlinear part, respectively. Note that (for sufficiently small t) the exact solution of the subproblem iu (t ) = μu(t ) is given by u(0) u(t ) = . 1 + i μtu(0) 123 Found Comput Math (2018) 18:731–755 751 Thus, the associated Lie splitting scheme (of classical order one) reads 2 u n+1 iτ∂ L u = e (43) L n 1 + i μτ u and the associated Strang splitting scheme (of classical order two) is given by τ 2 n+1/2 i ∂ n u = e u S S n+1/2 (44) τ 2 u n+1 i ∂ S u = e . n+1/2 1 + i μτ u In Examples 5.1 and 5.2, we compare the above splitting methods as well as the classical first-order exponential integration scheme n+1 iτ∂ n 2 n u = e u − i μτ ϕ (iτ∂ ) u (45) E x E with the newly derived exponential-type integrator (35)for smooth and non-smooth solutions, respectively, where we set μ = 1 and integrate up to T = 1. The numerical experiments in particular indicate that the exponential-type integrator (35)ispreferable over the commonly used splitting methods for low regularity solutions. This is due to the fact that within this exponential integration scheme all the stiff parts (i.e., the terms involving the differential operator ∂ ) are solved exactly. In Example 5.3 we consider a quadratic Schrödinger equation of type (27) with a small nonlinearity on a longer time interval T = 10. The numerical experiment indicates the favorable behavior of the exponential integrator (35) even in the smooth case. This is due to the fact that its error constant is triggered by the nonlinearity (and not by the differential operator ∂ ). Example 5.1 (Smooth solutions) We choose the smooth initial value u(0, x ) = sin x cos x ∈ C (T) (46) 2 2 normalized in L . The error at time T = 1 measured in a discrete L norm of the exponential-type integrator (35), the Lie splitting scheme (43), the Strang splitting scheme (44) and the classical exponential integrator (45) is illustrated in Fig. 6.As a reference solution, we take the Strang splitting scheme (44) with a very small time step size. All the methods show their classical orders of convergence, as expected. Example 5.2 (Non-smooth solutions) We choose the non-smooth initial value (41) 2 2 normalized in L . The error at time T = 1 measured in a discrete L norm of the exponential-type integrator (35), the Lie splitting scheme (43), the Strang splitting scheme (44) and the classical exponential integrator (45) is illustrated in Fig. 7.Asa reference solution, we take the schemes themselves with a very small time step size. For this example, the classical exponential integrator is not convergent. Lie and Strang splitting suffer from a strong order reduction (down to order one half). Here, the new 123 752 Found Comput Math (2018) 18:731–755 Fig. 6 (quadratic NLS) Error of -1 the Lie splitting (43)(blue, dotted), Strang splitting (44) -2 (magenta, solid line), classical exponential integrator (45) -3 (green, dashed)and exponential-type integration scheme (35)(red, squares)for a -4 smooth solution with initial value (46)and μ = 1. The slope -5 of the continuous lines is one and two, respectively (Color figure online) -6 -2 -1 10 10 Fig. 7 (quadratic NLS) Error of 0 the Lie splitting (43)(blue, dotted), Strang splitting (44) (magenta, solid line), classical exponential integrator (45) -1 (green, dashed)and exponential-type integration scheme (35)(red, squares)for a non-smooth solution with initial -2 value (41)and μ = 1. The slope of the continuous lines is one half and one, respectively (Color figure online) -3 -2 -1 10 10 exponential-type integrator is by far the best method. It is first-order convergent, as predicted by Theorem 4.3. Example 5.3 (Small nonlinearity) We consider the quadratic Schrödinger equation of type (27) with small nonlinearity, i.e., 2 2 i ∂ u =−∂ u + μu ,μ = 0.01, T = 10. (47) The error at time T = 10 measured in a discrete L norm of the exponential-type integrator (35), the Lie splitting scheme (43), the Strang splitting scheme (44) and the classical exponential integrator (45) for smooth and non-smooth solutions is illustrated in Fig. 8 left and right, respectively. As a reference solution, we take the Strang splitting scheme for smooth and the schemes themselves for non-smooth solutions with a very small time step size. For smooth initial data, all methods show their classical orders of convergence. The new exponential-type integrator, however, is by far the most accurate scheme (for the considered range of step sizes). For non-smooth initial data, the classical exponential integrator does not converge; Lie and Strang splitting suffer 123 Found Comput Math (2018) 18:731–755 753 -3 -2 -4 -3 -5 10 -4 -6 10 -5 -7 -6 -8 -2 -1 -2 -1 10 10 10 10 Fig. 8 (quadratic NLS) Error of the Lie splitting (43)(blue, dotted), Strang splitting (44)(magenta, solid line), classical exponential integrator (45)(green, dashed) and exponential-type integration scheme (35) (red, squares)fora small nonlinearity (47). Left picture Smooth initial value (46). The slope of the continuous lines is one and two, respectively. Right picture Non-smooth initial value (41). The slope of the continuous lines is one half, and one, respectively (Color figure online) 0 0 10 10 -1 -1 -2 -2 -3 -3 -4 10 -4 -5 -5 -6 -2 -1 -2 -1 10 10 10 10 0 0 10 10 -1 10 -1 -2 -2 -3 -3 -4 -4 -2 -1 -2 -1 10 10 10 10 3 1 1 Fig. 9 Exponent p = 1(upper row, left)and p = (upper row, right), p = (lower row, left)and p = 4 2 4 (lower row, right). Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39) (magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type integration scheme (25)(red, dash dotted) with smooth initial value (48). The slope of the continuous and dotted line is one and two, respectively (Color figure online) 123 754 Found Comput Math (2018) 18:731–755 again from a strong order reduction down to order one half. Only the new exponential- type integrator is very accurate and first-order convergent, as expected. 5.3 Numerical Experiments for Semilinear Schrödinger Equations with Non-integer Exponents p In this section, we numerically investigate the convergence rates in the case of semi- linear Schrödinger equations of type (6) with non-integer exponents p.InFig. 9 we plot the error at time T = 1 measured in a discrete H norm of the exponential-type integrator (25), the Lie splitting scheme (38), the Strang splitting scheme (39)aswell as the classical exponential integrator (40) for different values of p > 0. In all the simulations, we take the smooth initial value u(0, x ) = sin x ∈ C (T) (48) and μ = 1. Note, however, that the nonlinearity is only smooth for p > . For p = 1, all methods show their classical orders of convergence, as expected. For smaller values of p, however, the convergence behavior of Strang and Lie splitting deteriorates. Surprisingly, both exponential integrators retain their first-order conver- gence down to very small values of p. This is the subject of future research. Acknowledgements Open access funding provided by University of Innsbruck and Medical University of Innsbruck. K. Schratz gratefully acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG) through CRC 1173. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. I. Bejenaru, T. Tao, Sharp well-posedness and ill-posedness results for a quadratic non-linear Schrödinger equation. J. Funct. Anal. 233:228–259 (2006). 2. C. Besse, B. Bidégaray, S. Descombes, Order estimates in time of splitting methods for the nonlinear Schrödinger equation, SIAM J. Numer. Anal. 40:26–40 (2002). 3. J. Bourgain, Fourier transform restriction phenomena for certain lattice subsets and applications to nonlinear evolution equations. Part I: Schrödinger equations. Geom. Funct. Anal. 3:209–262 (1993). 4. B. Cano, A. González-Pachón, Exponential time integration of solitary waves of cubic Schrödinger equation. Appl. Numer. Math. 91:26–45 (2015). 5. T. Cazenave, F.B. Weissler, The Cauchy problem for the critical nonlinear Schrödinger equation. Nonlinear Anal., 14:807–836 (1990). 6. E. Celledoni, D. Cohen, B. Owren, Symmetric exponential integrators with an application to the cubic Schrödinger equation. Found. Comput. Math. 8:303–317 (2008). 7. T.Chen,N.Pavlovic, ´ The quintic NLS as the mean field limit of a boson gas with three-body interactions. J. Funct. Anal. 260:959–997 (2011). 8. D. Cohen, L. Gauckler, One-stage exponential integrators for nonlinear Schrödinger equations over long times. BIT 52:877–903 (2012). 9. G. Dujardin, Exponential Runge-Kutta methods for the Schrödinger equation. Appl. Numer. Math. 59:1839–1857 (2009). 123 Found Comput Math (2018) 18:731–755 755 10. J. Eilinghoff, R. Schnaubelt, K. Schratz, Fractional error estimates of splitting schemes for the non- linear Schrödinger equation. J. Math. Anal. Appl. 442:740–760 (2016). 11. E. Faou, Geometric Numerical Integration and Schrödinger Equations. European Math. Soc. Publish- ing House, Zürich 2012. 12. L. Gauckler, Convergence of a split-step Hermite method for the Gross–Pitaevskii equation. IMA J. Numer. Anal. 31:396–415 (2011). 13. P. Germain, N. Masmoudi, J. Shatah, Global solutions for 3D quadratic Schrödinger equations. Int. Math. Res. Notices 3:414–432 (2009). 14. E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Second edition, Springer, Berlin 2006. 15. E. Hairer, S. P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems. Second edition. Springer, Berlin 1993. 16. M. Hochbruck, A. Ostermann, Exponential integrators. Acta Numer. 19:209–286 (2010). 17. M. Hofmanova, K. Schratz, An exponential-type integrator for the KdV equation. Numer. Math. (2016). doi:10.1007/s00211-016-0859-1 18. H. Holden, K. H. Karlsen, K.-A. Lie, N. H. Risebro, Splitting for Partial Differential Equations with Rough Solutions. European Math. Soc. Publishing House, Zürich 2010. 19. L. I. Ignat, A splitting method for the nonlinear Schrödinger equation. J. Differential Equations 250:3022–3046 (2011). 20. A.-K. Kassam, L. N. Trefethen, Fourth-order time-stepping for stiff PDEs. SIAM J. Sci. Comput. 26:1214–1233 (2005). 21. C. Kenig, G. Ponce, L. Vega, Quadratic forms for the 1-D semilinear Schrödinger equation. Trans. Amer. Math. Soc. 346:3323–3353 (1996). 22. C. Kenig, G. Ponce, L. Vega, On the ill-posedness of some canonical dispersive equations. Duke Math. J. 106:617–633 (2001). 23. N. Kishimoto, Low-regularity bilinear estimates for a quadratic nonlinear Schrödinger equation. J. Differential Equations 247:1397–1439 (2009). 24. J. D. Lawson, Generalized Runge–Kutta processes for stable systems with large Lipschitz constants. SIAM J. Numer. Anal. 4:372–380 (1967). 25. C. Lubich, On splitting methods for Schrödinger–Poisson and cubic nonlinear Schrödinger equations. Math. Comp. 77:2141–2153 (2008). 26. R.I. McLachlan, G.R.W. Quispel, Splitting methods, Acta Numer. 11:341–434 (2002). 27. K. Nakanishi, H. Takaoka, Y. Tsutsumi, Counterexamples to bilinear estimates related with the KdV equation and the nonlinear Schrödinger equation. Methods Appl. Anal. 8:569–578 (2001). 28. K. Schratz (joint with A. Ostermann), Derivation of a low regularity exponential-type integrator for semilinear Schrödinger equations with polynomial nonlinearities. Oberwolfach Reports 18:928–931 (2016). 29. T. Tao, Multilinear weighted convolution of L functions, and applications to nonlinear dispersive equations. Amer. J. Math, 123:839–908 (2001). 30. T. Tao, Nonlinear Dispersive Equations. Local and Global Analysis. Amer. Math. Soc., Providence 31. M. Thalhammer, Convergence analysis of high-order time-splitting pseudo-spectral methods for non- linear Schrödinger equations. SIAM J. Numer. Anal. 50:3231–3258 (2012). 32. Y. Tsutsumi, L solutions for nonlinear Schrödinger equations and nonlinear groups. Funkcial. Ekvac., 30:115–125 (1987). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Foundations of Computational Mathematics Springer Journals

Low Regularity Exponential-Type Integrators for Semilinear Schrödinger Equations

Free
25 pages

Loading next page...
 
/lp/springer_journal/low-regularity-exponential-type-integrators-for-semilinear-schr-dinger-gnyIxCDZtO
Publisher
Springer Journals
Copyright
Copyright © 2017 by The Author(s)
Subject
Mathematics; Numerical Analysis; Economics, general; Applications of Mathematics; Linear and Multilinear Algebras, Matrix Theory; Math Applications in Computer Science; Computer Science, general
ISSN
1615-3375
eISSN
1615-3383
D.O.I.
10.1007/s10208-017-9352-1
Publisher site
See Article on Publisher Site

Abstract

Found Comput Math (2018) 18:731–755 https://doi.org/10.1007/s10208-017-9352-1 Low Regularity Exponential-Type Integrators for Semilinear Schrödinger Equations 1 2 Alexander Ostermann · Katharina Schratz Received: 24 March 2016 / Revised: 22 March 2017 / Accepted: 3 May 2017 / Published online: 26 May 2017 © The Author(s) 2017. This article is an open access publication Abstract We introduce low regularity exponential-type integrators for nonlinear Schrödinger equations for which first-order convergence only requires the bound- edness of one additional derivative of the solution. More precisely, we will prove r r +1 first-order convergence in H for solutions in H (with r > d/2) of the derived schemes. This allows us lower regularity assumptions on the data than for instance required for classical splitting or exponential integration schemes. For one-dimensional quadratic Schrödinger equations, we can even prove first-order convergence without any loss of regularity. Numerical experiments underline the favorable error behav- ior of the newly introduced exponential-type integrators for low regularity solutions compared to classical splitting and exponential integration schemes. Keywords Nonlinear Schrödinger equations · Exponential-type time integrator · Low regularity · Convergence Mathematics Subject Classification 65M12 · 65M70 · 35Q41 Communicated by Arieh Iserles. Alexander Ostermann alexander.ostermann@uibk.ac.at Katharina Schratz katharina.schratz@kit.edu Department of Mathematics, University of Innsbruck, Technikerstr. 13, 6020 Innsbruck, Austria Fakultät für Mathematik, Karlsruhe Institute of Technology, Englerstr. 2, 76131 Karlsruhe, Germany 123 732 Found Comput Math (2018) 18:731–755 1 Introduction Semilinear Schrödinger equations, in particular those of type 2 p i ∂ u =−u + μ|u| u, p ∈ N (1) with μ ∈ R are nowadays extensively studied numerically. In this context, split- ting methods (where the right-hand side is split into the kinetic and nonlinear part, respectively) as well as exponential integrators (including Lawson-type Runge–Kutta methods) contribute particularly attractive classes of integration schemes. For an extensive overview on splitting and exponential integration methods, we refer to [14,16,18,26], and for their rigorous convergence analysis in the context of semi- linear Schrödinger equations, we refer to [2,4,6,8,9,11,12,25,31] and the references therein. However, within the construction of all these numerical methods the stiff part (i.e., the terms involving the differential operator ) is approximated in one way or another which generally requires the boundedness of two additional spatial derivatives of the exact solution. In particular, convergence of a certain order only holds under sufficient additional regularity assumptions on the solution. In the following, we will illustrate the local error behavior of classical splitting and exponential integration methods and the thereby introduced smoothness requirements. Splitting schemes The Strang splitting for the cubic Schrödinger equation (i.e., p = 1) n+1/2 i  n u = e u , n+1/2 n+1/2 −iτμ u n+1/2 u = e u , + − n+1/2 n+1 i u = e u , where the right-hand side is split into the kinetic T (u) = i u and nonlinear part V (u) =−i μ|u| u is rigorously analyzed in [25]. In particular, its second-order con- r r +4 r vergence in H for solutions in H and its first-order convergence in H for solutions r +2 in H holds for all r ≥ 0, see [25] and [10], respectively. The result follows from the fact that the local error of the splitting scheme can be expressed through the double Lie commutator [T , [T , V ]](u) and the Lie commutator [T , V ](u) for second- and first-order methods, respectively. The latter reads [T , V ](u) = (∇u ·∇u) u + (∇uu) ·∇u + (u∇u) ·∇u + uu u, 2μ see [25, Section 4.2]. Due to the appearance of u in the local error, the boundedness of at least two additional derivatives of the exact solution is required. Based on [25], fractional error estimates for the Strang splitting were established in [10] which require the boundedness of 2 + 2γ additional derivatives for convergence of order 1 + γ , with 0 <γ < 1. Furthermore, in [19] first-order convergence of a filtered Lie splitting method for Schrödinger equations of type (1) with 1 ≤ 2 p < 4/d was shown in 2 d 2 d L (R ) for solutions in H (R ). 123 Found Comput Math (2018) 18:731–755 733 Classical exponential integrators The classical first-order exponential integrator in the cubic case reads e − 1 n+1 iτ n n 2 n u = e u − i μτ ϕ (iτ) |u | u with ϕ (z) = . (2) 1 1 Its construction is based on Duhamel’s formula iτ i (τ −s) 2 p u(t + τ) = e u(t ) − i μ e |u(t + s)| u(t + s)ds n n n n by applying the approximation u(t + s) ≈ u(t ) (3) n n in the integral terms and solving the remaining integral over the free Schrödinger i (τ −s) group e exactly. Due to the approximation (3), the above scheme introduces a local error of the form 2 2 2 4 2 2 ∂ i |u| u = i ∂ uu + 2i |u| ∂ u = μ|u| u − 2|u| u + u u, t t t see [16]. Hence, first-order convergence also requires the boundedness of at least two derivatives of the exact solution due to the appearance of u and its complex conjugate counterpart. The main novelty in this work lies in the development and analysis of a first-order exponential-type integrator for Schrödinger equations of type (1) e − 1 p+1 p n+1 iτ n n u = e u − i μτ u ϕ (−2iτ) u with ϕ (z) = (4) 1 1 which only requires the boundedness of one additional derivative of the exact solution. The construction is based on Duhamel’s formula looking at the twisted variable v(t ) = −it  2 p e u(t ) and treating the dominant term triggered by u in an exact way. This idea of twisting the variable is widely used in the analysis of partial differential equations in low regularity spaces (see for instance [3]) and also well known in the context of numer- ical analysis, see [24] for the introduction of Lawson-type Runge–Kutta methods. However, instead of approximating the appearing integrals with a Runge–Kutta method (see for instance [20]) we integrate the dominant second-order stiff parts exactly. A similar approach has been successfully applied to the one-dimensional KdV equa- tion, see [17]. However, due to the Burgers-type nonlinearity, additional smoothness assumptions on the exact solution are necessary. More precisely, first-order conver- 1 3 gence in H is only guaranteed for solutions in H in the KdV setting. The introduced exponential integrators for nonlinear Schrödinger equations are in contrast first-order r r +1 convergent in H for solutions in H (with r > d/2). Furthermore, we show 2 2 that our approach yields methods for quadratic nonlinearities of type u and |u| in r r one dimension which are convergent in H for solutions in H , i.e., no additional 123 734 Found Comput Math (2018) 18:731–755 smoothness of the solution is required. Such generous convergence results for non- linear Schrödinger equations are up to our knowledge not yet know in the numerical analysis literature. For practical implementation issues, we impose periodic boundary conditions and refer to [3] for local wellposedness (LWP) results of nonlinear Schrödinger equations r d in H (T ). Note that in the case of a cubic nonlinearity we obtain LWP for r > 1/2 if d ≤ 3. In the case of a quadratic nonlinearity of type u (or u ), LWP can be pushed down to r > −1/2 on the one-dimensional torus, see [21]. In the following let r r d r > d/2. We denote by · the standard H = H (T ) Sobolev norm. In particular, we exploit the well-known bilinear estimate fg ≤ c  f  g (5) r r,d r r which holds for some constant c > 0. Furthermore, we denote by c a generic r,d constant which may depend on r, d, p and μ. 2 Low Regularity Exponential Integrators We consider nonlinear Schrödinger equations (NLS) of type 2 p i ∂ u(t, x ) =−u(t, x ) + μ|u(t, x )| u(t, x ), 0 d u(0, x ) = u (x ), (t, x ) ∈ R × T , p ∈ N (6) d 2 with μ ∈ R and where we have set  = ∂ . Furthermore, we employ the k=1 x so-called ϕ function defined as e − 1 ϕ (z) = . In the following, we derive exponential-type integrators for Schrödinger equations of −it type (6) based on iterating Duhamel’s formula in the twisted variable v(t ) = e u(t ). Note that the twisted variable satisfies −it  it  2 p it  0 i ∂ v(t ) = μe |e v(t )| e v(t ) ,v(0) = v (7) with mild solution given by −i (t +s) i (t +s) 2 p i (t +s) n n n v(t + τ) = v(t ) − i μ e |e v(t + s)| e v(t + s) ds, n n n n (8) where t = nτ . In order to derive our numerical scheme, we proceed as follows. As it  r e is a linear isometry on H for all t ∈ R and the bilinear estimate (5) holds we obtain that 123 Found Comput Math (2018) 18:731–755 735 2 p+1 v(t + s) − v(t ) ≤|μ| v(t + ξ) dξ n n r n r 2 p+1 ≤ s|μ| sup v(t + ξ) , r > d/2. (9) n r 0≤ξ ≤s In this sense, we have for |s|≤ τ v(t + s) ≈ v(t ) (10) n n for a small time step τ . Remark 2.1 Note that the twisted variable allows the formal expansion 2 p v(t + s) − v(t ) = O s|v| v n n such that the approximation (10) holds without additional regularity assumptions (see (9) for the rigorous estimate). The classical exponential integrator (2), in contrast, is based on the approximation (3) for the original solution u. This approximation error is small only under additional smoothness assumptions, since is 2 p u(t + s) − u(t ) = e − 1 u(t ) + O s|u| u . n n n Plugging (10)into(8) yields the approximation −i (t +s) i (t +s) 2 p i (t +s) n n n v(t + τ) ≈ v(t ) − i μ e |e v(t )| e v(t ) ds (11) n n n n which is the basis of our numerical scheme. Hence, we are left with deriving a numerical approximation to the integral τ −i (t +s) i (t +s) 2 p i (t +s) n n n I (w, t ) := e |e w| e w ds. To illustrate the idea, we will first consider the cubic case p = 1 in Sect. 2.1 below. In Sect. 2.2 we will deal with general nonlinearities p ≥ 1. The special case of quadratic 2 2 nonlinearities of type u and |u| , respectively, will be treated in Sect. 4. 2.1 Cubic Nonlinearities p = 1 For notational simplicity, we first illustrate the idea in one dimension. In Sect. 2.1.2 we will give the generalization to arbitrary dimensions d ≥ 1. 123 736 Found Comput Math (2018) 18:731–755 2.1.1 Cubic Nonlinearities p = 1 in Dimension d = 1 2 ikx Let f ∈ L (T). Then, we will denote its Fourier expansion by f (x ) = f e . k∈Z −1 Furthermore, we define a regularization of ∂ through its action in Fourier space by −1 (ik) if k = 0 −1 −1 −1 ikx (∂ ) := , i.e., ∂ f (x ) = (ik) f e , k k x x 0if k = 0 k∈Z\{0} and, by continuity, 2 2 2 it ∂ −itk it ∂ x x e − 1 e − 1 e − 1 = , = 1. (12) 2 2 2 it ∂ −itk it ∂ x x k=0 k=0 In the case of a cubic nonlinearity [i.e., p = 1in (6)] in one spatial dimension, the integral in (11) can be expressed in terms of the Fourier expansion as follows 2 2 2 τ −i (t +s)∂ −i (t +s)∂ i (t +s)∂ n n n x x x I (w, t ) = e e w e w ds 2 2 2 2 i (t +s)  +k −k −k i x 1 2 3 = e ds w ˆ w ˆ w ˆ e k k k 1 2 3 ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 (13) 2 2 2 2 it  +k −k −k i x 1 2 3 = e w ˆ w ˆ w ˆ e k k k 1 2 3 ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 is 2k −2k (k +k )+2k k 1 2 3 2 3 e ds. Naturally, the next step would be the exact integration of the appearing integral is 2k −2k (k +k )+2k k 1 2 3 2 3 e ds. (14) However, the obtained convolution cannot be converted into the physical space straightforwardly and therefore only yields a practical scheme in dimension one as the iteration needs to be carried out fully in Fourier space. In order to obtain an efficient and practical low regularity implementation, we therefore only treat the dominant quadratic term 2k exactly. For this purpose, we write the integrand in (14)as iβ 2 e − 1 is 2k −2k (k +k )+2k k 1 2 3 2 3 2isk γ 1 1 (15) e = e 1 +|β| |β| 123 Found Comput Math (2018) 18:731–755 737 with β = s(−2k (k + k ) + 2k k and some 0 ≤ γ ≤ 1. Inserting (15)into(13) 1 2 3 2 3 and integrating the first term 2isk 2 e ds = τϕ (2i τ k ) yields that 2 2 2 2 τ it  +k −k −k 2 i x 1 2 3 I (w, t ) = τ e ϕ (2i τ k )w ˆ w ˆ w ˆ e n 1 k k k 1 1 1 2 3 ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 (16) + R (w, t ) 2 2 2 −it ∂ it ∂ 2 −it ∂ τ n n n x x x = τ e e w ϕ (−2iτ∂ )e w + R (w, t ), 1 n x 1 where 2 2 2 2 τ 2 2r it  +k −k −k 1 2 3 R (w, t ) = (1 +||) e 2k (k + k ) − 2k k n 1 2 3 2 3 1 r ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 ⎛ ⎞ is −2k (k +k )+2k k τ 1 2 3 2 3 2 e − 1 i x 2isk γ 1 ⎝ ⎠ w ˆ w ˆ w ˆ e e s ds . k k k 1 2 3 γ s |2k (k + k ) − 2k k 0 1 2 3 2 3 −γ iβ Since |β| (e − 1) is uniformly bounded for β ∈ R and 0 ≤ γ ≤ 1, the integral in 1+γ the above formula is of the order τ . This shows that ⎛ ⎞ ⎜   ⎟ τ 2 2+2γ 2r ⎜ ⎟ R (w, t ) ≤ cτ (1 +||) |k k |+|k k |+|k k | |ˆ w ||w ˆ ||w ˆ | n 1 2 1 3 2 3 k k k 1 r 1 2 3 ⎝ ⎠ ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 ⎛ ⎞ ⎜ ⎟ γ γ γ 2+2γ 2r ⎜ ⎟ ≤ cτ (1 +||) 1 +|k | 1 +|k | 1 +|k | |ˆ w ||w ˆ ||w ˆ | 1 2 3 k k k 1 2 3 ⎝ ⎠ ∈Z k ,k ,k ∈Z 1 2 3 =−k +k +k 1 2 3 (17) for some constant c > 0 and 0 ≤ γ ≤ 1. ikx Next we define the auxiliary function g(x ) = g ˆ e through its Fourier k∈Z coefficients g ˆ = (1 +|k|) |ˆ w | (18) k k which allows us to express the bound on the remainder given in (17)asfollows τ 1+γ 3 R (w, t ) ≤ cτ g  . n r r 123 738 Found Comput Math (2018) 18:731–755 With the aid of the bilinear estimate (5) we can thus conclude that for r > 1/2 3/2 τ 1+γ 3 1+γ 2(r +γ) 2 R (w, t ) ≤ cτ g = cτ (1 +|k|) |ˆ w | , n r k 1 r k∈Z where we have used the definition of the Fourier coefficients in (18). Hence, we can conclude the following bound on the remainder τ 1+γ 3 R (w, t ) ≤ cτ w for 0 ≤ γ ≤ 1, r > 1/2 (19) n r r +γ for some constant c > 0. Next, we replace the integral in (11) with v in place of v(t ) by the term (16) τ n without remainder R and w = v . Twisting the solution back again, i.e., setting n it ∂ n u := e v yields the following exponential-type integration scheme for the cubic NLS n+1 iτ∂ n n 2 2 u = e u − iτμ(u ) ϕ (−2iτ∂ )u . (20) Next we generalize the above scheme to arbitrary dimensions d ≥ 1. 2.1.2 Cubic Nonlinearities p = 1 in Dimension d ≥ 1 In the following we use the notation j := ( j ,..., j ) ∈ Z , j · x := j x + ··· + j x . 1 d 1 1 d d In the case of a cubic nonlinearity p = 1 in dimension d, the integral in (11) can be expressed in terms of the Fourier expansion as follows τ −i (t +s) −i (t +s) i (t +s) n n n I (w, t ) = e e w e w ds 2 2 2 2 i (t +s)  +j −k −l n i ·x = e ds w ˆ w ˆ w ˆ e . j k l d d ∈Z j,k,l∈Z =−j+k+l Note that we have 2 2 2 2 2 + j − k − l = 2j − 2j · (k + l) + 2k · l. In order to obtain an efficient and practical low regularity implementation, we again treat the dominant quadratic term 2j exactly. 123 Found Comput Math (2018) 18:731–755 739 This yields that 2 2 2 2 2 τ it  +j −k −l i ·x 2isj τ I (w, t ) = e w ˆ w ˆ w ˆ e e ds + R (w, t ) n j k l n 1 1 d d ∈Z j,k,l∈Z =−j+k+l 2 2 2 2 it  +j −k −l n i ·x 2 τ = e w ˆ w ˆ w ˆ e ϕ (2i τ j ) + R (w, t ), j k l 1 n d d ∈Z j,k,l∈Z =−j+k+l (21) where similarly to above [see (19)] we have that τ 1+γ 3 R (w, t ) ≤ cτ w for 0 ≤ γ ≤ 1, r > d/2 n r 1 r +γ for some constant c > 0. Inserting the approximation (21)into(11) and twisting the solution back again, i.e., setting n it  n u := e v yields the generalization of the exponential-type integrator (20) to dimensions d ≥ 1 n+1 iτ n n 2 u = e u − iτμ(u ) ϕ (−2iτ)u . Next we consider the general case of (6) with p ∈ N. 2.2 Nonlinearities with Integer p In the case of a general nonlinearity p ∈ N in (6), the integral in (11) can be expressed in terms of the Fourier expansion as follows τ −i (t +s) i (t +s) 2 p i (t +s) n n n I (w, t ) = e |e w| e w ds. Similar observations to above (see also [28]) yield the following numerical scheme p+1 p n+1 n −it  it  n −it n n n n v = v − i μτ e e v ϕ (−2iτ) e v (22) which satisfies n+1 n τ n τ n v = v − i μI (v , t ) − i μR (v , t ), (23) n n p p with 2 p+1 τ n 1+γ n R (v , t ) ≤ cτ v  for 0 ≤ γ ≤ 1, r > d/2(24) n r p r +γ for some constant c > 0. 123 740 Found Comput Math (2018) 18:731–755 In the original variable u, the exponential-type integration scheme for the nonlinear Schrödinger equation (6) reads p+1 p n+1 iτ n n u = e u − i μτ u ϕ (−2iτ) u . (25) In the following section, we give an error analysis for the above scheme. 3 Error Analysis In this section, we carry out the error analysis of the exponential-type integrator (22). In the following, we set p+1 p τ −it  it  −it (w) := w − i μτ e e w ϕ (−2iτ) e w k+1 τ k such that in particular v = (v ). Lemma 3.1 (Stability) Let r > d/2 and f, g ∈ H . Then, for all t ∈ R we have τ τ ( f ) − (g) ≤ (1 + τ |μ|L) f − g , r r t t where L depends on  f  and g . r r Proof Note that for all t ∈ R and w ∈ H it holds by (12) that ϕ (2it )w ≤ cw . 1 r it  r r Thus, as e is a linear isometry on H and H is an algebra for r > d/2 we obtain that p+1 p τ τ it  −it ( f ) − (g) ≤ f − g + τ |μ| e f ϕ (−2iτ) e f r r 1 t t p+1 p it  −it − e g ϕ (−2iτ) e g 2 p 2 p−l ≤ f − g + τ c|μ| f − g  f  g . r r r l=0 This proves the assertion. Lemma 3.2 (Local error) Let r > d/2 and 0 ≤ γ ≤ 1. Assume that v(t + t ) = t r +γ t φ (v(t )) ∈ H for 0 ≤ t ≤ τ , where φ denotes the exact flow of (7). Then, τ τ 1+γ φ (v(t )) − (v(t )) ≤ cτ , k k r where c depends on sup φ (v(t )) . k r +γ 0≤t ≤τ 123 Found Comput Math (2018) 18:731–755 741 Proof Note that by (23)wehavethat τ −i (t +s) i (t +s) 2 p i (t +s) k k k (v(t )) = v(t ) − i μ e |e v(t )| e v(t ) ds k k k k + i μR (v(t ), t ). k k Furthermore, (8) implies that τ −i (t +s) i (t +s) 2 p i (t +s) k k k φ (v(t )) = v(t ) − i μ e |e v(t + s)| e v(t + s) ds. k k k k it  r r Thus as e is a linear isometry on H and H is an algebra for r > d/2 we obtain that 2 p τ τ φ (v(t )) − (v(t )) ≤|μ| sup v(t + s) v(t + s) − v(t ) ds k k r k r k k r 0≤s≤τ 0 +|μ|R (v(t ), t ) . k k r Together with the estimate on the difference v(t + s) − v(t ) given in (9) and the k k estimate on the remainder R (v, t ) given in (24) this yields for r > d/2 and 0 ≤ γ ≤ that 4 p+1 2 p+1 τ τ 2 2 1+γ φ (v(t )) − (v(t )) ≤ τ c|μ| sup v(t + s) + τ c|μ|v(t ) k k r k r k t r +γ 0≤s≤τ which concludes the stated local error bound. These two lemmata allow us to prove the following convergence bound. Theorem 3.3 Let r > d/2 and 0 <γ ≤ 1. Assume that the exact solution of (7) r +γ satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (22) 0 n is bounded by n γ v(t ) − v  ≤ cτ , n r where c depends on sup v(t ) . r +γ 0≤t ≤T Proof The triangle inequality yields that k+1 τ τ k v(t ) − v  =φ (v(t )) − (v ) k+1 r k r τ τ τ τ k ≤φ (v(t )) − (v(t )) + (v(t )) − (v ) . k k r k r t t t k k k k r Thus, we obtain with the aid of Lemmata 3.2 and 3.1 that as long as v ∈ H we have that k+1 1+γ k (26) v(t ) − v  ≤ cτ + (1 + τ |μ|L) v(t ) − v  , k+1 r k r 123 742 Found Comput Math (2018) 18:731–755 where c depends on sup v(t ) and L depends on v(t ) as well as on r +γ k r t ≤t ≤t k k+1 k r +γ k r v  . As long as v(t ) ∈ H and v ∈ H for 0 ≤ k ≤ n we therefore obtain by r k iterating the estimate (26) that n+1 1+γ n v(t ) − v  ≤ cτ + (1 + τ |μ|L) v(t ) − v n+1 r n r 1+γ τ |μ|L 1+γ n−1 ≤ cτ + e cτ + (1 + τ |μ|L) v(t ) − v n−1 r 1+γ nτ |μ|L γ t |μ|L ≤ cnτ e ≤ ct τ e . The assertion then follows by induction, respectively, a Lady Windermere’s fan argu- ment (see, for example [11,15,25]). The above theorem immediately yields a convergence result for the exponential- type integration scheme (22) approximating the solution of the nonlinear Schrödinger equation (6). Corollary 3.4 Let r > d/2 and 0 <γ ≤ 1. Assume that the exact solution of (6) r +γ satisfies u(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (25) 0 n is bounded by n γ u(t ) − u  ≤ cτ , n r where c depends on sup u(t ) . r +γ 0≤t ≤T it  r Proof As e is a linear isometry on H for all t ∈ R we have that n it  n n u(t ) − u  =e v(t ) − v  =v(t ) − v  . n r n r n r Thus, the assertion follows from Theorem 3.3. 4 Quadratic Schrödinger Equations in One Space Dimension In this section, we focus on quadratic nonlinear Schrödinger equations of type (see, e.g., [1]) 2 2 0 i ∂ u =−∂ u + μu , u(0, x ) = u (x ), (t, x ) ∈ R × T,μ ∈ R (27) as well as of type (see, e.g, [21]) 2 2 0 i ∂ u =−∂ u + μ|u| , u(0, x ) = u (x ), (t, x ) ∈ R × T,μ ∈ R, (28) which have recently gained a lot of attention in the literature, see for instance [1,3,5, 13,21–23,27,29,32] and the references therein. Here, the key relations 2 2 2 (k + k ) − k − k = 2k k (29) 1 2 1 2 1 2 123 Found Comput Math (2018) 18:731–755 743 and 2 2 2 (k − k ) − k + k =−2k (k − k ) (30) 1 2 2 1 2 1 2 allow us to derive first-order exponential-type integrators for Schrödinger equations of type (27) and (28) without imposing any additional regularity assumptions on the solution. 4.1 Quadratic Schrödinger Equations of Type (27) −it ∂ Set v(t ) = e u(t ), where u is the solution of the Schrödinger equation (27). Then the twisted variable v(t ) satisfies 2 2 −it ∂ it ∂ x x i ∂ v = μe e v(t ) (31) which by Duhamel’s formula yields that 2 2 −i (t +s)∂ i (t +s)∂ n n x x v(t + τ) = v(t ) − i μ e e v(t + s) ds. (32) n n n To construct our numerical scheme, we proceed as above. We start from the following approximation 2 2 −i (t +s)∂ i (t +s)∂ τ n n x x v(t + τ) ≈ v(t ) − i μ e e v(t ) ds =: (v(t )) (33) n n n n and rewrite the above integral in Fourier space as follows 2 2 τ −i (t +s)∂ i (t +s)∂ n n x x I (w, t ) = e e w ds 2 2 2 −i (t +s)∂ i (t +s)∂ ik x i (t +s)∂ ik x n n 1 n 2 x x x = e e w ˆ e e w ˆ e ds k k 1 2 k k 1 2 2 2 2 i (t +s) (k +k ) −k −k n 1 2 i (k +k )x 1 2 1 2 = e w ˆ w ˆ e ds. k k 1 2 k ,k 1 2 The key relation (29) now shows that 2 2 2 i τ (k +k ) −k −k 1 2 1 2 2 2 2 e − 1 it (k +k ) −k −k τ n 1 2 i (k +k )x 1 2 1 2 I (w, t ) = e w ˆ w ˆ e n k k 1 2 2ik k 1 2 k ,k 1 2 k ,k =0 1 2 ik x 2 + 2τ w ˆ w ˆ e + τ w ˆ 0 k 1 0 k =0 123 744 Found Comput Math (2018) 18:731–755 2 2 2 i τ (k +k ) −k −k 1 2 1 2 2 2 2 i e − 1 it (k +k ) −k −k i (k +k )x n 1 2 1 2 1 2 = e w ˆ w ˆ e k k 1 2 2 (ik )(ik ) 1 2 k ,k 1 2 k ,k =0 1 2 ik x 2 + 2τ w ˆ w ˆ e − τ w ˆ 0 k 1 0 k ∈Z 2 2 i 2 2 2 2 −it ∂ −iτ∂ i (t +τ)∂ −1 it ∂ −1 n n n x x x x = e e e ∂ w − e ∂ w x x + 2τ w ˆ w − τ w ˆ . Together with the approximation in (33) this yields the following integration scheme: n+1 τ n v = (v ) n n n 2 = 1 − 2i μτ v ˆ v + i μτ (v ˆ ) 0 0 (34) 2 2 μ 2 2 2 2 −it ∂ −iτ∂ i (t +τ)∂ −1 n it ∂ −1 n n n n x x x x + e e e ∂ v − e ∂ v . x x In order to obtain an approximation to the original solution u(t ) of the quadratic Schrödinger equation (27) at time t = nτ , we twist the variable back again, i.e., we n it ∂ n set u := e v . This yields the following exponential-type integrator 2 2 2 μ 2 μ 2 n+1 n iτ∂ n n 2 iτ∂ −1 n iτ∂ −1 n x x x u = 1 − 2i μτ u ˆ e u + i μτ (u ˆ ) + e ∂ u − e ∂ u . x x 0 0 2 2 (35) 4.2 Quadratic Schrödinger Equations of Type (28) For the quadratic Schrödinger equation (28), the equation in the twisted variable v(t ) = −it ∂ e u(t ) reads 2 2 −it ∂ it ∂ x x i ∂ v = μe e v(t ) . To construct our numerical scheme, we proceed as above. We start from the following approximation 2 2 −i (t +s)∂ i (t +s)∂ n n x x v(t + τ) ≈ v(t ) − i μ e e v(t ) ds (36) n n n and rewrite the above integral in Fourier space as follows 2 2 τ −i (t +s)∂ i (t +s)∂ n n x x I (w, t ) = e e w ds 2 2 2 i (t +s) (k −k ) −k +k n 1 2 i (k −k )x 1 2 1 2 = e w ˆ w ˆ e ds. k k 1 2 k ,k 1 2 123 Found Comput Math (2018) 18:731–755 745 The key relation (30) now shows that 2 2 2 i τ (k −k ) −k +k 1 2 1 2 2 2 2 e − 1 it (k −k ) −k +k τ n 1 2 i (k −k )x 1 2 1 2 I (w, t ) = e w ˆ w ˆ e n k k 1 2 −2ik (k − k ) 2 1 2 k ,k 1 2 k =k ,k =0 1 2 2 ik x 2 +τ w ˆ w ˆ e + τ w 0 k 2 i 2 2 2 −1 −i (t +τ)∂ i (t +τ)∂ −i (t +τ)∂ −1 n n n x x x = ∂ e e w e ∂ w x x 2 2 2 −1 −it ∂ it ∂ −it ∂ −1 2 n n n x x x − ∂ e e w e ∂ w + τ w ˆ w + τ w . 0 2 x x Together with the approximation in (36) this yields (by twisting the variable back it ∂ again) the following integration scheme in the original solution u(t ) = e v(t ) n 2 n+1 iτ∂ n n 2 u = 1 − i μτ u ˆ e u − i μτ u (37) μ 2 2 2 −1 iτ∂ n −iτ∂ −1 n iτ∂ n −1 n x x x + ∂ e u e ∂ u − e u ∂ u . x x x 4.3 Error Analysis In this section, we carry out the error analysis of the exponential-type integrators (35) and (37). In the following let r > 1/2. We commence with the quadratic Schrödinger equation of type (27). Let φ denote the exact flow of (31), i.e., v(t ) = φ (v(0)). The following lemma provides a local error bound of the scheme (34). t r Lemma 4.1 Let r > 1/2 and assume that v(t + t ) = φ (v(t )) ∈ H for 0 ≤ t ≤ τ , k k where φ denotes the exact flow of (27). Then τ τ 2 φ (v(t )) − (v(t )) ≤ cτ , k k r where c depends on sup φ (v(t )) . k r 0≤t ≤τ it ∂ r Proof As e is a linear isometry on H for all t ∈ R we obtain by taking the difference of (32) with the approximation (33) and using the bilinear estimate (5) that 2 2 2 2 τ τ i (t +s)∂ i (t +s)∂ k k x x φ (v(t )) − (v(t )) ≤|μ|  e v(t + s) − e v(t )  ds k k r k k 2 2 i (t +s)∂ i (t +s)∂ k k x x ≤ c|μ|τ sup e (v(t + s) − v(t ))  e (v(t + s) + v(t )) k k r k k r 0≤s≤τ ≤ c|μ|τ sup v(t + s) − v(t ) k k r 0≤s≤τ 2 i (t +ξ)∂ ≤ c|μ| τ sup  e v(t + ξ)  dξ, 0≤s≤τ 0 where c depends on sup φ (v(t )) . This yields the assertion. k r 0≤s≤τ 123 746 Found Comput Math (2018) 18:731–755 Next we state the stability estimate. Lemma 4.2 Let r > 1/2 and f, g ∈ H . Then, for all t ∈ R we have τ τ ( f ) − (g) ≤ (1 + τ |μ|L)  f − g , r r t t where L depends on  f + g . Proof The assertion follows from the representation of the numerical flow given in (33) together with the bilinear estimate (5). The above lemmata allow us to prove the following global convergence result. Theorem 4.3 Let r > 1/2. Assume that the exact solution of (31) satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (34) is bounded by 0 n v(t ) − v  ≤ cτ, n r where c depends on sup v(t ) . 0≤t ≤T Proof The proof follows the line of argumentation to the proof of Theorem 3.3 thanks to the local error bound in Lemma 4.1 and the stability estimate of Lemma 4.2. The above theorem immediately yields first-order convergence of the exponential- type integrator (35)(in H ) without any additional smoothness assumptions on the exact solution of (27). Corollary 4.4 Let r > 1/2. Assume that the exact solution of (27) satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (35) is bounded by 0 n u(t ) − u  ≤ cτ, n r where c depends on sup u(t ) . 0≤t ≤T it ∂ Proof The bound follows directly from Theorem 4.3 as e is a linear isometry on H for all t ∈ R. Similarly to above we obtain the following first-order convergence result for the exponential-type integrator (37) approximating the solution of the quadratic Schrödinger equation (28). Corollary 4.5 Let r > 1/2. Assume that the exact solution of (28) satisfies v(t ) ∈ H for 0 ≤ t ≤ T . Then, there exists a constant τ > 0 such that for all step sizes 0 <τ ≤ τ and t = nτ ≤ T we have that the global error of (37) is bounded by 0 n u(t ) − u  ≤ cτ, n r where c depends on sup u(t ) . 0≤t ≤T Proof The proof is similar to that of Corollary 4.4 and will therefore be omitted here. 123 Found Comput Math (2018) 18:731–755 747 5 Numerical Experiments In this section, we compare the newly developed exponential-type integrator (25) applied to Schrödinger equations with low regularity solutions in the energy space with the classical Lie splitting, the classical Strang splitting and the classical first- order exponential integrator. In the numerical experiments, we use a standard Fourier pseudospectral method for the space discretization where we choose the largest Fourier mode K = 2 (i.e., the spatial mesh size x = 0.0061). The numerical experiments show that the new integrator (25) is preferable over the classical splitting and expo- nential integration methods for low regularity solutions. In particular, the experiments indicate that the proposed integrator (25) is convergent of order one even for rough solutions whereas the classical methods suffer from order reduction. Although the new integrator (25) was derived for nonlinearities with integer p ≥ 1, it can be applied to problem (6) with non-integer p as well. Preliminary numerical results for small non-integer values of p are presented in Sect. 5.3. 5.1 Numerical Experiments for Schrödinger Equations of Type (6) In this section, we consider nonlinear Schrödinger equations of type (6). The associated Lie splitting (of classical order one) reads n+1/2 iτ n u = e u , (38) 2 p n+1/2 n+1 −iτμ u n+1/2 u = e u , L L and the associated Strang splitting (of classical order two) is given by n+1/2 i  n u = e u , 2 p n+1/2 n+1/2 −iτμ u n+1/2 − (39) u = e u , + − n+1 n+1/2 u = e u . Furthermore, the classical first-order exponential integrator (of classical order one)is defined through n+1 iτ n n 2 p n u = e u − i μτ ϕ (iτ) |u | u . (40) E E E In our numerical experiments, we choose μ = 1 (unless otherwise noted) and set K 2K U := [u(0, x ), u(0, x ), . . . , u(0, x )]= rand(2K , 1) + i rand(2K , 1) ∈ C , 1 2 2K (41) where x = (−1 + )π and rand(2K , 1) return 2K uniformly distributed random numbers between 0 and 1. We compare the above schemes with the newly derived exponential-type integrator (25) for initial values 123 748 Found Comput Math (2018) 18:731–755 0.17 0.17 0.16 0.16 0.15 0.15 Fig. 1 Initial values (42) normalized in L for two different values of ϑ. Left ϑ = . Right ϑ = 5 -2 -2 -3 -3 -4 -4 -5 -2 -1 -2 -1 10 10 10 10 Fig. 2 (cubic NLS) Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39) (magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type integration 3/2 2 scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. The slope of the continuous lines is one half and one, respectively (Color figure online) −ϑ |k| if k = 0, K −ϑ K −ϑ u =|∂ | U , |∂ | := (42) x ,K x ,K 0if k = 0 for different values of ϑ normalized in L . (For typical initial values, see Fig. 1.) All experiments are carried out with constant step sizes τ = for j = 1,..., 512. 5.1.1 Cubic Schrödinger Equations In this section, we consider the classical cubic Schrödinger equation, i.e., p = 1in (6). The error at time T = 1 measured in a discrete H norm of the exponential-type integrator (25), the Lie splitting scheme (38), the Strang splitting scheme (39)aswell as the classical exponential integrator (40) for the initial value (42) normalized in L is illustrated in Fig. 2 for ϑ = 3/2, respectively, ϑ = 2 and in Fig. 3 for ϑ = 3, respectively, ϑ = 5. As a reference solution, we take the schemes themselves with a very small time step size if ϑ = 5. For ϑ = 5 we take the Strang splitting scheme (39) with a very small time step size as the reference solution for all schemes. The experiments show that Strang splitting is second-order convergent for ϑ ≥ 5, and 123 Found Comput Math (2018) 18:731–755 749 -2 -2 10 10 -3 -3 10 10 -4 -4 10 10 -5 -5 10 10 -6 -6 10 10 -7 -7 10 10 -8 -8 10 10 -9 -9 10 10 -2 -1 -2 -1 10 10 10 10 Fig. 3 (cubic NLS) Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39)(magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type inte- 3 5 gration scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. The slope of the continuous lines is one and two, respectively (Color figure online) that Lie splitting is first-order convergent for ϑ ≥ 3. These observations are in line with the convergence proofs from the literature [10,25], see also our discussion in the introduction. For smaller values of ϑ, both Lie and Strang splitting show a zigzag behavior, depending on whether the local errors accumulate or happy error cancelation occurs. The classical exponential integrator of order one is first-order convergent for ϑ ≥ 3 as expected (see the discussion in the introduction). For smaller values of ϑ, the convergence behavior of the exponential integrator gets less regular and order reduction occurs. On the other hand, the new exponential-type integrator introduced in this paper is first-order convergent for all ϑ ≥ 2. This behavior is in line with Theorem 3.3. The integrator even shows first-order convergence and small errors for ϑ = . 5.1.2 Quintic Schrödinger Equations In this section, we consider the quintic Schrödinger equation, i.e., p = 2in (6), which appears for instance as the mean field limit of a Boson gas with three-body interactions, see [7]. The error at time T = 1 measured in a discrete H norm of the exponential- type integrator (25), the Lie splitting scheme (38), the Strang splitting scheme (39)as well as the classical exponential integrator (40) for the initial value (42) normalized in L is illustrated in Fig. 4 for ϑ = 3/2, respectively, ϑ = 2 and in Fig. 5 for ϑ = 3, respectively, ϑ = 5. As a reference solution we take the schemes themselves with a very small time step size if ϑ = 5. For ϑ = 5 we take the Strang splitting scheme (39) with a very small time step size as the reference solution for all schemes. The outcome of these numerical experiments is almost the same as for the cubic nonlinear Schrödinger equation. The only notable difference is the exponential-type integrator for ϑ = , where the convergence behavior is now less regular. Note, however, that the errors are still much smaller compared to the other methods. 123 750 Found Comput Math (2018) 18:731–755 -2 -3 -3 -4 -4 -5 -5 -6 -6 10 10 -2 -1 -2 -1 10 10 10 10 Fig. 4 (quintic NLS) Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39)(magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type inte- 3/2 2 gration scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. Left picture The slope of the continuous lines is one quarter and one, respectively. Right picture The slope of the continuous lines is one half and one, respectively (Color figure online) -4 -4 -5 -5 10 -6 -6 -7 -7 -8 -8 -9 -9 -2 -1 -2 -1 10 10 10 10 Fig. 5 (quintic NLS) in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39) (magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type integra- 3 5 tion scheme (25)(red, squares). Left picture H solutions. Right picture H solutions. The slope of the continuous lines is one and two, respectively (Color figure online) 5.2 Numerical Experiments for Quadratic Schrödinger Equations In this section, we consider quadratic Schrödinger equations of type (27). In order to derive the splitting methods, we split the right-hand side of (27) into the linear and nonlinear part, respectively. Note that (for sufficiently small t) the exact solution of the subproblem iu (t ) = μu(t ) is given by u(0) u(t ) = . 1 + i μtu(0) 123 Found Comput Math (2018) 18:731–755 751 Thus, the associated Lie splitting scheme (of classical order one) reads 2 u n+1 iτ∂ L u = e (43) L n 1 + i μτ u and the associated Strang splitting scheme (of classical order two) is given by τ 2 n+1/2 i ∂ n u = e u S S n+1/2 (44) τ 2 u n+1 i ∂ S u = e . n+1/2 1 + i μτ u In Examples 5.1 and 5.2, we compare the above splitting methods as well as the classical first-order exponential integration scheme n+1 iτ∂ n 2 n u = e u − i μτ ϕ (iτ∂ ) u (45) E x E with the newly derived exponential-type integrator (35)for smooth and non-smooth solutions, respectively, where we set μ = 1 and integrate up to T = 1. The numerical experiments in particular indicate that the exponential-type integrator (35)ispreferable over the commonly used splitting methods for low regularity solutions. This is due to the fact that within this exponential integration scheme all the stiff parts (i.e., the terms involving the differential operator ∂ ) are solved exactly. In Example 5.3 we consider a quadratic Schrödinger equation of type (27) with a small nonlinearity on a longer time interval T = 10. The numerical experiment indicates the favorable behavior of the exponential integrator (35) even in the smooth case. This is due to the fact that its error constant is triggered by the nonlinearity (and not by the differential operator ∂ ). Example 5.1 (Smooth solutions) We choose the smooth initial value u(0, x ) = sin x cos x ∈ C (T) (46) 2 2 normalized in L . The error at time T = 1 measured in a discrete L norm of the exponential-type integrator (35), the Lie splitting scheme (43), the Strang splitting scheme (44) and the classical exponential integrator (45) is illustrated in Fig. 6.As a reference solution, we take the Strang splitting scheme (44) with a very small time step size. All the methods show their classical orders of convergence, as expected. Example 5.2 (Non-smooth solutions) We choose the non-smooth initial value (41) 2 2 normalized in L . The error at time T = 1 measured in a discrete L norm of the exponential-type integrator (35), the Lie splitting scheme (43), the Strang splitting scheme (44) and the classical exponential integrator (45) is illustrated in Fig. 7.Asa reference solution, we take the schemes themselves with a very small time step size. For this example, the classical exponential integrator is not convergent. Lie and Strang splitting suffer from a strong order reduction (down to order one half). Here, the new 123 752 Found Comput Math (2018) 18:731–755 Fig. 6 (quadratic NLS) Error of -1 the Lie splitting (43)(blue, dotted), Strang splitting (44) -2 (magenta, solid line), classical exponential integrator (45) -3 (green, dashed)and exponential-type integration scheme (35)(red, squares)for a -4 smooth solution with initial value (46)and μ = 1. The slope -5 of the continuous lines is one and two, respectively (Color figure online) -6 -2 -1 10 10 Fig. 7 (quadratic NLS) Error of 0 the Lie splitting (43)(blue, dotted), Strang splitting (44) (magenta, solid line), classical exponential integrator (45) -1 (green, dashed)and exponential-type integration scheme (35)(red, squares)for a non-smooth solution with initial -2 value (41)and μ = 1. The slope of the continuous lines is one half and one, respectively (Color figure online) -3 -2 -1 10 10 exponential-type integrator is by far the best method. It is first-order convergent, as predicted by Theorem 4.3. Example 5.3 (Small nonlinearity) We consider the quadratic Schrödinger equation of type (27) with small nonlinearity, i.e., 2 2 i ∂ u =−∂ u + μu ,μ = 0.01, T = 10. (47) The error at time T = 10 measured in a discrete L norm of the exponential-type integrator (35), the Lie splitting scheme (43), the Strang splitting scheme (44) and the classical exponential integrator (45) for smooth and non-smooth solutions is illustrated in Fig. 8 left and right, respectively. As a reference solution, we take the Strang splitting scheme for smooth and the schemes themselves for non-smooth solutions with a very small time step size. For smooth initial data, all methods show their classical orders of convergence. The new exponential-type integrator, however, is by far the most accurate scheme (for the considered range of step sizes). For non-smooth initial data, the classical exponential integrator does not converge; Lie and Strang splitting suffer 123 Found Comput Math (2018) 18:731–755 753 -3 -2 -4 -3 -5 10 -4 -6 10 -5 -7 -6 -8 -2 -1 -2 -1 10 10 10 10 Fig. 8 (quadratic NLS) Error of the Lie splitting (43)(blue, dotted), Strang splitting (44)(magenta, solid line), classical exponential integrator (45)(green, dashed) and exponential-type integration scheme (35) (red, squares)fora small nonlinearity (47). Left picture Smooth initial value (46). The slope of the continuous lines is one and two, respectively. Right picture Non-smooth initial value (41). The slope of the continuous lines is one half, and one, respectively (Color figure online) 0 0 10 10 -1 -1 -2 -2 -3 -3 -4 10 -4 -5 -5 -6 -2 -1 -2 -1 10 10 10 10 0 0 10 10 -1 10 -1 -2 -2 -3 -3 -4 -4 -2 -1 -2 -1 10 10 10 10 3 1 1 Fig. 9 Exponent p = 1(upper row, left)and p = (upper row, right), p = (lower row, left)and p = 4 2 4 (lower row, right). Error in the energy space H of the Lie splitting (38)(blue, dotted), Strang splitting (39) (magenta, solid line), classical exponential integrator (40)(green, dashed) and exponential-type integration scheme (25)(red, dash dotted) with smooth initial value (48). The slope of the continuous and dotted line is one and two, respectively (Color figure online) 123 754 Found Comput Math (2018) 18:731–755 again from a strong order reduction down to order one half. Only the new exponential- type integrator is very accurate and first-order convergent, as expected. 5.3 Numerical Experiments for Semilinear Schrödinger Equations with Non-integer Exponents p In this section, we numerically investigate the convergence rates in the case of semi- linear Schrödinger equations of type (6) with non-integer exponents p.InFig. 9 we plot the error at time T = 1 measured in a discrete H norm of the exponential-type integrator (25), the Lie splitting scheme (38), the Strang splitting scheme (39)aswell as the classical exponential integrator (40) for different values of p > 0. In all the simulations, we take the smooth initial value u(0, x ) = sin x ∈ C (T) (48) and μ = 1. Note, however, that the nonlinearity is only smooth for p > . For p = 1, all methods show their classical orders of convergence, as expected. For smaller values of p, however, the convergence behavior of Strang and Lie splitting deteriorates. Surprisingly, both exponential integrators retain their first-order conver- gence down to very small values of p. This is the subject of future research. Acknowledgements Open access funding provided by University of Innsbruck and Medical University of Innsbruck. K. Schratz gratefully acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG) through CRC 1173. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References 1. I. Bejenaru, T. Tao, Sharp well-posedness and ill-posedness results for a quadratic non-linear Schrödinger equation. J. Funct. Anal. 233:228–259 (2006). 2. C. Besse, B. Bidégaray, S. Descombes, Order estimates in time of splitting methods for the nonlinear Schrödinger equation, SIAM J. Numer. Anal. 40:26–40 (2002). 3. J. Bourgain, Fourier transform restriction phenomena for certain lattice subsets and applications to nonlinear evolution equations. Part I: Schrödinger equations. Geom. Funct. Anal. 3:209–262 (1993). 4. B. Cano, A. González-Pachón, Exponential time integration of solitary waves of cubic Schrödinger equation. Appl. Numer. Math. 91:26–45 (2015). 5. T. Cazenave, F.B. Weissler, The Cauchy problem for the critical nonlinear Schrödinger equation. Nonlinear Anal., 14:807–836 (1990). 6. E. Celledoni, D. Cohen, B. Owren, Symmetric exponential integrators with an application to the cubic Schrödinger equation. Found. Comput. Math. 8:303–317 (2008). 7. T.Chen,N.Pavlovic, ´ The quintic NLS as the mean field limit of a boson gas with three-body interactions. J. Funct. Anal. 260:959–997 (2011). 8. D. Cohen, L. Gauckler, One-stage exponential integrators for nonlinear Schrödinger equations over long times. BIT 52:877–903 (2012). 9. G. Dujardin, Exponential Runge-Kutta methods for the Schrödinger equation. Appl. Numer. Math. 59:1839–1857 (2009). 123 Found Comput Math (2018) 18:731–755 755 10. J. Eilinghoff, R. Schnaubelt, K. Schratz, Fractional error estimates of splitting schemes for the non- linear Schrödinger equation. J. Math. Anal. Appl. 442:740–760 (2016). 11. E. Faou, Geometric Numerical Integration and Schrödinger Equations. European Math. Soc. Publish- ing House, Zürich 2012. 12. L. Gauckler, Convergence of a split-step Hermite method for the Gross–Pitaevskii equation. IMA J. Numer. Anal. 31:396–415 (2011). 13. P. Germain, N. Masmoudi, J. Shatah, Global solutions for 3D quadratic Schrödinger equations. Int. Math. Res. Notices 3:414–432 (2009). 14. E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. Second edition, Springer, Berlin 2006. 15. E. Hairer, S. P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems. Second edition. Springer, Berlin 1993. 16. M. Hochbruck, A. Ostermann, Exponential integrators. Acta Numer. 19:209–286 (2010). 17. M. Hofmanova, K. Schratz, An exponential-type integrator for the KdV equation. Numer. Math. (2016). doi:10.1007/s00211-016-0859-1 18. H. Holden, K. H. Karlsen, K.-A. Lie, N. H. Risebro, Splitting for Partial Differential Equations with Rough Solutions. European Math. Soc. Publishing House, Zürich 2010. 19. L. I. Ignat, A splitting method for the nonlinear Schrödinger equation. J. Differential Equations 250:3022–3046 (2011). 20. A.-K. Kassam, L. N. Trefethen, Fourth-order time-stepping for stiff PDEs. SIAM J. Sci. Comput. 26:1214–1233 (2005). 21. C. Kenig, G. Ponce, L. Vega, Quadratic forms for the 1-D semilinear Schrödinger equation. Trans. Amer. Math. Soc. 346:3323–3353 (1996). 22. C. Kenig, G. Ponce, L. Vega, On the ill-posedness of some canonical dispersive equations. Duke Math. J. 106:617–633 (2001). 23. N. Kishimoto, Low-regularity bilinear estimates for a quadratic nonlinear Schrödinger equation. J. Differential Equations 247:1397–1439 (2009). 24. J. D. Lawson, Generalized Runge–Kutta processes for stable systems with large Lipschitz constants. SIAM J. Numer. Anal. 4:372–380 (1967). 25. C. Lubich, On splitting methods for Schrödinger–Poisson and cubic nonlinear Schrödinger equations. Math. Comp. 77:2141–2153 (2008). 26. R.I. McLachlan, G.R.W. Quispel, Splitting methods, Acta Numer. 11:341–434 (2002). 27. K. Nakanishi, H. Takaoka, Y. Tsutsumi, Counterexamples to bilinear estimates related with the KdV equation and the nonlinear Schrödinger equation. Methods Appl. Anal. 8:569–578 (2001). 28. K. Schratz (joint with A. Ostermann), Derivation of a low regularity exponential-type integrator for semilinear Schrödinger equations with polynomial nonlinearities. Oberwolfach Reports 18:928–931 (2016). 29. T. Tao, Multilinear weighted convolution of L functions, and applications to nonlinear dispersive equations. Amer. J. Math, 123:839–908 (2001). 30. T. Tao, Nonlinear Dispersive Equations. Local and Global Analysis. Amer. Math. Soc., Providence 31. M. Thalhammer, Convergence analysis of high-order time-splitting pseudo-spectral methods for non- linear Schrödinger equations. SIAM J. Numer. Anal. 50:3231–3258 (2012). 32. Y. Tsutsumi, L solutions for nonlinear Schrödinger equations and nonlinear groups. Funkcial. Ekvac., 30:115–125 (1987).

Journal

Foundations of Computational MathematicsSpringer Journals

Published: May 26, 2017

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off