A linearized energy--conservative finite element method for the nonlinear Schr\"{o}dinger equation with wave operator

TL;DR

This paper introduces a linearized finite element method for the nonlinear Schrödinger equation with wave operator that conserves energy, offers higher stability than non-conservative schemes, and provides optimal error estimates validated by numerical experiments.

Contribution

The paper presents a novel energy-conservative, linearized FEM for the nonlinear Schrödinger equation with wave operator, including a modified leap-frog scheme and rigorous error analysis.

Findings

Method conserves energy in discrete norm

Achieves higher stability compared to non-conservative schemes

Numerical results confirm optimal convergence and stability

Abstract

In this paper, we propose a linearized finite element method (FEM) for solving the cubic nonlinear Schr\"{o}dinger equation with wave operator. In this method, a modified leap-frog scheme is applied for time discretization and a Galerkin finite element method is applied for spatial discretization. We prove that the proposed method keeps the energy conservation in the given discrete norm. Comparing with non-conservative schemes, our algorithm keeps higher stability. Meanwhile, an optimal error estimate for the proposed scheme is given by an error splitting technique. That is, we split the error into two parts, one from temporal discretization and the other from spatial discretization. First, by introducing a time-discrete system, we prove the uniform boundedness for the solution of this time-discrete system in some strong norms and obtain error estimates in temporal direction. With the help of the preliminary temporal estimates, we then prove the pointwise uniform boundedness of the finite element solution, and obtain the optimal $L^2$-norm error estimates in the sense that the time step size is not related to spatial mesh size. Finally, numerical examples are provided to validate the convergence-order, unconditional stability and energy conservation.

Tables

Table 1. Table 1: Convergent analysis in spatial direction for example 1

	P1 element	P2 element
$1/h=16$	2.509E$-$01	9.096E$-$03
$1/h=32$	5.246E$-$02	1.098E$-$03
$1/h=64$	1.430E$-$02	1.486E$-$04
$orde{r}_{avg}$	2.06	2.97

Table 2. Table 2: Convergent analysis in temporal direction for example 1.

	P1 element	P2 element
$1/\tau =16$	2.456E$-$01	2.457E$-$01
$1/\tau =32$	6.312E$-$02	6.322E$-$02
$1/\tau =64$	1.598E$-$02	1.606E$-$03
$orde{r}_{avg}$	1.97	1.97

Table 3. Table 3: Convergent analysis in spatial direction for example 2.

	P1 element	P2 element
$1/h=16$	1.136E$-$02	6.159E$-$05
$1/h=32$	2.843E$-$03	7.511E$-$06
$1/h=64$	7.126E$-$04	9.321E$-$07
$orde{r}_{avg}$	2.00	3.02

Table 4. Table 4: Convergent analysis in temporal direction for example 2.

	P1 element	P2 element
$1/\tau =16$	2.367E$-$02	2.365E$-$02
$1/\tau =32$	6.196E$-$03	6.083E$-$03
$1/\tau =64$	1.572E$-$03	1.503E$-$03
$orde{r}_{avg}$	1.96	1.99

Equations

$$u_{tt}-\Delta u+\mathbf{i}{u}_t+|u{|}^{2}u+w(\mathbf{x})u=0,$$

$$u(\mathbf{x},t)=0,$$

$$u(\mathbf{x},0)=u_0(\mathbf{x}),u_t(\mathbf{x},0)=u_1(\mathbf{x}),$$

$$\frac{d}{dt}{\int }_{\Omega }\left(|u_t{|}^2+|\nabla u{|}^2+\frac{|u{|}^4}{2}+w(\mathbf{x})|u{|}^2\right)dx=0.$$

$$\Vert R_h{u}^n-U_h^n{\Vert }_{L^{\infty }}\le C{h}^{-\frac{d}{2}}\Vert R_h{u}^n-U_h^n{\Vert }_{L^2}\le C{h}^{-\frac{d}{2}}({\tau }^p+h^{r+1}),$$

$$(u,v)={\int }_{\Omega }u(x)(v(x){)}^{\ast }dx,$$

$$(\nabla (v-R_{h}v),\nabla \omega )=0,\forall \omega \in V_h.$$

$$\Vert v-R_{h}v{\Vert }_{L^2}+h\Vert \nabla (v-R_{h}v){\Vert }_{L^2}\le C{h}^{r+1}\Vert v{\Vert }_{H^{r+1}},$$

$$\Vert R_{h}v{\Vert }_{W^{1,p}}\le C\Vert v{\Vert }_{W^{1,p}},p>1,$$

$$\Vert v{\Vert }_{L^{\infty }}\le C{h}^{-\frac{d}{2}}\Vert v{\Vert }_{L^2},$$

