A note on optimal H^1-error estimates for Crank-Nicolson approximations to the nonlinear Schrödinger equation

In this paper we consider a mass- and energy--conserving Crank-Nicolson time discretization for a general class of nonlinear Schrödinger equations. This scheme, which enjoys popularity in the physics community due to its conservation properties, was already subject to several analytical and numerical studies. However, a proof of optimal L^∞(H^1)-error estimates is still open, both in the semi-discrete Hilbert space setting, as well as in fully-discrete finite element settings. This paper aims at closing this gap in the literature.

Authors

• 8 publications
• 2 publications
06/09/2020

High-order mass- and energy-conserving SAV-Gauss collocation finite element methods for the nonlinear Schrödinger equation

A family of arbitrarily high-order fully discrete space-time finite elem...
03/08/2021

Unconditionally optimal convergence of an energy-conserving and linearly implicit scheme for nonlinear wave equations

In this paper, we present and analyze an energy-conserving and linearly ...
08/31/2020

Discrete conservation laws for finite element discretisations of multisymplectic PDEs

In this work we propose a new, arbitrary order space-time finite element...
06/28/2021

Finite Element Approximations of a Class of Nonlinear Stochastic Wave Equation with Multiplicative Noise

Wave propagation problems have many applications in physics and engineer...
07/02/2021

Error estimates for semi-discrete finite element approximations for a moving boundary problem capturing the penetration of diffusants into rubber

We study a semi-discrete finite element approximation of weak solutions ...
05/04/2021

Unconditional energy dissipation and error estimates of the SAV Fourier spectral method for nonlinear fractional generalized wave equation

