Convergence of locally divergence-free discontinuous-Galerkin ...

2 downloads 0 Views 325KB Size Report
Proof. See [2]. D. Now we define the local L2 projection operator πh from S m ... Rk is entirely determined in a triangle K by its two sets of moments: Me(v), MK(v).
ESAIM: Mathematical Modelling and Numerical Analysis

ESAIM: M2AN Vol. 39, No 6, 2005, pp. 1177–1202 DOI: 10.1051/m2an:2005051

CONVERGENCE OF LOCALLY DIVERGENCE-FREE DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS OF THE 2D-MHD SYSTEM ∗

¨ ner 2 Nicolas Besse 1 and Dietmar Kro Abstract. We present the convergence analysis of locally divergence-free discontinuous Galerkin methods for the induction equations which appear in the ideal magnetohydrodynamic system. When we use a second order Runge Kutta time discretization, under the CFL condition ∆t ∼ h4/3 , we obtain error estimates in L2 of order O(∆t2 + hm+1/2 ) where m is the degree of the local polynomials. Mathematics Subject Classification. 76W05, 65J10. Received: January 17, 2005. Revised: July 4, 2005.

1. Introduction The magnetohydrodynamic (MHD) equations modelize electrically conducting fluid flow in which the electromagnetic forces can be of the same order or even greater than hydrodynamic ones. The ideal MHD system which combines the equations of gas dynamics with Maxwell equations in which relativistic, viscous and resistive effects are neglected can be written in the following three dimensional conservative form ∂t ρ + ∇ · (ρu) = 0     1 2 ∂t (ρu) + ∇ · ρu ⊗ u + p + 2 |B| I − B ⊗ B = 0 ∂t (ρe) + ∇ ·

∂t B + ∇ · (u ⊗ B − B ⊗ u) = 0   ρe + p + 12 |B|2 u − B(u · B) = 0



∇·B

= 0

(mass conservation) (momentum conservation) (induction equation) (energy conservation) (divergence constraint)

where ρ is the density, u the velocity field, B the magnetic field, p the pressure, e the total energy and I identity matrix. If the initial data are divergence free, that is to say that ∇ · B0 = 0, then an exact solution will satisfy this constraint for all the times. For smooth solution this is obvious because the induction equation can be Keywords and phrases. Magnetohydrodynamics, discontinuous-Galerkin methods, convergence analysis. ∗ The authors acknowledge financial support from the HYKE project No HPRN-CT-2002-00282 on “Hyperbolic and Kinetic Equations: Asymptotics, Numerics, Applications”, funded by European Union. 1 Institut de Recherche Mathematique Avanc´ ee, Universit´e Louis Pasteur, CNRS, 7 rue Ren´e Descartes, 67084 Strasbourg Cedex, France. [email protected] 2 Institut f¨ ur Angewandte Mathematik Universit¨ at Freiburg, Hermann-Herder-Str. 10, 79104 Freiburg i. Br., Germany. [email protected] c EDP Sciences, SMAI 2005 

Article published by EDP Sciences and available at http://www.edpsciences.org/m2an or http://dx.doi.org/10.1051/m2an:2005051

1178

¨ N. BESSE AND D. KRONER

rewritten as ∂t B + curl(B × u) = 0 and we have the identity ∇ · (∇ × ·) = 0. The numerical preservation of the ∇ · B = 0 condition is an important and much debated problem for numerical MHD codes [3, 4, 10–12, 20, 22]. Because ∇ · B errors arise in numerical simulations and may increase in time, numerical instabilities can appear and lead to unphysical behaviour of the system. For example numerically incorrect magnetic field topologies lead to unphysical plasma transport orthogonal to the magnetic field. The non enforcement of the ∇ · B constraint leads to the loss of momentum and energy conservation and allows fictitious forces to develop parallel to the magnetic field. These effects are discussed in [3, 4]. In this paper we present the convergence analysis and error estimates of locally divergence-free Runge Kutta discontinuous Galerkin schemes for smooth solutions of the two dimensional induction equation which arises in the MHD system ∂B + curl(B × u) = 0 (1) ∂t where B = (Bx (t, x), By (t, x)) is the magnetic field, and u = (ux (t, x), uy (t, x)) is the velocity field, with the notations x = (x, y). We assume for this paper that the velocity field u is given. The construction and the convergence analysis of this discontinuous Galerkin method rest on three ingredients. First we rewrite the induction equation as a Friedrichs system. Then we write a discontinuous Galerkin formulation of this new system in which we choose upwind flux as definition of the flux at the cell interface. In order to have a locally divergence free scheme we use piecewise solenoidal functions that are totally discontinuous across interface cells but which are pointwise divergence free on each element. This basis functions were developed in [2, 16] in the context of nonconforming finite element method for the stationary Navier-Stokes equations and were also used in [17]. We also use the Nedelec finite element in H(curl) [19] in two dimensions which is obtained by rotating the Raviart-Thomas element [21] by π/2. Then we used the framework developed in [6, 24] to show the convergence and obtained error estimates for our scheme. In [7] the authors develop a locally divergence free discontinuous Galerkin scheme for numerically solving the Maxwell equation. They also show an error estimate of the form O(hk+1/2 ) where k is the degree of the local polynomials for a scheme semi-discretized in space only.

2. The numerical scheme 2.1. The Friedrichs formulation Thanks to the divergence constraint, ∇ · B = 0, the induction equation (1) can be rewritten as the following Friedrichs system: ∂B ∂(Ax B) ∂(Ay B) + + + CB = 0 in Ω × [0, T ], (2) ∂t ∂x ∂y where     ux (t, x) uy (t, x) 0 0 Ax (t, x) = , Ay (t, x) = , 0 ux (t, x) 0 uy (t, x) and   ∂x ux (t, x) ∂y ux (t, x) C(t, x) = − . ∂x uy (t, x) ∂y uy (t, x) Sometimes we will use the notation (v 1 , v 2 ) instead of (v x , v y ) where v is a vector or matrix fields. For suitable boundary conditions on ∂Ω Friedrichs [14] has proved the existence and uniqueness in L2 (Ω) of a weak and strong solution to (2), under appropriate smoothness conditions and the additional assumption C + CT +

∂Ay ∂Ax + ≥ α I, ∂x ∂y

in Ω,

(3)