$${\delta }_{\tau }^2{U}^n=\frac{U^{n+1}-2U^n+U^{n-1}}{{\tau }^2},D_{\tau }{U}^n=\frac{U^{n+1}-U^{n-1}}{2\tau },U^{\overline{n}}=\frac{U^{n-1}+U^{n+1}}{2}.$$

$$({\delta }_{\tau }^2{U}_h^n,v_h)+(\nabla U_h^{\overline{n}},\nabla v_h)+\mathbf{i}(D_{\tau }{U}_h^n,v_h)+\left(|U_h^n{|}^2{U}_h^{\overline{n}},v_h\right)+(w{U}_h^{\overline{n}},v_h)=0,n=1,2,\cdots ,N-1,$$

$$U_h^0=R_h{u}_0,U_h^1=R_h\left(u_0+\tau u_1+\frac{{\tau }^2}{2}(\Delta u_0-\mathbf{i}{u}_1-|u_0{|}^2{u}_0-w{u}_0)\right).$$

$$E^{N-1}=E^{N-2}=\cdots =E^1=E^0,$$

$$E^n\triangleq {\Vert \frac{U_h^{n+1}-U_h^n}{\tau }\Vert }_{L^2}^2+\frac{1}{2}(\Vert \nabla U_h^{n+1}{\Vert }_{L^2}^2+\Vert \nabla U_h^n{\Vert }_{L^2}^2)+\frac{\Vert U_h^n{U}_h^{n+1}{\Vert }_{L^2}^2}{2}+\frac{(w{U}_h^n,U_h^n)+(w{U}_h^{n+1},U_h^{n+1})}{2}.$$

$$\text{Re}({\delta }_{\tau }^2{U}_h^n,U_h^{n+1}-U_h^{n-1})+\text{Re}(\nabla U_h^{\overline{n}},\nabla (U_h^{n+1}-U_h^{n-1}))+\text{Re}(\mathbf{i}{D}_{\tau }{U}_h^n,U_h^{n+1}-U_h^{n-1})$$

$$+\text{Re}(|U_h^n{|}^2{U}_h^{\overline{n}},U_h^{n+1}-U_h^{n-1})+\text{Re}(w{U}_h^{\overline{n}},U_h^{n+1}-U_h^{n-1})=0.$$

$$I_1={\Vert \frac{U_h^{n+1}-U_h^n}{\tau }\Vert }_{L^2}^2-{\Vert \frac{U_h^n-U_h^{n-1}}{\tau }\Vert }_{L^2}^2,$$

$$I_2=\frac{1}{2}\left(\Vert \nabla U_h^{n+1}{\Vert }_{L^2}^2+\Vert \nabla U_h^n{\Vert }_{L^2}^2\right)-\frac{1}{2}\left(\Vert \nabla U_h^n{\Vert }_{L^2}^2+\Vert \nabla U_h^{n-1}{\Vert }_{L^2}^2\right),$$

$$I_3=0,$$

$$I_4=\frac{\Vert U_h^{n+1}{U}_h^n{\Vert }_{L^2}^2}{2}-\frac{\Vert U_h^n{U}_h^{n-1}{\Vert }_{L^2}^2}{2},$$

$$I_5=\frac{(w{U}_h^{n+1},U_h^{n+1})+(w{U}_h^n,U_h^n)}{2}-\frac{(w{U}_h^n,U_h^n)+(w{U}_h^{n-1},U_h^{n-1})}{2}.$$

$$E^n=E^{n-1},n=1,\cdots ,N-1.$$

$$\Vert u_0{\Vert }_{H^{r+3}}+\Vert (u_0{)}_{ttt}{\Vert }_{H^2}+\Vert u{\Vert }_{L^{\infty }([0,T];H^{r+1})}+\Vert u{\Vert }_{L^{\infty }([0,T];W^{1,\infty })}+\Vert u_t{\Vert }_{L^{\infty }([0,T];H^{r+1})}$$

$$+\Vert u_{tt}{\Vert }_{L^{\infty }([0,T];H^{r+1})}+\Vert u_{tttt}{\Vert }_{L^{\infty }([0,T];H^1)}\le M_0,$$

$${\delta }_{\tau }^2{U}^n-\Delta U^{\overline{n}}+\mathbf{i}{D}_{\tau }{U}^n+|U^n{|}^2{U}^{\overline{n}}+w{U}^{\overline{n}}=0,n=1,2,\cdots ,N-1,$$

$$U^n(\mathbf{x})=0,\text{x∈∂Ω, n=2,3,⋯,N,}$$

