A note on optimal $H^1$-error estimates for Crank-Nicolson approximations to the nonlinear Schr\"odinger equation

TL;DR

This paper establishes optimal $H^1$-error estimates for a mass- and energy-conserving Crank-Nicolson scheme applied to nonlinear Schr"odinger equations, filling a gap in theoretical analysis and providing practical implementation insights.

Contribution

The paper proves optimal $H^1$-error bounds for the Crank-Nicolson method in both semi-discrete and fully-discrete settings, and introduces an efficient fixed point iteration for solving the nonlinear system.

Findings

Proved optimal $L^{ty}(H^1)$-error estimates for the scheme.

Demonstrated the efficiency of the fixed point iteration.

Validated theoretical results with numerical experiments.

Abstract

In this paper we consider a mass- and energy--conserving Crank-Nicolson time discretization for a general class of nonlinear Schr\"odinger 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^{\infty}(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. We also suggest a fixed point iteration to solve the arising nonlinear system of equations that makes the method easy to implement and efficient. This is illustrated by numerical experiments.

Tables

Table 1. Table 1 : H 1 superscript 𝐻 1 H^{1} -errors for harmonic potential test case ( 24 ). The reference solution u ref subscript 𝑢 ref u_{\text{ref}} is calculated with h = 0.0125 ℎ 0.0125 h=0.0125 and τ = 2 − 13 𝜏 superscript 2 13 \tau=2^{-13} . The energy is E ​ [ u h , τ n ] = 3.86874 𝐸 delimited-[] subscript superscript 𝑢 𝑛 ℎ 𝜏 3.86874 E[u^{n}_{h,\tau}]=3.86874 .

	$\tau =2^{-5}$	$\tau =2^{-6}$	$\tau =2^{-7}$	$\tau =2^{-8}$	$\tau =2^{-9}$	$\tau =2^{-10}$	$\tau =2^{-11}$
$h=0.2$	0.829	0.400	0.485	0.526	0.538	0.542	0.543
$h=0.1$	1.048	0.396	0.144	0.179	0.191	0.194	0.194
$h=0.05$	1.109	0.498	0.129	0.054	0.062	0.066	0.066
$h=0.025$	1.124	0.526	0.155	0.037	0.022	0.023	0.023
$h=0.0125$	1.128	0.533	0.163	0.039	0.009	0.002	0.0005

Table 2. Table 2 : Errors and computational times for the CN-FEM. For a relative comparison of the errors we recall ‖ | u τ , h | 2 ‖ L 1 = 1 subscript norm superscript subscript 𝑢 𝜏 ℎ 2 superscript 𝐿 1 1 \|\ |u_{\tau,h}|^{2}\ \|_{L^{1}}=1 and note E ​ [ u τ , h ] = 25.539397 𝐸 delimited-[] subscript 𝑢 𝜏 ℎ 25.539397 E[u_{\tau,h}]=25.539397 .

CN-FEM, $h=0.025$, $N_{DoF}={400}^2$
	${\Vert {\left|u_{\tau ,h}\right|}^2-{\left|u_{\text{ref}}\right|}^2\Vert }_{L^1}$	${\Vert \nabla \left(u_{h,\tau }-u_{\text{ref}}\right)\Vert }_{L^2}$	CPU [h]
$\tau =2^{-8}$	0.40	4.91	0.4
$\tau =2^{-9}$	0.18	3.20	0.8
$\tau =2^{-10}$	0.12	1.43	1.4

Table 3. Table 3 : Errors and computational times for the SP2. For a relative comparison of the errors we recall ‖ | u τ , h | 2 ‖ L 1 = 1 subscript norm superscript subscript 𝑢 𝜏 ℎ 2 superscript 𝐿 1 1 \|\ |u_{\tau,h}|^{2}\ \|_{L^{1}}=1 and note E ​ [ u τ , h ] = 25.539397 𝐸 delimited-[] subscript 𝑢 𝜏 ℎ 25.539397 E[u_{\tau,h}]=25.539397 .

SP2, $N_{DoF}={800}^2$
	${\Vert {\left|u_{\tau ,h}\right|}^2-{\left|u_{\text{ref}}\right|}^2\Vert }_{L^1}$	${\Vert \nabla \left(u_{h,\tau }-u_{\text{ref}}\right)\Vert }_{L^2}$	CPU [h]
$\tau =2^{-12}$	0.28	4.29	0.24
$\tau =2^{-13}$	0.28	2.42	0.49
$\tau =2^{-14}$	0.28	2.35	0.96

Equations

$$\mathrm{i}{\partial }_{t}u=-△u+Vu+\gamma (|u{|}^2)u$$

$$\gamma :[0,\infty )\to [0,\infty )$$

$$|\gamma (|v{|}^2)v-\gamma (|w{|}^2)w|\le L(K)|v-w|\text{mbox}forallv,w\in \mathbb{C}\text{mbox}with|v|,|w|\le K$$

