numerical analysis of a second-order pure lagrange-galerkin method ...

NUMERICAL ANALYSIS OF A SECOND-ORDER PURE LAGRANGE-GALERKIN METHOD FOR CONVECTION-DIFFUSION PROBLEMS. PART II: FULLY DISCRETIZED SCHEME AND NUMERICAL RESULTS∗ † ´ MARTA BEN´ITEZ† AND ALFREDO BERMUDEZ

Abstract. We analyze a second order pure Lagrange-Galerkin method for variable coefficient convection-(possibly degenerate) diffusion equations with mixed Dirichlet-Robin boundary conditions. In a previous paper the proposed second order pure Lagrangian time discretization scheme has been introduced and analyzed for the same problem. Moreover, the l∞ (H 1 ) stability and l∞ (H 1 ) error estimates of order O(∆t2 ) has been obtained. In the present paper l∞ (H 1 ) error estimates of order O(∆t2 ) + O(hk ) are obtained for the fully discretized pure Lagrange-Galerkin method. To prove these results we use some properties obtained in the previous paper. Finally, numerical tests are presented that confirm the theoretical results. Key words. convection-diffusion equation, Lagrangian method, characteristics method, Galerkin discretization, stability, error estimates, second order schemes AMS subject classifications. 65M12, 65M15, 65M25, 65M60

1. Introduction. For convection-diffusion problems with dominant convection, methods of characteristics for time discretization are extensively used (see the review paper [21]). These methods are based on time discretization of the material time derivative. For space discretization, they has been combined with finite differences [19], finite elements ([26], [9], [11], [24], [32], [31], [27]), spectral finite elements ([33], [1]), discontinuous finite elements ([3], [2], [4]), and so on. When combined with finite elements they are also called Lagrange-Galerkin methods. In particular, when the characteristics methods are formulated in Lagrangian coordinates (respectively, Eulerian coordinates) they are called pure Lagrangian methods (respectively, semi-Lagrangian methods). The Lagrange-Galerkin method has been mathematically analyzed and applied to different problems by several authors, primarily in the semi-Lagrangian version. Numerical solution of convection-diffusion partial differential equations by this kind of methods is addressed in ([19], [26], [32], [18], [5], [17], [13]) among others. In the present paper we will consider the combination of the pure Lagrangian method proposed and analyzed in [7] with a spatial discretization by using finite element spaces. There exists an extensive literature studying the classical first order characteristic method combined with finite elements applied to convection-diffusion equations. More precisely, if △t denotes the time step, h the mesh-size and k the degree of the finite elements space, estimates of the form O(hk ) + O(△t) in the l∞ (L2 (Rd ))-norm are shown in [32] (d denotes the dimension of the spatial domain). In [26] error estimates of the form O(hk ) + O(△t) + O(hk+1 /△t) in the l∞ (L2 (Ω))-norm are obtained under the assumption that the normal velocity vanishes on the boundary of Ω. All of these estimates involve constants depending on solution norms. For linear finite elements and for a velocity field vanishing on the boundary, convergence of order ∗ This work was supported by Xunta de Galicia under research project INCITE09 207 047 PR, and by Ministerio de Ciencia e Innovaci´ on (Spain) under research projects Consolider MATHEMATICA CSD2006-00032 and MTM2008-02483. † Dep. de Matem´ atica Aplicada, Universidade de Santiago de Compostela, Campus Vida, 15782 Santiago de Compostela, Spain ([email protected], [email protected]). The first author was supported by Ministerio de Educaci´ on.

1

2

´ M. BEN´ITEZ AND A. BERMUDEZ