$$U^0(\mathbf{x})=u_0(\mathbf{x}),U^1(\mathbf{x})=u_0(\mathbf{x})+\tau u_1(\mathbf{x})+\frac{{\tau }^2}{2}\left(\Delta u_0-\mathbf{i}{u}_1-|u_0{|}^2{u}_0-w{u}_0\right),\text{x∈Ω}.$$

$$\Vert u^n-U_h^n{\Vert }_{L^2}\le \Vert u^n-U^n{\Vert }_{L^2}+\Vert U^n-U_h^n{\Vert }_{L^2}.$$

$$a_j+\tau \sum _{k=0}^j{b}_k\le \tau \sum _{k=0}^j{\gamma }_k{a}_k+\tau \sum _{k=0}^j{c}_k+B,\mathrm{for}j\ge 0,$$

$$a_j+\tau \sum _{k=0}^j{b}_k\le \exp (\tau \sum _{k=0}^j{\gamma }_k{\sigma }_k)(\tau \sum _{k=0}^j{c}_k+B),\mathrm{for}j\ge 0.$$

$${\delta }_{\tau }^2{u}^n-\Delta u^{\overline{n}}+\mathbf{i}{D}_{\tau }{u}^n+|u^n{|}^2{u}^{\overline{n}}+w{u}^{\overline{n}}=P^n,n=1,2,\cdots ,N-1,$$

$$P^n=({\delta }_{\tau }^2{u}^n-u_{tt}^n)-\Delta (u^{\overline{n}}-u^n)+\mathbf{i}(D_{\tau }{u}^n-u_t^n)+\left(|u^n{|}^2+w\right)\left(u^{\overline{n}}-u^n\right).$$

$${\left(\tau \sum _{k=1}^n\Vert P^k{\Vert }_{L^2}^2\right)}^{\frac{1}{2}}+{\left(\tau \sum _{k=1}^n\Vert \nabla P^k{\Vert }_{L^2}^2\right)}^{\frac{1}{2}}\le C{\tau }^2,n=1,2,\cdots ,N-1.$$

$$\Vert U^m{\Vert }_{H^2}\le M+1,$$

Full text

AMS subject classifications. 65N12, 65N30, 35K61.

1 Introduction

In this paper, we study the following cubic nonlinear Schrödinger equation with wave operator

[TABLE]

for $\mathbf{x}\in\Omega$ and $t\in[0,T]$, where $\Omega$ is a bounded convex domain in $\mathbb{R}^{d}$ ($d=2,\ 3$), $\mathbf{i}=\sqrt{-1}$ is the complex unit, $w(\mathbf{x})$ is the real–valued potential function. Meanwhile, the initial and boundary conditions are defined by

[TABLE]

where $u_{0},u_{1}$ are two given functions. An important property of equation (1) is the energy conservation. Computing the inner product of Eq. (1) with $u_{t}$ in $L^{2}(\Omega)$, and taking the real parts, one can obtain the following energy conservative identity

[TABLE]

The nonlinear Schödinger equation is one of the most important equations in mathematical physics, which is originated from quantum mechanics. It has been widely used to model various nonlinear physical phenomena, such as underwater acoustics [1], nonlinear optics [2, 3], quantum condensates [4] and other nonlinear phenomena [5]. The cubic nonlinear Schödinger equation is one of the most important models, and it is also known as Gross–Pitaevskii equation (GPE), which plays a fundamental role in modeling the hydrodynamics of Bose–Einstein condensate [6, 7, 8].

Recently, nonlinear Schrödinger–type equations have also been widely studied. The nonlinear Schrödinger equation with wave operator is one of most important nonlinear Schrödinger–type equations, it has been derived from many physical areas. For example, the nonrelativistic limit of the Kelin–Gordon equation [9, 10, 11], the Langmuir wave envelope approximation in plasma [12] and the modulated planar pulse approximation of the sine–Gordon equation for light bullets [13, 14].

