High-order mass- and energy-conserving SAV-Gauss collocation finite element methods for the nonlinear Schrödinger equation
HHigh-order mass- and energy-conserving SAV–Gauss collocationfinite element methods for the nonlinear Schr¨odinger equation
Xiaobing Feng ∗ Buyang Li † Shu Ma † Abstract
A family of arbitrarily high-order fully discrete space-time finite element methods areproposed for the nonlinear Schr¨odinger equation based on the scalar auxiliary variable for-mulation, which consists of a Gauss collocation temporal discretization and the finite elementspatial discretization. The proposed methods are proved to be well-posed and conservingboth mass and energy at the discrete level. An error bound of the form O ( h p + τ k +1 ) in the L ∞ (0 , T ; H )-norm is established, where h and τ denote the spatial and temporal mesh sizes,respectively, and ( p, k ) is the degree of the space-time finite elements. Numerical experimentsare provided to validate the theoretical results on the convergence rates and conservationproperties. The effectiveness of the proposed methods in preserving the shape of a solitonwave is also demonstrated by numerical results. Key words:
Nonlinear Schr¨odinger equation, mass- and energy-conservation, high-orderconserving schemes, SAV-Gauss collocation finite element method, error estimates.
This paper is concerned with the development and analysis of high-order fully discrete numericalmethods for the following initial-boundary value problem of the nonlinear Schr¨odinger (NLS)equation: i ∂ t u − ∆ u − f ( | u | ) u = 0 in Ω × (0 , T ] , (1.1a) u = 0 on ∂Ω × (0 , T ] , (1.1b) u = u in Ω × { } , (1.1c)where Ω ⊂ R d is a polygonal or polyhedral domain with boundary ∂Ω , and u : Ω → C is acomplex-valued function, with i = √−
1, and f : R + → R is the derivative of some function F : R + → R . The best known examples are f ( s ) = ± s q − and F ( s ) = ± q + 1 s q +12 , with q > , (1.2)where the “ − ” and “+” cases are often referred to as defocusing and focusing models, respec-tively. In the focusing case, the solution will blow up in L ∞ ( Ω ) within finite time when the ∗ Department of Mathematics, The University of Tennessee, Knoxville, TN 37996, U.S.A. Email address:[email protected]). The work of this author was partially supported by the NSF-grant DMS-1620168. † Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, HongKong. Email address: [email protected], [email protected]. The work of B. Li was partiallysupported by a Hong Kong RGC grant (project ID: P0031035 ZZKQ) and the work of S. Ma was partiallysupported by a Hong Kong RGC grant (Project No. 15300817). a r X i v : . [ m a t h . NA ] J un nitial energy is negative; see [8, 36]. The NLS equation (1.1) arises from many applicationsin physics and engineering, and is one of the fundamental equations in mathematical physics[8, 36, 44, 27, 30].It is well known that the solutions of (1.1) conserve the mass and energy in the sense thatfor all t ≥ t (cid:90) Ω | u | d x = 0 , (mass conservation) (1.3)dd t (cid:90) Ω (cid:16) |∇ u | − F ( | u | ) (cid:17) d x = 0 . (energy conservation) (1.4)The development of numerical methods that can retain these conservation properties in numer-ical solutions is important for long-time numerical simulation, and therefore has been one of theresearch focuses in numerical approximation to the NLS equation.There exists a large amount of literature on numerical solutions and numerical analysis ofthe NLS equation. Delfour et al. [11] proposed a second-order modified Crank–Nicolson schemefor the NLS equation to preserve both mass and energy conservations. The construction of thismethod was motivated by a similar scheme for the Klein–Gordon equation in [34]. Sanz-Serna[28] extended the modified Crank–Nicolson time-stepping scheme to the NLS equation with moregeneral nonlinear term, and established an optimal order error estimate for the fully discretefinite element discretization. For the cubic NLS equation, more delicate results on the existence,uniqueness and convergence of numerical solutions of the modified Crank–Nicolson scheme wereproved by Akrivis et al. in [1]. The modified Crank–Nicolson scheme has been widely usedin computation and was combined with spatial finite difference methods to solve the NLS andGross–Pitaevskii equations in [2, 4, 5, 6], and with various spatial discretization methods toapproximate the NLS in [37, 21, 25, 39, 14]. Besides the modified Crank–Nicolson scheme, asecond order explicit leapfrog scheme for the NLS equation was proposed by Sanz-Serna andManoranjan in [29]. This scheme was proved to preserve mass conservation at the discrete level.Optimal rates of convergence of the fully discrete leapfrog finite element method was establishedin [28]. To avoid solving nonlinear systems, a linearly implicit leapfrog scheme was proposed byFei et al. in [45] for the NLS equation, and a linearly implicit relaxation scheme was proposedby Besse in [7]. Both schemes preserve mass and energy conservations at the discrete level.More recently, Feng et al. [13] constructed a class of second-order mass- and energy-conservingschemes for the NLS equation, including the modified Crank–Nicolson method, the implicitleapfrog method and a class of modified backward differentiation formulae as special cases.To the best of our knowledge, all the existing mass- and energy-conserving methods haveonly second-order accuracy in time. No higher-order time-stepping schemes, which conserve bothmass and energy, have been reported in the literature. Moreover, the existing error estimates fornonlinearly implicit schemes for the NLS equation generally require certain grid-ratio conditions.The standard grid-ratio conditions in the literature are τ = o ( h d ) for the cubic NLS equationand τ = o ( h d ) for general nonlinearity, where h and τ denote the spatial and temporal meshsizes. Karakashian and Makridakis [22, 23] proposed some continuous and discontinuous space-time Galerkin finite element methods for the cubic NLS equation and proved optimal-orderconvergence under a weaker grid-ratio condition τ k − | ln h | → k ≥ L ∞ -norm error estimate between thefully discrete and the semidiscrete-in-time numerical solutions. The error splitting techniqueallows to avoid grid-ratio conditions in using the inverse inequality.The objective of this paper is to develop a family of arbitrarily higher-order mass- andenergy-conserving fully discrete space-time finite element methods based on the scalar auxiliaryvariable (SAV) formulation of the NLS equation, and to establish the existence, uniquenessand optimal order convergence of numerical solutions without grid-ratio condition. Two keyideas are utilized in our construction of the method. First, the SAV reformulation of the NLSequation is used. This approach was introduced in [32, 31] as an enhanced version of the invariantenergy quadratization (IEQ) approach [40, 41, 42, 43], for developing energy-decay methods fordissipative (gradient flow) systems. Here we adapt the SAV approach to the dispersive NLSequation, and the SAV reformulation is essential to enable our methods to maintain the energyconservation property at the discrete level. Second, the Gauss collocation method is used fortime discretization in the SAV formulation of the NLS equation. The method can be viewed asan efficient implementation of the space-time finite element methods for the SAV formulationwith Gauss quadrature in time. The Gauss collocation method was combined with IEQ andSAV to preserve energy decay in solving phase field equations in [3, 18, 19]. We adopt thismethod here to preserve mass conservation without affecting the energy conservation structureof the SAV formulation.The SAV formulation introduces new difficulties to error analysis for the NLS equation due tothe presence of ∂ t u in the equation of r , see equation (2.6b), which leads to a consistency error ofsub-optimal order in time and introduces new difficulty in obtaining the stability estimate. As faras we know, rigorous analysis for convergence of numerical methods based on SAV formulationshas not been done for any wave equation so far. These difficulties are overcome by combiningthree techniques. First, inspired by the error analysis of Karakashian and Makridakis [23], ourproof makes use of properties of the Legendre polynomials on each interval I n , rewriting theGauss collocation method into a space-time Galerkin finite element method, which makes iteasier to choose suitable test functions in the error estimation. Second, we introduce a temporalRitz projection and use a super-approximation result of the temporal local L projection toeliminate the sub-optimal temporal consistency error caused by ∂ t u in the equation of r ; seeRemark 3.1 (hence, the proof of optimal order in Theorem 3.3 is one of our main contributions).Third, we estimate the time derivative of the error in H − ( Ω ) with a duality argument followingan H -norm error estimate. As a result, we obtain an optimal-order H -norm error estimatein the end. We prove the existence, uniqueness and optimal-order convergence of numericalsolutions based on Schaefer’s fixed point theorem in an L ∞ -neighborhood of the exact solution.This avoids grid-ratio conditions for the NLS equation with general nonlinearity.The rest of this paper is organized as follows. In Section 2, we present the SAV reformulationof the NLS equation and introduce our SAV space-time Gauss collocation finite element method.In Section 3, we first present an integral reformulation of the proposed numerical method andthen establish its mass and energy conservation properties. We also derive a consistency er-ror estimate for the proposed method, which is vitally used to prove an error estimate in thesubsequent section. In Section 4, we first establish the well-posedness of the numerical method3nd then prove an error bound of the form O ( h p + τ k +1 ) in the energy norm, where τ and h denote the temporal and spatial mesh sizes, respectively, with ( p, k ) denoting the degree ofpolynomials in the space-time finite element method. Finally, in Section 5, we present a fewnumerical experiments to validate the theoretical results, and to demonstrate the effectivenessof the proposed method in preserving the shape of a soliton wave.Throughout this paper, unless stated otherwise, C will be used to denote a generic positiveconstant which is independent of τ , h , n and N , but may depend on T and the regularity ofsolution. In this section, we construct a Gauss collocation finite element method based on the SAVreformulation of the NLS equation, and prove some desired properties of the numerical method,including the mass and energy conservation.
Let H k ( Ω ), k ≥
0, be the conventional complex-valued Sobolev space of functions on Ω , anddenote L ( Ω ) = H ( Ω ) and H ( Ω ) = { v ∈ H ( Ω ) : u = 0 on ∂Ω } . We denote by ( · , · ) and (cid:107) · (cid:107) the inner product and norm of the complex-valued Hilbert space L ( Ω ), respectively, defined by( u, v ) := (cid:90) Ω u v d x and (cid:107) u (cid:107) := (cid:112) ( u, u ) . For m, s ≥ ≤ p ≤ ∞ , the notation W m,p (0 , T ; H s ( Ω )) stands for the space-time Sobolevspace of functions which are W m,p in time and H s in space; see [12, Chapter 5.9]. We abbreviatethe norms of H s ( Ω ) and W m,p (0 , T ; H s ( Ω )) as (cid:107) · (cid:107) H k and (cid:107) · (cid:107) W m,p ( I n ; H s ) , respectively, omittingthe dependence on Ω in the subscripts. (1.1) The SAV formulation of the NLS equation (cf. [31]) introduces a scalar auxiliary variable r = (cid:113)(cid:82) Ω F ( | u | )d x + c with g ( u ) = f ( | u | ) (cid:113)(cid:82) Ω F ( | u | )d x + c , (2.5)with a positive c (which guarantees that the function r has a positive lower bound), andreformulate (1.1) as i∂ t u − ∆ u − rg ( u ) u = 0 in Ω × (0 , T ] , (2.6a)d r d t = Re (cid:0) g ( u ) u, ∂ t u (cid:1) in Ω × (0 , T ] , (2.6b) u = 0 on ∂Ω × (0 , T ] , (2.6c) u = u , r = r in Ω × { } , (2.6d)4here r = (cid:113)(cid:82) Ω F ( | u | )d x + c . The mass and energy conservation in the SAV formulationare dd t (cid:90) Ω | u | d x = 0 , and dd t (cid:18) (cid:90) Ω |∇ u | d x − r (cid:19) = 0 . (2.7) Let T h be a shape-regular and quasi-uniform triangulation of Ω with mesh size h ∈ (0 ,
1) and { t n } Nn =0 be a uniform partition of [0 , T ] with the time step size τ ∈ (0 , N is a positiveinteger and hence τ = TN . For an integer p ≥ Q p the space of complex-valuedpolynomials of degree ≤ p in space, and we denote by S h the complex-valued Lagrange finiteelement space subject to the triangulation of Ω , defined by S h = (cid:8) v ∈ C ( Ω ) : v | K ∈ Q p for all K ∈ T h , v = 0 on ∂Ω (cid:9) , where C ( Ω ) denotes the space of complex-valued uniformly continuous functions on Ω . Then S h is a complex Hilbert spaces with the inner product ( · , · ) and norm (cid:107) · (cid:107) .For an integer k ≥
1, let P k denote the space of real-valued polynomials of degree ≤ k in t .For a Banach space X , such as X = L ( Ω ) or X = S h , we define the following tensor-productspace: P k ⊗ X := span (cid:110) p ( t ) φ ( x ) : p ∈ P k , φ ∈ X (cid:111) = (cid:110)(cid:80) kj =0 t j φ j : φ j ∈ X (cid:111) . (2.8)Moreover, let P h : L ( Ω ) → S h denote the L projection operator defined by (cid:0) w − P h w, v h (cid:1) = 0 ∀ v h ∈ S h , ∀ w ∈ L ( Ω ) . The following stability properties are well-known (cf. [9]): (cid:107) P h w (cid:107) ≤ (cid:107) w (cid:107) ∀ w ∈ L ( Ω ) , (2.9a) (cid:107) P h w (cid:107) H ≤ C (cid:107) w (cid:107) H ∀ w ∈ H ( Ω ) , (2.9b)where C depends only on the shape-regularity and quasi-uniformity of the mesh.We also introduce the global space-time finite element spaces X τ,h = { v h ∈ C ([0 , T ]; S h ) : v h | I n ∈ P k ⊗ S h for n = 1 , . . . , N } , (2.10) Y τ,h = { q h ∈ C ([0 , T ]) : q h | I n ∈ P k for n = 1 , . . . , N } . (2.11) Let c j and w j , j = 1 , . . . , k , be the nodes and weights of the k -point Gauss quadrature rule inthe interval [ − ,
1] (see [33, Table 3.1]), and let t nj = t n − + (1 + c j ) τ / j = 1 , . . . , k denotethe Gauss points in the interval I n = [ t n − , t n ]. We define the following Gauss collocation finiteelement method for (2.6). Main Algorithm
Step 1:
Set u h := I h u and r h := r , where I h is the Lagrange interpolation operator ontothe finite element space. Determine ( u h , r h ) ∈ X τ,h × Y τ,h by the following two steps.5 tep 2: For n = 1 , , · · · , N , define (cid:8) ( u h ( t nj ) , r h ( t nj )) (cid:9) kj =1 ⊂ S h × R by solving recursively(in n ) the following nonlinear (algebraic) system:i (cid:0) ∂ t u h ( t nj ) , v h (cid:1) + (cid:0) ∇ u h ( t nj ) , ∇ v h (cid:1) (2.12a) − (cid:0) r h ( t nj ) g ( u h ( t nj )) u h ( t nj ) , v h (cid:1) = 0 , ∀ v h ∈ S h ,∂ t r h ( t nj ) = 12 Re (cid:0) g ( u h ( t nj )) , ∂ t u h ( t nj ) (cid:1) , (2.12b) u h ( t n − ) = u n − h and r h ( t n − ) = r n − h . (2.12c) Step 3:
Set u nh := u h ( t n ) and r nh := r h ( t n ). Remark 2.1 (a) We note that in (2.12a) and (2.12b), ∂ t u h ( t nj ) = ∂ t u h ( t ) | t = t nj and ∂ t r h ( t nj ) = ∂ t r h ( t ) | t = t nj . Main Algorithm really computes (cid:8) ( u h ( t nj ) , r h ( t nj )) (cid:9) kj =1 for each n ≥
1, however,since any k th order polynomial on I n is uniquely determined by its initial value at t n − and itsvalues at the k Gauss points t nj , j = 1 , . . . , k , then the Gauss-point values generated by MainAlgorithm uniquely determine the pair ( u h , r h ) ∈ X τ,h × Y τ,h .(b) Each of (2.12a) and (2.12b) consists of nonlinear algebraic equations, note that the testfunctions v h and q h can be different for different j , and one “side/initial condition” is prescribedfor each of u h and r h for each n . The number of equations imposed is the same as the degree offreedoms which equals the dimension of the space P k ⊗ S h for each n .(c) Main Algorithm can be obtained by applying the Gauss quadrature rule (in time) to a(continuous) space-time finite element method for (2.6). See Section 3.1. (2.12a) – (2.12b) In this subsection, we present several integral identities, including a reformulation of MainAlgorithm, and inequalities related to the proposed numerical method. These identities andinequalities will be used in the subsequent analysis of existence, uniqueness and convergence ofnumerical solutions.Consider the interval I n = [ t n − , t n ], then we define P nτ : L ( I n ; L ( Ω )) → P k − ⊗ L ( Ω ) tobe the L projection defined by (cid:90) I n ( u − P nτ u ) v d t = 0 ∀ v ∈ P k − ⊗ L ( Ω ) . (3.13)Thus u − P nτ u is orthogonal to all temporal polynomials of degree ≤ k −
1, which means that if u ∈ P k ⊗ L ( Ω ) then u − P nτ u = φ n − L k , (3.14)where φ n − ∈ L ( Ω ) and L k ( t ) := (cid:98) L k (cid:18) t − t n − − t n τ (cid:19) (3.15)6s the shifted Legendre polynomial (orthogonal to polynomials of lower degree on I n ). Thetemporal L projection operator P nτ has the following approximation property (cf. [10]):max t ∈ I n (cid:107) v − P nτ v (cid:107) X ≤ Cτ m max t ∈ I n (cid:107) ∂ mt v (cid:107) X , ≤ m ≤ k, (3.16)for all v ∈ C k ([0 , T ]; X ), where X = R or X = H s ( Ω ) for some s ∈ R .Since the k -point Gauss quadrature holds exactly for polynomials of degree 2 k − t nj , j = 1 , . . . , k , are the roots of the Legendre polynomial L k ( t )(cf. [24, p. 33]), it follows that the following two identities hold: (cid:90) I n v ( t )d t = τ k (cid:88) j =1 v ( t nj ) w j ∀ v ∈ P k − ⊗ S h , (3.17) v ( t nj ) = P nτ v ( t nj ) ∀ v ∈ P k ⊗ S h . (3.18)By choosing v h = τ v h ( t nj ) w j in (2.12a) and summing up the results for j = 1 , . . . , k , andusing (3.17)–(3.18) in the first two terms, we obtain the following integral identity: (cid:90) I n i (cid:0) ∂ t u h , v h (cid:1) d t + (cid:90) I n ( ∇ P nτ u h , ∇ v h )d t (3.19) − τ k (cid:88) j =1 w j ( r h ( t nj ) g ( u h ( t nj )) u h ( t nj ) , v h ( t nj )) = 0 ∀ v h ∈ P k ⊗ S h . Similarly, multiplying (2.12b) by τ q h ( t nj ) w j and summing up the results for j = 1 , . . . , k , andusing (3.17) in the first term, we have (cid:90) I n ∂ t r h q h d t = τ k (cid:88) j =1 w j (cid:0) g ( u h ( t nj )) u h ( t nj ) , ∂ t u h ( t nj ) q h ( t nj ) (cid:1) ∀ q h ∈ P k . (3.20)(3.19)–(3.20) provides a reformulation of Main Algorithm. It will be crucially used to showmass and energy conservations, as well as existence, uniqueness and convergence of numericalsolutions.From (3.14) we get (cid:107) φ n − (cid:107) = 1 | L k ( t n − ) | (cid:107) u h ( t n − ) − P nτ u h ( t n − ) (cid:107)≤ C (cid:107) u h ( t n − ) (cid:107) + C (cid:18) τ (cid:90) I n (cid:107) P nτ u h ( t ) (cid:107) d t (cid:19) , where we have used the inverse inequality in time. Thus, by using (3.14) again, we obtain thefollowing inequality: (cid:90) I n (cid:107) u h (cid:107) d t ≤ C (cid:90) I n (cid:107) P nτ u h (cid:107) d t + Cτ (cid:107) u h ( t n − ) (cid:107) ∀ u h ∈ P k ⊗ S h . (3.21)By using the two identities (3.17)–(3.18), one can also prove the following inequality: τ k (cid:88) j =1 w j (cid:107) v h ( t nj ) (cid:107) = (cid:90) I n (cid:107) P nτ v h ( t ) (cid:107) d t ≤ (cid:90) I n (cid:107) v h ( t ) (cid:107) d t ∀ v h ∈ P k ⊗ S h . (3.22)The inequalities (3.21)–(3.22) will be frequently used in the subsequent error analysis.7 .2 Mass and energy conservation properties In this subsection, we prove the following conservation properties of the numerical solution,which comprise of the first main theorem of this paper.
Theorem 3.1
Let ( u h , r h ) ∈ X τ,h × Y τ,h be a solution of Main Algorithm, then the followingmass and energy conservations hold:12 (cid:107) u h ( t n ) (cid:107) = 12 (cid:107) u h ( t ) (cid:107) for n ≥ , (cid:107)∇ u h ( t n ) (cid:107) − | r h ( t n ) | = 12 (cid:107)∇ u h ( t ) (cid:107) − | r h ( t ) | for n ≥ . Proof.
Setting v h = u h ∈ P k ⊗ S h in (3.19) and taking the imaginary part yieldIm (cid:90) I n i (cid:0) ∂ t u h , u h (cid:1) d t = − Im (cid:90) I n ( ∇ P nτ u h , ∇ u h )d t (3.23)+ Im (cid:20) τ k (cid:88) j =1 w j ( r h ( t nj ) g ( u h ( t nj )) , | u h ( t nj ) | ) (cid:21) = 0 , where we have used the definition of the projection operator P nτ , which impliesIm (cid:90) I n ( ∇ P nτ u h , ∇ u h )d t = Im (cid:90) I n ( ∇ P nτ u h , ∇ P nτ u h )d t = 0 . Then the mass conservation follows from (3.23) and the identityIm (cid:90) I n i (cid:0) ∂ t u h , u h (cid:1) d t = 12 (cid:107) u h ( t n ) (cid:107) − (cid:107) u h ( t n − ) (cid:107) . Alternatively, setting v h = ∂ t u h and q h = 2 r h in (3.19) and (3.20), respectively, and takingthe real parts yieldRe (cid:90) I n ( ∇ P nτ u h , ∇ ∂ t u h )d t = τ k (cid:88) j =1 w j ( r h ( t nj ) g ( u h ( t nj )) u h ( t nj ) , ∂ t u h ( t nj )) (3.24) | r h ( t n ) | − | r h ( t n − ) | = τ k (cid:88) j =1 w j (cid:0) r h ( t nj ) g ( u h ( t nj )) u h ( t nj ) , ∂ t u h ( t nj ) (cid:1) . (3.25)Since Re (cid:90) I n ( ∇ P nτ u h , ∇ ∂ t u h )d t = Re (cid:90) I n ( P nτ ∇ u h , ∇ ∂ t u h )d t = Re (cid:90) I n ( ∇ u h , ∇ ∂ t u h )d t = 12 (cid:107)∇ u h ( t n ) (cid:107) − (cid:107)∇ u h ( t n − ) (cid:107) , it follows that 12 (cid:107)∇ u h ( t n ) (cid:107) − (cid:107)∇ u h ( t n − ) (cid:107) (3.26)= τ k (cid:88) j =1 w j (cid:0) r h ( t nj ) g ( u h ( t nj )) u h ( t nj ) , ∂ t u h ( t nj ) (cid:1) . (cid:107)∇ u h ( t n ) (cid:107) − | r h ( t n ) | = 12 (cid:107)∇ u h ( t n − ) (cid:107) − | r h ( t n − ) | for n ≥ . (3.27)Thus, the energy conservation holds. The proof is complete. (cid:3) In this subsection, we prove that the average mass of numerical solutions at internal stages hasan upper bound unconditionally (independent of the regularity of solutions). This propertyfurthermore strengthens the stability of numerical solutions when the exact solution is notsmooth (for example, close to blow up).
Theorem 3.2
Let ( u h , r h ) ∈ X τ,h × Y τ,h be a solution of Main Algorithm, then the followinginequalities hold: max ≤ n ≤ N τ (cid:90) I n (cid:107) P nτ u h (cid:107) d t ≤ (cid:107) u h (0) (cid:107) , (3.28a)max ≤ n ≤ N max ≤ j ≤ k (cid:107) u h ( t nj ) (cid:107) ≤ C (cid:107) u h (0) (cid:107) , (3.28b)where C is a constant independent of τ , h and the regularity of the solution. Proof.
By the definition of the temporal L projection P nτ , we get (cid:90) I n (cid:107) P nτ u h ( t ) (cid:107) d t = Re (cid:90) I n ( u h ( t ) , P nτ u h ( t ))d t (3.29)= Re (cid:90) I n ( u h ( t n − ) , P nτ u h ( t n − ))d t + Re (cid:90) I n (cid:90) tt n − [( ∂ s u h ( s ) , P nτ u h ( s )) + ( u h ( s ) , ∂ s P nτ u h ( s ))]d s d t = Re ( u h ( t n − ) , P nτ u h ( t n − )) τ + Re (cid:90) I n ( ∂ t u h ( t ) , ( t n − t ) P nτ u h ( t ))d t + Re (cid:90) I n ( u h ( t ) , ( t n − t ) ∂ t P nτ u h ( t ))d t =: J + J + J , where we have interchanged the order of integration in deriving the second to last equality.Using H¨older’s and Young’s inequalities, we have J = Re ( u h ( t n − ) , P nτ u h ( t n − )) τ ≤ (cid:107) u h ( t n − ) (cid:107)(cid:107) P nτ u h ( t n − ) (cid:107) τ ≤ τ (cid:107) u h ( t n − ) (cid:107) + τ (cid:107) P nτ u h ( t n − ) (cid:107) . Setting v h = ( t n − t ) P nτ u h in (3.19) and taking the imaginary part yield J = Re (cid:90) I n ( ∂ t u h ( t ) , ( t n − t ) P nτ u h ( t ))d t = Im (cid:90) I n i (cid:0) ∂ t u h ( t ) , ( t n − t ) P nτ u h ( t ) (cid:1) d t = Im τ k (cid:88) j =1 w j (cid:0) r h ( t nj ) g ( u h ( t nj )) , | u h ( t nj ) | (cid:1) ( t n − t nj )9 Im (cid:90) I n (cid:107)∇ P nτ u h (cid:107) ( t n − t )d t = 0 , where we have used (3.18) in deriving the first term on the right-hand side. Since ( t n − t ) ∂ t P nτ u h ( t ) is a polynomial of degree k − t , it follows that (cid:90) I n ( u h ( t ) , ( t n − t ) ∂ t P nτ u h ( t ))d t = (cid:90) I n ( P nτ u h ( t ) , ( t n − t ) ∂ t P nτ u h ( t ))d t. Hence, J = Re (cid:90) I n ( P nτ u h ( t ) , ( t n − t ) ∂ t P nτ u h ( t ))d t = (cid:90) I n
12 dd t (cid:107) P nτ u h ( t ) (cid:107) ( t n − t )d t = − τ (cid:107) P nτ u h ( t n − ) (cid:107) + (cid:90) I n (cid:107) P nτ u h ( t ) (cid:107) d t. Substituting the estimates of J , J and J into (3.29), we obtain (cid:90) I n (cid:107) P nτ u h (cid:107) d t ≤ τ (cid:107) u h ( t n − ) (cid:107) = τ (cid:107) u h (0) (cid:107) , (3.30)where the last equality follows from mass conservation proved in Theorem 3.1. This proves(3.28a) holds.Substituting (3.30) into (3.21) and using the mass conservation again, we obtain (cid:82) I n (cid:107) u h (cid:107) d t ≤ Cτ (cid:107) u h (0) (cid:107) . Then, by using the inverse inequality, we obtainmax t ∈ I n (cid:107) u h ( t ) (cid:107) ≤ C (cid:107) u h (0) (cid:107) , which proves (3.28b). The proof is complete. (cid:3) Let I nτ u and I nτ r be the temporal Lagrange interpolation polynomials of u and r , respectively,interpolated at the k + 1 points t n − and t nj , j = 1 , . . . , k . Both I nτ u and I nτ r are temporalpolynomials of degree ≤ k . The one-dimensional temporal Lagrange interpolation operator I nτ has the following approximation property (cf. [10]):max t ∈ I n (cid:0) (cid:107) v − I nτ v (cid:107) X + τ (cid:107) ∂ t ( v − I nτ v ) (cid:107) X (cid:1) ≤ Cτ m +1 max t ∈ I n (cid:107) ∂ m +1 t v (cid:107) X (3.31)for all v ∈ C m +1 ([0 , T ]; X ) , ≤ m ≤ k , and X = R or X = H s ( Ω ) for some s ∈ R .To analyze the error of numerical solutions due to temporal discretization, we define a tempo-ral Ritz projection operator R nτ : W , ∞ ( I n ; L ( Ω )) → P k ⊗ L ( Ω ) by the following two conditions: (cid:90) I n ( ∂ t ( u − R nτ u ) , v )d t = 0 ∀ v ∈ P k − ⊗ L ( Ω ) , (3.32) u ( t n − ) − R nτ u ( t n − ) = 0 . (3.33)Clearly, ∂ t R nτ u is the temporal L projection of ∂ t u onto P k − ⊗ L ( Ω ). By using this prop-erty and the shifted Legendre polynomials defined in (3.15), we can express the temporal Ritzprojection as R nτ u ( t ) = u ( t n − ) + k − (cid:88) j =0 (cid:82) I n L j ( s ) ∂ s u ( s )d s (cid:82) I n | L j ( s ) | d s (cid:90) tt n − L j ( s )d s. (3.34)10his expression implies that if X ⊃ L ( Ω ) is a Banach space and u ∈ W , ∞ ( I n ; X ), then R nτ u is automatically in P k ⊗ X . Meanwhile, this temporal Ritz projection has the followingapproximation property. Lemma 3.1
Let X = R or H s ( Ω ) for some s ≥
0. For u ∈ W m +1 , ∞ ( I n ; X ), with 0 ≤ m ≤ k ,the following approximation property holds: (cid:107) u − R nτ u (cid:107) L ∞ ( I n ; X ) + τ (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; X ) ≤ Cτ m +1 (cid:107) u (cid:107) W m +1 , ∞ ( I n ; X ) . Proof.
We prove the result for the case X = H s ( Ω ) with s ≥
0. To this end, we denote by H s ( Ω ) (cid:48) the dual space of H s ( Ω ). Then, by the Riesz representation theorem, there exists acontinuous linear bijection J : H s ( Ω ) (cid:48) → H s ( Ω ) such that( u, J v ) H s = (cid:104) u, v (cid:105) ∀ u ∈ H s ( Ω ) and v ∈ H s ( Ω ) (cid:48) , (3.35)where ( · , · ) H s is the inner product of H s ( Ω ), and (cid:104)· , ·(cid:105) denotes the pairing between H s ( Ω ) andits dual space H s ( Ω ) (cid:48) , satisfying (cid:104) w, v (cid:105) = ( w, v ) ∀ w ∈ H s ( Ω ) , v ∈ L ( Ω ) (cid:44) → H s ( Ω ) (cid:48) . Then (3.32) implies that (cid:90) I n ( ∂ t ( u − R nτ u ) , J v ) H s d t = 0 ∀ v ∈ P k − ⊗ L ( Ω ) . (3.36)Since L ( Ω ) is dense in H s ( Ω ) (cid:48) , it follows that (3.36) actually holds for all v ∈ P k − ⊗ H s ( Ω ) (cid:48) .Since for any w ∈ P k − ⊗ H s ( Ω ) there exists v ∈ P k − ⊗ H s ( Ω ) (cid:48) satisfying J v = w , it followsthat (3.36) can be equivalently written as (cid:90) I n ( ∂ t ( u − R nτ u ) , w ) H s d t = 0 ∀ w ∈ P k − ⊗ H s ( Ω ) . (3.37)By using the inverse inequality in time and the property (3.37), we have (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; H s ) ≤ Cτ − (cid:90) I n (cid:0) ∂ t ( u − R nτ u ) , ∂ t ( u − R nτ u ) (cid:1) H s d t (3.38)= Cτ − (cid:90) I n (cid:0) ∂ t ( u − R nτ u ) , ∂ t ( u − g ) (cid:1) H s d t ≤ C (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; H s ) (cid:107) ∂ t ( u − g ) (cid:107) L ∞ ( I n ; H s ) , where g can be an arbitrary function in P k ⊗ H s ( Ω ), which implies ∂ t ( g − R nτ u ) ∈ P k − ⊗ H s ( Ω ).The inequality above implies (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; H s ) ≤ C inf g ∈ P k ⊗ H s ( Ω ) (cid:107) ∂ t ( u − g ) (cid:107) L ∞ ( I n ; H s ) (3.39) ≤ C (cid:107) ∂ t ( u − I nτ u ) (cid:107) L ∞ ( I n ; H s ) . This and the approximation property (3.31) together imply the following inequality: (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; H s ) ≤ Cτ m (cid:107) u (cid:107) W m +1 , ∞ ( I n ; H s ) . (3.40)11o estimate (cid:107) u − R nτ u (cid:107) L ∞ ( I n ; H s ) , we use a duality argument in time. Let w ∈ W , ∞ ( I n ; H s ( Ω ))be the solution of the following IVP: ∂ t w ( t ) = u ( t ) − ( R nτ u )( t ) , and w ( t n ) = 0 . (3.41)Then, testing (3.41) by J − ( u − R nτ u ), with the operator J defined in (3.35), we obtain (cid:90) I n (cid:107) u − R nτ u (cid:107) H s d t = (cid:90) I n (cid:0) ∂ t w, u − R nτ u (cid:1) H s d t = − (cid:90) I n (cid:0) w, ∂ t ( u − R nτ u ) (cid:1) H s d t = − (cid:90) I n (cid:0) w − I nτ w, ∂ t ( u − R nτ u ) (cid:1) H s d t (here (3.37) is used) ≤ Cτ (cid:107) w − I nτ w (cid:107) L ∞ ( I n ; H s ) (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; H s ) ≤ Cτ m +2 (cid:107) ∂ t w (cid:107) L ∞ ( I n ; H s ) (cid:107) u (cid:107) W m +1 , ∞ ( I n ; H s ) , where we have used (3.40) in the last inequality. By applying the temporal inverse inequalityand (3.41), we have (cid:107) u − R nτ u (cid:107) L ∞ ( I n ; H s ) ≤ Cτ − (cid:107) u − R nτ u (cid:107) L ( I n ; H s ) (3.42) ≤ Cτ m +1 (cid:107) ∂ t w (cid:107) L ∞ ( I n ; H s ) (cid:107) u (cid:107) W m +1 , ∞ ( I n ; H s ) ≤ Cτ m +1 (cid:107) u − R nτ u (cid:107) L ∞ ( I n ; H s ) (cid:107) u (cid:107) W m +1 , ∞ ( I n ; H s ) , which leads to (cid:107) u − R nτ u (cid:107) L ∞ ( I n ; H s ) ≤ Cτ m +1 (cid:107) u (cid:107) W m +1 , ∞ ( I n ; H s ) . (3.43)This completes the proof of Lemma 3.1. (cid:3) In addition to the above optimal-order approximation result, we also have the followingsuper-convergence result.
Lemma 3.2 (A super-approximation property) Let X = R or H s ( Ω ) for some s ≥
0. If w ∈ W k, ∞ ( I n ; W s, ∞ ( Ω )) and v ∈ P k − ⊗ X , then (cid:107) wv − P nτ ( wv ) (cid:107) L ( I n ; X ) ≤ Cτ (cid:107) v (cid:107) L ( I n ; X ) . Proof.
We only give a proof for the case X = H s ( Ω ) because the other cases are similar. Byapplying (3.16) with m = k , we have (cid:107) wv − P nτ ( wv ) (cid:107) L ( I n ; H s ) ≤ Cτ (cid:107) wv − P nτ ( wv ) (cid:107) L ∞ ( I n ; H s ) ≤ Cτ k + (cid:107) ∂ kt ( wv ) (cid:107) L ∞ ( I n ; H s ) ≤ C k − (cid:88) m =0 τ k + (cid:107) ∂ k − mt w∂ mt v (cid:107) L ∞ ( I n ; H s ) (since ∂ kt v = 0) ≤ C k − (cid:88) m =0 τ k + (cid:107) ∂ k − mt w (cid:107) L ∞ ( I n ; W s, ∞ ) (cid:107) ∂ mt v (cid:107) L ∞ ( I n ; H s ) C k − (cid:88) m =0 τ k + − m (cid:107) v (cid:107) L ∞ ( I n ; H s ) ≤ Cτ (cid:107) v (cid:107) L ∞ ( I n ; H s ) ≤ Cτ (cid:107) v (cid:107) L ( I n ; H s ) . here we have used the inverse inequality in time twice above. The proof is complete. (cid:3) Finally, we also recall the (spatial) Ritz projection operator R h : H ( Ω ) → S h defined by (cid:0) ∇ ( w − R h w ) , ∇ v h (cid:1) = 0 ∀ v h ∈ S h , ∀ w ∈ H ( Ω ) , and the discrete Laplacian operator ∆ h : S h → S h defined by(∆ h φ h , χ h ) := − ( ∇ φ h , ∇ χ h ) ∀ φ h , χ h ∈ S h . (3.44)It is known [9] that there hold the following identities: P h ∆ v = ∆ h R h v ∀ v ∈ H ( Ω ) , (3.45a) R nτ R h v = R h R nτ v ∀ v ∈ W , ∞ ( I n ; H ( Ω )) , (3.45b) R nτ ∆ h v h = ∆ h R nτ v h ∀ v ∈ W , ∞ ( I n ; S h ) . (3.45c)Moreover, there holds the following approximation property (cf. [9]): (cid:107) v − R h v (cid:107) H ≤ Ch p (cid:107) v (cid:107) H p +1 ∀ v ∈ H ( Ω ) ∩ H p +1 ( Ω ) . (3.46) (2.12a) – (2.12b) We define a pair of intermediate solutions (for comparison with the numerical solutions) u ∗ h = R nτ R h u and r ∗ h = R nτ r. Then, testing (2.6a) and (2.6b) by P nτ v h and P nτ q h , respectively, we obtain the following equa-tions for u ∗ h and r ∗ h : (cid:90) I n i (cid:0) ∂ t u ∗ h , P nτ v h (cid:1) d t + (cid:90) I n (cid:0) ∇ u ∗ h , ∇ P nτ v h (cid:1) d t (3.47) − τ k (cid:88) j =1 w j (cid:0) r ∗ h ( t nj ) g ( u ∗ h ( t nj )) u ∗ h ( t nj ) , P nτ v h ( t nj ) (cid:1) = (cid:90) I n ( d nu , P nτ v h )d t, (cid:90) I n ∂ t r ∗ h P nτ q h d t = τ k (cid:88) j =1 w j Re (cid:0) P nτ q h ( t nj ) g ( u ∗ h ( t nj )) u ∗ h ( t nj ) , ∂ t u ∗ h ( t nj ) (cid:1) (3.48)+ (cid:90) I n d nr P nτ q h d t, where d nu and d nr are consistency errors of the numerical method. By using the identities (3.17)and (3.32), we can express these consistency errors by d nu =i ∂ t R nτ ( R h u − u ) + ∆ h R h ( u − R nτ u ) + rg ( u ) u − I nτ [ r ∗ h g ( u ∗ h ) u ∗ h ] , (3.49)13 nr = 12 Re (cid:2)(cid:0) g ( u ) u, ∂ t u (cid:1) − I nτ (cid:0) g ( u ∗ h ) u ∗ h , ∂ t u ∗ h (cid:1)(cid:3) . (3.50)After using (3.13) and (3.18), equations (3.47)–(3.48) can be rewritten as (cid:90) I n i (cid:0) ∂ t u ∗ h , v h (cid:1) d t + (cid:90) I n (cid:0) ∇ P nτ u ∗ h , ∇ v h (cid:1) d t (3.51) − τ k (cid:88) j =1 w j (cid:0) r ∗ h ( t nj ) g ( u ∗ h ( t nj )) u ∗ h ( t nj ) , v h ( t nj ) (cid:1) = (cid:90) I n ( P nτ d nu , v h )d t, (cid:90) I n ∂ t r ∗ h q h d t = τ k (cid:88) j =1 w j Re (cid:0) q h ( t nj ) g ( u ∗ h ( t nj )) u ∗ h ( t nj ) , ∂ t u ∗ h ( t nj ) (cid:1) (3.52)+ (cid:90) I n P nτ d nr q h d t. Theorem 3.3
Suppose that the solution of (1.1) is sufficiently smooth, then d nu ∈ C ( I n ; H ( Ω ))and there hold sup t ∈ I n (cid:107) d nu (cid:107) H ≤ C ( h p + τ k +1 ) and sup t ∈ I n | P nτ d nr | ≤ C ( h p + τ k +1 ) . (3.53) Remark 3.1
The key observation for the consistency errors is that, although (3.50) containsan O ( h p + τ k ) error from ∂ t u − ∂ t u ∗ h , the temporal L projection operator P nτ acting on d nr furthermore reduces this error to O ( h p + τ k +1 ). This is proved by using the super-approximationresult in Lemma 3.2. Proof.
Since the spatial Ritz projection R h maps H ( Ω ) into S h ⊂ H ( Ω ), and the temporalRitz projection R nτ maps W , ∞ ( I n ; H ( Ω )) into P k ⊗ H ( Ω ), it follows that every term in (3.49)is in C ( I n ; H ( Ω )). This implies d nu ∈ C ( I n ; H ( Ω )).By using the triangle inequality, from (3.49) we getmax t ∈ I n (cid:107) d nu (cid:107) H ≤ max t ∈ I n (cid:0) (cid:107) ∂ t R nτ ( R h u − u ) (cid:107) H + (cid:107) ∆ h R h ( u − R nτ u ) (cid:107) H (cid:1) (3.54)+ max t ∈ I n (cid:0) (cid:107) rg ( u ) u − I nτ [ rg ( u ) u ] (cid:107) H + (cid:107) rg ( u ) u − r ∗ h g ( u ∗ h ) u ∗ h (cid:107) H (cid:1) = : D u + D u + D u + D u . Choosing m = 0 in Lemma 3.1, we obtain the following stability result: (cid:107) R nτ u (cid:107) W , ∞ ( I n ; H s ) ≤ C (cid:107) u (cid:107) W , ∞ ( I n ; H s ) . (3.55)Using (3.55) and (3.46), we can estimate D u as follows: D u = max t ∈ I n (cid:107) ∂ t R nτ ( R h u − u ) (cid:107) H ≤ (cid:107) R h u − u (cid:107) W , ∞ ( I n ; H ) ≤ Ch p (cid:107) R h u − u (cid:107) W , ∞ ( I n ; H p +1 ) . Similarly, using identity (3.45) and Lemma 3.1, we have D u = max t ∈ I n (cid:107) ∆ h R h ( u − R nτ u ) (cid:107) H = max t ∈ I n (cid:107) P h ∆( u − R nτ u ) (cid:107) H max t ∈ I n (cid:107) u − R nτ u (cid:107) H ≤ Cτ k +1 (cid:107) u (cid:107) W k +1 , ∞ ( I n ; H ) , and D u = max t ∈ I n (cid:107) rg ( u ) u − I nτ [ rg ( u ) u ] (cid:107) H ≤ Cτ k +1 . By using the triangle inequality, we decompose D u into two parts, D u ≤ max t ∈ I n (cid:0) (cid:107) rg ( u ) u − rg ( R h u ) R h u (cid:107) H + (cid:107) rg ( R h u ) R h u − R nτ rg ( R nτ R h u ) R nτ R h u (cid:107) H (cid:1) ≤ Ch p + Cτ k +1 . Then, substituting the estimates of D uj , j = 1 , , ,
4, into (3.54), we obtain the desired estimatefor (cid:107) d nu (cid:107) H .To estimate | P nτ d nr | , we rewrite (3.50) as d nr = 12 Re (cid:104)(cid:0) g ( u ) u, ∂ t ( u − u ∗ h ) (cid:1) + (cid:0) g ( u ) u − g ( u ∗ h ) u ∗ h , ∂ t u ∗ h (cid:1)(cid:105) + 12 Re (cid:104)(cid:0) g ( u ∗ h ) u ∗ h , ∂ t u ∗ h (cid:1) − I nτ (cid:0) g ( u ∗ h ) u ∗ h , ∂ t u ∗ h (cid:1)(cid:105) and test this expression by P nτ v in the time interval I n , with v ∈ P k . This yields (cid:90) I n P nτ d nr v d t = (cid:90) I n d nr P nτ v d t (3.56) ≤
12 Re (cid:90) I n (cid:0) g ( u ) u, ∂ t ( u − u ∗ h ) (cid:1) P nτ v d t + Cτ (cid:107) ( g ( u ) u − g ( u ∗ h ) u ∗ h , ∂ t u ∗ h ) (cid:107) L ∞ ( I n ) (cid:107) v (cid:107) L ( I n ) + Cτ k + (cid:107) ∂ k +1 t (cid:0) g ( u ∗ h ) u ∗ h , ∂ t u ∗ h (cid:1) (cid:107) L ∞ ( I n ) (cid:107) v (cid:107) L ( I n ) ≤
12 Re (cid:90) I n (cid:0) g ( u ) u, ∂ t ( u − u ∗ h ) (cid:1) P nτ v d t + Cτ ( h p + τ k +1 ) (cid:107) v (cid:107) L ( I n ) . The first term on the right-hand side of (3.56) can be estimated as follows.12 Re (cid:90) I n (cid:0) g ( u ) u, ∂ t ( u − u ∗ h ) (cid:1) P nτ v d t = (cid:90) I n ( g ( u ) u, ∂ t ( u − R nτ u )) P nτ v d t (3.57)+ (cid:90) I n ( g ( u ) u, ∂ t R nτ ( u − R h u )) P nτ v d t = : D r + D r , where D r = (cid:90) I n (cid:0) g ( u ) uP nτ v, ∂ t ( u − R nτ u ) (cid:1) d t = (cid:90) I n (cid:0) g ( u ) uP nτ v − P nτ ( g ( u ) uP nτ v ) , ∂ t ( u − R nτ u ) (cid:1) d t ≤ Cτ (cid:107) g ( u ) uP nτ v − P nτ ( g ( u ) uP nτ v ) (cid:107) L ( I n ; L ) (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; L ) Cτ (cid:107) P nτ v (cid:107) L ( I n ; L ) (cid:107) ∂ t ( u − R nτ u ) (cid:107) L ∞ ( I n ; L ) (we have used Lemma 3.2) ≤ Cτ k + (cid:107) v (cid:107) L ( I n ) (cid:107) ∂ k +1 t u (cid:107) L ∞ ( I n ; L ) (we have used Lemma 3.1) ,D r ≤ Cτ (cid:107) g ( u ) u (cid:107) L ∞ ( I n ; L ) (cid:107) u − R h u (cid:107) W , ∞ ( I n ; L ) (cid:107) v (cid:107) L ( I n ) ≤ Cτ h p (cid:107) v (cid:107) L ( I n ) (cid:107) u (cid:107) W , ∞ ( I n ; H p +1 ) . Substituting these estimates into (3.56), we obtain (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) I n P nτ d nr v d t (cid:12)(cid:12)(cid:12)(cid:12) ≤ Cτ ( h p + τ k +1 ) (cid:107) v (cid:107) L ( I n ) . Since this inequality holds for arbitrary v ∈ L ( I n ), it follows that (cid:107) P nτ d nr (cid:107) L ( I n ) ≤ Cτ ( h p + τ k +1 ) . Then, using the inverse inequality in time, we obtain the desired estimate for | P nτ d nr | . (cid:3) We define the error functions e uh = u h − u ∗ h and e rh = r h − r ∗ h , with the following abbreviations: e unj = e uh ( t nj ) and e rnj = e rh ( t nj ) ,u nj = u h ( t nj ) and r nj = r h ( t nj ) ,u ∗ nj = u ∗ h ( t nj ) and r ∗ nj = r ∗ h ( t nj ) ,v nj = v h ( t nj ) and q nj = q h ( t nj ) . Subtracting (3.51)–(3.52) from (3.19)–(3.20), we obtain the following error equations: i (cid:90) I n (cid:0) ∂ t e uh , v h (cid:1) d t = − (cid:90) I n (cid:0) ∇ P nτ e uh , ∇ v h (cid:1) d t + τ k (cid:88) j =1 w j (cid:16) e rnj g ( u nj ) u nj , v nj (cid:17) + τ k (cid:88) j =1 w j (cid:16) r ∗ nj (cid:2) g ( u nj ) u nj − g ( u ∗ nj ) u ∗ nj (cid:3) , v nj (cid:17) − (cid:90) I n ( P nτ d nu , v h )d t, (4.58a) (cid:90) I n ∂ t e rh q h d t = τ k (cid:88) j =1 w j Re (cid:0) q nj (cid:0) g ( u nj ) u nj − g ( u ∗ nj ) u ∗ nj (cid:1) , ∂ t u ∗ h ( t nj ) (cid:1) + τ k (cid:88) j =1 w j Re (cid:0) q nj g ( u nj ) u nj , ∂ t e uh ( t nj ) (cid:1) − (cid:90) I n P nτ d nr q h d t, (4.58b)which hold for all test functions v h ∈ P k ⊗ S h and q h ∈ P k . Remark 4.1
If (4.58) has a solution ( e uh , e rh ) ∈ X τ,h × Y τ,h with u h = u ∗ h + e uh and r h = r ∗ h + e rh ,then ( u h , r h ) is a solution of the numerical scheme (2.12). In the following, we prove existenceof a solution ( e uh , e rh ) to (4.58) with u h = u ∗ h + e uh and r h = r ∗ h + e rh .In this section, we prove existence and uniqueness of solutions to (4.58a)–(4.58b) by usingSchaefer’s Fixed Point Theorem, which is quoted below.16 heorem 4.1 (Schaefer’s Fixed Point Theorem [12, Chapter 9.2, Theorem 4]) Let B be aBanach space and let M : B → B be a continuous and compact mapping (possibly nonlinear).If the set (cid:8) φ ∈ B : ∃ θ ∈ [0 ,
1] such that φ = θM ( φ ) (cid:9) (4.59)is bounded in B , then the mapping M has at least one fixed point.We define X ∗ τ,h = (cid:110) v h ∈ X τ,h : max ≤ n ≤ N max ≤ j ≤ k (cid:107) v h ( t nj ) − u ∗ h ( t nj ) (cid:107) L ∞ ∩ H ≤ (cid:111) , (4.60) Y ∗ τ,h = (cid:110) q h ∈ Y τ,h : max ≤ n ≤ N max ≤ j ≤ k | q h ( t nj ) − r ∗ h ( t nj ) | ≤ (cid:111) , (4.61)where the norm (cid:107) · (cid:107) L ∞ ∩ H is defined as (cid:107) φ h (cid:107) L ∞ ∩ H := max (cid:0) (cid:107) φ h (cid:107) L ∞ , (cid:107) φ h (cid:107) H (cid:1) . For any element ( φ h , ϕ h ) ∈ X τ,h × Y τ,h , we define two associated numbers ρ [ φ h ] := min (cid:18) ≤ n ≤ N max ≤ j ≤ k (cid:107) φ h ( t nj ) (cid:107) L ∞ ∩ H , (cid:19) , (4.62a) ρ [ ϕ h ] := min (cid:18) ≤ n ≤ N max ≤ j ≤ k | ϕ h ( t nj ) | , (cid:19) , (4.62b)which are continuous with respect to ( φ h , ϕ h ) (because all norms are equivalent in the finite-dimensional space X τ,h × Y τ,h ). Furthermore, the two numbers defined above satisfy the followingestimates: max ≤ n ≤ N max ≤ j ≤ k (cid:107) ρ [ φ h ] φ h ( t nj ) (cid:107) L ∞ ∩ H ≤ , (4.63)max ≤ n ≤ N max ≤ j ≤ k | ρ [ ϕ h ] ϕ h ( t nj ) | ≤ . (4.64)Then we define u φ := u ∗ h + ρ [ φ h ] φ h and r ϕ := r ∗ h + ρ [ ϕ h ] ϕ h , (4.65)with the following abbreviations: u φnj = u φh ( t nj ) and ϕ nj = ϕ h ( t nj ) , and define ( e uh , e rh ) ∈ X τ,h × Y τ,h to be the solution of the following linear equations:i (cid:90) I n (cid:0) ∂ t e uh , v h (cid:1) d t + (cid:90) I n (cid:0) ∇ P nτ e uh , ∇ v h (cid:1) d t = τ k (cid:88) j =1 w j (cid:16) ϕ nj g ( u φnj ) u φnj , v nj (cid:17) (4.66)+ τ k (cid:88) j =1 w j (cid:16) r ∗ nj (cid:2) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:3) , v nj (cid:17) − (cid:90) I n ( P nτ d u , v h )d t (cid:90) I n ∂ t e rh q h d t = τ k (cid:88) j =1 w j Re (cid:0) q nj (cid:0) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:1) , ∂ t u ∗ h ( t nj ) (cid:1) (4.67)+ τ k (cid:88) j =1 w j Re (cid:0) q nj g ( u φnj ) u φnj , ∂ t φ h ( t nj ) (cid:1) − (cid:90) I n P nτ d r q h d t for all v h ∈ P k ⊗ S h and q h ∈ P k , n = 1 , . . . , N . We denote by M : X τ,h × Y τ,h → X τ,h × Y τ,h themapping from ( φ h , ϕ h ) to ( e uh , e rh ), and define the set B = (cid:8) ( φ h , ϕ h ) ∈ X τ,h × Y τ,h : ∃ θ ∈ [0 ,
1] such that ( φ h , ϕ h ) = θM ( φ h , ϕ h ) (cid:9) , (4.68)and the norm following norm on X τ,h × Y τ,h : for any ( φ h , ϕ h ) ∈ X τ,h × Y τ,h (cid:107) ( φ h , ϕ h ) (cid:107) X τ,h × Y τ,h := (cid:107) φ h (cid:107) L ∞ (0 ,T ; H ) + (cid:107) ϕ h (cid:107) L ∞ (0 ,T ) . (4.69) Lemma 4.1
The mapping M : X τ,h × Y τ,h → X τ,h × Y τ,h is well defined, continuous andcompact. Proof.
Since the right-hand sides of (4.66)–(4.67) are given, the linear equations (4.66)–(4.67)have a unique solution ( e uh , e rh ) ∈ X τ,h × Y τ,h for any given ( φ h , ϕ h ) ∈ X τ,h × Y τ,h . Thus themapping is well defined. Let (cid:96) u ( φ h , ϕ h ; v h ) and (cid:96) r ( φ h , ϕ h ; q h ) denote the right-hand sides of(4.66) and (4.67), respectively. Since all norms are equivalent in the finite-dimensional space X τ,h × Y τ,h , it follows that | (cid:96) u ( φ h , ϕ h ; v h ) − (cid:96) u ( ˆ φ h , ˆ ϕ h ; v h ) | ≤ o (1) (cid:107) v h (cid:107) L (0 ,T ; L ) , ∀ v h ∈ L (0 , T ; L ) , | (cid:96) r ( φ h , ϕ h ; q h ) − (cid:96) r ( ˆ φ h , ˆ ϕ h ; q h ) | ≤ o (1) (cid:107) q h (cid:107) L (0 ,T ) ∀ q h ∈ L (0 , T ; L ) , as ( φ h , ϕ h ) → ( ˆ φ h , ˆ ϕ h ) in X τ,h × Y τ,h . Where o (1) represents some quantity tending to zero.Using this property, it is easy to verify that ( e uh , e rh ) is continuous with respect to ( φ h , ϕ h ).Since X τ,h × Y τ,h is a finite-dimensional space, a continuous mapping is automatically com-pact. The proof is complete. (cid:3) We are now ready to state and prove the following key lemma.
Lemma 4.2
Let 1 ≤ d ≤ τ and h such that when τ ≤ τ and h ≤ h ,the following statement holds: If ( φ h , ϕ h ) ∈ B and ( e uh , e rh ) = M ( φ h , ϕ h ), then (cid:107) e uh (cid:107) L ∞ (0 ,T ; H ) + (cid:107) e rh (cid:107) L ∞ (0 ,T ) ≤ (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) , (4.70)max ≤ n ≤ N max ≤ j ≤ k (cid:107) e uh ( t nj ) (cid:107) L ∞ ∩ H ≤
12 and max ≤ n ≤ N max ≤ j ≤ k | e rh ( t nj ) | ≤ , (4.71) ρ [ φ h ] = 1 , ρ [ ϕ h ] = 1 . (4.72) Proof.
If ( φ h , ϕ h ) ∈ B and ( e uh , e rh ) = M ( φ h , ϕ h ), then( φ h , ϕ h ) = θM ( φ h , ϕ h ) = ( θe uh , θe rh ) , φ h = θe uh and ϕ h = θe rh . In this case, (4.66)–(4.67) can be rewritten asi (cid:90) I n (cid:0) ∂ t e uh , v h (cid:1) d t = − (cid:90) I n (cid:0) ∇ P nτ e uh , ∇ v h (cid:1) d t + θτ k (cid:88) j =1 w j (cid:16) e rnj g ( u φnj ) u φnj , v nj (cid:17) (4.73)+ τ k (cid:88) j =1 w j (cid:16) r ∗ nj (cid:2) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:3) , v nj (cid:17) − (cid:90) I n ( P nτ d nu , v h )d t, (cid:90) I n ∂ t e rh q h d t = τ k (cid:88) j =1 w j Re (cid:0) q nj (cid:0) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:1) , ∂ t u ∗ h ( t nj ) (cid:1) (4.74)+ θτ k (cid:88) j =1 w j Re (cid:0) q nj g ( u φnj ) u φnj , ∂ t e uh ( t nj ) (cid:1) − (cid:90) I n P nτ d nr q h d t, which hold for all v h ∈ P k ⊗ S h and q h ∈ P k , n = 1 , . . . , N . In the following, we derive estimatesfor e uh and e rh based on equations (4.73)–(4.74).From (4.63)–(4.64) and definition (4.65) we getmax ≤ n ≤ N max ≤ j ≤ k (cid:107) u φ ( t nj ) (cid:107) L ∞ ∩ H + max ≤ n ≤ N max ≤ j ≤ k | r ϕ ( t nj ) | (4.75) ≤ max ≤ n ≤ N max ≤ j ≤ k (cid:107) u ∗ h ( t nj ) (cid:107) L ∞ ∩ H + max ≤ n ≤ N max ≤ j ≤ k | r ∗ h ( t nj ) | + max ≤ n ≤ N max ≤ j ≤ k (cid:107) ρ [ φ h ] φ h ( t nj ) (cid:107) L ∞ ∩ H + max ≤ n ≤ N max ≤ j ≤ k | ρ [ ϕ h ] ϕ h ( t nj ) |≤ (cid:107) u ∗ h (cid:107) L ∞ (0 ,T ; L ∞ ∩ H ) + (cid:107) r ∗ h (cid:107) L ∞ (0 ,T ) + 2 . Thus (cid:107) u φ ( t nj ) (cid:107) L ∞ ∩ H and | r ϕ ( t nj ) | are bounded uniformly with respect to τ and h .Since the remainder of the proof is long, we divide it into four steps. Step 1: Estimation of (cid:107) e uh (cid:107) L ( I n ; H ) . Note that (cid:90) I n (cid:107)∇ P nτ e uh ( t ) (cid:107) d t = (cid:90) I n (cid:0) ∇ e uh ( t ) , ∇ P nτ e uh ( t ) (cid:1) d t (4.76)= (cid:90) I n (cid:20)(cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1) + (cid:90) tt n − ∂ s (cid:0) ∇ e uh ( s ) , ∇ P nτ e uh ( s ) (cid:1) d s (cid:21) d t = τ (cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1) + (cid:90) I n ∂ s (cid:0) ∇ e uh ( s ) , ∇ P nτ e uh ( s ) (cid:1) ( t n − s )d s = τ (cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1) + Re (cid:90) I n (cid:16) ∂ t ∇ e uh ( t ) , ∇ P nτ e uh ( t )( t n − t ) (cid:17) d t + Re (cid:90) I n (cid:16) ∇ e uh ( t ) , ( t n − t ) ∂ t ∇ P nτ e uh ( t ) (cid:17) d t = τ (cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1)
19 Re (cid:90) I n (cid:16) ∂ t ∇ e uh ( t ) , P nτ (cid:2) ∇ P nτ e uh ( t )( t n − t ) (cid:3)(cid:17) d t + Re (cid:90) I n (cid:16) ∇ P nτ e uh ( t ) , ( t n − t ) ∂ t ∇ P nτ e uh ( t ) (cid:17) d t =: τ (cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1) + J u + J u . In the second to last equality, we have used the identity (cid:90) I n (cid:16) ∇ e uh ( t ) , ( t n − t ) ∂ t ∇ P nτ e uh ( t ) (cid:17) d t = (cid:90) I n (cid:16) ∇ P nτ e uh ( t ) , ( t n − t ) ∂ t ∇ P nτ e uh ( t ) (cid:17) d t, which holds because ( t n − t ) ∂ t ∇ P nτ e uh ( t ) is a polynomial of degree k − J u = (cid:90) I n
12 dd t (cid:107)∇ P nτ e uh ( t ) (cid:107) ( t n − t )d t (4.77)= − (cid:107)∇ P nτ e uh ( t n − ) (cid:107) τ + (cid:90) I n (cid:107)∇ P nτ e uh ( t ) (cid:107) d t. Then substituting this into (4.76) yields (cid:90) I n (cid:107)∇ P nτ e uh ( t ) (cid:107) d t ≤ τ (cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1) + 2 J u . (4.78)Setting v h = ( − ∆ h ) P nτ (cid:2) P nτ e uh ( t )( t n − t ) (cid:3) in (4.73) and taking the imaginary part yield J u = Im (cid:90) I n (cid:16) i ∂ t e uh ( t ) , ( − ∆ h ) P nτ (cid:2) P nτ e uh ( t )( t n − t ) (cid:3)(cid:17) d t (4.79)= − Im (cid:90) I n (cid:0) ∇ P nτ e uh ( t ) , ∇ ( − ∆ h ) P nτ (cid:2) P nτ e uh ( t )( t n − t ) (cid:3)(cid:1) d t + θτ k (cid:88) j =1 w j Im (cid:16) e rnj g ( u φnj ) u φnj , ( − ∆ h ) P nτ (cid:2) P nτ e unj ( t n − t nj ) (cid:3)(cid:17) + τ k (cid:88) j =1 w j Im (cid:16) r ∗ nj (cid:0) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:1) , ( − ∆ h ) P nτ (cid:2) P nτ e unj ( t n − t nj ) (cid:3)(cid:17) − (cid:90) I n Im (cid:0) P nτ d nu , ( − ∆ h ) P nτ [ P nτ e u ( t )( t n − t )] (cid:1) d t =: (cid:88) m =1 J mu , where J u = − Im (cid:90) I n (cid:0) ∆ h P nτ e uh ( t ) , ∆ h P nτ (cid:2) P nτ e uh ( t )( t n − t ) (cid:3)(cid:1) d t (4.80)= − Im (cid:90) I n (cid:107) ∆ h P nτ e uh ( t ) (cid:107) ( t n − t ) d t = 0 . From identity (3.17) and inequality (3.22), we have J u = τ θ k (cid:88) j =1 w j Im (cid:0) ∇ P h [ e rnj g ( u φnj ) u φnj ] , P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:1) τ k (cid:88) j =1 w j | e rnj |(cid:107)∇ P h [ g ( u φnj ) u φnj ] (cid:107)(cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107)≤ τ k (cid:88) j =1 w j | e rnj |(cid:107)∇ [ g ( u φnj ) u φnj ] (cid:107)(cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) (here (2.9b) is used) ≤ Cτ k (cid:88) j =1 w j | e rnj |(cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) (here (4.75) is used) ≤ k (cid:88) j =1 τ τ w j (cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) + Cτ k (cid:88) j =1 w j | e rnj | ≤ τ (cid:90) I n (cid:107) P nτ [ ∇ P nτ e uh ( t )( t n − t )] (cid:107) d t + 2 Cτ (cid:90) I n | e rh | d t ≤ (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t + 2 Cτ (cid:90) I n | e rh | d t,J u = τ k (cid:88) j =1 w j Im (cid:0) r ∗ nj ∇ P h (cid:2) ( g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj ) (cid:3) , P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:1) ≤ τ k (cid:88) j =1 w j | r ∗ nj |(cid:107)∇ P h (cid:2) ( g ( P nτ u φnj ) P nτ u φnj − g ( P nτ u ∗ nj ) P nτ u ∗ nj ) (cid:3) (cid:107)× (cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) (here u φnj = P nτ u φnj is used) ≤ Cτ k (cid:88) j =1 w j (cid:107)∇ P nτ e unj (cid:107)(cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) (here (4.75) is used) ≤ Cτ (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t + Cτ − (cid:90) I n (cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) d t ≤ Cτ (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t,J u = − (cid:90) I n Im (cid:0) ∇ P h P nτ d nu , P nτ [ ∇ P nτ e u ( t )( t n − t )] (cid:1) d t ≤ Cτ (cid:90) I n (cid:107)∇ P h P nτ d nu (cid:107) d t + 14 τ (cid:90) I n (cid:107) P nτ [ ∇ P nτ e unj ( t n − t nj )] (cid:107) d t ≤ Cτ max t ∈ I n (cid:107) d nu (cid:107) H + 14 (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t. Substituting (4.79) and the estimates of J mu , m = 1 , , ,
4, into (4.78), for sufficiently small stepsize τ we obtain (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t ≤ Cτ (cid:0) ∇ e uh ( t n − ) , ∇ P nτ e uh ( t n − ) (cid:1) + Cτ (cid:90) I n | e rh | d t + Cτ max t ∈ I n (cid:107) d nu (cid:107) H ≤ Cτ (cid:107)∇ e uh ( t n − ) (cid:107)(cid:107)∇ P nτ e uh ( t n − ) (cid:107) + Cτ (cid:90) I n | e rh | d t + Cτ max t ∈ I n (cid:107) d nu (cid:107) H Cτ (cid:107)∇ e uh ( t n − ) (cid:107) (cid:18) τ (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t (cid:19) + Cτ (cid:90) I n | e rh | d t + Cτ max t ∈ I n (cid:107) d nu (cid:107) H ≤ Cτ (cid:107)∇ e uh ( t n − ) (cid:107) + 12 (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t + Cτ (cid:90) I n | e rh | d t + Cτ max t ∈ I n (cid:107) d nu (cid:107) H , which then implies (cid:90) I n (cid:107)∇ P nτ e uh (cid:107) d t ≤ Cτ (cid:107)∇ e uh ( t n − ) (cid:107) + Cτ (cid:90) I n | e rh | d t + Cτ max t ∈ I n (cid:107) d nu (cid:107) H . By using inequality (3.21), we can remove the operator P nτ in the above inequality (meanwhilereplace the H seminorm by the full norm), i.e., (cid:90) I n (cid:107) e uh (cid:107) H d t ≤ Cτ (cid:107) e uh ( t n − ) (cid:107) H + Cτ (cid:90) I n | e rh | d t + Cτ max t ∈ I n (cid:107) d nu (cid:107) H . (4.81) Step 2: Estimation of (cid:107) e rh (cid:107) L ( I n ) . To estimate the second term on the right-hand side of(4.81), we proceed similarly as (4.76), that is, (cid:90) I n | P nτ e rh | d t = (cid:90) I n e rh ( t ) P nτ e rh ( t ) d t (4.82)= (cid:90) I n (cid:20) e rh ( t n − ) P nτ e rh ( t n − ) + (cid:90) tt n − ∂ s (cid:0) e rh ( s ) P nτ e rh ( s ) (cid:1) d s (cid:21) d t = τ e rh ( t n − ) P nτ e rh ( t n − ) + (cid:90) I n ∂ t e rh ( t ) P nτ (cid:2) P nτ e rh ( t )( t n − t ) (cid:3) d t + (cid:90) I n e rh ( t )( t n − t ) ∂ t P nτ e rh ( t )d t = τ e rh ( t n − ) P nτ e rh ( t n − ) + (cid:90) I n ∂ t e rh ( t ) P nτ (cid:2) P nτ e rh ( t )( t n − t ) (cid:3) d t + (cid:90) I n P nτ e rh ( t )( t n − t ) ∂ t P nτ e rh ( t )d t =: τ e rh ( t n − ) P nτ e rh ( t n − ) + J r + J r , where we have changed the order of integration in the third equality, and used the the followingidentity in the second to last equality: (cid:90) I n e rh ( t )( t n − t ) ∂ t P nτ e rh ( t )d t = (cid:90) I n P nτ e rh ( t )( t n − t ) ∂ t P nτ e rh ( t )d t, which holds because ( t n − t ) ∂ t P nτ e rh ( t ) is a polynomial of degree k − J r = (cid:90) I n
12 dd t | P nτ e rh ( t ) | ( t n − t )d t (4.83)= − | P nτ e rh ( t n − ) | τ + (cid:90) I n | P nτ e rh ( t ) | d t ≤ (cid:90) I n | P nτ e rh ( t ) | d t. Then substituting (4.83) into (4.82) yields (cid:90) I n | P nτ e rh | d t ≤ τ (cid:0) e rh ( t n − ) , P nτ e rh ( t n − ) (cid:1) + 2 J r Cτ | e rh ( t n − ) |(cid:107) P nτ e rh (cid:107) L ∞ ( I n ) + 2 J r ≤ Cτ | e rh ( t n − ) | (cid:18) τ (cid:90) I n | P nτ e rh | d t (cid:19) + 2 J r ≤ Cτ | e rh ( t n − ) | + 12 (cid:90) I n | P nτ e rh | d t + 2 J r , which then implies (cid:90) I n | P nτ e rh | d t ≤ Cτ | e rh ( t n − ) | + 4 J r , (4.84)In order to estimate J r , we choose q h = P nτ (cid:2) P nτ e rh ( t )( t n − t ) (cid:3) in (4.74), which yields thefollowing identity: J r = (cid:90) I n ∂ t e rh ( t ) P nτ (cid:2) P nτ e rh ( t )( t n − t ) (cid:3) d t (4.85)= − τ k (cid:88) j =1 w j Re (cid:16) P nτ [ P nτ e rnj ( t n − t nj )] (cid:0) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:1) , ∂ t u ∗ h ( t nj ) (cid:17) + τ θ k (cid:88) j =1 w j Re (cid:16) P nτ [ P nτ e rnj ( t n − t nj )] g ( u φnj ) u φnj , ∂ t e uh ( t nj ) (cid:17) − (cid:90) I n P nτ d nr P nτ [ P nτ e r ( t )( t n − t )]d t =: (cid:88) m =1 J mr . By using (3.22), we have J r ≤ τ k (cid:88) j =1 w j | P nτ [ P nτ e rnj ( t n − t nj )] |(cid:107) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:107)(cid:107) ∂ t u ∗ h ( t nj ) (cid:107) (4.86) ≤ Cτ k (cid:88) j =1 w j | P nτ [ P nτ e rnj ( t n − t nj )] |(cid:107) e unj (cid:107)≤ τ τ k (cid:88) j =1 w j | P nτ [ P nτ e rnj ( t n − t nj )] | + Cτ τ k (cid:88) j =1 w j (cid:107) e unj (cid:107) ≤ τ (cid:90) I n | P nτ [ P nτ e rh ( t )( t n − t )] | d t + Cτ (cid:90) I n (cid:107) P nτ e uh (cid:107) d t ≤ (cid:90) I n | P nτ e rh | d t + Cτ (cid:90) I n (cid:107) e uh (cid:107) d t,J r ≤ τ k (cid:88) j =1 w j | P nτ [ P nτ e rnj ( t n − t nj )] |(cid:107) g ( u φnj ) u φnj (cid:107) L ∞ (cid:107) ∂ t e uh ( t nj ) (cid:107) (4.87) ≤ τ τ k (cid:88) j =1 w j | P nτ [ P nτ e rnj ( t n − t nj )] | + Cτ τ k (cid:88) j =1 w j (cid:107) ∂ t e uh ( t nj ) (cid:107) τ (cid:90) I n | P nτ [ P nτ e rh ( t )( t n − t )] | d t + Cτ (cid:90) I n (cid:107) ∂ t e uh (cid:107) d t ≤ (cid:90) I n | P nτ e rh | d t + C (cid:90) I n (cid:107) e uh (cid:107) d t,J r = − (cid:90) I n P nτ d nr P nτ [ P nτ e r ( t )( t n − t )]d t (4.88) ≤ Cτ (cid:90) I n | P nτ d nr | d t + 18 τ (cid:90) I n | P nτ [ P nτ e rh ( t )( t n − t )] | d t ≤ Cτ max t ∈ I n | P nτ d nr | + 18 (cid:90) I n | P nτ e rh | d t. Substituting (4.86)–(4.88) into (4.84)–(4.85), and using (3.21), we get (cid:90) I n | e rh | d t ≤ Cτ | e rh ( t n − ) | + C (cid:90) I n (cid:107) e uh (cid:107) d t + Cτ max t ∈ I n | P nτ d nr | . (4.89)Then, combining (4.81) and (4.89), we obtain (for sufficiently small τ ) (cid:90) I n (cid:107) e uh (cid:107) H d t ≤ Cτ (cid:2) (cid:107) e uh ( t n − ) (cid:107) H + | e rh ( t n − ) | (4.90)+ τ max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:3) , (cid:90) I n | e rh | d t ≤ Cτ (cid:2) (cid:107) e uh ( t n − ) (cid:107) H + | e rh ( t n − ) | + τ max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:3) . (4.91) Step 3: Estimation of (cid:107)∇ e uh ( t n ) (cid:107) and | e rh ( t n ) | . Setting v h = ∂ t e uh in (4.73) and taking the realpart, we get12 (cid:107)∇ e uh ( t n ) (cid:107) − (cid:107)∇ e uh ( t n − ) (cid:107) = θτ k (cid:88) j =1 w j Re (cid:0) e rnj g ( u φnj ) u φnj , ∂ t e uh ( t nj ) (cid:1) (4.92)+ τ k (cid:88) j =1 w j Re (cid:16) r ∗ nj (cid:2) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:3) , ∂ t e uh ( t nj ) (cid:17) − (cid:90) I n Re ( d nu , ∂ t e uh )d t ≤ Cτ k (cid:88) j =1 w j (cid:0) | e rnj | + (cid:107)∇ e unj (cid:107) (cid:1) (cid:107) ∂ t e uh ( t nj ) (cid:107) H − + C (cid:90) I n (cid:107) d nu (cid:107) H (cid:107) ∂ t e uh (cid:107) H − d t ≤ C (cid:90) I n ( (cid:107) e uh (cid:107) H + | e rnj | + (cid:107) d nu (cid:107) H )d t + 12 (cid:90) I n (cid:107) ∂ t e uh (cid:107) H − d t. In order to estimate the last term (cid:82) I n (cid:107) ∂ t e uh (cid:107) H − d t above, we consider (4.73), from which we canderive the following estimate for any test function v ∈ L ( I n ; H ): (cid:12)(cid:12)(cid:12)(cid:90) I n ( ∂ t e uh , v )d t (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) i (cid:90) I n (cid:0) ∇ e uh , ∇ P h P nτ v (cid:1) d t − i θτ k (cid:88) j =1 w j (cid:16) e rnj g ( u φnj ) u φnj , P h P nτ v ( t nj ) (cid:17) − i τ k (cid:88) j =1 w j (cid:16) r ∗ nj (cid:2) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:3) , P h P nτ v ( t nj ) (cid:17)
24 i (cid:90) I n ( d nu , P h P nτ v )d t (cid:12)(cid:12)(cid:12) ≤ C (cid:0) (cid:107) e uh (cid:107) L ( I n ; H ) + | e rh | L ( I n ) + (cid:107) d nu (cid:107) L ( I n ; H − ) (cid:1) (cid:107) P h P nτ v (cid:107) L ( I n ; H ) ≤ C (cid:0) (cid:107) e uh (cid:107) L ( I n ; H ) + | e rh | L ( I n ) + (cid:107) d nu (cid:107) L ( I n ; H − ) (cid:1) (cid:107) v (cid:107) L ( I n ; H ) . By the duality between L ( I n ; H − ) and L ( I n ; H ), we obtain (cid:90) I n (cid:107) ∂ t e uh (cid:107) H − d t ≤ C (cid:90) I n ( (cid:107) e uh (cid:107) H + | e rh | + (cid:107) d nu (cid:107) H − )d t. (4.93)Then, summing up (4.92) and (4.93), we have (cid:107)∇ e uh ( t n ) (cid:107) − (cid:107)∇ e uh ( t n − ) (cid:107) + (cid:90) I n (cid:107) ∂ t e uh (cid:107) H − d t (4.94) ≤ C (cid:90) I n (cid:0) (cid:107) e uh (cid:107) H + | e rh | + (cid:107) d nu (cid:107) H (cid:1) d t. Setting q h = 2 e rh in (4.74) yields | e rh ( t n ) | − | e rh ( t n − ) | = τ k (cid:88) j =1 w j Re (cid:0) e rnj (cid:0) g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj (cid:1) , ∂ t u ∗ h ( t nj ) (cid:1) + θτ k (cid:88) j =1 w j Re (cid:0) e rnj g ( u φnj ) u φnj , ∂ t e uh ( t nj ) (cid:1) − (cid:90) I n P nτ d nr e r d t (4.95) ≤ Cτ k (cid:88) j =1 w j | e rnj |(cid:107) e unj (cid:107) + Cτ k (cid:88) j =1 w j | e rnj |(cid:107) ∂ t e unj (cid:107) H − + C (cid:90) I n | P nτ d nr || e r | d t ≤ C (cid:90) I n (cid:0) (cid:107) e uh (cid:107) + | e rh | + (cid:107) ∂ t e uh (cid:107) H − (cid:1) d t + C (cid:90) I n | P nτ d nr | d t ≤ C (cid:90) I n (cid:0) (cid:107) e uh (cid:107) H + | e rh | (cid:1) d t + C (cid:90) I n ( (cid:107) d nu (cid:107) H + | P nτ d nr | )d t, where we have used (4.93) to obtain the last inequality. Step 4: Completion of the proof.
Summing up (4.94) and (4.95) yields (cid:107)∇ e uh ( t n ) (cid:107) + | e rh ( t n ) | − (cid:107)∇ e uh ( t n − ) (cid:107) − | e rh ( t n − ) | + (cid:90) I n (cid:107) ∂ t e uh (cid:107) H − d t (4.96) ≤ C (cid:90) I n (cid:0) (cid:107) e uh (cid:107) H + | e rh | (cid:1) d t + C (cid:90) I n ( (cid:107) d nu (cid:107) H + | P nτ d nr | )d t. Then, substituting (4.90)–(4.91) into the inequality above, we obtain (cid:0) (cid:107)∇ e uh ( t n ) (cid:107) + | e rh ( t n ) | (cid:1) − (cid:0) (cid:107)∇ e uh ( t n − ) (cid:107) + | e rh ( t n − ) | (cid:1) + (cid:90) I n (cid:107) ∂ t e uh (cid:107) H − d t (4.97)25 Cτ (cid:0) (cid:107)∇ e uh ( t n − ) (cid:107) + | e rh ( t n − ) | (cid:1) + C (cid:90) I n ( (cid:107) d nu (cid:107) H + | P nτ d nr | )d t. It follows from Gronwall’s inequality thatmax ≤ n ≤ N (cid:0) (cid:107)∇ e uh ( t n ) (cid:107) + | e rh ( t n ) | (cid:1) + C (cid:90) T (cid:107) ∂ t e uh (cid:107) H − d t (4.98) ≤ C (cid:0) (cid:107) e uh (0) (cid:107) H + | e rh (0) | (cid:1) + C N (cid:88) n =1 (cid:90) I n ( (cid:107) d nu (cid:107) H + | P nτ d nr | )d t. Then, substituting this inequality into (4.90)–(4.91) and using temporal inverse inequlaity, weobtain max t ∈ [0 ,T ] (cid:0) (cid:107) e uh ( t ) (cid:107) H + | e rh ( t ) | (cid:1) ≤ C (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) . (4.99)Hence, (4.70) holds.When τ and h are sufficiently small, inequality (4.99) implies thatmax t ∈ [0 ,T ] (cid:107) e uh ( t ) (cid:107) H ≤
12 and max t ∈ [0 ,T ] | e rh ( t ) | ≤ . (4.100)On the one hand, by the inverse inequality, we havemax t ∈ [0 ,T ] (cid:107) e uh ( t ) (cid:107) L ∞ ≤ C(cid:96) h max t ∈ [0 ,T ] (cid:107) e uh ( t ) (cid:107) H ≤ C(cid:96) h (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) , (4.101)where (cid:96) h = d = 1 , ln(2 + 1 /h ) if d = 2 ,h − if d = 3 . On the other hand, by choosing a test function v in (4.73) satisfying the properties v ( t nj ) = 1and v ( t ni ) = 0 for i (cid:54) = j , and using property (3.18), we obtain (cid:107) ∆ h e unj (cid:107) = (cid:13)(cid:13)(cid:13) i ∂ t e unj − θP h (cid:2) e rnj g ( u φnj ) u φnj (cid:3) + P h d unj (4.102) − P h (cid:2) r ∗ nj ( g ( u φnj ) u φnj − g ( u ∗ nj ) u ∗ nj ) (cid:3)(cid:13)(cid:13)(cid:13) ≤ Cτ − (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) , where we have used (4.98)–(4.99) and an inverse inequality in time in estimating ∂ t e unj . By thediscrete Sobolev embedding inequality, for 1 ≤ d ≤ (cid:107) e unj (cid:107) L ∞ ≤ C (cid:107) e unj (cid:107) H (cid:107) ∆ h e unj (cid:107) ≤ Cτ − max ≤ n ≤ N (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) , (4.103)26here we have used (4.99) and (4.102) in the last inequality. Then, combining (4.101) and(4.103) yieldsmax ≤ n ≤ N max ≤ j ≤ k (cid:107) e u ( t nj ) (cid:107) L ∞ ≤ C min( (cid:96) h , τ − ) (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) ≤ C (cid:0) h p − + τ k + (cid:1) , where we have used the consistency estimate from Theorem 3.3. When τ and h are sufficientlysmall, the inequality above impliesmax ≤ n ≤ N max ≤ j ≤ k (cid:107) e u ( t nj ) (cid:107) L ∞ ≤ . (4.104)This together with (4.100) gives (4.71).Furthermore, since φ h = θe uh and ϕ h = θe rh , it follows thatmax ≤ n ≤ N max ≤ j ≤ k (cid:107) φ h ( t nj ) (cid:107) L ∞ ∩ H ≤
12 and max ≤ n ≤ N max ≤ j ≤ k | ϕ h ( t nj ) | ≤ , which imply ρ [ φ h ] = ρ [ ϕ h ] = 1 in view of the definition in (4.62). This proves (4.72). (cid:3) We now are ready to state and prove existence, uniqueness and convergence of numericalsolutions, which comprise of the second main theorem of this paper.
Theorem 4.2
Let 1 ≤ d ≤ τ and h such that when τ ≤ τ and h ≤ h , the numerical method (2.12) has a unique solution ( u h , r h ) ∈ X ∗ τ,h × Y ∗ τ,h . Moreover,this solution satisfies the following error estimate:max t ∈ [0 ,T ] (cid:16) (cid:107) u h ( t ) − u ∗ h ( t ) (cid:107) H + | r h ( t ) − r ∗ h ( t ) | (cid:17) ≤ C ( h p + τ k +1 ) . (4.105) Proof. Step 1: Existence.
By the definition of B , if ( φ h , ϕ h ) ∈ B and ( e uh , e rh ) = M ( φ h , ϕ h )then φ h = θe uh and ϕ h = θe rh . Thus (4.70) implies (cid:107) ( φ h , ϕ h ) (cid:107) X τ,h × Y τ,h = (cid:107) φ h (cid:107) L ∞ (0 ,T ; H ) + (cid:107) ϕ h (cid:107) L ∞ (0 ,T ) ≤ C, (4.106)which together with Schaefer’s fixed point theorem imply the existence of a fixed point ( φ h , ϕ h )for the mapping M (corresponding to θ = 1), with( e uh , e rh ) = ( φ h , ϕ h ) , u φ = u ∗ h + φ h and r φ = r ∗ h + ϕ h , satisfying (4.66)–(4.67), where we have used (4.72) in the expression (4.65). Consequently,( e uh , e rh ) is a solution of (4.58) with ( u h , r h ) = ( u φh , r ϕh ) = ( u ∗ h + e uh , r ∗ h + e rh ). Hence, in view ofthe discussions in Remark 4.1, ( u h , r h ) is a solution of the numerical scheme (2.12), and (4.71)implies ( u h , r h ) is in the set X ∗ τ,h × Y ∗ τ,h defined in (4.60)–(4.61). This proves existence of anumerical solution in X ∗ τ,h × Y ∗ τ,h . Step 2: Uniqueness.
Suppose that ( u h , r h ) and ( (cid:101) u h , (cid:101) r h ) in X ∗ τ,h × Y ∗ τ,h are two pairs of numer-ical solutions, and set e uh = u h − (cid:101) u h and e rh = r h − (cid:101) r h (abusing the notation). Subtracting the27orresponding equations satisfied by ( u h , r h ) and ( (cid:101) u h , (cid:101) r h ) shows that ( e uh , e rh ) satisfies equations(4.58) with d nu = d nr = 0. In the meantime, the definition in (4.60)–(4.61) implies (cid:107) e uh ( t nj ) (cid:107) L ∞ ∩ H ≤ | e rh ( t nj ) | ≤ . (4.107)Accordingly, ( e uh , e rh ) is a fixed point of the mapping M (corresponding to θ = 1 in B ) in thecase e uh (0) = e rh (0) = 0 and d nu = d nr = 0. Hence, an application of (4.70) yields (cid:107) e uh (cid:107) L ∞ (0 ,T ; H ) + (cid:107) e rh (cid:107) L ∞ (0 ,T ) ≤ C (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) = 0 . Thus, ( u h , r h ) = ( (cid:101) u h , (cid:101) r h ) and the uniqueness of the numerical solution is proved. Step 3: Error estimate.
Since the error functions e uh = u h − u ∗ h and e rh = r h − r ∗ h satisfy(4.58) and (4.107), it follows that ( e uh , e rh ) is a fixed point of the mapping M (corresponding to θ = 1 in B ). Hence, an application of (4.70) yields (cid:107) e uh (cid:107) L ∞ (0 ,T ; H ) + (cid:107) e rh (cid:107) L ∞ (0 ,T ) ≤ C (cid:104) (cid:107) e uh (0) (cid:107) H + | e rh (0) | + max ≤ n ≤ N max t ∈ I n (cid:0) (cid:107) d nu (cid:107) H + | P nτ d nr | (cid:1)(cid:105) . Substituting the consistency error estimates from Theorem 3.3 into the above inequality yieldsthe desired estimate (4.105). The proof is complete. (cid:3)
Remark 4.2
For the periodic and Neumann boundary conditions, the mass and energy con-servations in Theorem 3.1 and the error estimate in Theorem 4.2 can be proved similarly.
In this section, we present some one-dimensional numerical tests to validate the theoreticalresults proved in Theorems 3.1 and 4.2 about the mass and energy conservations, and the con-vergence rates of the proposed method. All the computations are performed using the softwarepackage FEniCS ( https://fenicsproject.org ).We consider the cubic nonlinear Schr¨odinger equationi ∂ t u − ∂ xx u − | u | u = 0 in ( − L, L ) × (0 , T ] ,u | t =0 = u in ( − L, L ) , with L = 20 , (5.108)subject to the periodic boundary condition. We choose u = sech( x ) exp(2i x ) so that the exactsolution is given by u ( x, t ) = sech( x + 4 t ) exp(i(2 x + 3 t )) . (5.109)This example contains a soliton wave and is often used as a benchmark for meansuring theeffectiveness of numerical methods for the NLS equation; see [35, 39, 26].28 .1 Convergence rates We solve problem (5.108) by the proposed method (2.12) and compare the numerical solutionswith the exact solution (5.109). Newton’s iteration is used to solve the nonlinear system. Theiteration is stopped when the error is below 10 − .The time discretization errors are presented in Table 1, where we have used finite elements ofdegree 3 with a sufficiently spatial mesh h = 2 L/ O ( τ k +1 ), which is consistent with the result proved in Theorem4.2.The spatial discretization errors are presented in Table 2, where we have chosen k = 3 witha sufficiently small time stepsize τ = 1 / O ( h p ) in the H norm. This is also consistent with the result proved in Theorem 4.2. We denote the mass and SAV energy of a numerical solution by M h ( t ) = (cid:90) Ω | u h ( t ) | d x and E h ( t ) = 12 (cid:90) Ω |∇ u h ( t ) | d x − r h ( t ) , (5.110)respectively. The evolution of mass and SAV energy of the numerical solutions is presented inFigure 1 with τ = 0 . h = 0 .
2. It is shown thatmass = 2 + O (10 − ) and SAV energy = − . O (10 − ) , which are much smaller than the error of the numerical solutions, as shown in Figure 2. Thisshows the effectiveness of the proposed method in preserving mass and energy (independent ofthe error of numerical solutions). The graph of | u ( x, t ) | is a soliton propagating towards left. Its shape remains unchanged forall t ≥ References [1] G. D. Akrivis, V. A. Dougalis and O. A. Karakashian,
On fully discrete Galerkin methods ofsecond-order temporal accuracy for the NLS equation , Numer. Math., 59(1991), pp. 31–53.[2] G. Akrivis,
Finite difference discretization of the cubic Schr¨odinger equation , IMA J. Numer.Anal. 13(1993), pp. 115–124.[3] G. Akrivis, B. Li and D. Li,
Energy-decaying extrapolated RK-SAV methods for the Allen–Cahn and Cahn–Hilliard equations , SIAM J. Sci. Comput., 41(2019), pp. A3703–A3727.29able 1: Time discretization errors of the proposed method, with h = L and T = 1. k τ p = 3 (cid:107) u ( x, t ) − u h ( x, t ) (cid:107) L ∞ (0 ,T ; H ) order2 1/60 3.7964E–05 –1/70 2.3429E–05 3.13121/80 1.5460E–05 3.11321/90 1.0733E–05 3.09851/100 7.7542E–06 3.08533 1/20 3.4019E–05 –1/20 1.3821E–05 4.03641/30 6.6322E–06 4.02751/35 3.5689E–06 4.02001/40 2.0886E–06 4.01234 1/8 1.2291E–04 –1/12 1.5120E–05 5.16811/14 6.8492E–06 5.13691/16 3.4634E–06 5.10671/20 1.1555E–06 4.9192 Table 2: Time discretization errors of the proposed method, with τ = and T = 1. p M k = 3 (cid:107) u ( x, t ) − u h ( x, t ) (cid:107) L ∞ (0 ,T ; H ) order1 1400 5.8670E–02 –1600 5.1134E–02 1.02951800 4.5330E–02 1.02292000 4.0719E–02 1.01832200 3.6964E–02 1.01492 240 1.9306E–02 –260 1.6438E–02 2.0094280 1.4167E–02 2.0062300 1.2338E–02 2.0041320 1.0842E–02 2.00273 90 1.6147E–02 –100 1.1661E–02 3.0894110 8.7112E–03 3.0599120 6.6844E–03 3.0436130 5.2435E–03 3.0334 M h ( t ) − M h (0) and SAV energy E h ( t ) − E h (0), with p = 3 and τ = h = 0 . p = 3 and τ = h = 0 . t ∈ [0 , | u ( · , t ) | .31 a) (b)(c) (d) Figure 4: Soliton propagation when t ∈ [0 , p = 1, M = 1200 and∆ t = 0 . (a) (b)(c) (d) Figure 5: Soliton propagation when t ∈ [0 , p = 1, M = 1200 and∆ t = 0 .
05. 324] X. Antoine, W. Bao and C. Besse,
Computational methods for the dynamics of the nonlinearSchr¨odinger/Gross–Pitaevskii equations , Comput. Phys. Commun., 184(2013), pp. 2621–2633.[5] W. Bao and Y. Cai,
Optimal error estimates of finite difference methods for the Gross-Pitaevskii equation with angular momentum rotation , Math. Comp., 82(2013), pp. 99–128.[6] W. Bao, Q. Tang and Z. Xu,
Numerical methods and comparison for computing dark andbright solitons in the NLS equation , J. Comput. Phys., 235(2013), pp. 423–445.[7] C. Besse,
A relaxation scheme for the nonlinear Schr¨odinger equation , SIAM J. Numer.Anal., 42(2004), pp. 934–952[8] J. Bourgain,
Global Solutions of Nonlinear Schr¨odinger Equations , vol. 46, American Math-ematical Society, 1999.[9] S. C. Brenner and L.R Scott,
The Mathematical Theory of FEMs , Third edition. Texts inApplied Mathematics, Vol. 15, Springer, New York, 2008.[10] C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang,
Spectral Methods: Fundamentalsin Single Domains , Springer, Berlin, 2007.[11] M. Delfour, M. Fortin and G. Payre,
Finite-difference solutions of a non-linear Schr¨odingerequation , J. Comput. Phys., 44(1981), pp.277–288.[12] L. C. Evans,
Partial Differential Equations , second edition, Graduate Studies in Mathe-matics 19, AMS, Providence, RI, 2010.[13] X. Feng, H. Liu and S. Ma,
Mass- and energy-conserved numerical schemes for nonlinearSchr¨odinger equations , Commun. Comput. Phys., 26(2019), pp. 1365–1396.[14] Z. Gao and S. Xie,
Fourth-order alternating direction implicit compact finite differenceschemes for two-dimensional Schr¨odinger equations , Appl. Numer. Math., 61(2011), pp.593–614.[15] D. Gilbarg and N. Trudinger,
Elliptic Partial Differential Equations of Second Order ,Springer, New York, 2001.[16] G. H. Golub and J. H. Welsch,
Calculation of Gauss quadrature rules , Math. Comput.,23(1969), pp. 221–230.[17] Y. Gong, Q. Wang and Y. Wang,
A conservative Fourier pseudo-spectral method for thenonlinear Schr¨odinger equation , J. Comput. Phys., 328(2017), pp. 354–370.[18] Y. Gong and J. Zhao,
Energy-stable RungeKutta schemes for gradient flow models usingthe energy quadratization approach , Appl. Math. Letters, 94(2019), pp. 224–231.[19] Y. Gong, J. Zhao and Q. Wang,
Arbitrarily high-order unconditionally energy stable schemesfor thermodynamically consistent gradient flow models , SIAM J. Sci. Comput. 42(2020), pp.B135–156. 3320] P. Henning and D. Peterseim,
CrankNicolson Galerkin approximations to nonlinearSchr¨odinger equations with rough potentials , Math. Models Meth. Appl. Sci., 27(2017),pp. 2147–2184.[21] J. Hong, Y. Liu, H. Munthe-Kaas and A. Zanna,
Globally conservative properties and errorestimation of a multi-symplectic scheme for Schr¨odinger equations with variable coefficients ,Appl. Numer.Math., 56(2006), pp. 814–843.[22] O. Karakashian and C. Makridakis,
A space-time finite element method for the nonlinearSchr¨odinger equation: the discontinuous Galerkin method , Math. Comp., 67(1998), pp.479–499.[23] O. Karakashian and C. Makridakis,
A space-time finite element method for the nonlinearSchrdinger equation: the continuous Galerkin method , SIAM J. Numer. Anal., 36(1999),pp. 1779–1807.[24] D. A. Kopriva,
Implementing spectral methods for partial differential equations: Algorithmsfor scientists and engineers , Springer Science & Business Media, 2009.[25] H. Liu, Y. Huang, W. Lu and N. Yi,
On accuracy of the mass-preserving DG method tomulti-dimensional Schr¨odinger equations , IMA J. Numer. Anal., 39(2019), pp. 760–791[26] W. Lu, and Y. Huang and H. Liu,
Mass preserving discontinuous Galerkin methods forSchr¨odinger equations , J. Comput. Phys., 282(2015), pp. 210–226.[27] D. E. Pelinovsky, V. V. Afanasjev and Y. S. Kivshar,
Nonlinear theory of oscillating, decay-ing, and collapsing solitons in the generalized nonlinear Schr¨odinger equation , Phys. Rev.E, 53(1996), pp. 1940–1953.[28] J. M. Sanz-Serna,
Methods for the numerical solution of the nonlinear Schr¨odinger equation ,Math. Comp., 43(1984), pp. 21–27[29] J. M. Sanz-Serna and V. S. Manoranjan,
A method for the integration in time of certainpartial differential equations , J. Comput. Phys., 52(1983), pp. 273–289[30] H. W. Sch¨urmann,
Traveling-wave solutions of the cubic-quintic nonlinear Schr¨odingerequation , Phys. Rev. E, 54 (1996), pp. 4312–4320.[31] J. Shen, J. Xu and J. Yang,
A new class of efficient and robust energy stable schemes forgradient flows , SIAM Rev., 61(2019), pp. 474–506.[32] J. Shen, J. Xu and J. Yang,
The scalar auxiliary variable (SAV) approach for gradient flows ,J. Comput. Phys., 353(2018), pp. 407–416.[33] J. Shen, T. Tang and L. L. Wang,
Spectral methods: algorithms, analysis and applications ,Springer Science & Business Media, 41(2011).[34] W. A. Strauss and L. Vazquez,
Numerical solution of a nonlinear Klein–Gordon equation ,J. Comput. Phys., 28(1978), pp.271–278.[35] N. Taghizadeh, M. Mirzazadeh and F. Farahrooz,
Exact solutions of the nonlinearSchr¨odinger equation by the first integral method , J. Math. Anal. Appl., 374(2011), pp.549–553. 3436] T. Tao,
Nonlinear Dispersive Equations: Local and Global Aalysis , American MathematicalSociety, 2006.[37] J. Wang,
A new error analysis of Crank–Nicolson Galerkin FEMs for a generalized nonlinearSchr¨odinger equation , J. Sci. Comput., 60(2014), pp. 390–407.[38] T. Wang, B. Guo and Q. Xu,
Fourth-order compact and energy conservative differenceschemes for the nonlinear Schrdinger equation in two dimensions , J. Comput. Phys.,243(2013), pp. 382–399.[39] Y. Xu and C. W. Shu,
Local discontinuous Galerkin methods for nonlinear Schr¨odingerequations , J. Comput. Phys., 205(2005), pp. 72–97.[40] X. Yang,
Linear, first and second-order, unconditionally energy stable numerical schemesfor the phase field model of homopolymer blends , J. Comput. Phys. 327(2016), pp. 294–316.[41] X. Yang and L. Ju,
Efficient linear schemes with unconditional energy stability for thephase field elastic bending energy model , Comput. Meth. Appl. Mech. Engrg., 315(2017),pp. 691–712.[42] X. Yang and L. Ju,
Linear and unconditionally energy stable schemes for the binary fluid-surfactant phase field model , Comput. Meth. Appl. Mech. Engrg., 318(2017), pp. 1005–1029.[43] X. Yang, J. Zhao, Q. Wang and J. Shen,
Numerical approximations for a three componentsCahn–Hilliard phase-field model based on the invariant energy quadratization method , Math.Models Methods Appl. Sci., 27(2017), pp. 1993–2030.[44] N. J. Zabusky and M. D. Kruskal,
Interaction of ”solitons” in a collisionless plasma andthe recurrence of initial states , Phys. Rev. Lett., 15(1965), pp. 240–243.[45] F. Zhang, V. M. P´erez-Garcia and L. V´azquez,