$$0\le L(s)\le C{s}^q\text{mbox}fors\ge 0\text{mbox}and\{\begin{array}{ll}q\in [0,\infty )& \text{mbox}ford=2,\\ q\in [0,4)& \text{mbox}ford=3.\end{array}$$

$$u\in L^{\infty }([0,T],H_0^1(\mathcal{D}))\text{mbox}and{\partial }_{t}u\in L^{\infty }([0,T],H^{-1}(\mathcal{D}))$$

$$\mathrm{i}{\partial }_{t}u=-\Delta u+Vu+\gamma (|u{|}^2)u$$

$$u_{ttt}\in L^2(0,T;H^1(\mathcal{D})),{\partial }_t^{(k)}u\in L^2(0,T;H^2(\mathcal{D}))\text{mbox}for0\le k\le 2$$

$$u\in L^{\infty }(0,T;W^{1,\infty }(\mathcal{D})),$$

$${\int }_{\mathcal{D}}|u(t,x){|}^{2}dx={\int }_{\mathcal{D}}|u^0(x){|}^{2}dx=1$$

$$E[u(t)]=E[u^0],\text{mbox}whereE[u]:=\frac{1}{2}{\int }_{\mathcal{D}}|\nabla u(x){|}^2+V(x)|u(x){|}^2+\Gamma (|u(x){|}^2)dx,$$

$$\mathrm{i}{D}_{\tau }{u}_{\tau }^n=-\Delta u_{\tau }^{n+\frac{1}{2}}+V{u}_{\tau }^{n+\frac{1}{2}}+\frac{\Gamma (|u_{\tau }^{n+1}{|}^2)-\Gamma (|u_{\tau }^n{|}^2)}{|u_{\tau }^{n+1}{|}^2-|u_{\tau }^n{|}^2}{u}_{\tau }^{n+\frac{1}{2}},$$

$${\int }_{\mathcal{D}}|u_{\tau }^n{|}^{2}dx={\int }_{\mathcal{D}}|u^0{|}^{2}dx\text{mbox}andE[u_{\tau }^n]=E[u^0]\text{mbox}foralln\ge 0.$$

$$\underset{0\le n\le N}{\sup }\left(\Vert u_{\tau }^n{\Vert }_{L^{\infty }(\mathcal{D})}+\Vert u_{\tau }^n{\Vert }_{H^2(\mathcal{D})}\right)\le C(u),$$

$$\underset{0\le n\le N}{\sup }\Vert u_{\tau }^n-u^n\Vert \lesssim {\tau }^2,$$

$$\underset{0\le n\le N}{\sup }\Vert u_{\tau }^n-u^n{\Vert }_{H^1(\mathcal{D})}\lesssim {\tau }^2.$$

$$\mathrm{i}{D}_{\tau }{e}^n+\Delta e^{n+1/2}-V{e}^{n+1/2}-e_{\gamma }^n=T^n,$$

$$T^n:=\mathrm{i}(D_{\tau }{u}^n-{\partial }_{t}u(t_{n+1/2}))+\Delta (u^{n+1/2}-u(t_{n+1/2}))-V(u^{n+1/2}-u(t_{n+1/2}))$$

$${\xi }^n(x)$$

$${\xi }_{\tau }^n(x)$$

$$\mathrm{i}{D}_{\tau }{u}^n+\Delta u^{n+1/2}-V{u}^{n+1/2}-\frac{\Gamma (|u^{n+1}{|}^2)-\Gamma (|u^n{|}^2)}{|u^{n+1}{|}^2-|u^n{|}^2}{u}^{n+1/2}=T^n.$$

$$\sum _{k=0}^N\Vert T^k{\Vert }^2\le C{\tau }^3$$

$$\mathrm{i}{D}_{\tau }{e}^n+\Delta e^{n+1/2}-V{e}^{n+1/2}-e_{\gamma }^n=T^n$$

$$e_{\gamma }^n=\frac{\Gamma (|u^{n+1}{|}^2)-\Gamma (|u^n{|}^2)}{|u^{n+1}{|}^2-|u^n{|}^2}{u}^{n+1/2}-\frac{\Gamma (|u_{\tau }^{n+1}{|}^2)-\Gamma (|u_{\tau }^n{|}^2)}{|u_{\tau }^{n+1}{|}^2-|u_{\tau }^n{|}^2}{u}_{\tau }^{n+1/2}.$$

$$\frac{\Gamma (|u^{n+1}{|}^2)-\Gamma (|u^n{|}^2)}{|u^{n+1}{|}^2-|u^n{|}^2}=\frac{1}{|u^{n+1}{|}^2-|u^n{|}^2}{\int }_{|u^n{|}^2}^{|u^{n+1}{|}^2}\gamma (r)dr=:\gamma ({\xi }^n),$$

$$\frac{\Gamma (|u_{\tau }^{n+1}{|}^2)-\Gamma (|u_{\tau }^n{|}^2)}{|u_{\tau }^{n+1}{|}^2-|u_{\tau }^n{|}^2}=\frac{1}{|u_{\tau }^{n+1}{|}^2-|u_{\tau }^n{|}^2}{\int }_{|u_{\tau }^n{|}^2}^{|u_{\tau }^{n+1}{|}^2}\gamma (r)dr=:\gamma ({\xi }_{\tau }^n).$$

$$e_{\gamma }^n:=\gamma ({\xi }^n)u^{n+1/2}-\gamma ({\xi }_{\tau }^n)u_{\tau }^{n+1/2},$$

$$\frac{\Vert \nabla e^{n+1}{\Vert }^2-\Vert \nabla e^n{\Vert }^2}{2\tau }=\underset{\text{mbox}\mathrm{I}}{\underbrace{\text{Re}(⟨e_{\gamma }^n,D_{\tau }{e}^n⟩)}}\underset{\text{mbox}\mathrm{II}}{\underbrace{-\text{Re}(⟨T^n,D_{\tau }{e}^n⟩).}}$$

$$\mathrm{i}⟨D_{\tau }{e}^n,e_{\gamma }^n⟩=⟨\nabla e^{n+1/2},\nabla e_{\gamma }^n⟩+⟨V{e}^{n+1/2},e_{\gamma }^n⟩+\Vert e_{\gamma }^n{\Vert }^2+⟨T^n,e_{\gamma }^n⟩$$

$$|\text{mbox}\mathrm{I}|=|\text{Re}(⟨D_{\tau }{e}^n,e_{\gamma }^n⟩)|$$

$$f(a,b)=\frac{1}{b-a}{\int }_a^b\gamma (r)dr.$$

$$e_{\gamma }^n=(\gamma ({\xi }^n)-\gamma ({\xi }_{\tau }^n))u^{n+1/2}+\gamma ({\xi }_{\tau }^n)(u^{n+1/2}-u_{\tau }^{n+1/2})$$

Full text

AMS subject classifications

35Q55, 65M60, 65M15, 81Q05

1 Introduction

In this paper we consider nonlinear Schrödinger equations (NLS) seeking a complex function $u(t,x)$ such that

[TABLE]

in a bounded domain $\mathcal{D}\subset\mathbb{R}^{d}$, with a homogenous Dirichlet boundary condition on $\partial\mathcal{D}$ and a given initial value. Here, $V(x)$ is a known real-valued potential and $\gamma:[0,\infty)\rightarrow\mathbb{R}$ is a smooth (and possibly nonlinear) function that depends on the unknown density $|u|^{2}$. Of particular interest are cubic nonlinearities of the form $\gamma(|u|^{2})u=\kappa|u|^{2}u$, for some $\kappa\in\mathbb{R}$. In this case, the equation is called Gross–Pitaevskii equation. It has applications in optics [1, 19], fluid dynamics [38, 39] and, most importantly, in quantum physics, where it models for example the dynamics of Bose-Einstein condensates in a magnetic trapping potential [17, 27, 31]. Another relevant class is saturated nonlinearities, such as $\gamma(|u|^{2})=\kappa|u|^{2}(1+\alpha|u|^{2})^{-1}$ for some $\alpha\geq 0$, which appear in the context of nonlinear optical wave propagation in layered metallic structures [15, 23] or the propagation of light beams in plasmas [29]. In order to discretize nonlinear Schrödinger equations in time, splitting methods and exponential integrators typically yield highly efficient solution schemes that can be easily combined with a spectral discretization in space (cf. [9, 12, 14, 16, 18, 28, 35, 34] and the references therein). If the exact solution to the NLS admits high regularity, such discretization schemes typically show a remarkably good performance. However, if the regularity of the solution is strongly reduced, either by rough potentials $V$ (e.g. disorder potentials or optical lattices) or by rough initial values $u(0)$ (e.g. when effects close to phase transitions are studied), then the performance of these methods can drop dramatically. Here we refer exemplarily to the recent numerical experiments reported in [30, 26, 22]. To overcome this issue, Ostermann and Schratz proposed new low-regularity time-integrators [30, 26] which improve the convergence in low regularity regimes significantly. However, the approach still relies on a Fourier discretization in space, which is not an optimal choice due to the loss of spectral convergence for non-smooth solutions. Practically, the usage of a (low order) finite element space discretization is often desirable in order to account for spatial low regularity. In the following we will only discuss approaches that can be easily combined with finite elements in space, meaning that we put ourselves into the situation that we assume that the solution to the NLS does not admit much smoothness.

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 [22, 33], where it was also found that the complexity of the physical setup (or low-regularity) can stress this issue even further.

For the subclass of power law nonlinearities of the form $\gamma(|u|^{2})=\sum_{k=1}^{K}\kappa_{k}|u|^{2\sigma_{k}}$ for $\sigma_{k}\geq 0$ and $\alpha_{k}\in\mathbb{R}$, a mass and energy conserving relaxation scheme was proposed and analyzed by Besse [8, 10]. Thanks to its properties, the scheme shows a very good performance in realistic physical setups [22]. Despite the large variety of different numerical approaches for solving the time-dependent NLS (cf. [2, 3, 5, 20, 24, 25, 28, 32, 35, 34, 36, 37, 40] and the references therein) the literature knows however of 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 [32] and which is long-known in the physics community, is a Crank–Nicolson-type (CN) approach where the nonlinearity is approximated by a suitable difference quotient involving the primitive integral of $\gamma$. This is also the time discretization that we shall consider in this paper. Here we note it was analytically and numerically demonstrated that this method is applicable and reliable in low-regularity regimes [22, 21].

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 $P1$ finite element discretization in space, the first a priori error estimates for the arising method were obtained by Sanz-Serna 1984 [32] for cubic nonlinearities. He considers the case $d=1$ and derives optimal $L^{\infty}(L^{2})$-error estimates under the coupling constraint $\tau\lesssim h$, where $\tau$ denotes the time step size of the Crank-Nicolson method and $h$ 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 $L^{\infty}(L^{2})$ in dimension $d=1,2,3$ and under the relaxed coupling constraint $\tau\lesssim h^{d/4}$. Finally, in 2017 [21], the $L^{\infty}(L^{2})$-error estimates could be improved yet another time by showing that the coupling constraint can be fully removed. Furthermore, general nonlinearities could be considered, the influence of potentials could be taken into account and even convergence under weak regularity assumptions could be proved (with reduced convergence rates). However, so far, optimal error estimates in $L^{\infty}(H^{1})$ for this particular CN-discretization are still open in the literature.

One reason for this absence of $H^{1}$-results could be related to the techniques used for the error analysis in previous works (cf. [2, 20, 24, 25, 32, 40]) 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 $L^{\infty}(L^{2})$- and/or $L^{\infty}(H^{1})$-error estimates. 3. Using inverse estimates in the finite element space to show that the truncated approximations are uniformly bounded in $L^{\infty}(L^{\infty})$ by a term of the form $C(1+h^{-s}(\tau^{2}+h^{p}))$, with appropriate powers $p>s>0$ that depend on the considered space discretization, regularity and space dimension. 4. Concluding that if $\tau$ and $h$ are coupled in an appropriate way, then the truncated approximations are all uniformly bounded by a constant $C$ 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 $L^{\infty}(H^{1})$-error estimates for the Crank–Nicolson FEM. This is because it requires a suitable truncation of the primitive integral of $\gamma$ that is on the one hand consistent with the energy conservation and on the other hand allows for uniform bounds of the approximations in $L^{\infty}(W^{1,\infty})$. However, thanks to the new techniques developed in [37] and the CN error analysis suggested in [21] in the context of $L^{\infty}(L^{2})$-error estimates, the truncation step is no longer necessary and the desired $L^{\infty}(L^{\infty})$-bounds can be derived with elliptic regularity theory. With this, it is now possible to obtain $L^{\infty}(H^{1})$ 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 [21, 37] to fill the gap in the literature and prove optimal $L^{\infty}(H^{1})$-error estimates for the energy-conservative 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 $L^{\infty}(L^{2})$. Furthermore, we present and prove the new error estimate in $L^{\infty}(H^{1})$. The paper continues 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 $L^{\infty}(H^{1})$-error estimates, which is the main result of this paper. The paper concludes with a note on how to efficiently implement the method and two numerical experiments to confirm the convergence rates and to illustrate a setting in which it makes computational sense to use the CN-FEM instead of for example a spectral method.

2 Notation and Assumptions

We start with introducing the analytical setting of this work. Throughout the paper we assume that ${\mathcal{D}}\subset\mathbb{R}^{d}$ (for $d=2,3$) is a convex bounded domain with polyhedral boundary. On $\mathcal{D}$, the Sobolev space of complex-valued, weakly differentiable functions with a zero trace on $\partial\mathcal{D}$ and $L^{2}$-integratable partial derivatives is as usual denoted by $H^{1}_{0}(\mathcal{D}):=H^{1}_{0}(\mathcal{D},\mathbb{C})$. The potential $V\in L^{\infty}(\mathcal{D};\mathbb{R})$ is assumed to be real and nonnegative. Indirectly, we also assume that $V$ is sufficiently smooth so that it is compatible with the regularity assumptions for $u$ listed below (see [21] for a discussion on this aspect). The (possibly nonlinear) function

[TABLE]

is assumed to be $C^{2}[0,\infty)$, fulfills $\gamma(0)=0$ and its growth can be characterized with

[TABLE]

where $L$ is a function with the following growth properties

[TABLE]

Note that in [21] the admissible growth condition in $3d$ requires $q\in[0,2)$, which is however a typo and should be, as above, $q\in[0,4)$ (cf. [13, Proposition 3.2.5 and Remark 3.2.7] for the original result). Examples of 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 $\gamma(\rho)=\kappa\rho^{q}$ for $\kappa\geq 0$ and $0\leq q<\infty$ in $2d$ and $0\leq q<4$ in $3d$. Other physically relevant nonlinearities that fulfill the conditions are saturated nonlinearities appearing in the modeling of optical wave propagation such as $\gamma(\rho)=\kappa\rho(1+\alpha\rho)^{-1}$ for $\alpha,\kappa\geq 0$.

The above assumptions cover the regime of so-called defocussing (positive) nonlinearities and guarantees that the NLS and its Crank-Nicolson discretization are well-posed. For focussing (negative) nonlinearities, i.e. $\gamma:[0,\infty)\mapsto(-\infty,0]$, the well-posedness (of both the continuous and discrete models) can no longer be guaranteed without making additional technical assumptions. Typically, effects such as finite time blow ups can occur in this regime. To avoid constantly having to invoke a saving clause we restrict our attention to the defocussing case. We do however point out that under the assumptions that the NLS and the Crank-Nicolson discretizations are well-posed (without blow-up in the time interval $[0,T]$) then all our error estimates hold without changes.

For the initial value we assume that $u^{0}\in H^{1}_{0}(\mathcal{D})\cap H^{2}(\mathcal{D})$ and, without loss of generality, that it has a normalized mass, i.e. $\int_{\mathcal{D}}|u^{0}(x)|^{2}\hskip 2.0ptdx=1$. With this, the considered nonlinear Schrödinger equation (NLS) reads as follows. For a maximum time $T>0$ and an initial value $u^{0}$, we seek

[TABLE]

such that $u(\cdot,0)=u^{0}$ and

[TABLE]

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

[TABLE]

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

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

[TABLE]

and the energy, i.e.

[TABLE]

with $\Gamma(\rho):=\int_{0}^{\rho}\gamma(r)\hskip 2.0ptdr$.

For brevity, we shall denote the $L^{2}$-norm of a function $v\in L^{2}(\mathcal{D}):=L^{2}(\mathcal{D},\mathbb{C})$ by $\|v\|$. The $L^{2}$-inner product is denoted by $\langle v,w\rangle=\int_{\mathcal{D}}v(x)\hskip 2.0pt\overline{w(x)}\hskip 2.0ptdx$. Here, $\overline{w}$ denotes the complex conjugate of $w$.

Throughout the paper we will use the notation $A\lesssim B$, to abbreviate $A\leq CB$, where $C$ is a constant that only depends on $u$, $T$, $d$, $\mathcal{D}$, $V$ and $\gamma$, but not on the discretization.

Remark 2.1**.**

In the analysis we restrict our attention to homogeneous Dirichlet boundary conditions. Typically these boundary conditions can be motivated by physical reasoning. For example in the context of Bose Einstein condensates, the magnetic potential $V$ is a trapping potential that becomes very quickly very large and hence traps the condensate in a bounded region. Mathematically this leads to an exponential decay of the solution $u$ to zero (in moderate distances from the origin of the coordinate system) and hence justifies to truncate the computational domain to a simple geometric object on which the problem is solved with zero boundary conditions. A typical alternative found in the literature are periodic boundary conditions which are e.g. favorable for spectral methods. Both the formulation of the Crank-Nicolson method and its error analysis can be easily generalized to that case.

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 $L^{\infty}(H^{1})$-error estimates in the Hilbert space setting. For that, let $T$ denote the final time of computation, $N$ the number of time-steps, and $\tau=T/N$ the time step size. By $t_{n}$ we shall mean $t_{n}=n\tau$. The exact solution at time $t_{n}$ shall be denoted by $u^{n}:=u(t_{n},\cdot)$. We also introduce a short hand notation for discrete time derivatives which is $D_{\tau}u^{n}:=(u^{n+1}-u^{n})/\tau$ and analogously $D_{\tau}u_{\tau}^{n}:=(u_{\tau}^{n+1}-u_{\tau}^{n})/\tau$.

3.1 Method formulation and main result

With the notation above, the semi-discrete Crank–Nicolson approximation $u_{\tau}^{n+1}\in H^{1}_{0}(\mathcal{D})$ to $u^{n+1}$ is given recursively as the solution (in the sense of distributions) to the equation

[TABLE]

where $u_{\tau}^{n+\frac{1}{2}}:=(u_{\tau}^{n}+u_{\tau}^{n+1})/2$. The initial value is selected as $u_{\tau}^{0}=u^{0}$. It is easily seen that the discretization conserves both mass and energy, i.e.

[TABLE]

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

Theorem 3.1**.**

Under the general assumptions of this paper, there exists a constant $C(u)>0$ and a solution $u_{\tau}^{n}\in H^{1}_{0}(\mathcal{D})$ to the semi-discrete Crank-Nicolson scheme (4) that is uniquely characterized by the property that

[TABLE]

and the a priori estimate for the $L^{2}$-error

[TABLE]

where $u$ is the (unique) exact solution with the regularity property (2).

Our main result on optimal error estimates in the $L^{\infty}(H^{1})$ reads as follows.

Theorem 3.2** (Optimal $H^{1}$-error estimates for the semi-discrete method).**

Consider the setting of Theorem 3.1, then the $L^{\infty}(H^{1})$-error converges with optimal order in $\tau$, i.e.

[TABLE]

The theorem is proved in Section 3.2 below.

3.2 Proof of Theorem 3.2

In this section we will prove Theorem 3.2. Let us introduce some notation that is used throughout the proofs. We recall $D_{\tau}e^{n}=(e^{n+1}-e^{n})/\tau$. Furthermore, we let $e^{n+1/2}:=(e^{n+1}+e^{n})/2$ and $u^{n+1/2}:=(u^{n+1}+u^{n})/2$. For time derivatives at fixed time $t^{n}$, we also write $\partial_{t}u^{n}:=\partial_{t}u(t^{n},\cdot)$.

We begin by establishing a differential equation for the time discrete error $e^{n}=u^{n}-u^{n}_{\tau}$. This is stated in the following lemma.

Lemma 3.3** (Consistency error).**

The error e^{n}=u^{n}-u_{\tau}^{n}\ fulfills the identity

[TABLE]

where the consistency error $T^{n}$ is given by

[TABLE]

Here, $e^{n}_{\gamma}:=\gamma(\xi^{n})u^{n+1/2}-\gamma(\xi^{n}_{\tau})u^{n+1/2}_{\tau}$ for some bounded functions $\xi^{n},\xi^{n}_{\tau}\in L^{\infty}(\mathcal{D})$ with the properties that

[TABLE]

for almost all $x\in\mathcal{D}$.

Proof.

It is easily verified that exact solution fulfills

[TABLE]

By the regularity assumptions we can apply Taylor expansion arguments to $T^{n}$ to see:

[TABLE]

The argument that proves (9) is elaborated in Appendix A, where it also becomes visible how the regularity assumptions enter explicitly in the estimate. Next, subtracting (8) from (4) we find that $e^{n}=u^{n}-u_{\tau}^{n}$ satisfies:

[TABLE]

where $e^{n}_{\gamma}$ denotes the error coming from the nonlinear term, defined by

[TABLE]

Recalling the definition of $\Gamma$ we have:

[TABLE]

likewise

[TABLE]

The expression for $e^{n}_{\gamma}$ is thus simplified to

[TABLE]

where $\xi^{n}$ is a function taking values between $|u^{n}|^{2}$ and $|u^{n+1}|^{2}$ and $\xi^{n}_{\tau}$ a function taking values between $|u^{n}_{\tau}|^{2}$ and $|u_{\tau}^{n+1}|^{2}$. ∎

The differential equation in Lemma 3.3 is now used to derive a recurrence formula for the $H^{1}$-norm of the error. Multiplying (6) by $D_{\tau}e^{n}$, integrating and taking the real part yields:

[TABLE]

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 (6) by $e^{n}_{\gamma}$ results in:

[TABLE]

and consequently

[TABLE]

In order to use Grönwall’s inequality we need to bound $\|e^{n}_{\gamma}\|$ and $\|\nabla e^{n}_{\gamma}\|$ in terms of $\|e^{n}\|$,$\|\nabla e^{n}\|$ and terms of $\mathcal{O}(\tau^{2})$. These bounds are formulated in the two following lemmas.

Lemma 3.4**.**

Given the optimal $L^{2}$-convergence of Theorem 3.1 and the uniform bounds (5), the error coming from the nonlinear term behaves as $\tau^{2}$, i.e. $\|e^{n}_{\gamma}\|\lesssim\tau^{2}$.

Proof.

We introduce the function $f$ to denote how $\gamma(\xi^{n})$ depends on $|u^{n}|^{2}$ and $|u^{n+1}|^{2}$

[TABLE]

The derivative of $f$ with respect to to its $i$-th input is denoted $f_{i}$, i.e. $f_{1}(a,b)=\partial_{a}f(a,b)$. A standard application of the mean value theorem yields that, for some function $\eta^{n}$ taking values between $|u^{n}|^{2}$ and $|u^{n}_{\tau}|^{2}$ a.e. and for some function $\eta^{n+1}$ likewise between $|u^{n+1}|^{2}$ and $|u^{n+1}_{\tau}|^{2}$, it holds:

[TABLE]

A quick sanity check shows that $f_{1}$ and $f_{2}$ are bounded by the derivative of $\gamma$:

[TABLE]

where $c,\theta^{-}$ and $\theta^{+}$ lie somewhere between $a$ and $b$. With the $L^{\infty}$-bounds on $u^{n}$ and $u^{n}_{\tau}$ it is now straightforward to show $\|e^{n}_{\gamma}\|\leq C_{\gamma^{\prime}}(\|e^{n}\|+\|e^{n+1}\|)$:

[TABLE]

As $\||u^{n}|^{2}-|u^{n}_{\tau}|^{2}\|\leq\||u^{n}|+|u^{n+1}|\|_{L^{\infty}}\|e^{n}\|$, it now follows that

[TABLE]

and the lemma is proved. ∎

Lemma 3.5**.**

Given Theorem 3.1, the gradient of the error coming from the nonlinear term is bounded as $\|\nabla e^{n}_{\gamma}\|\lesssim\|\nabla e^{n+1}\|+\|\nabla e^{n}\|+\tau^{2}.$

Proof.

The steps are much the same as in lemma (3.4), with the exception that we need to use $W^{1,\infty}$-bounds on $u^{n}$ (3), which are not available for $u^{n}_{\tau}$. We begin by splitting $\nabla e^{n}_{\gamma}$ into terms so that the previous lemma may be used.

[TABLE]

By the previous lemma and the $W^{1,\infty}$-bound on $u^{n}$ we may conclude,

[TABLE]

What is left to bound is the term $\|\nabla(\gamma(\xi^{n})-\gamma(\xi^{n}_{\tau}))\|$. We consider its dependence on $|u^{n}|^{2},|u^{n+1}|^{2},|u^{n}_{\tau}|^{2}$ and $|u^{n+1}_{\tau}|^{2}$:

[TABLE]

Where $f_{1}$ is to be read as $f_{1}(|u^{n}|^{2},|u^{n+1}|^{2})$ and $f_{1,\tau}$ as $f_{1}(|u^{n}_{\tau}|^{2},|u^{n+1}_{\tau}|^{2})$ and likewise for $f_{2}$ and $f_{2,\tau}$. Another application of the mean value theorem yields:

[TABLE]

for some $\theta^{n}$, $\vartheta^{n}$ between $|u^{n}|^{2}$ and $|u^{n}_{\tau}|^{2}$ and some $\theta^{n+1}$, $\vartheta^{n+1}$ between $|u^{n+1}|^{2}$ and $|u^{n+1}_{\tau}|^{2}$. The following quick calculations show that the partial derivatives of $f$ of order two are bounded by $\gamma^{\prime\prime}$.

[TABLE]

Where it was used that $c=(b+a)/2+C_{\gamma^{\prime\prime}}(a-b)^{2}$. It thus becomes clear that $\|f_{1,1}+f_{1,2}+f_{2,2}\|_{L^{\infty}}\leq C_{\gamma^{\prime\prime}}$. This gives us the following $L^{2}$-bound on (3.5)

[TABLE]

Continuing from (12), we now conclude that:

[TABLE]

It is noted that $\nabla(|u^{n}|^{2}-|u^{n}_{\tau}|^{2})$ may be written

[TABLE]

Using the $W^{1,\infty}$-bound available for $u^{n}$ we have that $\|\nabla(|u^{n}|^{2}-|u^{n}_{\tau}|^{2})\|\lesssim\|\nabla(u^{n}-u^{n}_{\tau})\|$. With this, eq. (14) becomes $\|\nabla e^{n}_{\gamma}\|\lesssim\|\nabla e^{n}\|+\|\nabla e^{n+1}\|+\tau^{2}$ and the lemma is proved. ∎

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 (11), we have the estimate

[TABLE]

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

[TABLE]

We start with estimating $\mbox{\rm IIa},\mbox{\rm IIc}$ and IId, which can be bounded in a similar way.

*Step 1, bounding IIa:

*By replacing $D_{\tau}e^{n}$ using (6) (i.e. time derivative is replaced by regularity in space) we have

[TABLE]

The term $\|D_{\tau}u^{n}-\partial_{t}u(t_{n+1/2})\|_{L^{2}}^{2}$ was absorbed in $\lesssim\|T^{n}\|^{2}$. Here we see how the assumption that $\partial_{ttt}u\in L^{2}(0,T;H^{1}(\mathcal{D}))$ will enter, as it will allow us to conclude that $\sum\|D_{\tau}u^{n}-\partial_{t}u(t_{n+1/2})\|_{H^{1}}^{2}\lesssim\tau^{3}$.

*Step 2, bounding IIc:

*We use the same idea as for IIa to get

[TABLE]

*Step 3, bounding IId:

*We start from

[TABLE]

and replace $D_{\tau}e^{n}$ again using (6). Furthermore, in virtue of the assumptions it holds that $\|\nabla(\gamma(\xi^{n})-\gamma(|u(t_{n+1/2})|^{2}))\|\lesssim\|u^{n+1/2}-u(t_{n+1/2})\|_{H^{1}}$, this is made explicit in the Appendix (A). We thus obtain

[TABLE]

*Step 4, bounding IIb:

*The previous technique does not work on this term since replacing the discrete time derivative with regularity in space would give rise to the term $\nabla\Delta(u^{n+1/2}-u(t_{n+1/2}))$, which we can not afford. Instead we use summation by parts in time to get the factor $D_{\tau}\Delta(u^{n+1/2}-u(t_{n+1/2}))$, which when integrated against $e^{n+1/2}$ can be handled. First we recall:

[TABLE]

Using this on term IIb yields:

[TABLE]

Collecting the estimates we have the following estimate for term II.

Lemma 3.7**.**

For term $\mbox{\rm II}=-\text{Re}(\langle T^{n},D_{\tau}e^{n}\rangle)$ it holds the estimate

[TABLE]

Here we note the importance of not estimating the absolute value of the first term since it is necessary to use the fact that $n$ of these terms cancel when summed up, i.e. $\sum_{k}D_{\tau}a^{k}=\frac{1}{\tau}(a^{n+1}-a^{0})$. We are now ready to finish the proof of the first main result.

Proof of Theorem 3.2.

We pick off where we left (10) and find by using Lemma 3.6 and 3.7:

[TABLE]

Summing this up and using $e^{0}=0$ gives

[TABLE]

and therefore, recalling (9),

[TABLE]

Young’s inequality with $\epsilon>0$ is used on the last term:

[TABLE]

Which holds since,

[TABLE]

where we have $\|\partial_{tt}u\|_{L^{\infty}(H^{1})}\lesssim\|\partial_{tt}u\|_{L^{2}(H^{1})}+\|\partial_{ttt}u\|_{L^{2}(H^{1})}$ by Sobolev embeddings. Finally we arrive at

[TABLE]

and for e.g. $\epsilon=1/2$ we can absorb $\epsilon\|\nabla e^{n+1}\|^{2}$ in the left hand side and conclude

[TABLE]

Grönwall’s inequality now yields:

[TABLE]

∎

4 Fully-discrete Crank-Nicolson scheme

We shall now consider the fully-discrete setting that is based on a finite element discretization in space. For that, we let $S_{h}\subset H^{1}_{0}(\mathcal{D})$ denote the space of P1 Lagrange finite elements on a quasi-uniform simplicial mesh on $\mathcal{D}$ with mesh size $h$. In this setting we have by standard finite element theory (cf. [11]) the following estimate for any $u\in H^{1}_{0}(\mathcal{D})\cap H^{2}(\mathcal{D})$:

[TABLE]

Here $P_{h}:H^{1}_{0}(\mathcal{D})\rightarrow S_{h}$ denotes (for example) the Ritz-projection into the finite element space. For a given discrete initial value $u_{h}^{0}\in S_{h}$ the CN-FEM approximation $u^{n}_{\tau,h}\in S_{h}$ to $u^{n}$ is given by the fully discrete equation

[TABLE]

for all $v\in S_{h}$. The initial value is selected as $u^{0}_{\tau,h}=P_{h}u^{0}$. As in the semi-discrete case, the discretization conserves both mass and energy, i.e.

[TABLE]

The scheme is well-posed and the corresponding approximations converge in the $L^{\infty}(L^{2})$-norm with optimal order in space and time to the exact solution. A proof of this statement can be easily extracted from [21, Theorem 3.1 and Lemma 5.3]. In particular, we have the following result.

Theorem 4.1**.**

Under the general assumptions of this paper, there exists a solution $u^{n}_{\tau,h}\in H^{1}_{0}(\mathcal{D})$ to the fully discrete Crank-Nicolson scheme (20) such that the following a priori error estimates hold

[TABLE]

With this we are ready to state our final theorem.

Theorem 4.2** (Optimal $H^{1}$-error estimates for the fully discrete method).**

Let $u_{\tau}^{n}\in H^{1}_{0}(\mathcal{D})$ denote the fully-discrete Crank-Nicolson approximations from Theorem 4.1, then it holds

[TABLE]

Proof.

First, we recall the inverse estimate on quasi-uniform meshes (cf. [11]), i.e. $\|\nabla v_{h}\|\leq Ch^{-1}\|v_{h}\|$ for all $v_{h}\in S_{h}$, which implies

[TABLE]

With this, the $H^{1}$ convergence result (18) together with Theorem 4.1 suffice to show optimal $H^{1}$-convergence rates for the fully discrete method. This is made clear by the following splitting.

[TABLE]

Here we have made use of the inequality (19), the uniform $H^{2}$-regularity of $u^{n}_{\tau}$, i.e. $\|u^{n}_{\tau}\|_{H^{2}}\leq C(u)$ (cf. (5)) and the optimal $L^{2}$-estimates. In virtue of Theorem 3.2 we may thus conclude:

[TABLE]

∎

Detailed numerical studies that confirm the optimal convergence rates stated in Theorem 4.1 and Theorem 4.2 are presented in [21, 22].

5 Implementation and Numerical Examples

In this section we will discuss how the Crank-Nicolson FEM discretization can be efficiently implemented and practically used. Afterwards, we present two numerical experiments. The first one is to confirm the theoretically predicted convergence rates in Theorem 4.1 and the second experiment demonstrates that our approach is fully competitive in low regularity regimes, where we compare it with a time-splitting spectral method.

5.1 Efficient implementation

The Crank Nicolson method (20), albeit popular, suffers from the drawback that it requires solving a fully nonlinear system of equations in each time step. Furthermore, this system of equations is often solved through a Newton step, the implementation of which can become complicated and expensive for general nonlinearities. We present here a competitive fixed point solver which makes the method perform on par in terms of computational time with linearized time-stepping methods such as the RE-FEM proposed by C. Besse [8] which was found to be best performing in [22].

To detail the proposed fixed-point iteration, let $U^{n}\in\mathbb{R}^{N}$ denote the vector of nodal values that belongs to the function $u^{n}_{\tau,h}\in S_{h}$. Introducing the following matrix notation:

[TABLE]

the equation (20) in matrix form becomes :

[TABLE]

Let $L_{1}=M+\mathrm{i}\tau/2(A+M_{V})$ and $L_{2}=M-\mathrm{i}\tau/2(A+M_{V})$. Our fixed point iteration takes the form:

[TABLE]

Here we note that matrix $L_{1}$ does not change with time. Hence, the above iteration can be done efficiently by precomputing the LU-factorization of $L_{1}$. After it is precomputed, each time step only involves matrix-vector multiplications, but no longer the solving of a linear system of equations. In fact, the main cost in each time step account for the assembly of the updated mass matrix $M_{\Gamma}(U^{n+1}_{i},U^{n})$ with the densities from the previous time step and the previous iteration. Compared to this, all other costs are essentially negligible. We find that typically, but dependent on $\tau$, 4-8 iterations are required to reach a tolerance of machine epsilon. To illustrate the efficiency we conclude with two numerical test problems.

5.2 Harmonic potential

First we consider a smooth potential to confirm the expected convergences rates for a type of nonlinearity that complements previous test cases [21, 22]. Here we seek $u(x,t)$ with

[TABLE]

where we consider the saturated nonlinearity $\gamma(r):=r/(1+r)$ (cf. [29, 15, 23]) . Furthermore, $\mathcal{D}=[-5,5]^{2}$ is the computational domain, the maximum time is selected as $T=1$ and the trapping potential $V(x,y)=(\nu_{x}x)^{2}+(\nu_{y}y)^{2}$. For the time-dependent problem we set the trapping frequencies to $\nu_{x}=2$ and $\nu_{y}=3$. The initial value $u_{0}$ is the unique positive ground state with $\int_{\mathcal{D}}|u_{0}|^{2}=1$ to the problem with $\nu_{x}=\nu_{y}=1$, i.e. it solves the eigenvalue problem

[TABLE]

with ground state eigenvalue (chemical potential) $\lambda_{0}>0$. The $H^{1}$-errors are presented in Table 1. The $\mathcal{O}(h)$-convergence is best seen in column $\tau=2^{-9}$, where initially the convergence is $\mathcal{O}(h^{1.5})$ but flattens out to $\mathcal{O}(h^{1.2})$ for the last data point. Since the reference solution is also computed with $h=0.0125$, it is expected that this last order of convergence is an overestimate. Using the values in row $h=0.0125$ we estimate the order of convergence with respect to $\tau$ to be $1.9$, hence, confirming the theoretically predicted rates from Theorem 4.1.

5.3 Discontinuous potential

This example illustrates a moderate setting where, due to reduced regularity of the exact solution, finite element based methods are preferable over spectral methods. In the following we compare the Crank-Nicolson approach with a Strang splitting spectral method of order 2 (SP2) [7] which is known to show a very good performance in smooth settings.

In this test problem we seek $u(x,t)$ with

[TABLE]

where we consider the saturated nonlinearity $\gamma(r):=10r/(1+r)$. For a fair comparison with the SP2, we consider our problem with periodic boundary conditions which are easier to handle by the spectral method. The generalization of the Crank-Nicolson method to periodic boundary conditions is straightforward. Furthermore, $\mathcal{D}=[-5,5]^{2}$ is again the computational domain and the maximum time is selected as $T=1$. The trapping potential $V(x,y)=(\nu_{x}x)^{2}+(\nu_{y}y)^{2}+100(\mathds{1}_{|x|\geq 1}(x)+\mathds{1}_{|y|\geq 1}(y))$, with trapping frequencies $\nu_{x}=1$ and $\nu_{y}=3$, is discontinuous and causes a slight loss of regularity. We stress that this is a moderate test case, as illustrated in Fig. 1 most of the dynamics take place within the unit cube where the potential is smooth. The initial value $u_{0}$ is the unique positive ground state with $\int_{\mathcal{D}}|u_{0}|^{2}=1$ and $V_{0}(x,y)=x^{2}+y^{2}$, i.e. it solves the eigenvalue problem

[TABLE]

The errors and the computational times of the CN-FEM are presented in Table 2 and the errors and computational times of the SP2 in Table 3. The reference solution, $u_{\text{ref}}$, is computed using the CN-FEM with $h=0.0125$ and $\tau=2^{-13}$. The implantation was done in Julia. It is important to keep in mind that the SP2 uses mostly inbuilt functions such as the fast Fourier transform from the C subroutine library (FFTW). These functions are heavily optimized and show an extremely good performance. In spite of this we see, comparing the errors of the CN-FEM for $h=0.025$ and $\tau=2^{-10}$ to those of the SP2 for $N_{DoF}=3.2\cdot 10^{6}$ and $\tau=2^{-13}$, that they are on par with respect to CPU time relative to accuracy, with a slight computational advantage for the CN-FEM. This advantage becomes clearer, the larger the region of reduced regularity (e.g. in the context of optical lattices or disorder potentials). This justifies the usage of CN-FEM in low regularity regimes. Furthermore, it is clearly seen that the space discretization dominates the error of the SP2, in order for the spectral method to catch up in terms of accuracy with the CN-FEM with $h=0.0125$, an estimated 10 to 40 million degrees of freedom would be needed and thus the memory cost would become an issue.

Acknowledgements. We thank the anonymous referees for their helpful and insightful comments that improved the contents of this paper.

Appendix A Explicit decomposition of the consistency error

Here we make the consistency error (7) and its estimate (9) explicit and highlight where the regularity assumptions come in. For this we use the following standard integral remainder of Taylor expansion (for Sobolev functions $v\in H^{n+1}(a,b)$):

[TABLE]

where $a$ is the point of expansion and $T_{n}^{a,b}(v)$ the Taylor polynomial of degree $n$ of $v$ evaluated in $a$ and $b$. Recalling that $u^{n+1/2}=(u^{n+1}+u^{n})/2$, the Taylor expansion implies in our case:

[TABLE]

Thus

[TABLE]

Likewise we have

[TABLE]

and

[TABLE]

For the estimate of the term coming from the nonlinearity, we set for the sake of brevity $a:=|u^{n}|^{2}$, $b:=|u^{n+1}|^{2}$ and $c:=|u(t_{n+1/2})|^{2}$ to obtain

[TABLE]

Where $c^{\prime}$ lies between $a$ and $b$. Thus

[TABLE]

which proves (9). Finally we use the above expression for $\gamma(\xi^{n})-\gamma(|u(t_{n+1/2})|^{2})$, to bound $\|\nabla(\gamma(\xi^{k})-\gamma(|u(t_{k+1/2})|^{2}))\|^{2}$,

[TABLE]