Due to the wide applications of the Schrödinger and Schrödinger–type equations, performing efficient and accurate numerical simulations plays an essential role in many real applications. It is remarkable that nonconservative schemes for Schrödinger–type equations may lead to numerical blow–up [15]. Therefore, conservative schemes become very important for Schrödinger and Schrödinger–type equations. In the last several decades, numerical simulations of both the nonlinear Schrödinger equation and the nonlinear Schrödinger–type equation have been studied extensively. For examples, finite difference methods [15, 17, 16, 22, 23, 18, 19, 20, 21], finite element methods [24, 27, 28, 25, 26, 29, 30, 31, 32] and Fourier spectral method [33]. In the field of finite difference methods, an implicit nonconservative difference scheme had been developed in [16] for solving nonlinear Schrödinger equation, the method needs lots of algebraic operators. Zhang et al. [15] had pointed out that the nonconservative schemes may easily lead to the numerical solution blow up. Thus, a conservative difference scheme was provided for solving the nonlinear Schrödinger in their work. Subsequently, conservative finite difference methods were developed to solve the Schrödinger equation with wave operator [22, 23, 18, 19, 20, 21]. In the field of finite element methods, Galerkin finite element methods were used to solve the generalized nonlinear Schrödinger equation [24, 25, 26], the coupled nonlinear Schrödinger equation [27] and the nonlinear Schrödinger–Helmholtz system [28]. In these works, optimal error estimates are achieved in the sense that the time step size is not related to spatial mesh size. Moreover, the local discontinuous Galerkin methods were used to simulate the 1D nonlinear Schrödinger equation [30, 31] and multi–dimensional nonlinear Schrödinger equation with wave operator [32]. However, in [32], optimal error estimates were only obtained for the linear equation at semi–discrete level. Therefore, in the context of numerical analysis, rigorous error estimates for the numerical scheme of nonlinear Schrödinger equation with wave operator are needed.

Previously, when using Galerkin FEMs for solving partial differential equations (PDEs)

[34, 35, 36], to obtain the error estimates of linearized explicit (or semi–implicit), the pointwise uniform boundedness of numerical solution in certain strong–norms was often required. Traditionally, the inverse inequality and mathematical induction were applied to obtain the pointwise boundness of the numerical solution,

[TABLE]

where $u^{n}$ and $U^{n}_{h}$ are the exact and numerical solutions at time level $t_{n}$, respectively. $R_{h}$ is Ritz projection operator and $d$ is the dimension, $\tau,h$ are temporal and spatial mesh sizes, $p,r$ are positive integers, referring to the convergence order of the FEM. But, to get the uniform boundedness of numerical solutions, the above inequality results in an unnecessary restriction between the time step size and spatial mesh size. Recently, a new technique [37, 38] was introduced to analyze error estimates of the linearized semi–implicit FEMs for time–dependent nonlinear PDEs. In Li and Sun’s works [37, 38], errors were split into two parts, one part was from the temporal discretization and the other part was from the spatial discretization. By analyzing the introduced time–discrete PDEs, the pointwise uniform boundedness of FEM solutions can be proved in the sense that there is no restriction between the time step size and spatial mesh size. Comparing to previous error analysis with conditional stability in [34, 35, 36], optimal error estimates were obtained unconditionally in [37, 38]. Recently, this new technique has been used to analyze linearized FEMs for nonlinear Schrödinger type equation [24, 27, 28, 25, 26] and many other PDEs [39, 40, 41, 42, 43, 44]. To our best knowledge, this new technique is mainly carried out for nonlinear parabolic type of equations and also possibly coupled with elliptic type of equations. However, when considering the nonlinear Schrödinger equation with wave operator, it will possess the combination of both parabolic and hyperbolic properties, and this could raise some complexity for the error analysis.

In this paper, a linearized FEM is proposed to solve the cubic nonlinear Schrödinger equation with wave operator (1) subject to initial boundary conditions (2)–(3). In this scheme, system (1)–(3) is discretized by a modified frog–leap scheme in time direction and the Galerkin finite element method in spatial direction. The proposed scheme is a semi-implicit linear method which only needs to solve a linear system at each time step. Thus, the proposed scheme is simpler and more efficient than implicit nonlinear schemes, which need to do iteration at each time step. More importantly, discrete energy of the proposed method is conserved so that the scheme will not yield blow up. Subsequently, we will apply the error splitting technique [37, 38] to study the proposed linearized energy–conservative FEM. By introducing a time–discrete system, we will prove the uniform boundedness of time–discrete solutions in certain strong norms, and give the error estimates of time–discrete solutions. Based on the uniform boundedness of time–discrete solutions and mathematical induction, we get the pointwise uniform boundedness of fully–discrete FEM solutions in the sense that there is no restriction between the time step size and spatial mesh size. With the help of the above pointwise uniform boundedness of the FEM solutions and the traditional error analysis method, we can obtain the optimal error estimates in $L^{2}$–norm. Our work in this paper can be regarded as a complementary to Guo and Xu’s work [32].

The rest of the paper is organized as follows. In section 2, a linearized Galerkin FEM for the nonlinear Schrödinger equation with wave operator is given, together with some assumptions and notations. Meanwhile, the energy conservative property of the fully–discrete system is presented and the time–discrete system is introduced. In section 3, the uniform boundedness of the time–discrete solutions is proved in some strong norms. Moreover, the error estimates of the time–discrete solutions are obtained. In section 4, based on the uniform boundedness of time–discrete solution in $H^{2}$–norm, we obtain the uniform boundedness of the fully discrete solutions in $L^{\infty}$–norm. In section 5, we derive the optimal error estimates in $L^{2}$–norm unconditionally. Section 6 shows the numerical results, which confirm well with the theoretical findings. Finally, a brief conclusion is provided in the final section.