In this paper, we consider a second-order scalar auxiliary variable (SAV...
12/03/2020

Carleman estimates and controllability results for fully-discrete approximations of 1-D parabolic equations

In this paper, we prove a Carleman estimate for fully-discrete approxima...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this paper we consider nonlinear Schrödinger equations (NLS) seeking a complex function such that

 i∂tu=−△u+Vu+γ(|u|2)u

in a bounded domain , with a homogenous Dirichlet boundary condition on and a given initial value. Here, is a known real-valued potential and is a smooth (and possibly nonlinear) function that depends on the unknown density . Of particular interest are cubic nonlinearities of the form , for some . In this case, the equation is called Gross–Pitaevskii equation. It has applications in optics [1, 13], fluid dynamics [30, 31] and, most importantly, in quantum physics, where it models for example the dynamics of Bose-Einstein condensates in a magnetic trapping potential [12, 20, 23]. Another relevant class are saturated nonlinearities, such as for some , which appear in the context of nonlinear optical wave propagation in layered metallic structures [11, 17] or the propagation of light beams in plasmas [22].

Nonlinear Schrödinger equations come with important physical invariants, where the mass and the energy are considered as two of the most crucial ones. When solving a NLS numerically it is therefore of great importance to also reproduce this conservation on the discrete level. This aspect was emphasized by various numerical studies [16, 25].

For the subclass of power law nonlinearities of the form for and , a mass and energy conserving relaxation scheme was proposed and analyzed by Besse [7, 8]. Thanks to its properties, the scheme shows a very good performance in realistic physical setups [16]. Despite the large variety of different numerical approaches for solving the time-dependent NLS (cf. [2, 3, 5, 14, 18, 19, 21, 24, 27, 26, 28, 29, 32] and the references therein) the literature knows however only one time discretization that conserves both mass and energy simultaneously for arbitrary (smooth) nonlinearities. This discretization, which was first mathematically studied by Sanz-Serna [24] and which is long-known in the physics community, is a Crank-Nicolson-type approach where the nonlinearity is approximated by a suitable difference quotient involving the primitive integral of . This is also the time discretization that we shall consider in this paper.

A combination of the method with a finite difference space discretization was proposed and analyzed by Bao and Cai [4, 6]. Combining the Crank–Nicolson time discretization with a finite element discretization in space, the first a priori error estimates for the arising method were obtained by Sanz-Serna 1984 [24] for cubic nonlinearities. He considers the case and derives optimal -error estimates under the coupling constraint , where denotes the time step size of the Crank-Nicolson method and the mesh size of the finite element discretization in space. In 1991, Akrivis et al. [2] improved this result by showing optimal convergence rates in in dimension and under the relaxed coupling constraint . Finally, in 2017 [15], the -error estimates could be improved yet another time by showing that the coupling constraint can be fully removed. Furthermore, general nonlinearities could be considered and the influence of potentials could be taken into account. However, so far, optimal error estimates in are still open in the literature.

One reason for this absence of -results could be related to the techniques used for the error analysis in previous works (cf. [2, 14, 18, 19, 24, 32]) which is based on the following steps: 1. Appropriate truncation of the nonlinearity to obtain a problem with bounded growth. 2. Analyzing the scheme with truncation in the FE space and deriving corresponding - and/or -error estimates. 3. Using inverse estimates in the finite element space to show that the truncated approximations are uniformly bounded in by a term of the form , with appropriate powers that depend on the considered space discretization, regularity and space dimension. 4. Concluding that if and are coupled in an appropriate way, then the truncated approximations are all uniformly bounded by a constant and hence coincide with a solution to the scheme without truncation.

This strategy does not only have the disadvantage that it produces unnecessary coupling conditions, but also that it becomes impractically technical when considering -error estimates for the Crank–Nicolson FEM. This is because it requires a suitable truncation of the primitive integral of that is on the one hand consistent with the energy conservation and on the other hand allows for uniform bounds of the approximations in . However, thanks to the new techniques developed in [29] and the CN error analysis suggested in [15] in the context of -error estimates, the truncation step is no longer necessary and the desired -bounds can be derived with elliptic regularity theory. With this, it is now possible to obtain estimates in a direct way, not only in the finite element setting, but also in the semi-discrete Hilbert space setting.

In this paper we will therefore build upon the results from [15, 29] to fill the gap in the literature and prove optimal -error estimates for the Crank–Nicolson approach without coupling constraints and for a general class of nonlinearities. The paper is structured as follows. In Section 2 we present the notation and the analytical assumptions on the problem. In Section 3 we present the time–discrete Crank–Nicolson method, we recall its well-posedness and optimal error estimates in . Furthermore, we present and prove the new error estimate in . The paper concludes with the fully-discrete setting presented in Section 4, where the time discretization is combined with a finite element discretization in space. We recall what is known about this discretization and finally prove corresponding -error estimates, which is the main result of this paper.

2 Notation and Assumptions

We start with introducing the analytical setting of this work. Throughout the paper we assume that (for ) is a convex bounded domain with polyhedral boundary. On , the Sobolev space of complex-valued, weakly differentiable functions with a zero trace on and -integratable partial derivatives is as usual denoted by . The potential is assumed to be real and nonnegative. Indirectly, we also assume that is sufficiently smooth so that it is compatible with the regularity assumptions for listed below (see [15] for a discussion on this aspect). The (possibly nonlinear) function

 γ:[0,∞)→[0,∞)

is assumed to be smooth, fulfills and its growth can be characterized with

 |γ(|v|2)v−γ(|w|2)w|≤L(K)|v−w|for % all v,w∈C with |v|,|w|≤K

where is a function with the following growth properties

 0≤L(s)≤Csqfor s≥0and {q∈[0,∞)for d=2,q∈[0,4)for d=3.

Note that in [15] the admissible growth condition in requires , which is however a typo and should be, as above, (cf. [10, Proposition 3.2.5 and Remark 3.2.7] for the original result). Examples for nonlinearities that fulfill these assumptions are mentioned in the introduction. The most common and physically relevant choices covered by our setting are power law nonlinearities for and in and in . Other physically relevant nonlinearities that fulfill the conditions are saturated nonlinearities appearing in the modeling of optical wave propagation such as for .

For the initial value we assume that and, without loss of generality, that it has a normalized mass, i.e. . With this, the considered nonlinear Schrödinger equation (NLS) reads as follows. For a maximum time and an initial value , we seek

 u∈L∞([0,T],H10(D))and∂tu∈L∞([0,T],H−1(D))

such that and

 i∂tu=−Δu+Vu+γ(|u|2)u (1)

in the sense of distributions. Problem (1) admits at least one solution, that is even unique for repulsive cubic nonlinearities in and (cf. [10] in general and [15, Remark 2.1] for precise references). We assume that the solution admits the following additional regularity, which is

 uttt∈L2((0,T),L2(D))and∂(k)tu∈L2(0,T;H2(D))for 0≤k≤2, (2)

where we note that any solution with such increased regularity must be unique (cf. [15, Lemma 3.1]). In the rest of the paper hence always refers to this uniquely characterized solution.

It is well known that solutions to the NLS (1) preserve the mass, i.e.

 ∫D|u(t,x)|2dx=∫D|u0(x)|2dx=1

and the energy, i.e.

 E[u(t)]=E[u0],where E[u]:=12∫D|∇u(x)|2+V(x)|u(x)|2+Γ(|u(x)|2)dx,

with .

For brevity, we shall denote the -norm of a function by . The -inner product is denoted by . Here, denotes the complex conjugate of .

Throughout the paper we will use the notation , to abbreviate , where is a constant that only depends on , , , , and , but not on the discretization.

3 Time-discrete Crank-Nicolson scheme

In this section we will state the semi-discrete Crank-Nicolson scheme, recall its well-posedness and available stability bounds, and then use these results to prove optimal -error estimates in the Hilbert space setting. For that, let denote the final time of computation, the number of time-steps, and the time step size. By we shall mean . The exact solution at time shall be denoted by . We also introduce a short hand notation for discrete time derivatives which is and analogously .

3.1 Method formulation and main result

With the notation above, the semi-discrete Crank–Nicolson approximation to is given recursively as the solution (in the sense of distributions) to the equation

 iDτunτ=−Δun+12τ+Vun+12τ+Γ(|un+1τ|2)−Γ(|unτ|2)|un+1τ|2−|unτ|2un+12τ, (3)

where . The initial value is selected as . It is easily seen that the discretization conserves both mass and energy, i.e.

The scheme (3) is well-posed and admits a set of a priori error estimates. The properties are summarized in the following theorem that is proved in [15, Theorem 4.1].

Theorem 3.1.

Under the general assumptions of this paper, there exists a constant and a solution to the semi-discrete Crank-Nicolson scheme (3) that is uniquely characterized by the property that

 sup0≤n≤N(∥unτ∥L∞(D)+∥unτ∥H2(D))≤C(u), (4)

and the a priori estimate for the -error

 sup0≤n≤N∥unτ−un∥≲τ2,

where is the (unqiue) exact solution with the regularity property (2).

Our main result on optimal error estimates in the reads as follows.

Theorem 3.2 (Optimal H1-error estimates for the semi-discrete method).

Consider the setting of Theorem 3.1, then the -error converges with optimal order in , i.e.

 sup0≤n≤N∥unτ−un∥H1(D)≲τ2.

The theorem is proved in Section 3.2 below.

3.2 Proof of Theorem 3.2

In this section we well prove Theorem 3.2. Let us introduce some notation that is used throughout the proofs. We recall . Furthermore, we let and . For time derivatives at fixed time , we also write .

We begin by establishing a differential equation for the time discrete error . This is stated in the following lemma.

Lemma 3.3.

The error fulfills the identity

 iDτen+Δen+1/2−Ven+1/2−enγ=Tn, (5)

where

 Tn=i(Dτun−∂tu(tn+1/2))+Δ(un+1/2−u(tn+1/2))−V(un+1/2−u(tn+1/2)) (6) −(Γ(|un+1|2)−Γ(|un|2)|un+1|2−|un|2−γ(|u(tn+1/2)|2)),

and for some bounded functions with the properties that

 ξn(x) ∈[min(|un|2,|un+1|2),max(|un|2,|un+1|2)]and ξnτ(x) ∈[min(|unτ|2,|un+1τ|2),max(|unτ|2,|un+1τ|2)]

for almost all .

Proof.

It is easily verified that exact solution fulfills

 iDτun+Δun+1/2−Vun+1/2−Γ(|un+1|2)−Γ(|un|2)|un+1|2−|un|2=Tn. (7)

By the regularity assumptions we can apply Taylor expansion arguments to to see:

 max0≤n≤N∥Tn∥≤Cτ2 (8)

Subtracting (7) from (3) we find that satisfies:

 iDτen+Δen+1/2−Ven+1/2−enγ=Tn

where denotes the error coming from the nonlinear term, defined by

 enγ=Γ(|un+1|2)−Γ(|un|2)|un+1|2−|un|2un+1/2−Γ(|un+1τ|2)−Γ(|unτ|2)|un+1τ|2−|unτ|2un+1/2τ.

Recalling the definition of we have:

 Γ(|un+1|2)−Γ(|un|2)|un+1|2−|un|2=1|un+1|2−|un|2∫|un+1|2|un|2γ(r)dr=:γ(ξn),

likewise

The expression for is thus simplified to

 enγ=γ(ξn)un+1/2−γ(ξnτ)un+1/2τ,

where is a function taking values between and and takes values between and . ∎

The differential equation in Lemma 3.3 is now used to derive a recurrence formula for the -norm of the error. Multiplying (5) by , integrating and taking the real part yields:

 ∥∇en+1∥2−∥∇en∥22τ=Re(⟨enγ,Dτen⟩)\rm I% −Re(⟨Tn,Dτen⟩).\rm II (9)

The idea is to bound the terms I and II in such a way that Grönwall’s inequality can be used. We proceed to bound term I. Multiplying the error PDE (5) by results in:

 i⟨Dτen,enγ⟩=⟨∇en+1/2,∇enγ⟩+⟨Ven+1/2,enγ⟩+∥enγ∥2+⟨Tn,enγ⟩

and consequently

 |\rm I|=|Re(⟨Dτen,enγ⟩)| ≤|Im(⟨∇en+1/2,∇enγ)⟩|+|Im⟨Ven+1/2,enγ⟩|+|Im(⟨Tn,enγ⟩)| ≲∥∇en+1/2∥2+∥∇enγ∥2+∥V∥∞(∥en+1/2∥2+∥enγ∥2)+τ4+∥enγ∥2 ≲∥∇en+1∥2+∥∇en∥2+∥∇enγ∥2+∥enγ∥+τ4. (10)

In order to use Grönwall’s inequality we need to bound and in terms of , and terms of . These bounds are formulated in the two following lemmas.

Lemma 3.4.

Given the optimal -convergence of Theorem 3.1 and the uniform bounds (4), the error coming from the nonlinear term behaves as , i.e. .

Proof.

We recall that and that is defined by

 γ(ξnτ)=1|un+1τ|2−|unτ|2∫|un+1τ|2|unτ|2γ(r)dr. (11)

Since is twice differentiable we have:

 γ(ξnτ)=γ(|un+1τ|2)+γ(|unτ|2)2+Rnτ

where . Multiplying (3) by and taking the imaginary part gives us together with Theorem 3.1 that

 ∥|un+1τ|2−|unτ|2∥ =τ∥Im(Δun+1/2τ¯¯¯un+1/2τ)∥τ (12)

and hence . The same argument, with the small difference that is bounded in virtue of the assumptions, can be applied to find

 γ(ξn)=γ(|un+1|2)+γ(|un|2)2+Rn. (13)

Using this and the standard mean value theorem we find that

 enγ=(γ(ξn)−γ(ξnτ))un+1/2+γ(ξnτ)(un+1/2−un+1/2τ) = (γ(|un+1|2)+γ(|un|2)2−γ(|un+1τ|2)+γ(|unτ|2)2)un+/2 +(Rn−Rnτ)un+1/2+γ(ξnτ)en+1/2 = Cγ′12(|un+1|2−|un+1τ|2)un+1/2+12Cγ′(|un|2−|unτ|2)un+1/2+γ(ξn)en+1/2+O(τ2).

Since it now follows that

Lemma 3.5.

Given Theorem 3.1, the gradient of the error coming from the nonlinear term is bounded as

Proof.

Starting from the equality

 ∇enγ=∇((γ(ξn)−γ(ξnτ))un+1/2+γ(ξnτ)en+1/2),

we recall that

 ∇1b(x)−a(x)∫b(x)a(x)γ(r)dr=∇bγ(b)−γ(ξ)b−a+∇aγ(ξ)−γ(a)b−a, (14)

where . Using formula (14) with and we find:

 (|un+1τ|2−|unτ|2)∇γ(ξnτ) = ∇|un+1τ|2(γ(|un+1τ|2)−γ(ξnτ))+∇|unτ|2(γ(ξnτ)−γ(|unτ|2)) = ∇|un+1τ|2(γ(|un+1τ|2)−γ(|un+1τ|2)+γ(|unτ|2)2−Rnτ) +∇|unτ|2(γ(|un+1τ|2)+γ(|unτ|2)2+Rnτ−γ(|unτ|2)) = 12(γ(|un+1τ|2)−γ(|unτ|2))(∇|un+1τ|2+∇|unτ|2)+Rnτ(∇(|un+1τ|2−|unτ|2)),

where the remainder is bounded as . It is straightforward to check that . Taking the gradient of (3), multiplying by and integrating ( is interpreted weakly, is therefore integrated by parts), gives us:

 ∥∇(un+1τ−unτ)(¯¯¯un+1τ+¯¯¯unτ)∥Cτ.

Recalling it becomes clear that

 ∥∇(|un+1τ|2−|unτ|2)∥≲τ. (15)

This means that . Thus,

 ∇γ(ξnτ)=12γ(|un+1τ|2)−γ(|unτ|2)|un+1τ|2−|unτ|2(∇|un+1τ|2+∇|unτ|2)+O(τ2). (16)

Using this and the mean value theorem we find:

 ∇γ(ξnτ)= 12γ′(ηnτ)(∇|un+1τ|2+∇|unτ|2)+O(τ2) (17) where γ′(ηnτ)=1|un+1τ|2−|unτ|2∫|un+1τ|2|unτ|2γ′(r)dr.

As the exact same steps are valid for , with the difference that is bounded in virtue of the assumptions, the term becomes:

 12(γ′(ηn)−γ′(ηnτ))∇(|un+1|2+|un|2)+12γ′(ηnτ)∇(|un+1|2−|un|2−|un+1τ|2+|unτ|2)+O(τ2).

Since is three times differentiable it follows that:

 γ′(ηn)−γ′(ηnτ)=γ′(|un+1|2+|un|22)−γ′(|un+1τ|2+|unτ|22)+Rnη

where the remainder is bounded as . Using the mean value theorem we find:

 γ′(ηn)−γ′(ηnτ)=Cγ′′(|un+1|2+|un|2−|un+1τ|2−|unτ|2)+Cγ′′′τ2.

Thus we are able to conclude:

 ||∇enγ||≲||∇en+1||2+||∇en||+τ2.

With Lemma 3.4 and 3.5 we now have the following bound on term I.

Lemma 3.6.

For term I which is given by (3.2), we have the estimate

 |\rm I|≲∥∇en+1∥2+∥∇en∥2+τ4. (18)

We can now proceed to bound term II. Here we explicate the Taylor term using (6) to see

 \rm II=−Re(⟨Tn,Dτen⟩)≤ |⟨(Dτun−∂tu(tn+1/2)),Dτen⟩|\rm IIa+Re⟨−Δ(un+1/2−u(tn+1/2)),Dτen⟩\rm IIb +|⟨V(un+1/2−u(tn+1/2)),Dτen⟩|\rm IIc +|⟨(Γ(|un+1|2)−Γ(|un|2)|un+1|2−|un|2−γ(|u(tn+1/2)|2),Dτen⟩)|\rm IId.

We start with estimating and IId, which can be bounded in a similar way.
Step 1, bounding IIa:
By replacing