with α a strictly positive constant. We suppose that u ∈ L∞ (0, T ; C ∞ (Ω)), ∇ · u ∈ L∞ (0, T ; L∞ (Ω)) and C ∈ L∞ (0, T ; C ∞(Ω)). We recall that if B0 ∈ Hs ( 2 ) with s ∈ , then the system (2) admits a unique weak  0 s 2 solution B ∈ C [0, +∞[; H ( ) ∩ C 1 [0, +∞[; Hs−1 ( 2 ) (see [1]). In particular, if s > 2 the weak solution  lies in C 1 [0, +∞[× 2 and is in fact a classical solution.

Ê

Ê

Ê

Ê

Ê

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1179

 x) = e−λt B(t, x), with λ > 0 instead of B. Remark 1. The condition (3) can be easily satisfied by using B(t,  satisfies Then B    ∂B ∂(Ax B) ∂(Ay B)  = 0 in Ω × [0, T ], + + + CB (4) ∂t ∂x ∂y where C = C + λ I. Since we assume ∇u ∈ L∞ (Ω) we can choose λ sufficiently large such that condition (3) is satisfied for (4).

2.2. The approximation spaces First we introduce the space ∞ (0, T ; X) defined by  ∞  (0, T ; X) := f : {t0 , ..., tNT } → X| ||f ||∞ (0,T ;X) =

n

max ||f (t )||X

1≤n≤NT

 1 are not divergence free spaces. Now we present a trace result and a continuous embedding which will be useful later. Let K a convex polygon with Lipschitz boundary in 2 . First we define the space V(K) and W(K) by W(K) = {v ∈ H(curl; K)∩ H(div; K), v × ν ∈ L2 (∂K)} equipped with the norm v W(K) = v L2 (K) + curl v L2 (K) + div v L2 (K) +  v × ν L2 (∂K) and V(K) = {v ∈ H(curl; K) ∩ H(div; K), v · ν ∈ L2 (∂K) equipped with the norm v V(K) = v L2 (K) + curl v L2 (K) + div v L2 (K) + v · ν L2 (∂K) . It is proved in [8,9] that V(K) = W(K) and then that the mapping v −→ γ0 v defined on D(K) has a unique linear continuous extension as an operator from W(K) onto L2 (∂K), where γ0 v is the boundary values of v on ∂K. Now let us define H 1/2 (e) ∀e ∈ ∂K equipped with the norm Q(K) = {v ∈ H(curl; K) ∩ H(div; K), v × ν e ∈ v Q(K) = v L2 (K) + curl v L2 (K) + div v L2 (K) + e∈∂K v × ν e H 1/2 (e) . Then we have the following theorem where a proof can be found in [13]

Ê

Theorem 3. Let K be a convex polygon in

Ê2 with Lipschitz boundary. Then we have the continuous embedding: Q(K) → H1 (K).

(13)

2.3. The discontinuous Galerkin method In this section we describe our discontinuous Galerkin method. If we take the scalar product of the equation (1) with a test function ϕ, integrate the scalar product on a cell K and use the following Green formula

K

v∂i wdx = −

K



w∂i vdx +

e∈∂K

e

i wvνe,K dΓ

(14)

y 1 2 x where ν e,K = (νe,K , νe,K )T = (νe,K , νe,K )T denotes the outward unit normal to the face e of K, we obtain

 ∂t

K

 

2 2   i B · ϕdx − Ai B · ∂i ϕdx + Ai B · ϕνe,K dΓ + CB · ϕdx = 0. i=1

K

i=1 e∈∂K

e

(15)

K

Now we replace respectively B and ϕ by Bh and ϕh in (15) where Bh , ϕh ∈ Vhk (Ω). Nevertheless the terms arising from the boundary of the cell K in (15) are not well defined or have no sense since Bh and ϕh are discontinuous on the boundary ∂K of the element K. Then we replace these terms by a numerical 2 i and upwind flux that we are going to define in the following lines. Let us define Ae,K (t, x) = i=1 (Ai )|e νe,K − + − + Ce,K (t, x) = −(Ae,K (t, x)) , De,K (t, x) = (Ae,K (t, x)) where A and A denote respectively the negative and the positive part of A. Then we define the upwind numerical flux g(νe,K , v, w) = −Ce,K w + De,K v. By noting that |A| = A+ − A− and A = A+ + A− we can rewrite g(νe,K , v, w) as the Lax-Friedrichs flux − |Ae,K | (w−v) , where v (resp. w) denotes the interior (resp. exterior) value of the g(νe,K , v, w) = Ae,K (v+w) 2 2 solution from the cell K. Now we obtain the following semi-discretized scheme in space  ∂t

K

 

2 2   Bh · ϕh dx − Ai Bh · ∂i ϕh dx + g(ν e,K , BK , BKe ) · ϕK dΓ + CBh · ϕh dx = 0 i=1

K

i=1 e∈∂K

e

K

for all ϕh ∈ Vhk where Ke is the neighbourhood of K along the face e. Let T be the final time and ∆t = T /NT the time step. Now we use a second order Runge-Kutta scheme for the discretization in time. Then the fully

1182

¨ N. BESSE AND D. KRONER

discretized scheme reads

K

K

· ϕh dx =

Bn+1 · ϕh dx = h

where n FK (W, ϕh )

Yhn

=

2  i=1

Ani Wn

K

n Bnh · ϕh dx + ∆tFK (Bnh , ϕh ) 

 1 n 1 n ∆t n+1 n Bh + Yh · ϕh dx + F (Yh , ϕh ) 2 2 2 K K

(16)

K

· ∂i ϕh dx +

2   i=1 e∈∂K

with

n

e

g (W)e · ϕK dΓ −

K

(17)

C n Wn · ϕh dx

n n n n (WK (WK + WK ) − WK ) e e − |Ane,K | 2 2 or W k . If we expand Bh and ϕh , that is to say

n n g n (W)e = g n (νe,K , WK , WK ) = Ane,K e

and Ane,K = Ae,K (tn , x), C n = C(tn , x). Let X k be V k Bnh (x)

=



dim(X k )

½

BnK (x) K (x),

K∈Th

BnK (x)

=

 i=1

i,n σK Ψi (x),

dim(X k ) n YK (x)

=

 i=1

(18)

i,n θK Ψi (x)

with Ψi ∈ X k , then we obtain the scheme  n  = [ΣK ]n + ∆tLnK ([Σh ]n )  [ΘK ] 1 n+1 n  = ([ΣK ]n + [ΘK ]n ) + ∆tLn+1  [ΣK ] K ([Θh ] ) 2 T T   k dim (X k ),n 1,n where [ΣK ]n = σK , . . . , σK , [Σh ]n = σ 1,n , . . . , σ card(Th )dim(X ),n and for all i ∈ {1, ..., dim(X k )} [LnK

n

([Σh ] )]i =

2 dim(X  

k

)

l=1 m,j=1

+

2   i=1 e∈∂K



m,j=1

with (MK )ij =

K

K

  Anl (x)Ψj (x) · M−T K m,i ∂l Ψm (x)dx

dim(X k )



e m,j=1

k

dim(X )



j,n σK

j,n σK

K

  j,n n j,n n [σK De,K (x(Γ)) − σK Ce,K (x(Γ))]Ψj (x(Γ)) · M−T K m,i Ψm (x(Γ))dΓ e

  C n (x)Ψj (x) · M−T K m,i Ψm (x)dx

Ψi (x) · Ψj (x)dx.

3. Analysis of the semi-discretized scheme in space In this section we show the L2 -stability, the convergence and present some error estimates for the semidiscretized scheme in space, continuous in time.   Theorem 4. Let u and B0 sufficiently, regular, typically we consider that u ∈ L∞ [0, +∞[; W 1,∞ ( 2 ) and B0 ∈ Hm+1 ( 2 ) then there exists a constant C = C( ∇u L∞ ([0,T ]×Ω) , T ) independent of h such that

Ê

Ê

Bh L∞ (0,T ;L2 (Ω)) ≤ C Bh (0) L2 (Ω) ,

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1183

where Bh (0) ∈ Vhm or Bh (0) ∈ Whm with m = 1. Let [Bh ]e be the jump of Bh across an edge e of Eh , the set of all the edges of the partition Th then, there exists a constant C = C( ∇u L∞ ([0,T ]×Ω) , T, B L∞ (0,T ;Hm+1 (Ê2 )) ) independent of h such that   T   |Ae |[Bh ]e · [Bh ]e dΓdt ≤ Chm+η B − Bh L∞ (0,T ;L2 (Ω)) +  0

where η =

1 2

e

e∈Eh

if Bh ∈ Vhm and η = − 21 if Bh ∈ Whm with m = 1.

Proof. The stability result is proved in the Appendix 5. We will show simultaneously the convergence of the scheme and derivate some error estimates. First we begin to derivate estimates which are valid as well for Bh ∈ Vhm than Bh ∈ Wh1 . Then we continue the proof by distinguishing the case where Bh ∈ Vhm and the case where Bh ∈ Wh1 . First by noting that Lh (Bh , ϕh ) = 0

∀ϕh ∈ Vhm (Ω) or ∀ϕh ∈ Wh1 (Ω)

and Lh (B, ϕh ) = 0

then

Lh (eh , ϕh ) = 0, ∀ϕh ∈ Vhm (Ω) or ∀ϕh ∈ Wh1 (Ω) (19) where eh = B − Bh . By noting that πh Bh = Bh and πh eh − eh = πh (B − Bh ) − (B − Bh ), then using (19) we obtain (20) Lh (πh eh , πh eh ) = Lh (πh eh − eh , πh eh ) = Lh (πh B − B, πh eh ). From the stability analysis done in Appendix 5 the left hand side of (20) is the expression

1 1 1 T  2 2 Lh (πh eh , πh eh ) = πh eh (T ) L2 (Ω) − πh eh (0) L2 (Ω) + |Ae |[πh eh ]e · [πh eh ]e dΓdt 2 2 2 0 e∈Eh e

T

1 T + Cπh eh · πh eh dxdt + ∇ · u πh eh 22 dxdt. 2 0 Ω 0 Ω

(21)

Now it remains to estimate the right hand side of (20) that is to say find an estimate for Lh (πh B − B, ϕh ),

∀ϕh ∈ Vhm (Ω) or ∀ϕh ∈ Wh1 (Ω).

Let us decompose Lh (πh B − B, ϕh ) in several terms which will be estimated in the sequel.

Lh (πh B − B, ϕh ) =

T

0

− −



∂t (πh B − B) · ϕh dxdt

2 

T

Ω i=1 0

T  0

e∈Eh T

+ 0



(22)

e

Ai (πh B − B) · ∂i ϕh dxdt

(23)

g(πh B − B)e · [ϕh ]e dΓdt

(24)

C(πh B − B) · ϕh dxdt.

(25)

Case where Bh ∈ Vhm From the definition of πh for the term (22) we have

0

T



(πh (∂t B) − ∂t B) · ϕh dxdt = 0.

(26)

1184

¨ N. BESSE AND D. KRONER

Now we define a new projection operator Ph by (Ph g)|K = PK g =

1 |K|

K

gdx, ∀g ∈ L1 (Ω) ∩ W 1,∞ (Ω),

and Ph =



PK

K∈Th

½K .

Then by Taylor expansion it can be easily shown that there exists a constant C independent of h such that Ph g − g L∞ (Ω) ≤ Ch g W 1,∞ (Ω) .

(27)

Then the i-th term of (23) can be recasted in

T

− 0

T

− 0





Ai (πh B − B) · ∂i ϕh dxdt =



T

0

(πh B − B) · Ph Ai ∂i ϕh dxdt =



0

Ω T

(Ph Ai − Ai )(πh B − B) · ∂i ϕh dxdt



(Ph Ai − Ai )(πh B − B) · ∂i ϕh dxdt

(28)

because Ai is symmetric, Ph Ai ∂i ϕh = Ph ui ∂i ϕh ∈ Vhm (Ω) and because of the definition of πh . Now we give an estimate for the term (28). Thanks the approximation properties (7) and (27), the inverse properties (8), and the Cauchy-Schwarz inequality we obtain    T    (Ph Ai − Ai )(πh B − B) · ∂i ϕh dxdt   0 Ω  1/2  1/2

T   ≤ πh B − B 22 dx (Ph Ai − Ai )∂i ϕh 22 dx dt 0

K

K∈Th

K

  ≤ C T, Ai L∞ (0,T ;W 1,∞ (Ω)) , |B|L∞ (0,T ;Hm+1 (Ω)) hm+1

0

T

ϕh L2 (Ω) dt.

(29)

For the term (25) we get      T   T       C(πh B − B) · ϕh dxdt ≤  C(πh B − B) · ϕh dxdt   0 Ω   0  K∈Th K 1/2

T   πh B − B 22 dx ϕh L2 (K) dt ≤ C L∞ (0,T ;L∞ (Ω)) 0

K∈Th

K

  ≤ C T, C L∞(0,T ;L∞ (Ω)) , |B|L∞ (0,T ;Hm+1 (Ω)) hm+1

0

T

ϕh L2 (Ω) dt.

Now it remains to estimate the term (24). First have g(πh B − B)e

= = ≤

1 Ae,K (πh B − B) − |Ae,K | [(πh B − B)]e 2 Ae,K − |Ae,K | Ae,K + |Ae,K | (πh B − B)Ke + (πh B − B)K 2 2 |Ae,K | (|πh B − B|Ke + |πh B − B|K ) .

(30)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1185

Since |Ae,K | is symmetric, using a Young inequality we obtain

0

T

 e

e∈Eh

g(πh B − B)e · [ϕh ]e dΓdt ≤

T

≤ 0

0

T

 e∈Eh

e

|Ae,K |−1/2 g(πh B − B)e · |Ae,K |1/2 [ϕh ]e dΓdt

2   1 T    |Ae |[ϕh ]e · [ϕh ]e dΓdt. |Ae,K |−1/2 g(πh B − B)e  dxdt + 4 0 2 e e

e∈Eh

(31)

e∈Eh

Now we want to estimate the first term of (31). Using (9) and the trace result γ0 v L2 (∂K) ≤ C v H1/2 (K) we get

T 

T   2   2 −1/2 ∞ ∞ | g(π B − B) dxdt ≤ 2 u πh B − B L2 (∂K) dt |Ae,K h e L (0,T ;L (Ω)) 0

2

e

e∈Eh

≤ 2 u L∞(0,T ;L∞ (Ω))

T

0

 K∈Th

0

K∈Th

  2 πh B − B H1/2 (K) dt ≤ C T, u L∞ (0,T ;L∞ (Ω)) , |B|L∞ (0,T ;Hm+1 (Ω)) h2m+1 . (32)

Finally, from (20)–(26), (29)–(32) we get

1 T  − + |Ae |[πh eh ]e · [πh eh ]e dΓdt ≤ C1 h2m+1 2 0 e e∈Eh

T

T   + C2 hm+1 πh eh (t) L2 (Ω) + 2 C L∞(Ω) + ∇ · u L∞ (Ω) πh eh (t) 2L2 (Ω) dt.

πh eh (T ) 2L2 (Ω)

πh eh (0) 2L2 (Ω)

0

(33)

0

Using a Young inequality, followed by a Gronwall inequality, from (33) we get, ∀t ≤ T ,   1 πh eh (t) L2 (Ω) ≤ C πh eh (0) L2 (Ω) + hm+1/2 e(CL∞ ([0,T ]×Ω) + 2 ∇·uL∞ ([0,T ]×Ω) +C2 /2)t ≤ Chm+1/2 . Finally we get, ∀t ≤ T , eh (t) L2 (Ω)



πh eh (t) L2 (Ω) + πh eh (t) − eh (t) L2 (Ω)



πh eh (t) L2 (Ω) + πh B(t) − B(t) L2 (Ω) ≤ Chm+1/2

which finishes the proof when Bh ∈ Vhm . Case where Bh ∈ Wh1 We have to give new estimates for the terms (22)–(25). For the term (22), using Cauchy-Schwarz inequality and the approximation property (10) we get

0

T



  (Πh (∂t B) − ∂t B) · ϕh dxdt ≤ C T, |∂t B|L∞ (0,T ;H1 (Ω)) h

0

T

ϕh L2 (Ω) dt.

(34)

Now let us estimate the term (24). The inequality (31) is still valid. It remains to give a new estimate for the first term of the right hand side of the inequality (31). Using the approximation properties (10)–(12), and the

1186

¨ N. BESSE AND D. KRONER

fact that divΠh B = divB = 0 then we get

0

T

2     |Ae,K |−1/2 g(Πh B − B)e  dxdt ≤ 2 u L∞ (0,T ;L∞ (Ω)) 2

e

e∈Eh

≤ 2 u L∞ (0,T ;L∞ (Ω)) 0

≤ 2 u L∞ (0,T ;L∞ (Ω))

0



T

T

 K∈Th

Πh B − B 2L2 (∂K) dt

2

K∈Th

Πh B − B W(K) dt

   2 2 Πh B − B H(curl;K) + (Πh B − B) × ν L2 (∂K) dt

T

0

K∈Th

  ≤ C T, u L∞ (0,T ;L∞ (Ω)) , B L∞ (0,T ;H2 (Ω)) h2 .

(35)

Now we give a new estimate for (23) and (25). Using the Green formula (14) we obtain −

2 

T

0

i=1





T

=

0

0

 2   K∈Th

T



Ai (Πh B − B) · ∂i ϕh dxdt +

K

0



C(Πh B − B) · ϕh dxdt 

∂i (Ai (Πh B − B)) + C(Πh B − B)

· ϕh dxdt

(36)

i=1

 

K∈Th e∈∂K

T

e

Ae,K (Πh B − B) · ϕh dΓdt.

(37)

First we look at the term (37). Using the approximation properties (10)–(12), inverse properties (8), the fact that div Πh B = div B = 0 on K and the Cauchy-Schwarz inequality we get

T

0

 

e

K∈Th e∈∂K

Ae,K (Πh B − B) · ϕh dΓdt ≤ 2 u L∞ (0,T ;L∞ (Ω))

≤ 2 u L∞ (0,T ;L∞ (Ω))



T

0

≤ 2 u L∞ (0,T ;L∞ (Ω))

K∈Th



T

0

K∈Th

0

T

 K∈Th

Πh B − B L2 (∂K) ϕh L2 (∂K) dt

h−1/2 Πh B − B W(K) ϕh L2 (K) dt   h−1/2 Πh B − B H(curl;K) + (Πh B − B) × ν L2 (∂K) ϕh L2 (K) dt





≤ C T, u L∞(0,T ;L∞ (Ω)) , B L∞ (0,T ;H2 (Ω)) h

1/2 0

T

ϕh L2 (Ω) dt.

(38)

As div Πh B = div B = 0 on K the term (36) can be recasted in

 2  

T

0

K∈Th T

= 0

=

0

T

K



K∈Th

K

K∈Th

K



 ∂i (Ai (Πh B − B)) + C(Πh B − B)

· ϕh dxdt

i=1

∇ × ((Πh B − B) × u) · ϕh dxdt {((Πh B − B) · ∇) u + (Πh B − B)div u + (u · ∇)(Πh B − B)} · ϕh dxdt.

(39)

1187

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

Now we have to estimates the three terms of (39). For the two first term of (39) we get



T

0

K

K∈Th

((Πh B − B) · ∇) u · ϕh dxdt ≤ 2 u L∞ (0,T ;W 1,∞ (Ω))

0



T

K∈Th

  ≤ C T, u L∞ (0,T ;W 1,∞ (Ω)) , |B|L∞ (0,T ;H1 (Ω)) h

T

0

Πh B − B L2 (K) ϕh L2 (K) dt ϕh L2 (Ω) dt

(40)

and



T

0

K∈Th

K

(Πh B − B)div u · ϕh dxdt ≤ divu L∞ (0,T ;L∞ (Ω))

T

0

 K∈Th

  ≤ C T, divu L∞ (0,T ;L∞ (Ω)) , |B|L∞ (0,T ;H1 (Ω)) h

0

Πh B − B L2 (K) ϕh L2 (K) dt

T

ϕh L2 (Ω) dt.

(41)

Using the approximation properties (10)–(12), the continuous imbedding (13), the fact that div Πh B = div B = 0 on K and the Cauchy-Schwarz inequality then we get

0

T

 K∈Th

K

(u · ∇)(Πh B − B) · ϕh dxdt

≤ 2 u L∞(0,T ;L∞ (Ω))

T 0

≤ 2 u L∞(0,T ;L∞ (Ω))

 K∈Th

T 0

 K∈Th

Πh B − B H1 (K) ϕh L2 (K) dt  Πh B − B H(curl;K) +

  ≤ C T, u L∞(0,T ;L∞ (Ω)) , B L∞ (0,T ;H2 (Ω)) h1/2

0

T

 e∈∂K

 (Πh B − B) × ν e H1/2 (e)

ϕh L2 (K) dt

ϕh L2 (Ω) dt.

(42)

Finally, by introducing the estimates (34)–(42) into the relations (20)–(25) and using a Young inequality followed by a Gronwall inequality we get Πh eh (t) L2 (Ω)

  1 ≤ C Πh eh (0) L2 (Ω) + h1/2 e(CL∞ ([0,T ]×Ω) + 2 ∇·uL∞ ([0,T ]×Ω) +C2 /2)t ≤ Ch1/2 ,

and ∀t ≤ T , we obtain eh (t) L2 (Ω)

≤ Πh eh (t) L2 (Ω) + Πh eh (t) − eh (t) L2 (Ω) , ≤ Ch1/2

which ends the proof when Bh ∈ Wh1 .



Remark 3. When Bh ∈ Wh1 we can not expect to obtain error estimate like O(h3/2 ) which is got when Bh ∈ Vh1 . The reason is that the space R1 does not reproduce all the polynomial fields of 1 which are divergence free. Then the approximation properties of R1 are the same as the ones of 0 . That’s the reason why we loose one order of convergence. Nevertheless we could loose one order more if there were not the approximation property (11). As the approximation in the space H(curl) and the space L2 is of the same order we can get convergence of the scheme whose the convergence rate is in O(h1/2 ), and that’s why the N´ed´elec space is interesting . Unfortunately it could not be generalized to the space Rk with k > 1 as div Rk = 0. We observe that 0 ⊂ R1 ⊂ 1 and that Wh1 ⊂ Vh1 with dim(Wh1 ) = 3 and dim(Vh1 ) = 5.









1188

¨ N. BESSE AND D. KRONER

4. Analysis of the fully discretized scheme In this section we present the stability, convergence and error analysis in L2 of the fully discretized scheme. The main result is   Theorem 5. Let u and B0 sufficiently regular, typically we consider that u ∈ W 1,∞ [0, +∞[×( 2 ) and B0 ∈ Hm+1 ( 2 ) and let us assume that there exists a constant β(α) depending on α such that the CFL condition 0 ∆t ≤ β(α)h4/3 holds. Moreover we assume the condition (3) and δB L2 (Ω) = B0 − πh B0 L2 (Ω) = O(hm+1 ). Then there exists a constant C = C( u W 1,∞ ([0,T ]×Ω) , T, α) independent of h such that

Ê

Ê

Bh l∞ (0,T ;L2 (Ω)) ≤ C Bh (0) L2 (Ω) ,  where Bh (0) ∈ Vhm or Bh (0) ∈ Whm with m = 1. Moreover there exists a constant C = C T, α, u W 1,∞ ([0,T ]×Ω) , B W 3,∞ (0,T ;Hm+1 (Ê2 )) independent of h such that   NT       |Ane |[Bnh − B(tn )]e · [Bnh − B(tn )]e dΓ ≤ C ∆t2 + hm+η B − Bh l∞ (0,T ;L2 (Ω)) + ∆t n=0 e∈Eh

where η =

1 2

e

if Bh ∈ Vhm and η = − 21 if Bh ∈ Whm with m = 1.

Remark 4. The CFL condition ∆t ≤ β(α)h4/3 is surprising, but it seems that it is linked to the time discretization with a Runge-Kutta scheme of order two. Indeed in [23] the author also obtains the convergence of a second order Runge-Kutta finite element scheme under the same CFL condition. The minimal regularity assumptions to obtain convergence are u ∈ W 1,∞ and B0 ∈ H1 when Bh ∈ Vh0 or B0 ∈ H1 (curl) when Bh ∈ Wh1 Proof. L2 -stability. We begin by proving the theorem when Bh ∈ Vhm and therefore we will modify some estimates to adapt the proof when Bh ∈ Wh1 . First we begin by proving the L2 -stability of the scheme. This proof remains true for the two choices of approximation space. If we take ϕh = Bnh in (16), and ϕh = Yhn in (17) then after summing over all elements K of the partition Th the sum (16)/2+(17) gives   2L2 (Ω) − Bnh 2L2 (Ω) − Bn+1 − Yhn 2L2 (Ω) = ∆t F n (Bnh , Bnh ) + F n+1 (Yhn , Yhn ) . Bn+1 h h

(43)

 n n where F n = K∈Th FK K . First we will give an estimate of terms of the form of F (ϕh , ϕh ) and in the second − Yhn L2 (Ω) . Using (129) and the Green formula (14) we moment we will give an estimate for the term Bn+1 h get

1  1 n n n F (ϕh , ϕh ) = − |Ae,K |[ϕh ]e · [ϕh ]e dΓ − C ϕh · ϕh dxdt − ∇ · un ϕh 22 dx. (44) 2 2 Ω e Ω

½

e∈Eh

From the condition (3), and the equations (44) and (43) we get Bn+1 2L2 (Ω) h

  ∆t ∆t ≤ 1−α − Yhn 2L2 (Ω) . Bnh 2L2 (Ω) − α Yhn 2L2 (Ω) + Bn+1 h 2 2

(45)

It remains to estimate the term Bn+1 − Yhn L2 (Ω) . To do this we remark that after summing over all elements h K of the partition Th the sum (17)−(16)/2 gives



(Bn+1 − Yhn ) · ϕh dx = h

 ∆t  n+1 n F (Yh , ϕh ) − F n (Bnh , ϕh ) . 2

(46)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1189

Let us estimate the right hand side of (46). We have 2      ∆t  n+1 n (Yh , ϕh ) − F n (Bnh , ϕh ) = ∆t An+1 Yhn − Ani Bnh · ∂i ϕh dx F i 2 K∈Th i=1 K

       − ∆t g n+1 (Yhn )e − g n (Bnh )e · ϕK dx − ∆t C n+1 Yhn − C n Bnh · ∂i ϕh dx. K∈Th e∈∂K

e

K∈Th

(47)

K

 † n+1 n By seeing that there exists a time t and t† such that Ane,K = An+1 e,K −∆t∂t Ae,K and |Ae,K | = |Ae,K |−∆t∂t |Ae,K |  n  † n †  † where t = t + θ ∆t and t = t + θ ∆t, with 0 ≤ θ ≤ 1 and 0 ≤ θ ≤ 1, then we get

g n+1 (Yhn )e − g n (Bnh )e

n

where g  (Bnh )e = ∂t Ae,K Bh −

n

= An+1 e,K Y h −

|An+1 e,K |

n

[Yhn ]e − Ane,K Bh +

2 = −g n+1 (Bnh − Yhn )e + ∆tg  (Bnh )e

∂t |A†e,K | [Bnh ]e . 2

|Ane,K | n [Bh ]e 2

Then (47) can be rewritten as

 ∆t ∆t  n+1 n (Yh , ϕh ) − F n (Bnh , ϕh ) = − F n+1 (Bnh − Yhn , ϕh ) + G(Bnh , ϕh ) F 2 2 where

G(Bnh , ϕh )

= ∆t

2

2  i=1



∂t A‡i i Bnh

· ∂i ϕh dx − ∆t

2

  K∈Th e∈∂K

g e



(Bnh )e

· ϕK dΓ − ∆t

2 Ω

∂t C  Bnh · ϕh dx (48)

with t‡i = tn + θ‡i ∆t and t = tn + θ ∆t such that Ani = An+1 − ∆t∂t A‡i i and C n = C n+1 − ∆t∂t C  . First we i estimate the term G(Bnh , ϕh ). From the Young and Cauchy-Schwarz inequalities and the inverse properties (8) we get  

   2   ‡i n ∂t Ai Bh · ∂i ϕh dx ∆t   K



K∈Th

    ∆t4  ε 2   2 n 2 C(ε) 2 ∂t A‡i i  ∞ B + ∂ ϕ h i h 2 dx h 2 h 3 L ([0,T ]×Ê2 ) K

K∈Th

≤ C(ε)

 

   2    n ∂t C Bh · ϕh dx ∆t   K K∈Th



∆t4 ε 2 2 Bnh L2 (Ω) + ϕh L2 (Ω) , h2 3

  ε 2 2 C ∂t C  L∞ ([0,T ]×Ê2 ) , ε ∆t4 Bnh L2 (Ω) + ϕh L2 (Ω) , 3

(49)

(50)

1190

¨ N. BESSE AND D. KRONER

and  

    2   g  (Bnh )e · ϕK dΓ ∆t   K∈Th e∈∂K e    

  ∂t Ae,K + ∂t |A†e,K | n ∂t Ae,K − ∂t |A†e,K | n  2    ≤ ∆t BK e + BK · ϕK dΓ   2 2 K∈Th e∈∂K e

     n   BK  ϕK + BnK ϕK dΓ ≤ C ∂t Ae,K L∞ ([0,T ]×Ê2 ) ∆t2 2 2 2 e 2 ≤ C(ε)

∆t h

K∈Th e∈∂K



4

K∈Th

(51)

e

ε  ∆t4 ε 2 2 2 2 Bnh L2 (∂K) + h ϕh L2 (∂K) ≤ C(ε) 2 Bnh L2 (Ω) + ϕh L2 (Ω) . 3 h 3 K∈Th

From (49)–(51), we get |G(Bnh , ϕh )|

 ≤ C(ε)

 ∆t4 2 2 4 + ∆t Bnh L2 (Ω) + ε ϕh L2 (Ω) h2

(52)

Now we will estimate the term ∆tF n+1 (Vh , ϕh ). From the Young and Cauchy-Schwarz inequalities and the inverse properties (8) we get      ∆t2 ε   An+1 Vh · ∂i ϕh dx ≤ C(ε) 2 Vh 2L2 (Ω) + ϕh 2L2 (Ω) , ∆t i   h 3 K

(53)

K∈Th

       C n+1 Vh · ϕh dx ∆t   K



K∈Th

and

C(ε)∆t2 Vh 2L2 (Ω) +

ε ϕh 2L2 (Ω) , 3

      ∆t2 ε   2 2 n+1 g (Vh )e · ϕK dΓ ≤ C(ε) 2 Vh L2 (Ω) + ϕh L2 (Ω) . ∆t   h 3 e

(54)

(55)

K∈Th e∈∂K

From (53)–(55), we get that |∆tF (Vh , ϕh )| ≤ C(ε)

∆t2 2 2 Vh L2 (Ω) + ε ϕh L2 (Ω) . h2

(56)

By taking ϕh = Bn+1 − Yhn in (46), and from (52) and (56) we obtain h     ≤ ∆tF (Bnh − Yhn , Bn+1 − Yhn ) + G(Bnh , Bn+1 − Yhn ) h h  4   n+1  ∆t ∆t2 2 n 2 n n 2 4  ≤ 2ε Bh − Yh L2 (Ω) + C(ε) 2 Bh − Yh L2 (Ω) + C(ε) + ∆t Bnh L2 (Ω) . h h2  n+1 2 B − Yhn  2 h

L (Ω)

(57)

Let us estimate Bnh − Yhn L2 (Ω) . By taking ϕh = Bnh − Yhn in (16) and from the estimate (56) we get 2

Bnh − Yhn L2 (Ω) = ∆tF n (Bnh , Bnh − Yhn ) ≤ ε Bnh − Yhn L2 (Ω) + C(ε) 2

2

∆t2 2 Bnh L2 (Ω) . h2

(58)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS 2

1191

n If ε is small enough, from (58) we deduce Bnh − Yhn L2 (Ω) ≤ C(ε) ∆t h2 Bh L2 (Ω) and using (57) we get 2

2

   n+1  ∆t4 ∆t4 2 n 2 4 B − Yh L2 (Ω) ≤ C(ε) ∆t + 4 + 2 Bnh L2 (Ω) . h h h

(59)

If we suppose ∆t ≤ C(ε)h4/3 , from (59) and (45) we get Bn+1 2L2 (Ω) ≤ (1 + C∆t) Bnh 2L2 (Ω) h

and Bn+1 L2 (Ω) ≤ eCT Bh (0) L2 (Ω) . h

Convergence and error estimates in L2 . Now we will prove the convergence of the scheme and show some error estimates. First of all we look at the case B ∈ Vhm and we will change some estimates to adapt the proof when B ∈ Wh1 . We set (60) Y(t, x) = B(t, x) + ∂t B(t, x)∆t. Using Taylor expansion in time we get B(t + ∆t, x) − B(t, x) − ∂t B(t, x)

∆t ∆t − ∂t B(t + ∆t, x) = O(∆t3 ). 2 2

(61)

Using Taylor expansion and from Equations (2) and (60) we obtain ∂t B(t + ∆t, x)

=

=

− −

2  i=1 2 

∂xi (Ai (t + ∆t, x)B(t + ∆t, x)) − C(t + ∆t, x)B(t + ∆t, x)    ∂xi Ai (t + ∆t, x) B(t, x) + ∂t B(t, x)∆t + O(∆t2 )

i=1

  −C(t + ∆t, x) B(t, x) + ∂t B(t, x)∆t + O(∆t2 ) =



2 

(62)

∂xi (Ai (t + ∆t, x)Y(t, x)) − C(t + ∆t, x)Y(t, x) + O(∆t2 ).

i=1

From (60) and (61) we get B(t + ∆t, x) =

1 ∆t 1 B(t, x) + Y(t, x) + ∂t B(t + ∆t, x) + O(∆t3 ). 2 2 2

(63)

From (62) and (63) we get 1 ∆t  1 B(t, x) + Y(t, x) − ∂x (Ai (t + ∆t, x)Y(t, x)) 2 2 2 i=1 i 2

B(t + ∆t, x)

=

∆t C(t + ∆t, x)Y(t, x) + O(∆t3 ). 2

(64)

∂xi (Ani (x)B(tn , x)) − ∆tC n (x)B(tn , x)

(65)

− Finally from (60) and (64) we obtain Y(tn , x)

= B(tn , x) − ∆t

2  i=1

B(tn+1 , x)

=

2   ∆t n+1 1 1 ∆t  B(tn , x) + Y(tn , x) − C(t ∂xi An+1 (x)Y(tn , x) − , x)Y(tn , x) + O(∆t3 ). i 2 2 2 i=1 2

(66)

1192

¨ N. BESSE AND D. KRONER

If we take the scalar product of (65) and (66) with a test function ϕh , integrate over a cell K and use the Green formula (14) then we get

n Yn · ϕh dx = Bn · ϕh dx + ∆tFK (Bn , ϕh ) (67) K K  

1 n 1 n ∆t n+1 n B + Y FK (Y , ϕh ) + Bn+1 · ϕh dx = (tn , x) · ϕh dx (68) · ϕh dx + 2 2 2 K K K where (tn , x) l∞ (0,T ;L2 (Ω)) = O(∆t3 ). n n enY = Yn − Yhn , IY = πh Yn − Yn , IB to (68) we obtain

n δY · ϕh dx K

n+1 δB · ϕh dx K

where n HK (ϕh ) =

and JKn+1 (ϕh ) =

K

K

n n We set the following notations δY = πh Yn − Yhn , δB = πh Bn − Bnh , = πh Bn − Bn , enB = Bn − Bnh . If we substract (16) to (67) and (17)

n n δB · ϕh dx + HK (ϕh ) K 

 1 n 1 n 1 δY + δB = · ϕh dx + JKn+1 (ϕh ) 2 2 2 K

=

n n n n (IY − IB ) · ϕh dx + ∆tFK (Bn , ϕh ) − ∆tFK (Bnh , ϕh )

n+1 n+1 n+1 n n (2IB − IY − IB + (tn , x)) · ϕh dx + ∆tFK (Yn , ϕh ) − ∆tFK (Yhn , ϕh ).

We set

Hn =

n If we take ϕh = δB in (69) and sum (69)/2+(70) gives

 n+1 2  2 δ B



K∈Th n ϕh = δY

L (Ω)

n HK

½K ,

and



J n+1 =

K∈Th

JKn+1

½K .

(69) (70)

(71) (72) (73)

in (70), after summing over all the element K of the partition Th the

 n+1  n 2 n 2 n n − δB L2 (Ω) − δB − δY = Hn (δB ) + J n+1 (δY ). L2 (Ω)

(74)

 n+1  n n n In the following we will give estimates for the three terms Hn (δB ), J n+1 (δY ) and δB − δY . L2 (Ω) n ). From equation (60) and the approximation properties (7) of πh we obtain Estimate of Hn (δB

  n+1 n n n I − IB + IY − IB L2 (Ω) ≤ C∆thm+1 |∂t B|L∞ (0,T ;Hm+1 (Ω)) . B L2 (Ω)

(75)

n Using (75), the Cauchy-Schwarz and Young inequalities we have for the first term of Hn (δB ) the estimate

 n n n n n n n 2 (IY − IB ) · δB dx ≤ IY − IB 2 δB 2 ≤ C∆th2m+2 + ∆t δB L2 (Ω) . (76) K∈Th

K



Now we estimate the term ∆tF n (Bn − πh Bn , ϕh ). First we notice the decomposition ∆t

 K∈Th

K

Ani (Bn − πh Bn ) · ∂i ϕh dx = ∆t − ∆t

 K∈Th

K

 K∈Th

K

(Ph Ani − Ani )(πh Bn − Bn ) · ∂i ϕh dx

(πh Bn − Bn ) · Ph Ani ∂i ϕh dx = ∆t

 K∈Th

K

(Ph Ani − Ani )(πh Bn − Bn ) · ∂i ϕh dx (77)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1193

because Ani is symmetric, Ph Ani ∂i ϕh = Ph uni ∂i ϕh ∈ Vhm (Ω) and because of the definition πh . Now we give an estimate for the term (77). Thanks the approximation properties (7) and (27), the inverse properties (8), the Cauchy-Schwarz and Young inequalities we obtain 

∆t

K∈Th



  K∈Th

K

K

(Ph Ani − Ani )(πh Bn − Bn ) · ∂i ϕh dx

(78)

 ∆t Ai 2L∞ (0,T ;W 1,∞ (Ω)) πh Bn − Bn 22 + ∆th2 ∂i ϕh 22 dx ≤ Ch2m+2 ∆t + ∆t ϕh 2L2 (Ω) .

Next we have ∆t

 K∈Th

n

K

n

n

C (πh B − B ) · ϕh dx ≤

  K∈Th

≤ Ch

K

2m+2

 ∆t C 2L∞ (0,T ;L∞ (Ω)) πh Bn − Bn 22 + ∆t ϕh 22 dx

∆t + ∆t ϕh 2L2 (Ω) .

(79)

As it has been done for the continuous case (see estimates (31)–(32)) we obtain ∆t

 e∈Eh

e

g n (πh Bn − Bn )e · [ϕh ]e dΓ ≤ C∆th2m+1 +

∆t  |Ane |[ϕh ]e · [ϕh ]e dΓ. 4 e

(80)

e∈Eh

From (78)–(80) we get ∆tF n (Bn − πh Bn , ϕh ) ≤ C∆th2m+1 + 3∆t ϕh 2L2 (Ω) +

∆t  |Ane |[ϕh ]e · [ϕh ]e dΓ. 4 e

(81)

e∈Eh

n n Now we want to estimate the term ∆tF n (Bn , δB ) − ∆tF n (Bnh , δB ). First we notice that n n n n n ∆tF n (Bn , δB ) − ∆tF n (Bnh , δB ) = ∆tF n (Bn − πh Bn , δB ) + ∆tF n (δB , δB ).

(82)

From (81) we obtain n n 2 ) ≤ C∆th2m+1 + 3∆t δB L2 (Ω) + ∆tF n (Bn − πh Bn , δB

∆t  n n |Ane |[δB ]e · [δB ]e dΓ. 4 e

(83)

e∈Eh

From (3) and (44) we get n n n 2 ∆tF n (δB , δB ) ≤ −α∆t δB L2 (Ω) −

∆t  n n |Ane |[δB ]e · [δB ]e dΓ. 2 e

(84)

e∈Eh

From (71), (73), (76), (82)–(84) we get n n 2 ) ≤ C∆t(h2m+1 + h2m+2 ) + (4 − α)∆t δB L2 (Ω) − Hn (δB

∆t  n n |Ane |[δB ]e · [δB ]e dΓ. 4 e e∈Eh

(85)

1194

¨ N. BESSE AND D. KRONER

n Estimate of J n+1 (δY ). Using (75), the Cauchy-Schwarz and Young inequalities we have for the first term of n+1 n J (δY ) the estimate

 K

K∈Th

n+1 n n n (2IB − IY − IB + (tn x)) · δY dx ≤

≤ Ω



n+1 n n n 2IB − IY − IB + (tn x)) 2 δY 2 dx

  n n+1 n n n 2(IB − IB ) 2 + IB − IY 2 + (tn x)) 2 δY 2 dx

n 2 L2 (Ω) . ≤ C(ε)∆t(∆t2 + hm+1 )2 + ε∆t δY

(86)

n n Now we estimate the term ∆tF n+1 (Yn , δY ) − ∆tF n+1 (Yhn , δY ). First we notice that n n n n n ∆tF n+1 (Yn , δY ) − ∆tF n+1 (Yhn , δY ) = ∆tF n+1 (Yn − πh Yn , δY ) + ∆tF n (δY , δY ).

(87)

The same proof for the estimate (81) leads to ∆tF

n+1

n

n

(Y − πh Y , ϕh ) ≤ C(ε)∆th

2m+1

+

ε∆t ϕh 2L2 (Ω)

∆t  + |Ane |[ϕh ]e · [ϕh ]e dΓ. 4 e

(88)

e∈Eh

From (88) we obtain n

n

∆tF (Y − πh Y

n

n , δY )

≤ C(ε)∆th

2m+1

+

n 2 ε∆t δY L2 (Ω)

∆t  n n + |Ane |[δY ]e · [δY ]e dΓ. 4 e

(89)

e∈Eh

From (3) and (44) we get ∆tF

n

n n (δY , δY )



n 2 −α∆t δY L2 (Ω)

∆t  n n − |Ane |[δY ]e · [δY ]e dΓ. 2 e

(90)

e∈Eh

From (72), (73), (86), (87), (89) and (90) we get

  ∆t  n n 2 n n ) ≤ C(ε)∆t (∆t2 + hm+1 )2 + h2m+1 ) + (4ε − α)∆t δY L2 (Ω) − |Ane |[δY ]e · [δY ]e dΓ. (91) J n+1 (δY 4 e e∈Eh

n+1 n Estimate of δB − δY L2 (Ω) . First we notice that after summing over all elements K of the partition Th the sum (70)−(69)/2 gives

 1  n+1 n+1 n J (δB − δY ) · ϕh dx = (ϕh ) − Hn (ϕh ) (92) 2 Ω where

 1 1  n+1 1 n+1 n n n J (ϕh ) − Hn (ϕh ) = (2IB − IY − IB + (tn x)) · ϕh dx − (I n − IB ) · ϕh dx + Ln (ϕh ) (93) 2 2 Ω 2 Ω Y

  with Ln (ϕh ) = ∆t F n+1 (Yn , ϕh ) − F n+1 (Yhn , ϕh ) − F n (Bn , ϕh ) + F n (Bnh , ϕh ) . From (75) and (76) we have n n |(IY − IB , ϕh )Ω | ≤  n+1   n n n  2I − IY − IB + (t x), ϕh Ω  ≤ B

C∆th2m+2 + ∆t ϕh 2L2 (Ω) 2

C∆t(∆t + h

m+1 2

) +

∆t ϕh 2L2 (Ω) .

(94) (95)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1195

Now let us estimate the term Ln (ϕh ). First we notice that Ln (ϕh ) =

 ∆t F n+1 (πh Yn − Yhn , ϕh ) + F n+1 (Yn − πh Yn , ϕh ) +F n (πh Bn − Bn , ϕh ) + F n (Bnh − πh Bn , ϕh )) .

(96)

From (56) and the approximation property (7) of πh we get |∆tF n (πh Bn − Bn , ϕh )| ≤   ∆tF n+1 (Yn − πh Yn , ϕh ) ≤

ε ϕ 2 2 + C(ε)∆t2 h2m B 2L∞ (0,T ;Hm+1 (Ω)) , 4 h L (Ω) ε ϕ 2 2 + C(ε)∆t2 h2m Y 2L∞ (0,T ;Hm+1 (Ω)) . 4 h L (Ω)

(97) (98)

n n , ϕh ) − ∆tF n+1 (δB , ϕh ). First we Now let us estimate the remaining term of Ln (ϕh ), that is to say ∆tF n+1 (δY notice that n n n n n ∆tF n+1 (δY , ϕh ) − ∆tF n+1 (δB , ϕh ) = −∆tF n+1 (δB − δY , ϕh ) + G(δB , ϕh ) (99)

where G(·, ϕh ) is defined by (48). From (56) we get   ∆t2 n n 2 ∆tF n+1 (δ n − δ n , ϕh ) ≤ ε ϕh 2 2 B Y L (Ω) + C(ε) 2 δB − δY L2 (Ω) 4 h

(100)

and from (52) we get n |G(δB , ϕh )| ≤ C(ε)



 ∆t4 ε 4 n 2 + ∆t L2 (Ω) + ϕh 2L2 (Ω) . δB 2 h 4

(101)

n n 2 n n − δY L2 (Ω) which appears in (100). If we take ϕh = δB − δY in (69) we get Now let us estimate the term δB n n 2 n n δB − δY L2 (Ω) = Hn (δY − δB )

(102)

where n n n n n n n n n n Hn (δY − δB ) = (IY − IB , δY − δB )Ω + ∆tF n (Bn − πh Bn , δY − δB ) + ∆tF n (πh Bn − Bnh , δY − δB ).

(103)

n n The estimate of the three terms in the right hand side of (103) gives for Hn (δY − δB ) the inequality

  ∆t2 n 2 n n n n 2 n n 2 − δB ) ≤ ∆t δY − δB L2 (Ω) + C(ε) ∆t2 h2m + ∆th2m+2 2 δB L2 (Ω) + 2ε δY − δB L2 (Ω) . Hn (δY h

(104)

If ε and ∆t is small enough, (102) and (104) leads to n n 2 δY − δB L2 (Ω)

  ∆t2 n 2 ≤ C(ε) ∆t2 h2m + ∆th2m+2 + 2 δB L2 (Ω) . h

(105)

n+1 n By taking ϕh = δB − δY in (92) and thanks to relation (93)–(101) and (105) we get

  n+1 n 2 − δY L2 (Ω) ≤ C(ε)∆t h2m+2 + (∆t2 + hm+1 )2 + ∆th2m + ∆t3 h2m−2 + ∆t2 h2m δB   ∆t4 ∆t4 n+1 n 2 n 2 + 2(∆t + ε) δB − δY L2 (Ω) + C(ε) ∆t4 + 2 + 4 δB L2 (Ω) . h h

(106)

1196

¨ N. BESSE AND D. KRONER

If ∆t and ε is small enough and ∆t ≤ C(ε)h4/3 then (106) leads to   n+1 n 2 δB − δY L2 (Ω) ≤ C(ε)∆t h2m+2 + (∆t2 + hm+1 )2 + h2m+4/3 + h2m+8/3     n 2 n 2 + C(ε) ∆t4 + ∆t5/2 + ∆t δB L2 (Ω) ≤ C(ε)∆t (∆t2 + hm+1 )2 + h2m+4/3 + C(ε)∆t δB L2 (Ω) .

(107)

Finally the estimates (74), (85) , (91) and (107) lead to

∆t  ∆t  n+1 2 n n n n n δB L2 (Ω) ≤ − |Ae |[δB ]e · [δB ]e dΓ − |Ane |[δY ]e · [δY ]e dΓ (108) 4 4 e e e∈Eh e∈Eh   n 2 n 2 + (1 + C(α, ε)∆t) δB L2 (Ω) + ∆t(4ε − α) δY L2 (Ω) + C(ε)∆t (∆t2 + hm+1 )2 + h2m+1 + h2m+4/3 . Using (105), (108) becomes   n+1 2 n 2 δB L2 (Ω) ≤ (1 + C(α, ε)∆t) δB L2 (Ω) + C(ε)∆t (∆t2 + hm+1 )2 + h2m+1 .

(109)

A discrete Gronwall inequality and (109) enable us to obtain 2  n+1 2 0 2 δB L2 (Ω) ≤ δB L2 (Ω) eC(α,ε)T + (∆t2 + hm+1 )2 + h2m+1 ≤ C ∆t2 + hm+1/2

(110)

0 2 because δB L2 (Ω) = B0 − πh B0 2L2 (Ω) = O(h2m+2 ). Finally, using (110) we get

1/2    n 2 Bnh (x) − B(tn , x) 2L2 (Ω) ≤ πh Bn (x) − B(tn , x) 2L2 (Ω) + δB L2 (Ω) ≤ C ∆t2 + hm+1/2 . To end the proof we notice that (108) implies that

2  ∆t  n n n 2 |Ane |[δB ]e · [δB ]e dΓ ≤ (1 + C(α, ε)∆t) δB L2 (Ω) + C∆t ∆t2 + hm+1/2 . 4 e

(111)

e∈Eh

If we make the sum on the index n from 0 to NT in (111) we obtain ∆t

NT   n=0 e∈Eh

e

 2 n n |Ane |[δB ]e · [δB ]e dΓ ≤ C ∆t2 + hm+1/2 .

n n ), J n+1 (δY ) and Now we suppose that Bh ∈ Wh1 . Then we see how change the estimates for Hn (δB n+1 n n n n δB − δY L2 (Ω) . First we give a new estimate for the term ∆tF (B − Πh B , ϕh ). Using the Green formula (14) we obtain

−∆t

2  i=1



Ani (Πh Bn − Bn ) · ∂i ϕh dx + ∆t = ∆t



 2   K∈Th

K

C n (Πh Bn − Bn ) · ϕh dx 

∂i (Ani (Πh Bn

i=1

− ∆t

n

n

n

n

− B )) + ∆tC (Πh B − B )

  K∈Th e∈∂K

e

· ϕh dx

(112)

Ane,K (Πh Bn − Bn ) · ϕh dΓ.

(113)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1197

First we look at the term (113). Using the approximation properties (10)–(12), inverse properties (8), the fact that div Πh Bn = div Bn = 0 on K, the Young and the Cauchy-Schwarz inequalities then we get

∆t

  K∈Th e∈∂K

e

Ane,K (Πh Bn − Bn ) · ϕh dΓ



2 u 2L∞ (0,T ;L∞ (Ω)) ∆th−1 Πh Bn − Bn 2L2 (∂K) + ∆th ϕh 2L2 (∂K) ≤ K∈Th





2 u 2L∞ (0,T ;L∞ (Ω)) ∆th−1 Πh Bn − Bn 2W(K) + ∆t ϕh 2L2 (K)

K∈Th



  

2 u 2L∞ (0,T ;L∞ (Ω)) ∆th−1 Πh Bn − Bn 2H(curl;K) + (Πh Bn − Bn ) × ν 2L2 (∂K) + ∆t ϕh 2L2 (K)

K∈Th

  2 ≤ C T, u L∞(0,T ;L∞ (Ω)) , B L∞ (0,T ;H2 (Ω)) ∆th + ∆t ϕh L2 (Ω) .

(114)

As div Πh Bn = div Bn = 0 on K the term (112) can be recasted in

∆t

 2   K∈Th

= ∆t

K

 ∂i (Ani (Πh Bn

n

n

n

n

− B )) + C (Πh B − B )

· ϕh dx

(115)

i=1



K

K∈Th

{((Πh Bn − Bn ) · ∇) un + (Πh Bn − Bn )div un + (un · ∇)(Πh Bn − Bn )} · ϕh dx.

Now we have to estimate the three terms of (115). For the two first term of (115) we get 

∆t

K∈Th





K∈Th

K

((Πh Bn − Bn ) · ∇) un · ϕh dx

4 u 2L∞ (0,T ;W 1,∞ (Ω)) ∆t Πh Bn − Bn L2 (K) + ∆t ϕh L2 (K) 2

2



  ≤ C T, u L∞ (0,T ;W 1,∞ (Ω)) , |B|L∞ (0,T ;H1 (Ω)) h2 + ∆t ϕh 2L2 (Ω) ,

(116)

and

∆t

 K∈Th

K

(Πh Bn − Bn )div un · ϕh dx



divu 2L∞ (0,T ;L∞ (Ω)) ∆t Πh Bn − Bn 2L2 (K) + ∆t ϕh 2L2 (K) ≤ K∈Th

  2 ≤ C T, divu L∞ (0,T ;L∞ (Ω)) , |B|L∞ (0,T ;H1 (Ω)) h2 + ∆t ϕh L2 (Ω) .

(117)

1198

¨ N. BESSE AND D. KRONER

Using the approximation properties (10)–(12), the continuous imbedding (13), the fact that div Πh B = div B = 0 on K and the Cauchy-Schwarz inequality then we get ∆t

 K∈Th





K∈Th





K

(un · ∇)(Πh Bn − Bn ) · ϕh dxdt

4 u 2L∞(0,T ;L∞ (Ω)) ∆t Πh Bn − Bn H1 (K) + ∆t ϕh L2 (K) 2



 C∆t Πh Bn −

K∈Th

Bn 2H(curl;K)



+

2



(Πh Bn − Bn ) ×

e∈∂K

ν e 2H1/2 (e)

 +

∆t ϕh 2L2 (K)

  2 ≤ C T, u L∞(0,T ;L∞ (Ω)) , B L∞ (0,T ;H2 (Ω)) h∆t + ∆t ϕh L2 (Ω) .

(118)

Using (112)–(118) and looking what it has been done in the case of the continuous case for the estimates (31) and (35) we obtain n

n

n

∆tF (B − Πh B , ϕh ) ≤ C∆th +

3∆t ϕh 2L2 (Ω)

∆t  + |Ane |[ϕh ]e · [ϕh ]e dΓ. 4 e

(119)

e∈Eh

By changing the approximation properties of πh by those of Πh , using (119), and following what it has been n+1 n n n ), J n+1 (δY ) and δB − δY L2 (Ω) are done in the case where B ∈ Vhm , the new estimates for Hn (δB n n 2 Hn (δB ) ≤ C∆t(h + h2 ) + (4 − α)∆t δB L2 (Ω) −

∆t  n n |Ane |[δB ]e · [δB ]e dΓ, 4 e

(120)

e∈Eh

  ∆t  n n 2 n n J n+1 (δY ) ≤ C(ε)∆t (∆t2 + h)2 + h) + (4ε − α)∆t δY L2 (Ω) − |Ane |[δY ]e · [δY ]e dΓ 4 e

(121)

e∈Eh

and

  n+1 n 2 n 2 δB − δY L2 (Ω) ≤ C(ε)∆t (∆t2 + h)2 + h4/3 + C(ε)∆t δB L2 (Ω) .

(122)

Therefore (120)–(122) and (74) leads to  Bnh (x)

n

− B(t , x) L2 (Ω) +

∆t

NT   n=0 e∈Eh

e

1/2 n |Ane |[δB ]e

·

n [δB ]e dΓ

  ≤ C ∆t2 + h1/2 . 

5. Conclusion Here we have proven the convergence of some locally divergence-free discontinuous Galerkin methods for the induction equation. The next step is to see how these schemes behave numerically and compare them with other known schemes. For this purpose, we can couple the discontinuous Galerkin scheme for the induction equation with a finite volume scheme for the hydrodynamic equations or use a discontinuous Galerkin scheme for both where the hydrodynamic variables can be approximated in a different finite dimensional space. This work will be the matter of a future collaboration.

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1199

Appendix A. Proof of theorem 2 Thanks to the theorem 1 there exists χ ∈ Vhm (Ω) such that |v|Hm+1 (K) , v − χ L2 (K) ≤ Chm+1 K

∀K ∈ Th .

(123)

Moreover we have πh v − v 2L2 (Ω) ≤

 K∈Th

πh v − v 2L2 (K) ≤

   πK (v − χ) 2L2 (K) + πK χ − χ 2L2 (K) + χ − v 2L2 (K) .

K∈Th

2

From (123) we have v − χ L2 (K) ≤ Ch2m+2 |v|2Hm+1 (K) . From (6), (123) and Cauchy-Schwarz inequality we K obtain that πK (v − χ) 2L2 (K) ≤ πK v − χ 2L2 (K) ≤ v − χ 2L2 (K) ≤ Ch2m+2 |v|2Hm+1 (K) . Using the fact K 2

that χ satisfies (6), if we choose φ = (πK χ − χ)|K , then we have πK χ − χ L2 (K) = 0. Then we obtain πh v − v 2L2 (Ω) ≤

 K∈Th

πh v − v 2L2 (K) ≤ C

 K∈Th

h2m+2 |v|2Hm+1 (K) ≤ Cσ 2 h2m+2 |v|2Hm+1 (Ω) . K

The proof of the inverse properties (8) which are due to the assumptions of the partition Th (see [2]) can be found in [5]. Following the proof of the approximation estimate (7), using (5) and the first inverse inequality of (8) we get πh v − v Hk (K) ≤ Chm+1−k |v|Hm+1 (K) , ∀K ∈ Th , v ∈ Sm+1 (Ω) with k ≤ m + 1. From this last approximation result and the interpolation theorem of Lions and Petree [18] or theorem 1.4 in [15] we obtain approximation result (9).

Appendix B. Proof of the stability result of theorem 4 Here we prove the L2 -stability of the semi-discretized scheme in space. This proof remains true for the two choices of approximation space. Let Lh (·, ·) the bilinear form defined by

Lh (Bh , ϕh ) =

T

0



∂t Bh · ϕh dxdt −

2  i=1



T

0

K

K∈Th

T

+ 0

Ai Bh · ∂i ϕh dxdt

 

K∈Th e∈∂K

e

g(Bh )e · ϕh dΓdt +

T

0



CBh · ϕh dxdt

(124)

  with ϕh ∈ Vhm or ϕh ∈ Wh1 . If we set [ϕh ]e = ϕKe − ϕK |e then we can rewritten (124) as

Lh (Bh , ϕh ) =

0

T



∂t Bh · ϕh dxdt −

2  i=1

0

T



Ai Bh · ∂i ϕh dxdt

T

− 0

 e∈Eh

e

g(Bh )e · [ϕh ]e dΓdt +

0

T



CBh · ϕh dxdt

(125)

1200

¨ N. BESSE AND D. KRONER

where Eh denotes the set of all the faces e of the partition Th . If we take ϕh = Bh in (125) we obtain Lh (Bh , ϕh )

1 B(T ) 2L2 (Ω) − 2 2 T   −

=

K∈Th i=1





T

0

e∈Eh

K

0

e

1 B(0) 2L2 (Ω) + 2

T

0

CBh · Bh dxdt Ω

Ai BK · ∂i BK dxdt

(126)

g(Bh )e · [Bh ]e dΓdt.

(127)

From the definition (18) of the upwind flux, the term (127) can be rewritten as follows:

T 0

 e∈Eh

e

T

g(Bh )e · [Bh ]e dΓdt = 0

 where Bh,e =

BK +BKe 2

 e∈Eh

1 Ae,K Bh,e · [Bh ]e dΓdt − 2 e

T

0

 e

e∈Eh

|Ae,K |[Bh ]e · [Bh ]e dΓdt (128)

 |e

. Now let us rewrite the term (126). By using the fact that

Ai BK · ∂i BK =

1 1 ∂i (Ai BK · BK ) − ∂i Ai BK · BK 2 2

(129)

and the Green formula (14), we obtain:



2  

1 2

K

0

K∈Th i=1

=

T

T

Ai BK · ∂i BK dxdt

0

∇ · u Bh 22 dxdt −



1 2

T

0

2    K∈Th e∈∂K i=1

e

i Ai BK · BK νe,K dΓdt.

(130)

From the fact that Ai is symmetric and hence Ae,K is, and the relation Ae,Ke = −Ae,K , the second term of (130) can rewritten as 1 − 2

2   

T 0

K∈Th e∈∂K i=1

=− 1 = 2

=

1 2

0

0 T

0 T



T

e∈Eh

 e∈Eh



e∈Eh

e

e

e

e

Ai BK ·

i BK νe,K dΓdt

1 =− 2

T

0

  K∈Th e∈∂K

e

Ae,K BK · BK dΓdt

(Ae,K BK · BK − Ae,K BKe · BKe )dΓdt

(BTK + BTKe )Ae,K (BKe − BK )dΓdt T

Bh,e Ae,K [Bh ]e dΓdt =

0

T

 e∈Eh

e

Ae,K Bh,e · [Bh ]e dΓdt.

(131)

ANALYSIS OF DISCONTINUOUS-GALERKIN METHODS FOR THE INDUCTION EQUATIONS

1201

From (128)–(131) we can write Lh (Bh , Bh )

1 1 1 T  2 2 Bh (T ) L2 (Ω) − Bh (0) L2 (Ω) + = |Ae,K |[Bh ]e · [Bh ]e dΓdt 2 2 2 0 e∈Eh e

T

1 T + CBh · Bh dxdt + ∇ · u Bh 22 dxdt. 2 0 0 Ω Ω

Since Lh (Bh , Bh ) = 0 then from the previous equation we get the inequality   Bh (t) 2L2 (Ω) ≤ Bh (0) 2L2 (Ω) + 2 C L∞ ([0,T ]×Ω) + ∇ · u L∞ ([0,T ]×Ω)

0

t

Bh (s) 2L2 (Ω) ; ds,

for all t ≤ T . Then using a Gronwall inequality we obtain 1

Bh (t) L2 (Ω) ≤ Bh (0) L2 (Ω) e(CL∞ ([0,T ]×Ω) + 2 ∇·uL∞ ([0,T ]×Ω) )t ,

∀t ≤ T.

References ´ [1] S. Alinhac and P. G´erard, Op´ erateurs pseudo-diff´ erentiels et th´ eor` eme de Nash-Moser, CNRS Editions (1991). [2] G.A. Baker, W.N. Jureidini and O.A. Karakashian, A piecewise solenoidal vector fields and the stokes problem. SIAM J. Numer. Anal. 27 (1990) 1466-1485. [3] D.S. Balsara and D.S. Spicer, A piecewise solenoidal vector fields and the stokes problem. J. Comput. Phys. 149 (1999) 270–292. [4] J.U. Brackbill and D.C. Barnes. The effect of nonzero ∇ · B on the numerical solution of the magnetohydrodynamic equations. J. Comput. Phys. 35 (1980) 426. [5] P.G. Ciarlet, Basic error estimates for elliptic problems, in Handbook of numerical analysis, P.G. Ciarlet and J.-L. Lions, Eds., North-Holland (1991) 17–351. [6] B. Cockburn, Discontinuous Galerkin methods for convection dominated problems, in High-order methods for computational physics, Springer, Berlin. Lect. Notes Comput. Sci. Eng. 9 (1999) 69–224. [7] B. Cockburn, F. Li and C.-W. Shu, Locally divergence-free discontinuous Galerkin-methods for the Maxwell equations. J. Comput. Phys. 194 (2004) 588–610. [8] M. Costabel, A remark on the regularity of solutions of Maxwell’s equations on Lipschitz domains. Math. Method Appl. Sci. 12 (1990) 365–368. [9] M. Costabel and M. Dauge, Un r´esultat de densit´e pour les ´equations de Maxwell r´egularis´ ees dans un domaine lipschitzien. C. R. Acad. Sci. Paris S´ er. I 327 849–854 (1998). [10] W. Dai and P.R. Woodward, On the divergence-free condition and conservation laws in numerical simulations for supersonic magnetohydrodynamic flow. Astrophys. J. 494 (1998) 317. [11] A. Dedner, F. Kemm, D. Kr¨ oner, C.-D. Munz, T. Schnitzer and M. Wesenberg, Hyperbolic Divergence cleaning for the MHD equations. J. Comput. Phys. 175 (2002) 645–673. [12] C.R. Evans and J.F. Hawley, Simulation of magnetohydrodynamic flows, a constrained transport method. Astrophys. J. 332 (1988) 659. [13] C. Foias and R. Temam, Remarques sur les ´equations de Navier-Stokes et les ph´enom´enes successifs de bifurcation. Ann. Sci. Norm. Sup. Pisa S´ er. IV 5 (1978) 29–63. [14] K.O. Friedrichs, Symmetric positive linear differential equations. Comm. Pure Appl. Math. XI (1958) 333–418. [15] V. Gilrault and P.-A. Raviart, Finite element methods for the Navier-Stokes equatons, Theory and algorithms. Springer Ser. Comput. Math. 5 (1986). [16] O.A. Karakashian and W.N. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations. SIAM J. Numer. Anal. 35 (1998) 93–120. [17] F. Li, C.-W. Shu, Locally divergence-free discontinuous Galerkin methods for MHD equations. SIAM J. Sci. Comput. 27 (2005) 413–442. [18] J.-L. Lions and J. Petree, Sur une classe d’espaces d’interpolation. Publ. I.H.E.S. 19 (1964) 5–68. [19] J.C. N´ed´ elec, Mixed finite element in 3 . Numer. Math. 35 (1980) 315–341.

Ê

1202

¨ N. BESSE AND D. KRONER

[20] K.G. Powell, An approximate Riemann solver for magnetohydrodynamics (that works in more than one dimension), ICASEReport 94-24 (NASA CR-194902), NASA Langley Research Center, Hampton, VA 23681-0001 (1994). [21] P.-A. Raviart and J.M. Thomas, A mixed finite element method for 2nd order elliptic problems, Mathematical aspects of finite elements methods, in Proc. of the conference held in Rome, 10–12 Dec. 1975, A. Dold, B. Eckmann, Eds., Springer, Berlin, Heidelberg, New York. Lect. Notes Math. 606 (1977). [22] G. T´ oth, The ∇ · B = 0 constraint in shock-capturing magnetohydrodynamics codes. J. Comput. Phys. 161 (2000) 605. [23] L. Ying, A second order explicit finite element scheme to multidimensional conservation laws and its convergence. Sci. China Ser. A 43 (2000) 945–957. [24] Q. Zhang and C.-W. Shu, Error Estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws. SIAM J. Numer. Anal. 42 (2004) 641–666.

To access this journal online: www.edpsciences.org