2 Linearized Galerkin FEM method

Let $\Omega$ be the convex polygon in $\mathbb{R}^{2}$ (or convex polyhedron in $\mathbb{R}^{3}$). Following the classical FEM theory, we define $\mathcal{T}_{h}$ be a quasi-uniform partition of $\Omega$ into triangles $\pi_{h}$ in $\mathbb{R}^{2}$ or tetrahedra in $\mathbb{R}^{3}$. Let $h=\max\limits_{\pi_{h}\in\mathcal{T}_{h}}{\rm diam}\,\pi_{h}$ denotes the meshsize. Let $V_{h}$ be the finite–dimensional subspace of $H^{1}_{0}(\Omega)$, which consists of continuous piecewise polynomials of degree $r$ ($r\geq 1$) on $\mathcal{T}_{h}$.

For any two complex functions $u,v$, the inner product is defined as

[TABLE]

where $v^{*}$ denotes the conjugate of $v$.

Let $R_{h}:H^{1}_{0}(\Omega)\rightarrow V_{h}$ be the Ritz projection operator defined by

[TABLE]

By the classical FEM theory [45], the following inequality is valid,

[TABLE]

for any $v\in H^{r+1}_{0}(\Omega)$.

In this paper, the following inverse inequality [45] is always used:

[TABLE]

for any $v\in V_{h}$ and $d=2,3.$

Let $\{t_{n}|t_{n}=n\tau;0\leq n\leq N\}$ be a uniform partition of $[0,T]$ with time step size $\tau=\frac{T}{N}$, and $U^{n}=u(\cdot,t_{n})$. We define

[TABLE]

With above notations, a linearized leap–frog Galerkin FEM is to seek $U^{n+1}_{h}\in V_{h}$ such that

[TABLE]

for any $v_{h}\in V_{h}.$ The initial and first step FEM solutions are defined by

[TABLE]

Next, we give a theorem to show that Eqs. (10)–(11) is an energy conservative scheme.

Theorem 1

The discrete energy of the FEM scheme Eqs. (10)–(11) is conservative, i.e.,

[TABLE]

where the discrete energy is defined as

[TABLE]

Proof. Putting $v_{h}=U^{n+1}_{h}-U^{n-1}_{h}$ in (10) and taking the real parts for the each term of the result equation, one has

[TABLE]

If we denote the terms of left side of above equation by $I_{i}$ ($i=1,\cdots,5$), then it is easy to see that

[TABLE]

Thus, from above results and (2), we have

[TABLE]

The proof is complete.  $\Box$

In this paper, we assume that the solution to the initial boundary problem (1)–(3) exists and satisfies,

[TABLE]

where $r\ (r\geq 1)$ is the order of the Galerkin FE space used in (10), $M_{0}$ is a positive constant depends only on $\Omega$. In addition, the potential function $w(\mathbf{x})$ is assumed to belong to $H^{r+1}(\Omega)$. Denote $M=\|u\|_{L^{\infty}([0,T];H^{2})}$, then by the above assumption, $M$ is a finite positive number which only depends on $\Omega$.

Now we introduce a corresponding time–discrete system for the Galerkin FEM scheme (10),

[TABLE]

with the following initial and boundary conditions

[TABLE]

With the solutions of time discrete system (14)–(16), we split the errors into two parts

[TABLE]

Under this splitting, we will prove that the first term in the right-hand side of above inequality is bounded by $O(\tau^{2})$ and the second term is bounded by $O(h^{2})$, with which, the classical inverse inequality and inductive assumption, we can obtain that the FE solutions are bounded in $L^{\infty}$–norm. For the simplicity of notations, we use $C$ to denote a generic positive constant and use $\epsilon$ to denote a generic small positive constant, where $C,\epsilon$ are independent of spatial and temporal meshsize.

Next, we give the following inequality, which will be used frequently.

Lemma 1

Discrete Gronwall’s inequality [46, 47] : Let $\tau$, $B$ and $a_{k}$, $b_{k}$, $c_{k}$, $\gamma_{k}$, for integers $k\geq 0$, be non–negative numbers such that

[TABLE]

suppose that $\tau\gamma_{k}<1$, for all $k$, and set $\sigma_{k}=(1-\tau\gamma_{k})^{-1}$. Then

[TABLE]

3 Temporal error estimates

In this subsection, we analyze the uniform boundedness of time–discrete solution $U^{n}$ in strong norms. Moreover, the error estimates of solution $U^{n}$ in certain norms are presented.