O(h2 ) + O(min(h, h2 /△t) + O(△t) in the l∞ (L2 (Ω))-norm is stated in [5], where the constants only depend on the data. In principle, the method of characteristics has been introduced for evolution equations but an adaption to solve convection-diffusion stationary problems has been proposed in [10]. In order to increase the order of time and space approximations, higher order schemes for the discretization of the material derivative and higher order finite element spaces would be used. In [30] a second order characteristics method for solving constant coefficient convection-diffusion equations with Dirichlet boundary conditions is studied. The Crank-Nicholson discretization has been used to approximate the material time derivative. For a divergence-free velocity field vanishing on the boundary and a smooth enough solution, stability and O(△t2 ) + O(hk ) error estimates in the l∞ (L2 (Ω))-norm are stated (see also [12] and [13] for further analysis). In [17], semi-Lagrangian and pure Lagrangian methods are proposed and analyzed for convection-diffusion equations. Error estimates for a Galerkin discretization of a pure Lagrangian formulation and for a discontinuous Galerkin discretization of a semi-Lagrangian formulation are obtained. The estimates are written in terms of the projections constructed in [15] and [16]. In the present paper, fully discretized pure Lagrange-Galerkin schemes are used for a more general problem. Specifically, we consider a (possibly degenerate) variable coefficient diffusive term instead of the simpler Laplacian one, a general mixed DirichletRobin boundary condition and a time dependent domain. Moreover, we analyze a scheme involving approximate characteristic curves. As in [7], the mathematical formalism of continuum mechanics (see for instance [23]) is used to introduce the schemes and to analyze the error. In most cases the exact characteristics curves are not easy to compute analytically, so, as in the first part of this work, our analysis include the case where the characteristics curves are approximated using a second order Runge-Kutta scheme. In [7] a l∞ (L2 ) stability inequality is stated and l∞ (L2 ) error estimates of order O(∆t2 ) are obtained; these estimates are uniform in the hyperbolic limit. Furthermore, stability and error estimates of order O(∆t2 ) are proved in the l∞ (H 1 )-norm. As a logical continuation of [7], fully discretized pure Lagrange-Galerkin scheme with a wide class of finite element spaces is analyzed in the present paper. More precisely, l∞ (L2 ) error estimates of order O(∆t2 )+O(hk ) are obtained; these estimates are bounded in the hyperbolic limit. Moreover, error estimates of order O(∆t2 ) + O(hk ) are proved in the l∞ (H 1 )-norm. Usually, the unconditional stability of characteristics methods is only proved under the assumption that the inner products in the Galerkin formulation are exactly calculated. This is rarely possible so in practice they are calculated using numerical quadrature. In general this adds some terms to the final error estimates and, in some cases, it produces the loss of unconditional stability. There are several papers in the literature analyzing the effect of numerical integration in Lagrange-Galerkin methods (see [24], [32], [28], [22], [34], [13]). In particular, in [24] Fourier analysis is developed for the classical Lagrange-Galerkin method involving piecewise linear finite elements, when it is applied to the one dimensional linear convection equation and combined with several quadrature formulas. Unconditional stability has been shown for the trapezoidal rule and unconditional instability has been proved when the mass matrix is exactly integrated and the term of characteristics is approximated by using the trapezoidal rule (Lemma 2.4 in [24]). In [8] an analogous approach is developed for the classical Lagrange-Galerkin method for piecewise linear finite element applied to the one dimensional linear convection equation. The term of characteristic is decomposed

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

3

into two parts; one of them is exactly integrated and the other one is approximated using the trapezoidal rule (see also [25] for more details). For this scheme conditional stability depending on the CFL number has been shown when the mass matrix is exactly integrated. Moreover, numerical results showing the influence of several quadrature formulas in the stability are presented. In the present paper, quadrature formulas leading to stable schemes are used for the practical implementation of the introduced methods. The paper is organized as follows. In Section 2 the convection-diffusion Cauchy problem is posed in a time dependent bounded domain, a weak formulation of this problem in Lagrangian coordinates is written and some notations and hypotheses are stated. In Section 3, we introduce the finite element spaces considered for spatial discretization, pose the corresponding fully discretized schemes, and state their stability properties. In Section 4, under suitable hypotheses on data and solution, l∞ (L2 ) and l∞ (H 1 ) error estimates of order O(∆t2 ) + O(hk ) for the solution of the fully discretized problem are derived. Finally, in Section 5 numerical examples showing the above theoretical results are presented. 2. Statement of the problem and weak formulation in Lagrangian coordinates. Let Ω be a bounded domain in Rd (d = 2, 3) with Lipschitz boundary Γ divided into two parts: Γ = ΓD ∪ ΓR , with ΓD ∩ ΓR = ∅. Let T be a positive constant and Xe : Ω × [0, T ] −→ Rd be a motion in the sense of Gurtin [23]. In particular, Xe ∈ C3 (Ω × [0, T ]) and for each fixed t ∈ [0, T ], Xe (·, t) is a one-to-one function satisfying (2.1)

det F (p, t) > 0 ∀p ∈ Ω,

being F (·, t) the Jacobian matrix of the deformation Xe (·, t). We call Ωt = Xe (Ω, t), D R R Γt = Xe (Γ, t), ΓD t = Xe (Γ , t) y Γt = Xe (Γ , t), for t ∈ [0, T ]. We assume that Ω0 = Ω. We will adopt the notation given in [7] for the trajectory of the motion (T ), the velocity (v) and the functional spaces involved (see §2 of [7] for more details). Let us introduce the set [ (2.2) O := Ωt . t∈[0,T ]

We denote by L the gradient of v. Let us consider the following initial-boundary value problem. (SP) STRONG PROBLEM. Find a function φ : T −→ R such that (2.3) ρ(x)

∂φ (x, t) + ρ(x)v(x, t) · grad φ(x, t) − div (A(x) grad φ(x, t)) = f (x, t), ∂t

for x ∈ Ωt and t ∈ (0, T ), subject to the boundary conditions (2.4) (2.5)

φ(·, t) = φD (·, t) on ΓD t , R αφ(·, t) + A(·) grad φ(·, t) · n(·, t) = g(·, t) on Γt ,

for t ∈ (0, T ), and the initial condition (2.6)

φ(x, 0) = φ0 (x) in Ω.

In the above equations, A : O −→ Sym denotes the diffusion tensor field, where Sym is the space of symmetric tensors in the d-dimensional space, ρ : O −→ R,

´ M. BEN´ITEZ AND A. BERMUDEZ

4

R f : T −→ R, φ0 : Ω −→ R, φD (·, t) : ΓD t −→ R and g(·, t) : Γt −→ R, t ∈ (0, T ) are given scalar functions and n(·, t) is the outward unit normal vector to Γt . For a Banach function space X and an integer m, spaces C m ([0, T ], X) and H m ((0, T ), X) will be abbreviated as C m (X) and H m (X), respectively, and endowed with norm   12   Z TX m ||ϕ||C m (X) := max max ||ϕ(j) (t)||X .||ϕ||H m (X) :=  ||ϕ(j) (t)||2X dt . t∈[0,T ]

j=0,...,m

0

j=0

In the above definitions, ϕ(j) denotes the j-th derivative of ϕ with respect to time. We define the material description Ψm of a spatial field Ψ by (2.7)

Ψm (p, t) = Ψ(Xe (p, t), t).

Similar definition is used for mappings Ψ defined in a subset of T or of O. Throughout this article, we use the notation em (p, t) := F −1 (p, t)Am (p, t)F −T (p, t) det F (p, t) ∀(p, t) ∈ Ω × [0, T ], A m(p, e t) := |F −T (p, t)m(p)| det F (p, t) ∀(p, t) ∈ Γ × [0, T ],

where m is the outward unit normal vector to Γ. We introduce the number of time steps, N , the time step ∆t = T /N , and the meshpoints, tn = n∆t for n = 0, 1/2, 1, . . . , N . In what follows we use the notation ψ n (y) := ψ(y, tn ) for a function ψ(y, t). In [7] the following Lagrangian formulation of the initial-boundary value problem (SP) has been deduced: (LSP) LAGRANGIAN STRONG PROBLEM. Find a function φm : Ω × [0, T ] −→ R such that h i em (p, t)∇φm (p, t) = fm (p, t) det F (p, t), (2.8) ρm (p, t)φ˙ m (p, t) det F (p, t) − Div A for (p, t) ∈ Ω × (0, T ), subject to the boundary conditions (2.9)

φm (p, t) = φD (Xe (p, t), t) on ΓD × (0, T ),

em (p, t)∇φm (p, t) · m(p) = m(p, αm(p, e t)φm (p, t) + A e t)g(Xe (p, t), t) on ΓR × (0, T ),

(2.10)

and the initial condition (2.11)

φm (p, 0) = φ0 (p) in Ω.

Moreover, the standard weak problem associated with this Lagrangian strong problem has been considered: (LWP) LAGRANGIAN WEAK PROBLEM. Find a function φm : Ω×[0, T ] −→ R such that Z Z ˙ em (p, t)∇φm (p, t) · ∇ψ(p) dp ρm (p, t)φm (p, t)ψ(p) det F (p, t) dp + A Ω Ω Z Z +α m(p, e t)φm (p, t)ψ(p) dAp = fm (p, t)ψ(p) det F (p, t) dp (2.12) ΓR ZΩ m(p, e t)gm (p, t)ψ(p) dAp , + ΓR

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

5

∀ψ ∈ HΓ1D (Ω) and t ∈ (0, T ). In [7] a second order pure Lagrangian time semidiscretized scheme have been proposed and analyzed. Stability and error estimates has been obtained under the following hypothesis on the problem data: Hypothesis 1. There exists a parameter δ > 0, such that the velocity field v is defined in [ δ T δ := (2.13) Ωt × {t}, t∈[0,T ]

being (2.14)

Ωδt :=

[

B(x, δ).

[

Ωt .

x∈Ωt

Moreover v ∈ C1 (T δ ). We recall some notations given in [7] (2.15) (2.16)

Oδ := TΓδR :=

δ

t∈[0,T ]

[

Gt × {t},

[

B(x, δ).

δ

t∈[0,T ]

being (2.17)

Gδt =

x∈ΓR t

Hypothesis 2. Function ρ is defined in Oδ and belongs to W 1,∞ (Oδ ), being Oδ the set introduced in (2.15). Moreover, (2.18)

a.e. x ∈ Oδ .

0 < γ ≤ ρ(x)

Let us denote ρ1,∞ = ||ρ||1,∞,Oδ . Hypothesis 3. The diffusion tensor, A, is defined in Oδ and belongs to W1,∞ (Oδ ). Moreover, A is symmetric and has the following form:   An1 Θ (2.19) A= , Θ Θ with An1 being a positive definite symmetric n1 × n1 tensor (n1 ≥ 1) and where Θ denotes appropriate zero mappings. Besides, there exists a strictly positive constant, Λ, which is a uniform lower bound for the eigenvalues of An1 . As a consequence of Hypothesis 3, there exists a unique positive definite symmetric n1 × n1 tensor function, Cn1 , such that An1 = (Cn1 )2 . Let us denote by C the symmetric and positive semidefinite d × d tensor   Cn1 Θ (2.20) C= . Θ Θ Notice that A = C 2 and C ∈ W1,∞ (Oδ ). Let us denote by G the matrix with coefficients Gij = | grad Cij |, 1 ≤ i, j ≤ d. At this point, let us introduce the constant (2.21)

cA = max{||G||2

∞, O

2 δ , ||C|| δ }, ∞, O

6

´ M. BEN´ITEZ AND A. BERMUDEZ

and the sequences of tensor functions (2.22) (2.23)

n n n n enRK := (FRK A )−1 A ◦ XRK (FRK )−T det FRK , p n n n −T n , e C det FRK RK := C ◦ XRK (FRK )

for 0 ≤ n ≤ N . Next, let us denote by B the d × d tensor   In1 Θ (2.24) B= , Θ Θ

where In1 is the n1 × n1 identity tensor. Clearly, under Hypothesis 3 we have (2.25)

Λ||Bw||2Ω ≤ hAw, wiΩ

∀w ∈ Rd .

Let us introduce the sequence of tensor functions p n n n eRK ∀n ∈ {0, 1, . . . , N }. B := B(FRK )−T det FRK

As far as the velocity field is defined in T δ (see Hypothesis 1), we can introduce the following assumption: Hypothesis 4. The velocity field satisfies (2.26)

(I − B)L(x, t)B = 0 ∀(x, t) ∈ T δ .

Remark 2.1. For any d × d tensor E of the form given in (2.19) it is easy to check that hEH T w1 , w2 i = hEH T Bw1 , Bw2 i, for any d × d tensor H satisfying (I − B)HB = 0, and vectors w1 , w2 ∈ Rd . It will be used below without explicitly stated. Hypothesis 5. Function f is defined in T δ and it is continuous with respect to the time variable, in space L2 . Hypothesis 6. Function g is defined in TΓδR and it is continuous with respect to the time variable, in space H 1 . Besides, coefficient α in boundary condition (2.5) is strictly positive. In what follows, cv denotes the positive constant (2.27)

cv := max ||v(·, t)||1,∞,Ωδt . t∈[0,T ]

Moreover, Cv (respectively, J and D) will denote a generic positive constant, related to the norm of the velocity field v (respectively, to the problem data), not necessarily the same at each occurrence. Let us introduce the notations  n+1 N −1 ψ − ψn d := {ψ n+1 + ψ n }N −1 , R\ S[ψ] [ψ] := , ∆t n=0 ∆t n=0

for a sequence ψb = {ψ n }N n=0 . Moreover, let us define the following sequence of functions of p. n n m e nRK = |(FRK )−T m| det FRK ,

for 0 ≤ n ≤ N . Hypothesis 7. Functions appearing in problem (2.3)-(2.6) satisfy:

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

7

• ρm ∈ C 2 (L∞ (Ω)), A ∈ W2,∞ (Oδ ), Am ∈ C 2 (W1,∞ (Ω)), • v ∈ C3 (T δ ), • fm ∈ C 2 (L2 (Ω)), f ∈ C 1 (T δ ), gm ∈ C 2 (L2 (ΓR )), g ∈ C 1 (TΓδR ) and α > 0. Hypothesis 8. Functions appearing in problem (2.3)-(2.6) satisfy: • ρm ∈ C 2 (L∞ (Ω)), A ∈ W2,∞ (Oδ ), Am ∈ C 3 (W1,∞ (Ω)), • v ∈ C3 (T δ ), • fm ∈ C 2 (L2 (Ω)), f ∈ C 1 (T δ ), gm ∈ C 3 (L2 (ΓR )), g ∈ C 2 (TΓδR ) and α > 0. 3. Space discretization. Finite element method. We propose a space discretization of the time semidiscretized problem introduced in [7] by using finite elements spaces Vhk , where h denotes the mesh-size and the positive integer k is the “approximation degree” in the following sense: Hypothesis 9. There exists an interpolation operator πh : C 0 (Ω) −→ Vhk satisfying ||πh ψ − ψ||s,2,Ω ≤ Qhr−s ||ψ||r,2,Ω

∀ψ ∈ C 0 (Ω) ∩ H r (Ω)

0 ≤ r ≤ k + 1, s = 0, 1,

for a positive constant Q independent of h. In order to obtain fully discrete schemes of the time semidiscretizated problem proposed in [7] (see §4 for more details), we use space Vhk to approximate space HΓ1D (Ω). Thus, we obtain the following fully discrete problem: (3.1)

(

 k N 0 k n N \ Given φ ∈ V , find φ = {φ } ∈ Vh such that m,∆t,h n=1 m,∆t,h h m,∆t,h E D E D n+ 12 n+ 21 \ k L∆t [φm,∆t,h ], ψh = F∆t , ψh ∀ψh ∈ Vh , for n = 0, . . . , N − 1. n+ 1

n+ 1

Mappings L∆t 2 [φ] ∈ (H 1 (Ω))′ and F∆t 2 ∈ (H 1 (Ω))′ are defined by * +
n+1 n+1 D E n n ρ ◦ XRK det FRK + ρ ◦ XRK det FRK φn+1 − φn n+ 21 L∆t [φ], ψ := ,ψ 2 ∆t Ω  *  en+1 +
n e ARK + ARK ∇φn+1 + ∇φn + , ∇ψ 2 2 Ω *
n+1
+ n n + φ m e n+1 + m e φ RK RK + α ,ψ , 2 2 ΓR   n+1 n+1 n+1 D E n n det FRK f ◦ XRK + det FRK f n ◦ XRK n+ 1 F∆t 2 , ψ := ,ψ 2 Ω  n+1 n+1  n+1 n m e RK g ◦ XRK + m e nRK g n ◦ XRK + ,ψ , 2 ΓR for φ ∈ C 0 (H 1 (Ω)) and ψ ∈ H 1 (Ω).

n+ 1

n+ 1

Remark 3.1. Regarding the definitions of L∆t 2 [φ] and F∆t 2 , only the values of function φ at discrete time steps {tn }N n=0 are required. Thus, the above definitions 1 N +1 can also be stated for a sequence of functions φb = {φn }N . n=0 ∈ [H (Ω)] By using the same procedures as the ones employed in [7] to get stability results of the semidiscretized scheme1 , we can obtain the following stability estimates for the fully discretized scheme. 1 By

replacing HΓ1 D (Ω) with Vhk .

´ M. BEN´ITEZ AND A. BERMUDEZ

8

Theorem 3.1. Let us assume Hypotheses 1, 2, 3, 5 and 6. Let φ\ m,∆t,h be the solution of (3.1) subject to the initial value φ0m,∆t,h ∈ Vhk . Then, there exist two positive constants J and D, independent of the diffusion tensor such that, for ∆t < D, we have r √ p Λ e \ \ γ det FRK φm,∆t,h + BRK S[∇φm,∆t,h ] 4 l∞ (L2 (Ω)) l2 (L2 (Ω)) r p α √ \ (3.2) + S [m e RK ]S[φm,∆t,h ] ≤J γ||φ0m,∆t,h ||Ω 8 l2 (L2 (ΓR ))  \ \ + det FRK f ◦ XRK 2 2 + m e RK g ◦ XRK 2 2 R . l (L (Ω))

l (L (Γ ))

Theorem 3.2. Let us assume Hypotheses 1, 2, 3, 4, 5 and 6, and let φ\ m,∆t,h be the solution of (3.1) subject to the initial value φ0m,∆t,h ∈ Vhk . Then there exist two positive constants J and D such that if ∆t < D then r r γ p Λ e \ \ B ∇φ S[ det F ]R [φ ] + RK m,∆t,h RK ∆t m,∆t,h 4 2 l∞ (L2 (Ω)) l2 (L2 (Ω)) r r α p \ Λ + m e RK φm,∆t,h ≤J B∇φ0m,∆t,h Ω 4 2 l∞ (L2 (ΓR )) (3.3)r α 0 \ + φm,∆t,h ΓR + det FRK f ◦ XRK 2 2 + m e RK\ g ◦ XRK ∞ 2 R 4 l (L (Ω)) l (L (Γ ))  + R∆t [m e\ g ◦ X ] . RK RK 2 2 R l (L (Γ ))

4. Error estimates for the fully discretized scheme. The aim of the present section is to estimate the difference between the discrete solution of (3.1), φ\ m,∆t,h := n N {φn }N , and the exact solution of the continuous problem, φc := {φ } . For m m,∆t,h n=0

m n=0

\ \ this, let us introduce the notations e\ m,∆t,h := φm,∆t,h − πh φm , ϑm,h := φm − πh φm . c \ [ [ Then, φm − φm,∆t,h = ϑm,h − e\ m,∆t,h and, since ϑm,h can be estimated by Hypothesis 9, the problem is reduced to establish a bound for e\ m,∆t,h . Notice that, according to solves the problem (2.12), for t 1 with 0 ≤ n ≤ N − 1, φc m n+ 2

(4.1)

D

E D E 1 n+ 12 Ln+ 2 [φc ,ψ m ], ψ = F

∀ψ ∈ HΓ1D (Ω),

1 1 ′ n+ 12 where Ln+ 2 [φc ∈ (H 1 (Ω))′ are defined by m ] ∈ (H (Ω)) and F D E   n+ 21  n+ 12 n+ 12 c n+ 12 L [φm ], ψ := ρ ◦ Xe det F φ˙ m ,ψ Ω D E D E 1 n+ 12 n+ 21 n+ 12 n+ 2 e + Am ∇φm , ∇ψ + α m e φm , ψ , Ω ΓR E D D E D E 1 1 1 1 1 1 n+ n+ 1 F n+ 2 , ψ := det F n+ 2 f n+ 2 ◦ Xe 2 , ψ + m e n+ 2 g n+ 2 ◦ Xe 2 , ψ

Ω

ΓR

,

∀ψ ∈ H 1 (Ω). In order to obtain error estimates in the l∞ (L2 (Ω))-norm let us state the following lemma.

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

9

Lemma 4.1. Assume Hypotheses 1, 2, 3, 9. Let φm ∈ C 1 (C 0 (Ω))∩C 0 (H k+1 (Ω))∩ H (H k (Ω)) be the solution of (4.1) and φ\ m,∆t,h be the solution of (3.1). Then there exist a positive constat c(v, T, δ) such that, for ∆t < c, the following inequality holds: 1

(4.2)

D

E n+ 1 [ n+1 n L∆t 2 [ϑ m,h ], em,∆t,h + em,∆t,h  2 1 en+1  n+1 n ≤ C ∇e + ∇e m,∆t,h RK m,∆t,h 8 Ω   1 en 2 n+1 n + CRK ∇em,∆t,h + ∇em,∆t,h 8 Ω q n o 2 α n+1 n+1 n n e RK + m e RK em,∆t,h + em,∆t,h + m 16 ΓR q 2 p 2 n+1 n+1 n en +γ det FRK em,∆t,h + γ det FRK m,∆t,h Ω Ω   2 1 ˙ n+1 2 n 2 2 2k + φm + ||φm ||k+1,2,Ω , +e cQ h φm 2 k+1,2,Ω ∆t L ((tn ,tn+1 ),H k (Ω))

being e c a positive constant, n ∈ {0, . . . , N − 1} and where α > 0 is the constant appearing in the Robin boundary condition (2.5). D E n+ 1 [ n+1 n n n n Proof. First, we decompose L∆t 2 [ϑ m,h ], em,∆t,h + em,∆t,h = I1 + I2 + I3 with +
n+1 n+1 n+1 n n ϑm,h − ϑnm,h n+1 ρ ◦ XRK det FRK + ρ ◦ XRK det FRK n = , em,∆t,h + em,∆t,h , 2 ∆t Ω   E 1 D en+1 n+1 n n en ARK + A ∇ϑn+1 + ∇ϑ , ∇e + ∇e , I2n = RK m,h m,∆t,h m,h m,∆t,h 4 D Ω   E
α n+1 n n I3n = m e n+1 e nRK ϑn+1 . RK + m m,h + ϑm,h , em,∆t,h + em,∆t,h 4 ΓR I1n

*

For I1n , aplying Cauchy-Schwarz inequality, Young’s inequality and Corollary 4.4 in [7], we first have 2 n+1 q 2 n ϑ p 2 − ϑ m,h m,h n+1 n+1 n n en I1 ≤ e c + γ det FRK em,∆t,h + γ det FRK m,∆t,h , ∆t Ω Ω Ω

(4.3)

where we have assumed that ∆t < K, being K the constant appearing in Corollary 4.4 in [7]. Here e c is a positive constant depending on v, T and ρ1,∞ /γ. Moreover, from Barrow’s rule, we have (4.4)

n ϑn+1 m,h (p) − ϑm,h (p)

∆t

1 = ∆t

Z

tn+1

ϑ˙ m,h (p, s)ds.

tn

Thus, by applying Cauchy-Schwarz inequality, we deduce (4.5)

Z tn+1  ϑn+1 (p) − ϑn (p) 2  12 1 m,h m,h ˙ ϑm,h (p, s) ds , ≤ √ ∆t ∆t tn

10

´ M. BEN´ITEZ AND A. BERMUDEZ

and then

(4.6)

n+1 2 Z Z tn+1  n ϑ 2 1 m,h − ϑm,h ϑ˙ m,h (p, s) ds dp ≤ ∆t ∆t Ω tn Z tn+1 Z Ω 2 2 1 1 ˙ ϑ˙ m,h (p, s) dp ds = . = ϑm,h 2 ∆t tn ∆t L ((tn ,tn+1 ),L2 (Ω)) Ω

Finally, by using Hypothesis 9 for s = 0 and r = k we obtain q 2 e cQ2 h2k ˙ 2 n+1 n+1 n + γ det FRK em,∆t,h I1 ≤ φm 2 ∆t L ((tn ,tn+1 ),H k (Ω)) Ω (4.7) p 2 n en +γ det FRK m,∆t,h . Ω

For I2n , we first apply Cauchy-Schwarz inequality and Young’s inequality, obtaining   2  2 1 en+1  n+1 n 2 I2n ≤ e c ∇ϑn+1 + ∇ϑ CRK ∇em,∆t,h + ∇enm,∆t,h m,h Ω + m,h 8 Ω Ω   2 1 en n+1 n + C RK ∇em,∆t,h + ∇em,∆t,h . 8 Ω

where we have used inequalities (4.11) and (4.15) from [7]. Here e c is a positive constant depending on v, T and cA and is bounded in the hyperbolic limit. From this inequality and by using Hypothesis 9 for s = 1 and r = k + 1 we obtain
2 n 2 I2n ≤ e cQ2 h2k ||φn+1 m ||k+1,2,Ω + ||φm ||k+1,2,Ω  2  2 1 e n+1  n+1 1 en  n+1 n n (4.8) + C ∇e + ∇e + C ∇e + ∇e . m,∆t,h RK m,∆t,h RK m,∆t,h m,∆t,h 8 8 Ω Ω

For I3n we apply Cauchy-Schwartz inequality, Young’s inequality and inequalities (4.11) and (4.15) from [7], getting   2 n 2 I3n ≤ e c ||ϑn+1 || + ||ϑ || R R m,h Γ m,h Γ (4.9) q n o 2 α n+1 n+1 n n + m e RK + m e RK em,∆t,h + em,∆t,h . 16 ΓR

Next, from the continuity of the trace mapping, there exist a positive constant cΩ such that ||ϑlm,h ||2ΓR ≤ cΩ ||ϑlm,h ||21,2,Ω , for l = n, n + 1. Then, by applying Hypothesis 9 for s = 1 and r = k + 1 in (4.9), we have
2 n 2 I3n ≤ e cQ2 h2k ||φn+1 m ||k+1,2,Ω + ||φm ||k+1,2,Ω (4.10) q n o 2 α n+1 n+1 n n + m e RK + m e RK em,∆t,h + em,∆t,h . 16 R Γ

Finally, summing up (4.7), (4.8) and (4.10) we get inequality (4.2). Theorem 4.2. Let us assume Hypotheses 1, 2, 3, 5, 6, 7 and 9, and Xe ∈ C5 (Ω × [0, T ]). Let (4.11)

φm ∈ C 3 (L2 (Ω)) ∩ C 1 (C 0 (Ω)) ∩ C 0 (H k+1 (Ω)) ∩ H 1 (H k (Ω)), ∇φm ∈ C 2 (H1 (Ω)), φm |ΓR ∈ C 2 (L2 (ΓR )),

11

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

be the solution of (4.1) and let φ\ m,∆t,h be the solution of (3.1) subject to the initial value φ0m,∆t,h = πh φ0m . Then, there exist two positive constants J and D, being the latter independent of the diffusion tensor, such that, if ∆t < D we have r γ p \ || det FRK (φm − φm,∆t,h )||l∞ (L2 (Ω)) 2 r Λ e \ B S [∇φ − ∇φ ] + RK m m,∆t,h 8 l2 (L2 (Ω)) r α p \ 2 + S [ m e ]S [φ − φ ] RK m m,∆t,h 2 2 R ≤ J ∆t ||φm ||C 3 (L2 (Ω)) (4.12) 8 l (L (Γ )) +||∇φm ||C 2 (H1 (Ω)) + ||∇φm · m||C 2 (L2 (ΓR )) + ||φm ||C 2 (L2 (ΓR ))

+|| det F fm ||C 2 (L2 (Ω)) + ||f ||C 1 (T δ ) + ||mg e m ||C 2 (L2 (ΓR )) + ||g||C 1 (T δR ) Γ   + ||φm ||C 0 (H k+1 (Ω)) . +Jhk φ˙ m 2 k



L (H (Ω))

k N +1 c [ \ Proof. First, recall that e\ . Then, by m,∆t,h = ϑm,h − φm + φm,∆t,h ∈ [Vh ] using (4.1) and (3.1) we have E D E D n+ 12 [ n+ 1 n+1 n+1 n n L∆t 2 [e\ m,∆t,h ], em,∆t,h + em,∆t,h = L∆t [ϑm,h ], em,∆t,h + em,∆t,h D  E D E 1 1 n+ 1 n+ 12 n+1 n n + Ln+ 2 − L∆t 2 [φc − F n+ 2 , en+1 + e m ], em,∆t,h + em,∆t,h + F∆t m,∆t,h , m,∆t,h

(4.13) for n ∈ {0, . . . , N − 1}. A lower bound for (4.13) is given by Lemma 4.8 in [7], namely, D E n+ 1 n+1 n (4.14) L∆t 2 [e\ m,∆t,h ], em,∆t,h + em,∆t,h q 2 2 1 1 p n+1 n+1 n+1 n det F n en ≥ ρ ◦ X det F e ρ ◦ XRK RK RK m,∆t,h − RK m,∆t,h ∆t ∆t Ω Ω     2 1 en+1 1 e n 2 n n ∇en+1 + C CRK ∇en+1 RK m,∆t,h + ∇em,∆t,h + m,∆t,h + ∇em,∆t,h 4 4 Ω Ω q 2   α n + m e n+1 e nRK en+1 RK + m m,∆t,h + em,∆t,h 4 ΓR ! q 2 2 p n+1 n+1 n n −b cγ det FRK em,∆t,h + det FRK em,∆t,h , Ω Ω

where b c = ρ1,∞ (cv + Cv ∆t)/γ and n ∈ {0, . . . , N − 1}. Next, by applying Lemmas 4.25 and 4.26 in [7] and Lemma 4.10 in [7] for the choices n+ 1

n+ 12

n n+1 ψ = en+1 = ξLΩ 2 , Gn+1 = ξLΓ m,∆t,h , ϕ = em,∆t,h , first for S n+ 12

−ξf

n+ 12

, Gn+1 = −ξg

and then for S n+1 =

, we have

(4.15) D  E D E 1 1 n+ 1 n+ 21 n+1 n n Ln+ 2 − L∆t 2 φc − F n+ 2 , en+1 m , em,∆t,h + em,∆t,h + F∆t m,∆t,h + em,∆t,h D E D E n+ 1 n+ 1 n+ 12 n+ 12 n+1 n n = ξLΩ 2 − ξf 2 , en+1 + e + ξ − ξ , e + e g m,∆t,h m,∆t,h m,∆t,h LΓ m,∆t,h Ω

ΓR

´ M. BEN´ITEZ AND A. BERMUDEZ

12

  4c   n+ 12 2 n+ 12 2 n+ 21 2 n+ 12 2 g ≤ cs ξLΩ + ξf + ξLΓ R + ξg R α Ω Γ Ω Γ q p n+1 n+1 2 n en +|| det FRK em,∆t,h ||2Ω + || det FRK m,∆t,h ||Ω q   2 α n+1 n+1 n n + m e RK + m e RK em,∆t,h + em,∆t,h , 16 ΓR

being cs and cg the positive constants appearing in Lemma 4.10 from [7]. By jointly considering the lower and upper bounds of (4.13) given in (4.14) and (4.15), respectively, and inequality (4.2) we deduce q 2 p 2 1 n+1 n+1 n+1 − 1 ρ ◦ X n det F n en ρ ◦ X det F e RK RK m,∆t,h RK RK m,∆t,h ∆t ∆t Ω Ω  2   2 1 e n+1  n+1 1 n+1 n en + CRK ∇em,∆t,h + ∇enm,∆t,h + C RK ∇em,∆t,h + ∇em,∆t,h 8 8 Ω Ω q 2   α n + m e n+1 e nRK en+1 RK + m m,∆t,h + em,∆t,h 8 ΓR 2 2  4c   (4.16)  1 1 n+ n+ n+ 12 2 n+ 1 2 g ≤ cs ξLΩ 2 + ξf 2 + ξLΓ R + ξg 2 R α Ω Ω Γ Γ ! 2 q 2 p n+1 n+1 n en em,∆t,h + det FRK +3b cγ det FRK m,∆t,h +e cQ2 h2k



Ω

Ω

 n+1 2 1 ˙ 2 n 2 + φm + ||φm ||k+1,2,Ω , φm 2 k+1,2,Ω ∆t L ((tn ,tn+1 ),H k (Ω))

where b c = max {1, 1/γ, ρ1,∞(cv + Cv ∆t)/γ} and e c is a positive constant. For n = 0, . . . N , let us introduce the notations p 2 n en θn1 := γ|| det FRK m,∆t,h ||Ω n−1   2 ΛX e s s+1 2 s θn := ∆t B RK ∇em,∆t,h + ∇em,∆t,h , 8 s=0 Ω q n−1   2 αX s+1 s s . θn := e s+1 + m e e + e ∆t m m,∆t,h RK RK m,∆t,h 8 s=0 ΓR

Now, for fixed q, 1 ≤ q ≤ N , let us sum (4.16) multiplied by ∆t from n = 0 to n = q − 1. We have, q−1 X

ρ1,∞ 1 θ γ 0 n=0 q  q   4c ∆t X  X n− 21 2 n− 12 2 n− 12 2 n− 12 2 g +cs ∆t + ξLΓ R + ξg R ξLΩ + ξf α n=1 Ω Ω Γ Γ n=1  q−1   X 2 2 2k n+1 2 n 2 ˙ +e cQ h + ∆t φm + ||φm ||k+1,2,Ω , φm k+1,2,Ω (1 −

(4.17)

n=0

3b c∆t)θq1

L2 ((tn ,tn+1 ),H k (Ω))

+

θq2

+ θq ≤ 6b c∆t

θn1 +

13

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

by using Hypotheses 2 and 3. Some of the terms on the right hand side of (4.17) can also be bounded. Firstly, we have q−1 X ˙ 2 φm 2

n=0

L ((tn ,tn+1 ),H k (Ω))

q−1  X n+1 2 φ + ∆t n=0

2 ≤ φ˙ m 2

2

+ ||φnm ||k+1,2,Ω k+1,2,Ω

m

L (H k (Ω))



+ 2T ||φm ||2C 0 (H k+1 (Ω)) .

Secondly, by using Lemmas 4.25 and 4.26 in [7] we get q  q   4c ∆t X  X n− 21 2 n− 12 2 n− 1 2 n− 12 2 g + ξLΩ + ξf ξLΓ R + ξg 2 R α n=1 Ω Ω Γ Γ n=1  ≤e c∆t4 ||φm ||2C 3 (L2 (Ω)) + ||∇φm ||2C 2 (H1 (Ω))

cs ∆t

+||∇φm · m||2C 2 (L2 (ΓR )) + ||φm ||2C 2 (L2 (ΓR )) + || det F fm ||2C 2 (L2 (Ω))  2 2 2 +||f ||C 1 (T δ ) + ||mg e m ||C 2 (L2 (ΓR )) + ||g||C 1 (T δ ) . ΓR

These estimates lead to

(1 −

3b c∆t)θq1

+

θq2

+ θq ≤ 6b c∆t

q−1 X

n=0

  e θn1 + e c θ01 + C

e contains the constant terms multiplied by h2k and ∆t4 . For ∆t small enough, where C we can apply the discrete Gronwall inequality (see, for instance, [29]) and take the \ [ maximun in q ∈ {1, . . . , N }. Then, noting that e0m,∆t,h = 0, φc m − φm,∆t,h = ϑm,h − e\ m,∆t,h , using Hypothesis 9, and bounds (4.11) and (4.15) from [7] the result follows. Remark 4.1. Notice that constant J appearing in the previous theorem is bounded in the limit when the diffusion tensor vanishes. In particular, Theorem 4.2 is also valid when A ≡ 0. Remark 4.2. In the particular case of pure convection problems, that is A ≡ 0, and assuming Dirichlet boundary conditions (ΓD ≡ Γ), an error estimate of the form O(hk+1 ) + O(∆t2 ) in the l∞ (L2 (Ω))-norm can be obtained by using analogous procedures to the ones in the previous theorem. In order to obtain error estimates in the l∞ (H 1 (Ω))-norm, let us state the following lemma. Lemma 4.3. Assume Hypotheses 1, 2, 3, 4 and 9. Let φm ∈ C 1 (C 1 (Ω)) ∩ C (H k+1 (Ω)) ∩ H 1 (H k+1 (Ω)) with ∇φm ∈ C 1 (C0 (Ω)) be the solution of (4.1), and φ\ m,∆t,h be the solution of (3.1). Then there exist a positive constant c(v, T, δ), such 0

´ M. BEN´ITEZ AND A. BERMUDEZ

14

that, for ∆t < c and q ∈ {1, . . . , N }, the following inequality holds: q−1 D X

n=0

n+ 1

n+1 n [ L∆t 2 [ϑ m,h ], em,∆t,h − em,∆t,h

E

q−1 q   2 1 X n+1 n+1 n+1 n n n ≤ ρ ◦ XRK det FRK + ρ ◦ XRK det FRK em,∆t,h − em,∆t,h 4∆t n=0 Ω q 2 Λ e q α 2 q + B m e qRK eqm,∆t,h RK ∇em,h + 4 8 Ω ΓR (4.18) q−1 q−1 q 2 2 X X e n e nRK enm,∆t,h R +b c∆tΛ cα∆t BRK ∇enm,∆t,h + b m Ω

n=1

n=1

Γ

 E
1 D e1 − ARK + A ∇ϑ1m,h + ∇ϑ0m,h , ∇e0m,∆t,h 4 Ω
1
0  α 1 0 − m e RK + 1 ϑm,h + ϑm,h , em,∆t,h ΓR 4   2 2 2 2k ˙ +e cQ h + ||φm ||C 0 (H k+1 (Ω)) , φm 2 k+1 L (H

(Ω))

where b c = max{Cv cA /Λ, Cv }, e c is a positive constant and α > 0 is the constant appearing in the Robin boundaryD condition (2.5). E n+ 1 [ n+1 n n n n Proof. First, we decompose L∆t 2 [ϑ m,h ], em,∆t,h − em,∆t,h = I1 + I2 + I3 , with * +
n+1 n+1 n+1 n n ϑm,h − ϑnm,h n+1 det FRK ρ ◦ XRK det FRK + ρ ◦ XRK n n I1 = , em,∆t,h − em,∆t,h , 2 ∆t Ω   E 1 D en+1 n+1 n n en I2n = ARK + A ∇ϑn+1 , RK m,h + ∇ϑm,h , ∇em,∆t,h − ∇em,∆t,h 4 D Ω   E
α n+1 n n I3n = m e n+1 e nRK ϑn+1 , RK + m m,h + ϑm,h , em,∆t,h − em,∆t,h 4 ΓR for 0 ≤ n ≤ N − 1. By applying Cauchy-Schwarz inequality and Young’s inequality we have q   2 1 n+1 n+1 n+1 n n n n I1 ≤ ρ ◦ XRK det FRK + ρ ◦ XRK det FRK em,∆t,h − em,∆t,h 4∆t Ω 2 e c n+1 n + ϑm,h − ϑm,h , ∆t Ω where we have used Corollary 4.4 in [7] and Hypothesis 2. Moreover, by using inequality (4.6) and Hypothesis 9 for s = 0 and r = k we obtain the following bound q−1 q q−1   2 X 1 X n+1 n+1 n+1 n n det F n I1n ≤ ρ ◦ X det F + ρ ◦ X e − e m,∆t,h RK RK RK RK m,∆t,h 4∆t n=0 Ω n=0 2 +e cQ2 h2k φ˙ m 2 k . L (H (Ω))

(4.19)

A bound for

q−1 X

I2n follows from the equality

n=0 q−1 X

q n q 1 0 hwn+1 , ∇en+1 m,∆t,h − ∇em,∆t,h iΩ = hw , ∇em,∆t,h iΩ − hw , ∇em,∆t,h iΩ

n=0

15

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

(4.20) −∆t

q−1  n+1 X w − wn

n=1

∆t

, ∇enm,∆t,h



, Ω

n+1 2 N being {wn }N = n=1   ∈ [L (Ω)] . More  precisely, by using equality (4.20) for w n+1 n+1 n n e e ARK + ARK ∇ϑm,h + ∇ϑm,h , Lemma 4.3, Corollaries 4.4 and 4.5, equality (4.48)

from [7], Cauchy-Schwarz inequality and Young’s inequality, the following estimate can be easily proved,  q−1  Λ 2 X q 2 q−1 2 e q q n + B I2 ≤ e c ∇ϑm,h + ∇ϑm,h RK ∇em,∆t,h 4 Ω Ω Ω n=0 q−1 2 X e n +Cv cA ∆t BRK ∇enm,∆t,h

Ω

n=1

(4.21)

 E
1 D e1 − ARK + A ∇ϑ1m,h + ∇ϑ0m,h , ∇e0m,∆t,h 4 Ω 2 q q−1 n+1 n X X ∇ϑ − ∇ϑ m,h m,h +e c∆t ||∇ϑnm,h ||2Ω + e c∆t . ∆t n=0

n=0

Ω

Some of the terms on the right hand side can be bounded. By using Hypothesis 9 for s = 1 and r = k + 1 we obtain q X q 2 q−1 2 ||∇ϑnm,h ||2Ω ∇ϑm,h + ∇ϑm,h + ∆t Ω

(4.22)

Ω

n=0

2 2k

≤ (2 + T )Q h

2

||φm ||C 0 (H k+1 (Ω)) .

Next, by using Barrow’s rule, Cauchy-Schwarz inequality and Hypothesis 9 for s = 1 and r = k + 1 we deduce 2 q−1 n+1 n 2 X ∇ϑm,h − ∇ϑm,h ∆t . ≤ Q2 h2k φ˙ m 2 k+1 ∆t L (H (Ω)) n=0

Ω

These estimates lead to q−1 q−1 2 2 X X 1 e q e n q n I2 ≤ CRK ∇em,∆t,h + Cv cA ∆t BRK ∇enm,∆t,h 4 Ω Ω n=0 n=1 D  E
1 (4.23) 0 e1 + A ∇ϑ1 + ∇ϑ0 − A RK m,h m,h , ∇em,∆t,h 4 Ω   2 2 2 2k ˙ +e cQ h + ||φm ||C 0 (H k+1 (Ω)) . φm L2 (H k+1 (Ω))

Now, we obtain an estimate for

q−1 X

I3n . For this purpose, we use the following equality

n=0

(4.24)

q−1 X

n=0

n hψ n+1 , en+1 m,∆t,h − em,∆t,h iΓR

=

hψ q , eqm,∆t,hiΓR − hψ 1 , e0m,∆t,h iΓR −∆t

q−1  n+1 X ψ − ψn

n=1

∆t

, enm,∆t,h



ΓR

,

´ M. BEN´ITEZ AND A. BERMUDEZ

16 for ψ n+1 = m e n+1 e nRK RK + m




 n ϑn+1 + ϑ m,h , as well as Lemma 4.3, Corollaries 4.4 and m,h

4.5 from [7], Cauchy-Schwarz inequality, and Young’s inequality. We get

2  2 2  α q q q I3n ≤ e c ϑqm,h R + ϑq−1 + m e e RK m,∆t,h m,h R 8 Γ Γ ΓR n=0 q−1 X

+Cv α∆t

(4.25)

q−1 q 2 X

 α 1 e nRK enm,∆t,h R − m e RK + 1 ϑ1m,h + ϑ0m,h , e0m,∆t,h ΓR m 4 Γ n=1 2 q q−1 n+1 n X X ϑm,h − ϑm,h n 2 +e c∆t ||ϑm,h ||ΓR + e c∆t . R ∆t n=0

n=0

Γ

We use the continuity of the trace operator and Hypothesis 9 for s = 1 and r = k + 1 to bound some terms on the right hand side of (4.25), getting q X q 2 q−1 2 ||ϑnm,h ||2ΓR ≤ (2 + T )cΩ Q2 h2k ||φm ||2C 0 (H k+1 (Ω)) ϑm,h R + ϑm,h R + ∆t Γ

Γ

n=0

and 2 q−1 n+1 n 2 X ϑm,h − ϑm,h ∆t . ≤ cΩ Q2 h2k φ˙ m 2 k+1 R ∆t L (H (Ω)) n=0

Γ

To obtain the last inequality, we also have used Barrow’s rule and Cauchy-Schwarz inequality. We substitute the preceding estimates into (4.25) to obtain q−1 X

n=0

(4.26)

I3n ≤

q 2 q−1 q 2 X α q q + Cv α∆t m e e e nRK enm,∆t,h R m RK m,∆t,h 8 Γ R Γ n=1



 α − m e 1RK + 1 ϑ1m,h + ϑ0m,h , e0m,∆t,h ΓR 4   2 2 2 2k ˙ +e cQ h + ||φm ||C 0 (H k+1 (Ω)) . φm 2 k+1 L (H

(Ω))

Finally, by summing up (4.19), (4.23) and (4.26) we get inequality (4.18).

Theorem 4.4. Let us assume Hypotheses 1, 2, 3, 4, 5, 6, 8 and 9, and Xe ∈ C5 (Ω × [0, T ]). Let (4.27)

φm ∈ C 3 (L2 (Ω)) ∩ C 1 (C 1 (Ω)) ∩ C 0 (H k+1 (Ω)) ∩ H 1 (H k+1 (Ω)) with ∇φm ∈ C 3 (H1 (Ω)) and φm |ΓR ∈ C 3 (L2 (ΓR )),

be the solution of (4.1) and φ\ m,∆t,h the solution of (3.1) subject to the initial value φ0m,∆t,h = πh φ0m . Then, there exist two positive constants J and D such that, if

17

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

∆t < D, we have r γ p \ S[ det FRK ]R∆t [φm − φm,∆t,h ] 8 l2 (L2 (Ω)) r Λ e \ BRK (∇φm − ∇φm,∆t,h ) + 4 l∞ (L2 (Ω)) r p α \ e RK (φm − φm,∆t,h ) ≤ J ∆t2 ||φm ||C 3 (L2 (Ω)) (4.28) + 4 m ∞ 2 R l

(L (Γ ))

+||∇φm ||C 2 (H1 (Ω)) + ||∇φm · m||C 3 (L2 (ΓR )) + ||φm ||C 3 (L2 (ΓR ))

+|| det F fm ||C 2 (L2 (Ω)) + ||f ||C 1 (T δ ) + ||mg e m ||C 3 (L2 (ΓR )) + ||g||C 2 (T δR ) Γ   k ˙ +J h φm 2 k+1 + ||φm ||C 0 (H k+1 (Ω)) . L (H



(Ω))

k N +1 [ c \ Proof. First, we recall that e\ . Then, by m,∆t,h = ϑm,h − φm + φm,∆t,h ∈ [Vh ] using (4.1) and (3.1) we deduce

D E D E n+ 1 n+ 12 [ n+1 n+1 n n L∆t 2 [e\ m,∆t,h ], em,∆t,h − em,∆t,h = L∆t [ϑm,h ], em,∆t,h − em,∆t,h D  E D E 1 1 n+ 1 n+ 12 n+1 n n + Ln+ 2 − L∆t 2 [φc − F n+ 2 , en+1 − e . m ], em,∆t,h − em,∆t,h + F∆t m,∆t,h m,∆t,h (4.29) A lower bound for (4.29) is given by Lemma 4.9 from [7], namely D

E n+ 1 n+1 n L∆t 2 [e\ m,∆t,h ], em,∆t,h − em,∆t,h q  2
 n+1 1 n+1 n+1 n n det F n ≥ ρ ◦ X det F + ρ ◦ X e − e m,∆t,h RK RK RK RK m,∆t,h 2∆t Ω 2 1 e n 1 en+1 n+1 2 n + CRK ∇em,∆t,h − CRK ∇em,∆t,h 2 2 Ω Ω (4.30) q 2 q 2 α n+1 − α m + m e n+1 e nRK enm,∆t,h R RK em,∆t,h 2 2 Γ R Γ   n n 2 e n+1 ∇en+1 ||2Ω + ||B eRK −b c∆tΛ ||B ∇e || m,∆t,h Ω RK m,∆t,h ! q 2 q 2 n+1 + m −b c∆tα m e n+1 e nRK enm,∆t,h , RK em,∆t,h ΓR

ΓR

where b c = max {cA Cv /Λ, Cv } and n ∈ {0, . . . , N − 1}. By applying Lemmas 4.25 and n 4.26 in [7], and Lemma 4.13 in [7] for the choices ψ = en+1 m,∆t,h , ϕ = em,∆t,h , first for n+ 1

n+ 12

S n+1 = ξLΩ 2 and then for S n+1 = −ξf

, an upper bound of (4.29) can be easily

´ M. BEN´ITEZ AND A. BERMUDEZ

18

obtained. By considering both estimates, the following inequality holds, q  2
 n+1 1 n+1 n+1 n n n ρ ◦ XRK det FRK + ρ ◦ XRK det FRK em,∆t,h − em,∆t,h 2∆t Ω 2 1 en+1 n+1 2 1 en n + CRK ∇em,∆t,h − CRK ∇em,∆t,h 2 2 Ω Ω 2 q q 2 α n+1 n en − α m + m e n+1 e e RK m,∆t,h RK m,∆t,h R 2 2 Γ R Γ D E 2c ∆t   n+ 12 n+ 1 n+ 1 s n+1 n [ ≤ L∆t [ϑ ||ξLΩ 2 ||2Ω + ||ξf 2 ||2Ω m,h ], em,∆t,h − em,∆t,h + γ q 2 (4.31) γ n+1 n+1 n n + det FRK + det FRK (em,∆t,h − em,∆t,h ) 8∆t Ω D E n+ 21 n+ 12 n+1 n + ξLΓ − ξg , em,∆t,h − em,∆t,h ΓR  2 2  e n+1 n+1 e n n +b c∆tΛ BRK ∇em,∆t,h + BRK ∇em,∆t,h Ω Ω ! q 2 q 2 n+1 n+1 n +b c∆tα m e RK em,∆t,h + m e nRK em,∆t,h , R Γ

ΓR

where cs is the constant appearing in Lemma 4.10 from [7]. For n = 0, . . . , N , let us introduce the notations n−1 q   2 γ X s+1 s+1 s s := det FRK + det FRK em,∆t,h − em,∆t,h 8∆t s=0 Ω q 2 2 α Λ e n n θn2 := B θn := m e nRK enm,∆t,h . RK ∇em,∆t,h , 4 4 ΓR Ω

θn1

Now, for fixed q ≤ 1, let us sum (4.31) from n = 0 to n = q − 1. Then, with the above notation we have θq1 + (1 − 4b c∆t)θq2 + (1 − 4b c∆t)θ q ≤ 12b c∆t

q−1 X

n=0

θn2 + 16b c∆t

q−1 X

θn

n=0

q  4c   2cs ∆t X  n− 12 2 n− 1 q− 1 q− 1 g + ||ξLΩ ||Ω + ||ξf 2 ||2Ω + ||ξLΓ 2 ||2ΓR + ||ξg 2 ||2ΓR γ n=1 α 2 n− 12 q−1 n+ 12  ∆tc X 1 ξ − ξ 1  12 2 g LΓ LΓ + ||ξLΓ ||ΓR + ||ξg2 ||2ΓR + 2α 2α n=1 ∆t

+



q−1 n+ 1 ∆tcg X ξg 2



ΓR

− ∆t

n− 1 ξg 2

2   2 2 2 2k ˙ cQ h + ||φm ||C 0 (H k+1 (Ω)) , + e φm 2 k+1 R L (H (Ω))

2α n=1 Γ (4.32) where we have used Hypotheses 2 and 3, inequality (4.18), e0m,∆t,h = 0, and Lemma n+ 1 n+1 4.14 in [7] for the choice ψb = e\ = ξLΓ 2 and then for Gn+1 = m,∆t,h , first for G n+ 21

−ξg

. Some of the terms on the right hand side of (4.32) can also be bounded as

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

19

follows: Firstly, by applying Lemmas 4.25 and 4.26 in [7] we deduce ∆t

q  X

n=1

(4.33)

n− 1

n− 12

||ξLΩ 2 ||2Ω + ||ξf

 q− 1 q− 1 ||2Ω + ||ξLΓ 2 ||2ΓR + ||ξg 2 ||2ΓR

 1 1 +||ξL2 Γ ||2ΓR + ||ξg2 ||2ΓR ≤ e c∆t4 ||φm ||2C 3 (L2 (Ω)) + ||∇φm ||2C 2 (H1 (Ω)) 2

2

+ ||φm ||C 2 (L2 (ΓR )) + ||mg e m ||C 2 (L2 (ΓR )) + ||∇φm · m||2C 2 (L2 (ΓR ))  +|| det F fm ||2C 2 (L2 (Ω)) + ||f ||2C 1 (T δ ) + ||g||2C 1 (T δ ) . ΓR

Secondly, the terms

\ R∆t [ξLΓ ] 2

l (L2 (ΓR ))

+ R\ ∆t [ξg ] 2

l (L2 (ΓR ))

are bounded as in Theorem 4.28 from [7]. We incorporate these estimates into (4.32) to get c∆t)θ q ≤ 12b θq1 + (1 − 4b c∆t)θq2 + (1 − 4b c∆t

q−1 X

n=0

θn2 + 16b c∆t

q−1 X

n=0

e θn + C,

e contains the constant terms multiplied by h2k and ∆t4 . For ∆t small enough, where C we can apply the discrete Gronwall inequality (see, for instance, [29]) and take the \ [ maximun in q ∈ {1, . . . , N }. Then, noting that φc m − φm,∆t,h = ϑm,h − e\ m,∆t,h , by using Hypothesis 9, bounds (4.11) and (4.15) in [7], and the following estimate (see, for instance, the proof of Lemma 4.1) p k ˙ S[ det F\ ]R [ϑ ] ≤ c e Qh , φm RK ∆t m,h L2 (H k (Ω))

l2 (L2 (Ω))

the result is concluded. Remark 4.3. In the particular case of diffusion tensor of the form A = ǫB with ǫ > 0, constants J and D appearing in the previous theorem are bounded as ǫ → 0. Approximate solution in Eulerian coordinates

In order to obtain an approximate solution of φn in Eulerian coordinates, we are going to calculate the spatial description of material field φnm,∆t,h . To do this, we distinguish two cases: b • Xe known. In this case, we calculate φ[ ∆t,h ∼ φ as follows (4.34)

φn∆t,h (x) := φnm,∆t,h (P (x, tn ))

∀x ∈ Ωtn , 0 ≤ n ≤ N,

where P is the reference map of motion Xe (see [7] for more details). • Xe unknown. In this case, we use accurate enough approximations of P preserving the error order of the method. More precisely, we use the second order Kunge-Kutta method considered to approximate the characteristics curves. Then, we calculate φ[ ∆t,h as follows (4.35)

n φn∆t,h (x) := φnm,∆t,h (PRK (x))

∀x ∈ Ωtn , 0 ≤ n ≤ N.

´ M. BEN´ITEZ AND A. BERMUDEZ

20

n being PRK the second order Kunge-Kutta approximation of P n . Notice that, n for a general velocity field, point PRK (x) can go out of the computational n domain. In this case, we approximate φn∆t,h (PRK (x)) by n φn∆t,h (PRK (x)) ≃ φnm,∆t,h (xf ),

(4.36)

n being xf the nearest point on the boundary to PRK (x). Notice that, if the velocity vanishes on the boundary of Ω and ∆t is small enough, then n PRK (Ω) = Ω (see Lemma 4.7 in [7]). In Example 2 below v satisfies this property. Remark 4.4. Notice that, from the estimates obtained in Lagrangian coordinates and by using appropriate changes of variable, we can deduce analogous estimates in Eulerian coordinates (see [6] for further details).

5. Numerical results. In order to assess the performance of the above numerical method and to check the convergence behavior predicted by the above theory, we solve two test problems in two space dimensions. The first one is the rotating Gaussian hill, for which we verify rates of convergence for the second order pure Lagrangian method described in the present paper and the analogous one of first order in time. The second example has a solution developing a steep layer and a velocity field which is not divergence-free. For this problem, we compare the numerical results obtained from the pure Lagrangian method proposed in this paper, with the analogous one of first order in time and with semi-Lagrangian methods. In Example 1, we calculate the error between discrete solution φh,∆t , given in (4.34), and exact solution φ. For this, we approximate the theoretical H 1 (Ωtn ) and L2 (Ωtn ) norms by using a quadrature formula exact for polynomials of degree 5. The functional spaces endowed with these norms are denoted by Hh1 (Ωtn ) and L2h (Ωtn ), respectively. Thus, we denote by l∞ (Hh1 (Ωtn )) and l∞ (L2h (Ωtn )) the spaces equipped with the norms N N b := max ||ψ n ||H 1 (Ωtn ) , ψb ∞ 2 := max ||ψ n ||L2 (Ωtn ) . ψ ∞ 1 l

(Hh (Ωtn ))

n=0

h

l

(Lh (Ωtn ))

n=0

h

Firstly, we show numerical results for the problem of the rotating Gaussian hill and then for the problem including a steep layer. Example 1

This is a convection-diffusion problem, see for instance [30] and [13], aiming to check the above properties of the proposed scheme and to compare them with the computed solution obtained by using a pure Lagrangian characteristics method of first order in time. Moreover, we compare the computed solution by using the standard first order characteristics method combined with piecewise linear finite elements, with the one obtained from the second order method proposed in this paper. The spatial domain is Ω = (−1., 1.) × (−1., 1.) and T = 2π. The diffusion tensor is A = σ1 I with σ1 given below. Moreover, v = (−x2 , x1 ), ρ = 1 and the right-hand side f = 0. We also impose appropriate Dirichlet boundary and initial conditions such that the solution of the problem is   σ2 (x(t) − xc )2 + (y(t) − yc )2 φ(x1 , x2 , t) = (5.1) exp − σ2 + 4σ1 t σ2 + 4σ1 t where x = x1 cos t + x2 sin t, y = −x1 sin t + x2 cos t,

21

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

(xc , yc ) = (0.25, 0), σ1 = 0.001, σ2 = 0.01. We solve this problem by using several pure Lagrangian methods. More precisely, let us denote by (LG)1 the method which arises from the Lagrangian weak problem (2.12), by approximating the material derivative at t = tn+1 by a first order backward formula and the characteristics by a first order Euler formula (see [7] for more details), combined with piecewise quadratic finite elements for space discretization. Similarly, we denote by (LG)3 the method which arises from replacing in (3.1) the second order Runge-Kutta approximation of Xe by a third order Runge-Kutta approximation. Finally, we denote by (LG)2 the second order scheme given by (3.1). We have also chosen for space discretization of problems (LG)2 and (LG)3 piecewise quadratic finite elements, that is k = 2. Moreover, we have also solved the pure convection problem (i.e. σ1 = 0) with the (LG)2 scheme. All these methods were combined with an exact quadrature formula for polynomials of degree 5 in all the terms. In Figure 5.1 we have ∞

∞

l (L2(Ω )) error curve h

2

l (H1(Ω )) error curve

t

h

n

2

10

t

n

10

0

10

0

10 −2

10

Relative error (%)

Relative error (%)

−2

−4

10

−6

10

(LG)

−8

10

(LG) −10

10

2

10

−4

10

(LG)

2

−6

10

(LG)3

3

3

3

y=C/N (LG)

10

y=C/N (LG)

−8

10

1

−12

1

y=C/N

y=C/N

−14

10

−10

1

2

10

3

4

10

5

10 10 10 N: number of time steps (∆ t=(2π)/N)

1

10

2

10

3

4

10 10 N: number of time steps (∆ t=(2π)/N)

10

Fig. 5.1. Example 1: computed l∞ (L2h (Ωtn )) (left) and l∞ (Hh1 (Ωtn )) (right) errors, in log-log scale, for σ1 = 0.001 versus the number of time steps, for a fixed spatial mesh of 133 × 133 vertices.


fixed a uniform spatial mesh of 133 × 133 vertices and shown the l∞ L2h (Ωtn ) and
l∞ Hh1 (Ωtn ) errors versus the number of time steps. These results show that, for this example, schemes (LG)2 and (LG)3 possess third-order accuracy in time and scheme (LG)1 has first-order accuracy in time. We notice that for the (LG)2 scheme we ∞

∞

l (H1(Ω )) error curve

l (L2(Ω )) error curve h

1

h

t

n

2

t

n

10

10

(LG)

(LG)

2

2

y=C/N2

y=C/N3

x

x

(LG)

1

0

(LG)

1

1

10 Relative error (%)

Relative error (%)

10

−1

10

−2

0

10

−1

10

10

−3

−2

10

10

2

10 N =N : number of degrees of freedom in each spatial direction (h=1/N ) x

y

x

2

10 N =N : number of degrees of freedom in each spatial direction (h=1/N ) x

y

x

Fig. 5.2. Example 1: computed l∞ (L2h (Ωtn )) (left) and l∞ (Hh1 (Ωtn )) (right) errors, in log-log scale, for σ1 = 0.001 versus 1/h, for ∆t = 2π/2000.

´ M. BEN´ITEZ AND A. BERMUDEZ

22 ∞

Error curves for the (LG) scheme, for ∆ t=(2π)/2000

l (L2(Ω )) error curve h

0

t

2

2

n

10

10

(LG)

2 2

1

y=C/N

−1

10

10

0

Relative error (%)

Relative error (%)

−2

10

−3

10

−4

10

10

−1

10

∞

1

l (Hh(Ω )) norm t

−2

10

n

2 x

y=C/N ∞

l

−3

−5

10

10

(L2(Ω h t

)) norm

n

3

y=C/Nx −4

−6

10

1

2

10

3

4

10 10 N: number of time steps (∆ t=(2π)/N)

10

2

10 N =N : number of degrees of freedom in each spatial direction (h=1/N )

10

x

y

x

Fig. 5.3. Example 1: computed errors for the (LG)2 scheme, in log-log scale, for σ1 = 0. On the left, the l∞ (L2h (Ωtn )) error versus the number of time steps, for a fixed spatial mesh of 265 × 265 vertices. On the right, the l∞ (L2h (Ωtn )) and l∞ (Hh1 (Ωtn )) errors versus 1/h, for ∆t = 2π/2000.

z

obtain a greater order than the one predicted by the theory (second-order). We can observe, for fixed h, that the error curves become horizontal as time step decreases below a threshold; this is because the term O(h2 )
dominates the global error. In
Figure 5.2 we represent, the computed l∞ L2h (Ωtn ) and l∞ Hh1 (Ωtn ) errors versus 1/h for a fixed small time step, namely ∆t = 2π/2000. We can observe that, as predicted by Theorems 4.2 and
4.4, the (LG)2 scheme possesses second-order accuracy ∞ 1 in space in the
l Hh (Ωtn ) -norm. Moreover, third-order accuracy in space in the ∞ 2 l Lh (Ωtn ) -norm is observed. In Figure 5.3 we represent the errors, obtained with the (LG)2 scheme for the pure convection problem (σ1 = 0). On the
left, we fix a uniform spatial mesh of 265 × 265 vertices, and show the l∞ L2h (Ωtn ) errors versus
∞ 2 the number of time steps. On the right, we represent the computed l L (Ω t n ) and h
l∞ Hh1 (Ωtn ) errors versus 1/h for a fixed small time step, ∆t = 2π/2000. Notice that, for the pure convection problem, the spatial error is dominant in the total error. These results show that, as predicted in Remark
4.2, the (LG)2 scheme possesses third-order accuracy in space, in the l∞ L2h (Ωtn ) -norm. Moreover, it is remarkable that even for the pure
convection problem, second-order accuracy in space is observed in the l∞ Hh1 (Ωtn ) -norm. In Figure 5.4 we can see the exact solution compared 0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 1

0 1 0.5

1 0.5

0

0

−0.5 y

0.5

1

−1

x

0

−0.5

−0.5 −1

0.5

0

x2

−0.5 −1

−1

x

1

Fig. 5.4. Exact (left) and computed (right) solution of Example 1 with σ1 = 0.001 at time T = 2π, with the classical first order scheme and mesh parameter h = 1/132 and ∆t = 2π/400.

with the solution computed by using the classical first order characteristics method

23

0.4

0.4

0.3

0.3

0.2

0.2

z

z

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

0.1

0.1

0 1

0 1 0.5

1

1 0.5

0

0

−0.5 y

0.5

0.5

0 −1

−1

0

−0.5

−0.5 y

x

−0.5 −1

−1

x

Fig. 5.5. Exact (left) and computed (right) solution of Example 1 with σ1 = 0.001 at time T = 2π, with the second order scheme (LG 2 ), mesh parameter h = 1/132 and ∆t = 2π/100.

combined with piecewise linear finite elements. In Figure 5.5 the exact solution is compared with the numerical solution obtained by using the second order method (LG)2 proposed in the present paper. In both cases a uniform spatial mesh of 133 × 133 vertices has
been used and we have chosen the number of time step minimizing the l∞ L2 (Ωtn ) error. Clearly, (LG)2 achieves better results than the corresponding classical first order method. Notice that, for this example, the problem (4.2) from [7] can be easily solved. Thus, the analytical expression for Xe is known, namely    cos(t) −sin(t) p1 Xe (p, t) = . sin(t) cos(t) p2 Example 2 We consider a second example to compare the numerical results obtained with semiLagrangian and pure Lagrangian methods. It has a solution developing a steep layer and a velocity field which is not divergence-free. This example has been solved in [17]. The spatial domain is Ω = (0, 1) × (0, 1), T = 1, and v = ∇ψ,

A = σ1 I,

f = 0,

ρ = 1,

being ψ = (1 − cos(2πx1 ))(1 − cos(2πx2 )),

σ1 = 0.001.

The initial data varies between φ0 (0, 0) = 0 and φ0 (1, 1) = 1 according to the following expression:  0 si ξ < 0,   1 0 (5.2) φ (x1 , x2 ) = (1 − cos(πξ)) si 0 ≤ ξ ≤ 1,   2 1 si 1 < ξ,

where ξ = x1 + x2 − 1/2. Notice that the velocity field is null on the boundary so Ωt = Ω ∀t ∈ [0, 1]. We impose Dirichlet boundary conditions given by the initial data, that is φD = φ0|Γ . In Figure 5.6 we plot the velocity field and the initial data. We solve this problem with the pure Lagrangian methods (LG)1 and (LG)2

´ M. BEN´ITEZ AND A. BERMUDEZ

24

1.000e+000 7.500e-001 5.000e-001 2.500e-001 0.000e+000

Fig. 5.6. Example 2: velocity field (left) and initial date (right).

and with two second order semi-Lagrangian methods. More precisely, we denote by (SLG)12 the semi-Lagrangian scheme analogous to (LG)2 , but re-initializing the transformation to the identity at the beginning of each time step (see [7] for more details), and by (SLG)22 a two-step second order semi-Lagrangian method. The latter has been proposed and analyzed for one-dimensional convection-diffusion equations in [20], and for the incompressible Navier-Stokes equations in [14]. In all cases we have chosen for space discretization piecewise quadratic finite elements. Moreover, an exact quadrature formula for polynomials of degree 2 is used to approximate all the integrals. For the (SLG)22 scheme, we use a first-order semi-Lagrangian method to calculate the numerical solution at the first time step. 1.2

1

0.8

z

0.6

0.4

0.2 Y 1.082e+000 Z

X

0

7.908e-001 5.000e-001 2.091e-001 -8.171e-002

−0.2

0

0.2

0.4

0.6

0.8

1

x

Fig. 5.7. Example 2: numerical solution contours at T = 1 (left) and the section x1 −→ 2 φN ∆t,h (x1 , 1/2) (right) for (SLG)2 semi-Lagrangian scheme, h = 1/16, ∆t = 1/60.

In Figures 5.7, 5.8, 5.9 and 5.10 we represent the numerical solution contours at final time T = 1 and the section x1 −→ φN ∆t,h (x1 , 1/2), computed by using the 2 1 (SLG)2 , (SLG)2 , (LG)1 and (LG)2 methods, respectively, and h = 1/16. The semiLagrangian methods present oscillations near the transition layer, so Gibbs phenomena is observed, while the pure Lagrangian methods are accurate even in the steep layer around the diagonal. These features can be observed on the plots of the sections. This problem has been also solved in [17] with a semi-Lagrangian method combined with a discontinuous Galerkin discretization, and also with a standard Galerkin scheme. The Gibbs phenomena is also observed for both methods even for very fine meshes, with h = 1/32. The oscillations produced by the standard Galerkin scheme are observed even far from the transition layer. Finally, approximate solution con-

25

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

1.2

1

0.8

z

0.6

0.4

0.2 Y

X

Z

0

1.010e+000 7.550e-001 5.000e-001 2.451e-001 -9.892e-003

−0.2

0

0.2

0.4

0.6

0.8

1

x

Fig. 5.8. Example 2: numerical solution contours at T = 1 (left) and the section x1 −→ 1 φN ∆t,h (x1 , 1/2) (right) for the (SLG)2 scheme, h = 1/16, ∆t = 1/60. 1.2

1

0.8

z

0.6

0.4

0.2 Y

Z

X

1.000e+000 7.501e-001 5.000e-001 2.499e-001 -2.173e-004

0

−0.2

0

0.2

0.4

0.6

0.8

1

x

Fig. 5.9. Example 2: numerical solution contours at T = 1 (left) and the section x1 −→ φN ∆t,h (x1 , 1/2) (right) for the (LG)1 scheme, h = 1/16, ∆t = 1/60.

tours in Lagrangian coordinates at T = 1, φN m,∆t,h , computed with the (LG)2 scheme are plotted in Figure 5.11.

´ M. BEN´ITEZ AND A. BERMUDEZ

26

1.2

1

0.8

z

0.6

0.4

0.2 Y

Z

1.000e+000 7.501e-001 5.000e-001 2.499e-001 -1.563e-004

X

0

−0.2

0

0.2

0.4

0.6

0.8

1

x

Fig. 5.10. Example 2: numerical solution contours at T = 1 (left) and the section x1 −→ φN ∆t,h (x1 , 1/2) (right) for the (LG)2 scheme, h = 1/16, ∆t = 1/60.

Y

Z

1.000e+000 7.500e-001 5.000e-001 2.500e-001 -4.127e-005

X

Fig. 5.11. Example 2: numerical solution contours at T = 1, φN m,∆t,h , for the (LG)2 scheme, h = 1/16, ∆t = 1/60.

6. Conclusions. We have performed the numerical analysis of a second-order pure Lagrange-Galerkin method for convection-diffusion equations with degenerate diffusion tensor and non-divergence-free velocity fields. Moreover, we have considered general Dirichlet-Robin boundary conditions. The method has been introduced and analyzed by using the formalism of continuum mechanics. In a previous paper the proposed second order pure Lagrangian time discretization scheme has been rigorously introduced and analyzed for the same problem. Although our analysis considers any velocity field and use approximate characteristic curves, error estimates of order O(∆t2 ) + O(hk ) have been obtained when smooth enough data and solutions are available. These results have been proved by using some properties obtained in the previous paper. Numerical tests have been presented to confirm the predicted behavior.

REFERENCES [1] M.D. Baker, E. S¨ uli, and A.F. Ware. Stability and convergence of the spectral LagrangeGalerkin method for mixed periodic/non-periodic convection-dominated diffusion problems. IMA J. Numer. Anal., 19:637–663, 1999. [2] J. Baranger, D. Esslaoui, and A. Machmoum. Error estimate for convection problem with characteristics method. Numer. Algorithms, 21 (1999):49–56. Numerical methods for partial differential equations (Marrakech, 1998).

HIGHER ORDER PURE LAGRANGE-GALERKIN METHODS

27

[3] J. Baranger and A. Machmoum. Une norme “naturelle” pour la m´ ethode des caract´ eristiques en ´ el´ ements finis discontinus: cas 1-D. RAIRO Mod´ el. Math. Anal. Num´ er., 30:549–574, 1996. [4] J. Baranger and A. Machmoum. A “natural” norm for the method of characteristics using discontinuous finite elements: 2D and 3D case. Math. Model. Numer. Anal., 33:1223–1240, 1999. [5] M. Bause and P. Knabner. Uniform error analysis for Lagrange-Galerkin approximations of convection-dominated problems. SIAM. J. Numer. Anal., 39:1954–1984, 2002. [6] M. Ben´ıtez. M´ etodos num´ ericos para problemas de convecci´ on difusi´ on. Aplicaci´ on a la convecci´ on natural, Ph.D. Thesis, Universidade de Santiago de Compostela. 2009. [7] M. Ben´ıtez and A. Berm´ udez. Numerical Analysis of a second-order pure Lagrange-Galerkin method for convection-diffusion problems. Part I: time discretization. Submitted. [8] M. Ben´ıtez and A. Berm´ udez. A second order characteristics finite element scheme for natural convection problems. Submitted. [9] M. Bercovier, O. Pironneau, and V. Sastri. Finite elements and characteristics for some parabolic-hyperbolic problems. Appl. Math. Modelling, 7:89–96, 1983. [10] A. Berm´ udez and J. Durany. La methode des caracteristiques pour des problemes de convectiondiffusion stationnaires. M2AN. Math. Mod. and Num. Anal., 21:7–26, 1987. [11] A. Berm´ udez and J. Durany. Numerical solution of steady-state flow through a porous dam. Comput. Methods Appl. Mech. Engrg., 68:55–65, 1988. [12] A. Berm´ udez, M.R. Nogueiras, and C. V´ azquez. Numerical solution of (degenerated) convection-diffusion-reaction problems with higher order characteristics/finite elements. Part I: Time discretization. SIAM. J. Numer. Anal., 44:1829–1853, 2006. [13] A. Berm´ udez, M.R. Nogueiras, and C. V´ azquez. Numerical solution of (degenerated) convection-diffusion-reaction problems with higher order characteristics/finite elements. Part II: Fully Discretized Scheme and Quadrature Formulas. SIAM. J. Numer. Anal., 44:1854–1876, 2006. [14] K. Boukir, Y. Maday, B. M´ etivet, and E. Razafindrakoto. A high-order characteristics/finite element method for the incompressible Navier-Stokes equations. Internat. J. Numer. Methods Fluids, 25:1421–1454, 1997. [15] K. Chrysafinos and N. J. Walkington. Error estimates for discontinuous Galerkin approximations of implicit parabolic equations. SIAM J. Numer. Anal., 43:2478–2499 (electronic), 2006. [16] K. Chrysafinos and N. J. Walkington. Error estimates for the discontinuous Galerkin methods for parabolic equations. SIAM J. Numer. Anal., 44:349–366 (electronic), 2006. [17] K. Chrysafinos and N. J. Walkington. Lagrangian and moving mesh methods for the convection diffusion equation. M2AN Math. Model. Numer. Anal., 42:25–55, 2008. [18] J. Douglas, Jr., C.-S. Huang, and F. Pereira. The modified method of characteristics with adjusted advection. Numer. Math., 83:353–369, 1999. [19] J. Douglas, Jr., and T.F. Russell. Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J. Numer. Anal., 19:871–885, 1982. [20] R.E. Ewing and T.F. Russel. Multistep galerkin methods along characteristics for convectiondiffusion problems. IMACS Publications. Advances in Computer Methods for Partial Differential Equations IV, pages 28–36, 1981. [21] R.E. Ewing and H. Wang. A summary of numerical methods for time-dependent advectiondominated partial differential equations. J. Comput. Appl. Math., 128:423–445, 2001. [22] G. Fourestey. Stabilit´ e des m´ ethodes de Lagrange-Galerkin du premier et du second ordre. tech. report, INRIA, Raport de recherche, 2002. [23] M.E. Gurtin. An Introduction to Continuum Mechanics. Mathematics in Science and Engineering, 158, Academic Press, San Diego, 1981. [24] K.W. Morton, A. Priestley, and E. S¨ uli. Stability of the Lagrange-Galerkin Method with nonexac integration. M2AN Math. Model. Numer. Anal., 22:625–653, 1998. [25] C. Par´ es. ´ estude math´ ematique et approximation num´ erique de quelques plobl` emes aux limites de la m´ ecanique des fluides incompressibles. Th` ese. Universit´ e Paris VI. 1992. [26] O. Pironneau. On the Transport-Diffusion Algorithm and Its Applications to the Navier-Stokes Equations. Numer. Math., 38:309–332, 1982. [27] O. Pironneau, J. Liou, and T. Tezduyar. Characteristic-Galerkin and Galerkin/least-squares space-time formulations for the advection-diffusion equations with time-dependent domains. Comput. Methods Appl. Mech. Engrg., 100:117–141, 1992. [28] A. Priestley. Exact projections and the Lagrange-Galerkin method: a realistic alternative to cuadrature. J. Comput. Phys., 112:316–333, 1994.

28

´ M. BEN´ITEZ AND A. BERMUDEZ

[29] A. Quarteroni and A. Valli. Numerical approximation of partial differential equations. Springer Series in Computational Mathematics, Springer-Verlag, 23, Berlin, 1994. [30] H. Rui and M. Tabata. A second order characteristic finite element scheme for convectiondiffusion problems. Numer. Math., 92:161–177, 2002. [31] E. S¨ uli. Convergence and nonlinear stability of the Lagrange-Galerkin method for the NavierStokes equations. Numer. Math., 53:459–483, 1988. [32] E. S¨ uli. Stability and convergence of the Lagrange-Galerkin method with non-exac integration. Academic Press, London. The mathematics of finite elements and applications, VI, pages 435–442, 1988. [33] E. S¨ uli and A. Ware. A spectral mehod of characteristics for hyperbolic problems. SIAM J. Numer. Anal., 28:423–445, 1991. [34] M. Tabata and S. Fujima. Robustness of a characteristic finite element scheme of second order in time increment. tech. report, MHF Preprint Series, 2004.