Let $u$ be the solution of the system (1)–(3), then $u^{n}$ satisfies

[TABLE]

where

[TABLE]

By using Taylor’s expansion and the assumption (2), one can obtain that

[TABLE]

Theorem 2

Suppose that the system (1)–(3) has a unique solution u satisfying (2). Then the time discrete system defined in (14)–(16) has a unique solution $U^{m},m=1,\cdots,N$. Moreover, there exists $\tau_{0}>0$, such that when $\tau<\tau_{0}$,

[TABLE]

where $e^{m}=u^{m}-U^{m}$, $C^{\prime}$ is a positive constant independent of $\tau$ and $h$.

Proof. Since the system (14) is a linear elliptic equation, we can get the existence and uniqueness of the solution $U^{m}$, $m=2,3,\cdots,N.$ First, we prove that there exists $\tau^{\prime}_{0}>0$ such that the estimate (20) holds for $m=1,\cdots,N.$

Obviously, when $m=0$, $e^{0}=0.$ From (2) and (16), it is easy to see that, when $\tau\leq\tau_{1}=\frac{C^{\prime}}{C_{1}}$,

[TABLE]

Hence, (20) holds for $m=1$.

Now suppose that (20) is valid for $m\leq k-1$. Then there exists a positive constant $\tau_{2}$ such that for $\tau\leq\tau_{2}=\frac{1}{\sqrt{C_{0}C^{\prime}}}$,

[TABLE]

Next, we prove (20) holds for $m=k$. Subtracting Eq. (14) from Eq. (17) gives

[TABLE]

where $R^{n}=(|u^{n}|^{2}+w)u^{\overline{n}}-(|U^{n}|^{2}+w)U^{\overline{n}}$.

Since

[TABLE]

one has

[TABLE]

where (2), (22) and the regularity of $w$ are used.

In addition,

[TABLE]

By using (2) and (22), one has

[TABLE]

where the analysis here is the same as (29).

[TABLE]

Adding the equations from (26) to (33), one has

[TABLE]

for $1\leq n\leq k-1$.

Computing the inner product with $D_{\tau}e^{n}$ on both sides of (23) and taking the real parts, one gets

[TABLE]

Multiplying above inequality both sides by $2\tau$ and summing up from $n=1$ to $n=k-1$, one obtains

[TABLE]

Now computing the inner product with $-\Delta D_{\tau}e^{n}$ on both sides of (23) and taking the real parts, one has

[TABLE]

Multiplying above inequality both sides by $2\tau$ and summing up from $n=1$ to $n=k-1$, one gets,

[TABLE]

Adding Eqs. (3) and (3), we get

[TABLE]

where we have used the Poincar$\acute{\rm{e}}$ inequality $\|e^{n}\|_{L^{2}}\leq C\|\nabla e^{n}\|_{L^{2}}$.

By using Gronwall’s inequality to (3), one obtains that, there exists $\tau_{3}>0$, such that for $\tau\leq\tau_{3}$,

[TABLE]

which implies that

[TABLE]

Now we take $C^{\prime}=C_{2}$. After $C^{\prime}$ is obtained, we let $\tau^{\prime}_{0}=\min\{\tau_{1},\tau_{2},\tau_{3}\}$, then (20) holds for $m=k$. The mathematical induction of (20) is finished. Therefore, (20) holds for $m=1,2,...,N.$ Thus, there exists a positive constant $\tau_{4}$ such that for $\tau\leq\tau_{4}=\frac{1}{\sqrt{C^{\prime}C_{0}}}$,

[TABLE]

for $m=1,2,...,N$. Let $\tau_{0}=\min\{\tau^{\prime}_{0},\tau_{4}\}$, we finish the proof of (19) and (20).\qed

4 Spatial error estimates

In this section, we will study the uniform boundedness of solution $U^{n}_{h}$ in strong norms.

Lemma 2

Suppose that the time–discrete system (14)–(16) has a unique solution $U^{n}$ ($n=0,\cdots,N$), then

[TABLE]

Proof. By using (19) and the Sobolev embedding theorem, we have

[TABLE]

Thus, the proof of this lemma is complete. \qed

The variation form of (14) can be defined by

[TABLE]

for $n=1,2,\cdots,N-1$.

Subtracting (10) from (44) and using (6), one gets

[TABLE]

where $Q^{n}=\left(|U^{n}|^{2}+w\right)U^{\overline{n}}-\left(|U^{n}_{h}|^{2}+w\right)U^{\overline{n}}_{h}$ and $e^{n}_{h}=R_{h}U^{n}-U^{n}_{h}$.

Theorem 3

Suppose that the system (1)–(3) has a unique solution $u$ satisfying (2). Then the finite element equation defined by (10) has a unique solution $U^{n}_{h},n=2,\cdots,N$. Moreover, there exists $\tau^{*}>0,h^{*}>0$, such that,

[TABLE]

when $\tau<\tau^{*},h<h^{*}$.

Proof. Firstly, we prove the existence and uniqueness of solution $U^{n}_{h}$ of system (10). If the solution $U^{n}_{h}$ is given for $n=0,1,2,\cdots,m-1$, then the system (10) has a unique solution if and only if the following homogeneous equation

[TABLE]

has only zero solution.

Let $v_{h}=W_{h}$ in above equation and take the imaginary part, one has

[TABLE]

which implies that $W_{h}=0$. Thus, equation (47) has only zero solution. It is natural to obtain that the uniqueness and existence for the solutions of (10).

Before the proof of (46) is given, we study the following result: there exists two constants $\tau^{*}>0$ and $h^{*}>0$, such that, when $\tau\leq\tau^{*}$ and $h\leq h^{*}$,

[TABLE]

As $k=0,1$, from (16) and (11), we can get $e^{0}_{h}=0$, $e^{1}_{h}=0$. Thus, (48) holds for $k=0,1$, and

[TABLE]

when $h\leq h_{1}=\left(\frac{1}{C_{3}}\right)^{\frac{6}{10-3d}}$ ($d=2,3$) ($h_{1}>0$), where the inverse inequality $\|e_{h}\|_{L^{\infty}}\leq C_{3}h^{-\frac{d}{2}}|e_{h}\|_{L^{2}}$ is used.

Now suppose that (48) is valid for $k\leq m-1$, we can get

[TABLE]

when $h\leq h_{1}$. Next, we want to show that (48) is also valid for $k=m$.

Putting $v_{h}=D_{\tau}e^{n}_{h}$($1\leq n\leq m-1$) on both sides of (45) and taking the real parts, one has

[TABLE]

And $\|Q^{n}\|_{L^{2}}$ can be further analyzed, thus

[TABLE]

Summing above inequality from $n=1$ to $n=m-1$, we obtain

[TABLE]

where $e^{0}_{h}=0$, $e^{1}_{h}=0$ and the Poincar$\acute{\rm{e}}$ inequality $\|e^{n}_{h}\|_{L^{2}}\leq C\|\nabla e^{n}_{h}\|_{L^{2}}$ are used.

Using the Gronwall’s inequality, there is a $\tau_{5}>0$, when $\tau\leq\tau_{5}$

[TABLE]

By the fact that $e^{n}_{h}\in H^{1}_{0}$, we have

[TABLE]

There exists a positive constant $h_{2}=\left(\frac{1}{C_{5}}\right)^{3}>0$, such that, when $h\leq h_{2}$,

[TABLE]

Thus, (48) holds for $k=m$. Let $\tau^{*}=\tau_{5}$, $h^{*}=\min\{h_{1},h_{2}\}$, the induction proof for (48) is closed.

Next, we prove the (46) holds for $n=2,\dots,N$. From (48), we can get

[TABLE]

when $h\leq h^{*}$.

The proof is complete.\qed

5 Convergent analysis for the full discrete scheme

In this section, we will provide the convergent analysis for the full discrete scheme (10)–(11).

Theorem 4

Suppose that the system (1)–(3) has a unique solution u satisfying (2). Then the finite element system defined in (10)–(11)* has a unique solution $U^{n}_{h},n=0,\cdots,N$ and there exists $\tau_{*}>0$, $h_{*}>0$ such that*

[TABLE]

when $\tau<\tau_{*},h<h_{*}$, where $C^{0}$ is a positive constant independent of $\tau$ and $h$.

Proof. As $n=0,1$, from (2) and (7), we have

[TABLE]

Thus, (56) holds for $n=0,1$.

Next, we prove (56) holds for $2\leq n\leq N$. The weak form of system (17) can be written as:

[TABLE]

where $1\leq n\leq N-1$.

Subtracting (58) from (10) and applying (6), we can get

[TABLE]

where $G^{n}=|u^{n}|^{2}u^{\overline{n}}-|U^{n}_{h}|^{2}U^{\overline{n}}_{h}$, $\sigma^{n}_{h}=R_{h}u^{n}-U^{n}_{h}$.

From (7), (2) and (46), we can get

[TABLE]

Putting $v_{h}=D_{\tau}\sigma^{n}_{h}$ in (5) and taking the real parts, we get

[TABLE]

Summing above inequality from $n=1$ to $n=N-1$, we have

[TABLE]

where the second inequality can be obtained by the fact that

[TABLE]

Applying Gronwall’s inequality and $\sigma^{0}_{h}=0$ to (5), there exists $\tau_{*}>0$, when $\tau\leq\tau_{*}$

[TABLE]

Due to $\sigma^{n}_{h}\in H^{1}_{0}$, we have

[TABLE]

Combining above result with (7), we have

[TABLE]

Thus, let $h_{*}=h^{*}$ and $C^{0}>C_{8}$ in above result, the proof of Theorem 4 is complete.\qed

6 Numerical results

In this section, two numerical examples are presented to verify our theoretical analysis, where a free software (Freefem++) is used to perform all computations.

Example 1

Consider the following cubic nonlinear Schrödinger equation with wave operator:

[TABLE]

where $w(x,y)=-x^{2}y^{2}$, and $\Omega=\{(x,y):(x-0.5)^{2}+(y-0.5)^{2}<0.5^{2}\}$. Moreover, the initial boundary conditions and the term $g$ in right-hand side are obtained by the exact solution

[TABLE]

A quasi–uniform triangulation is generated by FreeFEM++ with $M$ nodes distributed on the boundary of the circular domain $\Omega$. In the numerical implementation, we use both the linear and quadratic finite element when doing the spatial discretization, the final time is set to be $T=1$.

We first list the $L^{2}$–norm errors of numerical solutions in Table 1 with an extremely small time step size $\tau=2^{-14}$ so that the errors from temporal direction can be neglected. From Table 1, we can see that the numerical method reaches its optimal convergent order in spatial direction. Moreover, the convergence study for temporal direction using $L^{2}$–norm errors are also presented in Table 2 with the small spatial step size $h=1/512$. From Table 2, we can see that the numerical method is second–order accurate in temporal direction. To check unconditional stability, we take different spatial sizes $1/h=8,16,32,64,128$ with fixed time steps $\tau=0.1,0.05,0.01$. From Fig. 1, we can see that, for each fixed time step size, both $L^{2}$–norm errors from the linear and quadratic FEMs tend to be a constant when the mesh is refined gradually. It is shown that our numerical scheme is unconditionally stable. Thus, all numerical results completely match with our theoretical analysis.

Example 2

Consider the cubic nonlinear Schrödinger equation with wave operator:

[TABLE]

where $w=-(\sqrt{2}\pi+\sin^{2}(\pi x)\sin^{2}(\pi y))$ and $\Omega=[0,1]\times[0,1]$. Moreover, the initial boundary conditions are obtained by the exact solution

[TABLE]

Again, we solve the equation (68) by both linear and quadratic FEMs up to time $T=1$. Similar to the above example, the convergent study in spatial direction using $L^{2}$–norm errors are shown in table 3 by choosing extremely small time–step size $\tau=2^{-14}$. Meanwhile, the convergent study in temporal direction are presented with the small spatial step size $h=1/512$. Again, from Table 3 and Table 4, we can see that the proposed method reaches its optimal convergent order in both spatial and temporal directions.

Further, similar to the above example, we gradually refine the spatial mesh size for three fixed time–step sizes $\tau=0.1,0.05,0.01$ at $T=1$, From Fig. 2, we can see that the $L^{2}$–norm errors of both linear and quadratic FEMs tend to be a constant, which implies that our method is unconditionally stable.

Further, to illustrate the energy conservative property of our numerical scheme, we present the numerical results at different time stages $T=0,10,20,\cdots,100,$ by using linear element approximation. From Fig. 3, we can see that the discrete energy $E$ is conserved exactly (up to machine accuracy) with time evolution, which verifies our theoretical result in Theorem 1.

7 Conclusions

In this paper, an energy–conservative finite element method is present to solve nonlinear Schrödinger equation with wave operator. Comparing to previous works [25, 26], our scheme is proved to keep energy conservation in a certain discrete norm (Theorem 1). Thus, our scheme keeps higher stability. Moreover, we give unconditionally optimal error estimates of the modified leap–frog FEM for cubic nonlinear Schrödinger equation with wave operator. By introducing a time–discrete system, the error of numerical solutions is split into two parts: the temporal error and the spatial error. We present the uniform boundedness of time–discrete solutions in some strong norms and the error estimates in temporal direction. Based on these results, we get the $L^{2}$ optimal error estimates in the sense that the time step size is not related to spatial mesh size. At last, numerical examples are provided for verifying the convergence–order, unconditional stability and energy conservation of the proposed numerical method.

Acknowledgements

Dongdong He was supported by the president’s fund–research start–up fund from the Chinese University of Hong Kong, Shenzhen (No. PF01000857). Kejia Pan was supported by the National Natural Science Foundation of China (Nos. 41874086, 41474103), the Excellent Youth Foundation of Hunan Province of China (No. 2018JJ1042) and the Innovation–Driven Project of Central South University (No. 2018CX042).