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

3 downloads 10 Views 302KB Size Report
Apr 19, 2004 - We present the convergence analysis of locally divergence-free discontinuous .... For a subdomain D ⊂ Ω, Hm(D) denotes the usual Sobolev.

Convergence of locally divergence-free discontinuous-Galerkin methods for the induction equations of the MHD system ∗ Nicolas Besse



Dietmar Kr¨oner



April 19, 2004

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.

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 conservative form

∂t (ρu) + ∇ · ρu ⊗ u + p +

∂t (ρe) + ∇ ·

∂t ρ + ∇ · (ρu) = 0 (mass conservation)   I − B ⊗ B = 0 (momentum conservation)

1 2 2 |B|

∂t B + ∇ · (u ⊗ B − B ⊗ u) = 0 (induction equation)   ρe + p + 12 |B|2 u − B(u · B) = 0 (energy conservation) ∇ · B = 0 (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 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, 11, 12, 19, 21]. Because ∇ · B errors arise in numerical simulations and may increase in time, numerical instabilities can appear and lead to †

Institut de Recherche Mathematique Avanc´ee, Universit´e Louis Pasteur - CNRS, 7 rue Ren´e Descartes, 67084 Strasbourg Cedex, France ([email protected]) ‡ Institut f¨ ur Angewandte Mathematik Universit¨ at Freiburg Hermann-Herder-Str.10 D-79104 Freiburg i. Br., Germany ([email protected]) ∗ The authors acknoledge financial support from the HYKE project No HPRN-CT-2002-00282 on “Hyperbolic and Kinetic Equations: Asymptotics, Numerics, Applications”, funded by European Union. Mathematics Subject classifications: 76W05, 65J10. Keys words: magnetohydrodynamics, Discontinuous-Galerkin methods, convergence analysis.

1

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 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 is based on three things. 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 developped 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) [18] in two dimensions which is obtained by rotating the Raviart-Thomas element [20] by π/2. Then we used the framework developped in [6, 23] 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.

2

The numerical scheme

2.1

The Friedrichs formulation

Thanks to the divergence constraint, ∇ · B = 0, the induction equation (1) can be rewritten as ∂B ∂(Ax B) ∂(Ay B) + + + CB = 0 ∂t ∂x ∂y where

 Ax (t, x) =

ux (t, x) 0 0 ux (t, x)

and

 C(t, x) = −



in Ω × [0, T ], 

Ay (t, x) =

,

uy (t, x) 0 0 uy (t, x)

∂x ux (t, x) ∂y ux (t, x) ∂x uy (t, x) ∂y uy (t, x)

(2)  .

 .

Sometime 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 [13] has proved the existence and uniqueness in L2 (Ω) of a weak and strong solution to (2), under appropiate smoothness conditions and the additional assumption ∂Ax ∂Ay C + CT + + ≥ α I, in Ω, (3) ∂x ∂y with α a strictly positive constant. We suppose that u ∈ C ∞ (Ω), ∇ · u ∈ L∞ (Ω) and C ∈ C ∞ (Ω). We recall that if B0 ∈ Hs (R2 ) with s ∈ R, then the  system (2) admits a unique weak solution 0 s 2 1 s−1 2 B ∈ C [0, +∞[; H (R ) ∩ C  [0, +∞[; H (R ) (see [1]). In particular, if s > 2 the weak solution lies in C 1 [0, +∞[×R2 and is in fact a classical solution. 2

e x) = e−λt B(t, x), with λ > 0 Remark 1 The condition (3) can be easily satisfied by using B(t, e satisfies instead of B. Then B e e e ∂(Ay B) ∂(Ax B) ∂B e =0 + + + CeB ∂t ∂x ∂y

in Ω × [0, T ],

(4)

where Ce = 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 NT n ` (0, T ; X) := f : {t , ..., t } → X| ||f ||`∞ (0,T ;X) = max ||f (t )||X < ∞ 1≤n≤NT

where X denotes a functional space. For a subdomain D ⊂ Ω, H m (D) denotes the usual Sobolev space which is a Hilbert space with inner product X Z (u, v)m,D = ∂ α u∂ α vdx, u, v ∈ H m (D) |α|≤m D

and seminorm

1/2

 |v|H ( D) = 

X Z

|∂ α v|2 

.

|α|=m D

Then we define the space Hm (D) = (H m (D))2 with inner product (v, w)m,D =

2 X

(vi , wi )m,D ,

i=1

and the solenoidal vector fields Sm (D) = {v ∈ Hm (D) | ∇ · v = 0 in D} , m ≥ 1. Let Th a family of partitons of Ω that possesses properties described in [2, 16]. For integer k ≥ 0, Pk (D) will denote the linear space of polynomials in two variables, of degree less than or equal to k on D. Then we define  2 k (D) = Pk (D)

P

and

n V k (D) = v ∈

Pk (D) | ∇ · v = 0 in D

o

We have V k (D) ⊂ Sm (D) for all m ≥ 1. For k ≥ 0 we define Y V k (K) Vhk (Ω) = K∈Th

where K is a finite element (or cell) and N the number of finite elements. The way of constructing local bases for V k and hence for Vhk are described in [2, 16]. For examples, for k = 1 the bases functions can be constructed using the set of functions           1 0 y 0 x 1 Ξ = , , , , . 0 1 0 x −y 3

To obtain a basis function for k = 2 it suffices to augment the above set by  2     2    0 x −2xy y 2 , , , . Ξ = 0 x2 −2xy y2 Then we see that the basis functions for the approximation space Vhk can be constructed by adding suitable terms to Vhk−1 . Moreover these spaces possess optimal approximation properties in relation to the spaces Sm (Ω) which has been proved in [2]. We recall some of those. Theorem 1 Let m ≥ 0. If v ∈ Hm+1 (K), then there exists χ ∈

Pm(K) such that

kv − χkHj (K) ≤ Chm+1−j |v|Hm+1 (K) , 0 ≤ j ≤ m + 1 K

(5)

where hK is the diameter of the cell K. Let m ≥ 0. If v ∈ Sm+1 (K), then there exists χ ∈ Vm (K) such that kv − χkHj (K) ≤ Chm+1−j |v|Hm+1 (K) , 0 ≤ j ≤ m + 1. K

(6)

Proof. See [2].  Now we define the local L2 projection operator πh from Sm (Ω) onto Vhk (Ω), with m ≥ 1 by the equations Z (πK v − v) · φ = 0,

∀φ ∈ V k (K), ∀K ∈ Th .

(7)

πK 1K (x)

(8)

K

πh =

X K∈Th

where 1K (x) denotes the characteristic functions of K which is equal to 1 on K and 0 elsewhere. Thanks the properties of the partition Th (see [2]) there exists a constant σ such that hK ≤ σh, ∀K ∈ Th . For example we can take h = minK∈Th hK . Then we have the following approximation result. Theorem 2 Let Th be a partition of Ω satisfying the properties assumed in [2]. if v ∈ Sm+1 (Ω) then there exists a constant C such that kπh v − vkL2 (Ω) ≤ Chm+1 |v|Hm+1 (Ω) .

(9)

Proof. Thanks theorem 1 there exists χ ∈ Vhm (Ω) such that kv − χkL2 (K) ≤ Chm+1 K |v|Hm+1 (K) ,

∀K ∈ Th .

On the other hand we have kπh v − vk2L2 (Ω) ≤

X

kπh v − vk2L2 (K)

K∈Th



X 

 kπK (v − χ)k2L2 (K) + kπK χ − χk2L2 (K) + kχ − vk2L2 (K) .

K∈Th

From (10) we have kv − χk2L2 (K) ≤ Ch2m+2 |v|2Hm+1 (K) . K

4

(10)

From (7), (10) and Cauchy-Schwarz inequality we obtain that kπK (v − χ)k2L2 (K) ≤ kπK k kv − χk2L2 (K) ≤ kv − χk2L2 (K) ≤ Ch2m+2 |v|2Hm+1 (K) . K Finaly we have Z

∀φ ∈ V m (K), ∀K ∈ Th .

(πK χ − χ) · φ = 0, K

If we choose φ = (πK χ − χ)|K , then we have kπK χ − χk2L2 (K) = 0. Then we obtain X

kπh v − vk2L2 (Ω) ≤

kπh v − vk2L2 (K)

K∈Th

≤ C

X

2m+2 hK |v|2Hm+1 (K) ,

K∈Th

≤ Cσ 2 h2m+2 ≤ Ch

2m+2

X

|v|2Hm+1 (K)

K∈Th 2 |v|Hm+1 (Ω) .

 Moreover tkanks the assumptions of the partition Th (see [2]) we have the following inverse properties (see [5]). |vh |H1 (K) ≤

C kvh kL2 (K) , h

|vh |L2 (∂K) ≤

C kvh kL2 (K) , h1/2

∀vh ∈ V m (K)

(11)

Now we introduce an other discontinuous interpolation operator Πh based on the two-dimensional N´ed´elec element in H(curl) (see [18]). Let e k denote the space of homogeneous polynomials of degree k in R2 and consider the following subspace of k :

P

Rk =

P

Pk−1 ⊕ Sk

(12)

where Sk is defined by

P

n o Sk = p ∈ e k ; p(x) · x = 0, for all x = (x1 , x2 ) .

(13)

Then we define W k (K) by W k (K) = {v ∈ Rk (K)}

(14)

and Whk (Ω) =

Y

W k (K).

(15)

K∈Th

Following [14], if we define the two sets of moments of v ∈ Hs with s ≥ 1/2 on K: Z  m−1 Me (v) = (v × ν e ) q dΓ ∀q ∈ P (e) for all edges e of K

(16)

where ν e is the outward unit normal to the edges e of K, and Z  m−2 MK (v) = v · qdx, ∀q ∈ (K) ,

(17)

e

P

K

5

then a vector field v of Rk is entirely determined in a triangle K by its two sets of moments: Me (v), MK (v). Moreover the tangential components of v on a given edge e of K depend only the moments Me (v) defined on that edges. Besides we have the following approximation properties (see [14, 18]): For all v ∈ Hm+1 (Ω) we have kv − Πh vkL2 (K) ≤ Chm |v|Hm (K)

(18)

  kv − Πh vkH(curl;K) ≤ Chm |v|Hm (K) + |v|Hm+1 (K)

(19)

and k(v − Πh v) × ν e kHs (e) ≤ Chm−s |v|Hm (e) ∀e ∈ ∂K

s≤m

(20)

with Πh v|K = ΠK v where ΠK v ∈ W m (K) is determined by MK (ΠK v − v) = {0},

Me (ΠK v − v) = {0}.

(21)

P

If m = 1 then R1 = 0 ⊕S1 where S1 = span(x⊥ ) and then it is obvious that div R1 ≡ 0. Then the vector fields of R1 is divergence free. Unfortunately the spaces Rk with k > 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 R2 . First we define the space V(K) and W(K) by  W(K) = v ∈ H(curl; K) ∩ H(div; K), v × ν ∈ L2 (∂K) (22) equipped with the norm kvkW(K) = kvkL2 (K) + kcurl vkL2 (K) + kdiv vkL2 (K) + kv × νkL2 (∂K)

(23)

 V(K) = v ∈ H(curl; K) ∩ H(div; K), v · ν ∈ L2 (∂K)

(24)

and equipped with the norm kvkV(K) = kvkL2 (K) + kcurl vkL2 (K) + kdiv vkL2 (K) + kv · νkL2 (∂K) .

(25)

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 define n o Q(K) = v ∈ H(curl; K) ∩ H(div; K), v × ν e ∈ H 1/2 (e) ∀e ∈ ∂K

(26)

equipped with the norm kvkQ(K) = kvkL2 (K) + kcurl vkL2 (K) + kdiv vkL2 (K) +

X e∈∂K

Then we have the following theorem

6

kv × ν e kH 1/2 (e)

(27)

Theorem 3 Let K be a convex polygon in R2 with lipschitz boundary. Then we have the continuous embedding: Q(K) ,→ H1 (K) (28) Proof. Let K a convex polygon where {Γi }{i=1,...Ne } denote the edges of K. We suppose that the vector field u in R2 is such that u ∈ H1 (K), div u, curl u ∈ L2 (K) and u × ν|Γ ∈ H 1/2 (Γi ) i with ν the normal of ∂K. Let ψ the solution of the problem ∆ψ = curl u ∂ν ψ = u × ν Since K is convex and u × ν is piecewise in H 1/2 from a result in [15] or in [14] (theorem 1.10) there exists a constant C such that   X kψkH 2 (K) ≤ C kcurl ukL2 (K) + ku × νkH 1/2 (Γi )  . (29) Γi ∈∂K

Then we construct u0 such that u0 = ∇⊥ ψ where ∇⊥ = (∂2 , −∂1 )T then u0 × ν = u × ν on ∂K and from (29) we obtain   X ku0 kH 1 (K) ≤ C kcurl ukL2 (K) + ku × νkH 1/2 (Γi )  . (30) Γi ∈∂K

Now we set v = u − u0 with v × ν = 0 on ∂K. From the proposition 3.1 of [14] we get  kvkH 1 (K) ≤ C kvkL2 (K) + kcurl vkL2 (K) + kdiv vkL2 (K) .

(31)

Therefore, from (30) and (31) we obtain kukH 1 (K) ≤ kvkH 1 (K) + ku0 kH 1 (K) 

(32) 

≤ C kvkL2 (K) + kcurl ukL2 (K) + kdiv ukL2 (K) +

X

ku × νkH 1/2 (Γi )  , (33)

Γi ∈∂K

which ends the proof.

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 Z Z X Z i wvνe,K dΓ (34) v∂i wdx = − w∂i vdx + K

K

e∈∂K

e

1 , ν 2 )T = (ν x , ν y )T denotes the outward unit normal to the face e of K, we where ν e,K = (νe,K e,K e,K e,K obtain Z  X Z 2 Z 2 X Z X i CB · ϕdx = 0. (35) B · ϕdx − Ai B · ∂i ϕdx + Ai B · ϕνe,K dΓ + ∂t K

i=1

K

i=1 e∈∂K

7

e

K

Now we replace respectively B and ϕ by Bh and ϕh in (35) where Bh , ϕh ∈ Vhk (Ω). Nevertheless the terms arising from the boundary of the cell K in (35) are not well defined or have no sense since Bh and ϕh are discontinuous on the bondary ∂K of the element K. Then we replace these terms by a numerical upwind flux that we are going to define in the following lines. Let define 2 X i . Ae,K (t, x) = (Ai )|e νe,K i=1

and 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 (v + w) (w − v) g(νe,K , v, w) = Ae,K − |Ae,K | . 2 2 Now we obtain the following semi-discretized scheme in space Z  X Z 2 Z 2 X Z X ∂t Bh · ϕh dx − Ai Bh ·∂i ϕh dx+ g(ν e,K , BK , BKe )·ϕK dΓ+ CBh ·ϕh dx = 0. K

K

i=1

i=1 e∈∂K

e

K

for all ϕh ∈ Vhk where Ke is the neighbor 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 discretized scheme reads Z Z n n (36) (Bnh , ϕh ) Yh · ϕh dx = Bnh · ϕh dx + ∆tFK K K  Z Z  1 n 1 n ∆t n+1 n n+1 Bh · ϕh dx = Bh + Yh · ϕh dx + F (Yh , ϕh ) (37) 2 2 2 K K K where n FK (W, ϕh ) =

2 Z X i=1

Ani Wn · ∂i ϕh dx +

K

i=1 e∈∂K

with n

g (W)e = g

n

2 X Z X

n n (νe,K , WK , WK ) e

=

n n (WK Ae,K

g n (W)e · ϕK dΓ −

Z

e

C n Wn · ϕh dx

K

n − Wn ) n ) + WK (WK K n e e − |Ae,K | 2 2

and Ane,K = Ae,K (tn , x), C n = C(tn , x). Let X k be V k or W k . If we expand Bh and ϕh , that is to say dim(X k )

Bnh (x)

=

X K∈Th

1

BnK (x) K (x),

BnK (x)

=

X

dim(X k ) i,n σK Ψi (x),

n YK (x)

=

i=1

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 8

X i=1

i,n θK Ψi (x)

(38)

where [ΣK ]n =



 dim (X k ),n T

1,n σK , . . . , σK



, [Σh ]n =

σ 1,n , . . . , σ card(Th )dim(X

k ),n

T

and for all i ∈

{1, ..., dim(X k )} [LnK

n

([Σh ] )]i =

2 dim(X X X

k)

j,n σK

Z K

l=1 m,j=1

+

2 X Z X i=1 e∈∂K

X



m,j=1

j,n σK

X

∂l Ψm (x)dx

Z K

  j,n n j,n n −T [σK De,K (x(Γ)) − σK C (x(Γ))]Ψ (x(Γ)) · M j K e e,K

  C n (x)Ψj (x) · M−T K

with

m,i

m,i

Ψm (x(Γ))dΓ

Ψm (x)dx

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

(MK )ij =

3

m,i

dim(X k )

e m,j=1

dim(X k )

  Anl (x)Ψj (x) · M−T K

K

Analysis of the semi-discretized scheme in space

In this section we show the L2 -stabilty, the convergence and present some error estimates for the semi-discretized scheme in space, continuous in time.  Theorem 4 Let u and B0 sufficiently, regular, typically we consider that u ∈ C 0 [0, +∞[; Hm+1 (R2 ) and B0 ∈ Hm+1 (R2 ) then there exists a constant C = C(k∇ukL∞ (Ω) , T ) independent of h such that kBh kL∞ (0,T ;L2 (Ω)) ≤ CkBh (0)kL2 (Ω) , 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(k∇ukL∞ (Ω) , T, kBkL∞ (0,T ;Hm+1 (R2 )) ) independent of h such that v uZ u kB − Bh kL∞ (0,T ;L2 (Ω)) + t

T

0

where η =

1 2

XZ e∈Eh

|Ae |[Bh ]e · [Bh ]e dΓdt ≤ Chm+η

e

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

Proof. 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 -stabilty of the scheme. This proof remains true for the two choices of approximation space. Let Lh (·, ·) the bilinear form defined by Z

T

Z

Lh (Bh , ϕh ) =

∂t Bh · ϕh dxdt − 0

Z + 0

Ω T

i=1

X X Z K∈Th e∈∂K

2 Z X 0

T

X Z

Ai Bh · ∂i ϕh dxdt

K

K∈Th

Z

T

Z

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

e

CBh · ϕh dxdt. 0

9



(39)

with ϕh ∈ Vhm or ϕh ∈ Wh1 . If we set [ϕh ]e = ϕKe − ϕK Z

T

Z

Lh (Bh , ϕh ) =

∂t Bh · ϕh dxdt − 0

Z

Ω T

− 0



2 Z X i=1

XZ

then we can rewrite (39) as

|e

T

0

Ai Bh · ∂i ϕh dxdt Ω

Z

T

Z

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

e

e∈Eh

Z

CBh · ϕh dxdt 0

(40)



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

1 1 kB(T )k2L2 (Ω) − kB(0)k2L2 (Ω) 2 2 2 Z T Z X X Ai BK · ∂i BK dxdt K∈Th i=1 0 Z T XZ

(41)

K

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

0

Z

e∈Eh T

(42)

e

Z CBh · Bh dxdt.

+ 0



From the definition (38) of the upwind flux, the term (42) can be rewritten as following Z 0

T

XZ e∈Eh

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

e

0

T

XZ e∈Eh

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

T

Z 0

XZ e∈Eh

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

e

(43) where Bh,e =



BK +BKe 2

 |e

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

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

(44)

and the Green formula (34) then we obtain −

2 Z X X K∈Th i=1

1 = 2

T

Z Ai BK · ∂i BK dxdt

0

Z

K T

Z ∇·

0



ukBh k22 dxdt

1 − 2

Z 0

10

T

2 Z X X X K∈Th e∈∂K i=1

e

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

(45)

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

1 2

Z 0

T

2 Z X X X

i Ai BK · BK νe,K dΓdt = −

e

K∈Th e∈∂K i=1

1 2

Z

T

0

Z

X X Z K∈Th e∈∂K

Ae,K BK · BK dΓdt

e

T

XZ 1 = − (Ae,K BK · BK − Ae,K BKe · BKe )dΓdt 2 0 e∈Eh e Z Z 1 T X = (BTK + BTKe )Ae,K (BKe − BK )dΓdt 2 0 e e∈Eh Z T XZ T = Bh,e Ae,K [Bh ]e dΓdΓdt 0

Z

e∈Eh T

= 0

e

XZ e∈Eh

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

(46)

e

From (43)-(46) we can write Z Z 1 1 1 T X 2 2 Lh (Bh , Bh ) = kBh (T )kL2 (Ω) − kBh (0)kL2 (Ω) + |Ae,K |[Bh ]e · [Bh ]e dΓdt 2 2 2 0 e e∈Eh Z TZ Z TZ 1 + CBh · Bh dxdt + ∇ · ukBh k22 dxdt. 2 0 Ω 0 Ω Since Lh (Bh , Bh ) = 0 then from the previous equation we get the inequality kBh (t)k2L2 (Ω)



kBh (0)k2L2 (Ω)

+ 2kCkL∞ (Ω) + k∇ukL∞ (Ω)



Z 0

t

kBh (s)k2L2 (Ω) ds,

for all t ≤ T . Then using a Gronwall inequality we obtain 1 kBh (t)kL2 (Ω) ≤ kBh (0)kL2 (Ω) e(kCkL∞ (Ω) + 2 k∇ukL∞ (Ω) )t ,

∀t ≤ T.

Now we will show simultaneously the convergence of the scheme and derivate some error estimates. First we treat the case where Bh ∈ Vhm . First by noting that Lh (Bh , ϕh ) = 0

and Lh (B, ϕh ) = 0

∀ϕh ∈ Vhm (Ω)

then Lh (eh , ϕh ) = 0,

∀ϕh ∈ Vhm (Ω)

(47)

where eh = B − Bh . By noting that πh Bh = Bh and πh eh − eh = πh (B − Bh ) − (B − Bh ) then using (47) we obtain Lh (πh eh , πh eh ) = Lh (πh eh − eh , πh eh ) = Lh (πh B − B, πh eh ).

11

(48)

From the stabilty analysis done previously the left hand side of (48) is the expression Z Z 1 1 1 T X kπh eh (T )k2L2 (Ω) − kπh eh (0)k2L2 (Ω) + |Ae |[πh eh ]e · [πh eh ]e dΓdt 2 2 2 0 e e∈Eh Z TZ Z TZ 1 + Cπh eh · πh eh dxdt + ∇ · ukπh eh k22 dxdt. (49) 2 0 Ω 0 Ω

Lh (πh eh , πh eh ) =

Now it remains to estimate the right hand side of (48) that is to say find an estimate for ∀ϕh ∈ Vhm (Ω).

Lh (πh B − B, ϕh ), Let us write Lh (πh B − B, ϕh ). Z

T

Z

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

∂t (πh B − B) · ϕh dxdt 0 Ω 2 Z T X

Z

Ω i=1 0 Z T XZ 0

Z

e∈Eh T

(50)

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

(51)

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

(52)

e

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

+ 0

(53)



From the definition of πh for the term (50) we have Z

T

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

0

(54)



Now we define a new projection operator Ph by Z 1 (Ph g)|K = PK g = gdx, ∀g ∈ L1 (Ω) ∩ W 1,∞ (Ω), |K| K

and Ph =

X

P K 1K .

K∈Th

Then by Taylor expansion it can be easily shown that there exists a constant C independent of h such that kPh g − gkL∞ (Ω) ≤ ChkgkW 1,∞ (Ω) . (55) Then the i-th term of (51) can be recasted in Z − 0

T

Z TZ Ai (πh B − B) · ∂i ϕh dxdt = (Ph Ai − Ai )(πh B − B) · ∂i ϕh dxdt Ω 0 Ω Z TZ − (πh B − B) · Ph Ai ∂i ϕh dxdt 0 Ω Z TZ = (Ph Ai − Ai )(πh B − B) · ∂i ϕh dxdt

Z

0

(56)



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 term (56). Thanks the approximation properties (9) and (55), the inverse

12

properties (11), and the Cauchy-Schwarz inequality we obtain Z T Z (Ph Ai − Ai )(πh B − B) · ∂i ϕh dxdt 0



T

Z ≤ 0

X Z K∈Th

1/2 Z

kπh B −

Bk22 dx

k(Ph Ai −

K

Ai )∂i ϕh k22 dx

1/2 dt

K

Z

T

≤ kAi kL∞ (0,T ;W 1,∞ (Ω)) 0

Z ≤ kAi kL∞ (0,T ;W 1,∞ (Ω)) 0

X

Chm+2 |B|L∞ (0,T ;Hm+1 (K)) k∂i ϕh kL2 (K) dt

K∈Th T

X

Chm+1 |B|L∞ (0,T ;Hm+1 (K)) kϕh kL2 (K) dt

K∈Th

≤ CkAi kL∞ (0,T ;W 1,∞ (Ω)) hm+1

Z

T

1/2 

 X 

0

|B|2L∞ (0,T ;Hm+1 (K)) 

K∈Th

1/2 X



kϕh k2L2 (K) 

dt

K∈Th

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

Z

T

kϕh kL2 (Ω) dt Z  m+1 T ≤ C T, kAi kL∞ (0,T ;W 1,∞ (Ω)) , |B|L∞ (0,T ;Hm+1 (Ω)) h kϕh kL2 (Ω) dt. 0

(57)

0

For the term (53) we get Z

T

0

Z T X Z C(πh B − B) · ϕh dxdt ≤ C(πh B − B) · ϕh dxdt Ω 0 K∈Th K 1/2 Z T X Z 2 ≤ kCkL∞ (0,T ;L∞ (Ω)) kπh B − Bk2 dx kϕh kL2 (K) dt

Z

0

K∈Th

K

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

Z

T

kϕh kL2 (Ω) dt. 0

(58) Now it remains to estimate the term (52). First have 1 g(πh B − B)e = 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 ) . Since |Ae,K | is symmetric and using a Young inequality we obtain Z 0

T

XZ e∈Eh

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

e

0

Z ≤ 0

T

T

XZ e∈Eh

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

e

Z Z

2 XZ 1 T X

−1/2 |Ae |[ϕh ]e · [ϕh ]e dΓdt. (59) g(πh B − B)e dxdt +

|Ae,K | 4 0 2 e e e∈Eh

e∈Eh

13

Now we want to estimate the first term of (59). Z 0

T

2 XZ

|Ae,K |−1/2 g(πh B − B)e dxdt 2

e

e∈Eh

T

Z ≤ 2 max kAi kL∞ (0,T ;L∞ (Ω)) i=1,2

0 T

Z

X

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

≤ 2kukL∞ (0,T ;L∞ (Ω)) h−1

XZ  e∈Eh

kπh B − Bk2L2 (∂K) dt

K∈Th T X

Z

0

e

 k(πh B − B)Ke k22 + k(πh B − B)K k22 dxdt

kπh B − Bk2L2 (K) dt

K∈Th

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

(60)

Finally, from (48)-(54), (57), (58), (59) and (60) we get kπh eh (T )k2L2 (Ω) − kπh eh (0)k2L2 (Ω) + C2 hm+1

Z

1 2

Z

T

T

0

Z

XZ e∈Eh T Z

kπh eh (t)kL2 (Ω) + 0

0



|Ae |[πh eh ]e · [πh eh ]e dΓdt ≤ C1 h2m+1 +

e

 2kCkL∞ (Ω) + k∇ukL∞ (Ω) kπh eh (t)k2L2 (Ω) dt. (61)

Using a Young inequality, from (61) we get that ∀t ≤ T kπh eh (t)k2L2 (Ω) ≤ C1 h2m+1 + C2 h2m+2 + kπh eh (0)k2L2 (Ω) Z  t + 2kCkL∞ (Ω) + k∇ukL∞ (Ω) + C2 kπh eh (s)k2L2 (Ω) ds.

(62)

0

From (62) and a Gronwall inequality we get   1 kπh eh (t)kL2 (Ω) ≤ C kπh eh (0)kL2 (Ω) + hm+1/2 e(kCkL∞ (Ω) + 2 k∇ukL∞ (Ω) +C2 /2)t ≤ Chm+1/2 . Finaly we get that ∀t ≤ T keh (t)kL2 (Ω) ≤ kπh eh (t)kL2 (Ω) + kπh eh (t) − eh (t)kL2 (Ω) ≤ kπh eh (t)kL2 (Ω) + kπh B(t) − B(t)kL2 (Ω) ≤ C(hm+1 + hm+1/2 ) ≤ Chm+1/2 which finishes the proof when Bh ∈ Vhm . Now we suppose that Bh ∈ Wh1 . We have to give new estimates for the terms (50)-(53). For the term (50), using Cauchy-Schwarz inequality and the approximation property (18) we get Z

T

Z



Z

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



T

kϕh kL2 (Ω) dt.

(63)

0

Now let us estimate the term (52). The inequality (59) is still valid. It remains to give a new estimate for the first term of the right hand side of the inequality (59). Using the approximation 14

properties (18)-(20), and the fact that divΠh B = divB = 0 then we get Z 0

T

2 XZ

|Ae,K |−1/2 g(Πh B − B)e dxdt e∈Eh

2

e

Z

T

≤ 2 max kAi kL∞ (0,T ;L∞ (Ω)) i=1,2

0 T

Z

X

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

X

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

kΠh B − Bk2L2 (∂K) dt kΠh B − Bk2W(K) dt

K∈Th T

Z

e

e∈Eh

 k(Πh B − B)Ke k22 + k(πh B − B)K k22 dxdt

K∈Th T

Z

XZ 

X 

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

 kΠh B − Bk2H(curl;K) + k(Πh B − B) × νk2L2 (∂K) dt

K∈Th

 ≤ C T, kukL∞ (0,T ;L∞ (Ω)) , kBkL∞ (0,T ;H2 (Ω)) h2 .

(64)

Now we give a new estimate for (51) and (53). Using the Green formula (34) we obtain −

2 Z X

T

Z

0

i=1

Z

T

Z

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

Z = 0

T

Z − 0

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

2 X

X Z K

K∈Th



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

· ϕh dxdt

(65)

i=1

X X Z K∈Th e∈∂K

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

(66)

e

First we look at the term (66). Using the approximation properties (18)-(20), inverse properties (11), the fact that div Πh B = div B = 0 on K and the Cauchy-Schwarz inequality then we get Z 0

T

X X Z K∈Th e∈∂K

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

e

Z

T

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

Z

K∈Th T

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

Z

T

0

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

X

∂K

kΠh B − Bk2 kϕh k2 dΓdt

kΠh B − BkL2 (∂K) kϕh kL2 (∂K) dt

K∈Th

≤ 2kukL∞ (0,T ;L∞ (Ω)) Z

X Z

X

h−1/2 kΠh B − BkW(K) kϕh kL2 (K) dt

K∈Th T

X

  h−1/2 kΠh B − BkH(curl;K) + k(Πh B − B) × νkL2 (∂K) kϕh kL2 (K) dt

K∈Th

 ≤ C T, kukL∞ (0,T ;L∞ (Ω)) , kBkL∞ (0,T ;H2 (Ω)) h1/2

Z 0

15

T

kϕh kL2 (Ω) dt.

(67)

As div Πh B = div B = 0 on K the term (65) can be recast in ! Z T X Z 2 X ∂i (Ai (Πh B − B)) + C(Πh B − B) · ϕh dxdt 0

i=1 K∈Th K Z T X Z

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

=

0

Z

K

K∈Th T

= 0

X Z

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

(68)

K

K∈Th

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

Z 0

X Z K∈Th

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

K T

Z

X Z

≤ 2kukL∞ (0,T ;W 1,∞ (Ω)) 0 T

Z

K

K∈Th

X

≤ 2kukL∞ (0,T ;W 1,∞ (Ω)) 0

kΠh B − Bk2 kϕh k2 dxdt

kΠh B − BkL2 (K) kϕh kL2 (K) dt

K∈Th



T

Z

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

kϕh kL2 (Ω) dt.

(69)

and Z

T

0

X Z K∈Th

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

K

Z

T

X Z

≤ kdivukL∞ (0,T ;L∞ (Ω)) 0

Z

K∈Th T

X

≤ kdivukL∞ (0,T ;L∞ (Ω)) 0

K

kΠh B − Bk2 kϕh k2 dxdt

kΠh B − BkL2 (K) kϕh kL2 (K) dt

K∈Th

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

Z 0

T

kϕh kL2 (Ω) dt.

(70)

Using the approximation properties (18)-(20), the continuous imbedding (28), the fact that div Πh B = div B = 0 on K and the Cauchy-Schwarz inequality then we get Z 0

T

X Z K∈Th

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

K

Z

T

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

Z ≤ 2kukL∞ (0,T ;L∞ (Ω)) 0

X

kΠh B − BkH1 (K) kϕh kL2 (K) dt

K∈Th T

! X

kΠh B − BkH(curl;K) +

K∈Th

X

k(Πh B − B) × ν e kH1/2 (e)

kϕh kL2 (K) dt

e∈∂K



≤ C T, kukL∞ (0,T ;L∞ (Ω)) , kBkL∞ (0,T ;H2 (Ω)) h

1/2

Z

T

0

16

kϕh kL2 (Ω) dt.

(71)

Finally, from (48)-(53) and (63)-(71), we get kΠh eh (T )k2L2 (Ω) − kΠh eh (0)k2L2 (Ω) + C2 (h + h1/2 )

Z

1 2

Z

T

0

T

XZ e∈Eh e T Z

Z kΠh eh (t)kL2 (Ω) +

0

0



|Ae |[Πh eh ]e · [Πh eh ]e dΓdt ≤ C1 h2 +  2kCkL∞ (Ω) + k∇ukL∞ (Ω) kΠh eh (t)k2L2 (Ω) dt. (72)

Using a Young inequality, from (72) we get that ∀t ≤ T kΠh eh (t)k2L2 (Ω) ≤ C1 h2 + C2 (h2 + h)kΠh eh (0)k2L2 (Ω) Z  t kπh eh (s)k2L2 (Ω) ds. + 2kCkL∞ (Ω) + k∇ukL∞ (Ω) + C2

(73)

0

From (73) and a Gronwall inequality we get   1 kΠh eh (t)kL2 (Ω) ≤ C kΠh eh (0)kL2 (Ω) + h1/2 e(kCkL∞ (Ω) + 2 k∇ukL∞ (Ω) +C2 /2)t ≤ Ch1/2 . Finaly we get that ∀t ≤ T keh (t)kL2 (Ω) ≤ kΠh eh (t)kL2 (Ω) + kΠh eh (t) − eh (t)kL2 (Ω) ≤ kΠh eh (t)kL2 (Ω) + kΠh B(t) − B(t)kL2 (Ω) ≤ C(h + h1/2 ) ≤ Ch1/2 which ends the proof when Bh ∈ Wh1 .



Remark 2 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 (19). 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 6= 0. We observe that 0 ⊂ R1 ⊂ 1 and that Wh1 ⊂ Vh1 with dim(Wh1 ) = 3 and dim(Vh1 ) = 5.

P P

P

4

P

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 ∈ C 1 [0, +∞[; Hm+1 (R2 ) and B0 ∈ Hm+1 (R2 ) and let us assume that there exists a constant β(α) depending on α such that 0k the CFL condition ∆t ≤ β(α)h4/3 holds. Moreover we assume the condition (3) and kδB L2 (Ω) =

17

kB0 − πh B0 kL(Ω) = O(hm+1 ). Then there exists a constant C = C(kukW 1,∞ ([0,T ]×Ω) , T, α) independent of h such that kBh kl∞ (0,T ;L2 (Ω)) ≤ CkBh (0)kL2 (Ω) , where Bh (0) ∈ Vhm or Bh (0) ∈ Whm with  m = 1. Moreover there exists a constant C = C (T, α, kukW 1,∞ ([0,T ]×Ω) , kBkW 3,∞ (0,T ;Hm+1 (R2 )) independent of h such that v u NT X Z u X  kB − Bh kl∞ (0,T ;L2 (Ω)) + t∆t |Ane |[Bnh − B(tn )]e · [Bnh − B(tn )]e dΓ ≤ C ∆t2 + hm+η n=0 e∈Eh

where η =

1 2

e

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

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 -stabilty of the scheme. This proof remains true for the two choices of approximation space. If we take ϕh = Bnh in (36), and ϕh = Yhn in (37) then after summing over all elements K of the partition Th the sum (36)/2+(37) gives  n+1 2 n 2 kBn+1 − Yhn k2L2 (Ω) = ∆t F n (Bnh , Bnh ) + F n+1 (Yhn , Yhn ) . (74) h kL2 (Ω) − kBh kL2 (Ω) − kBh P n 1 . First we will give an estimate of terms of the form of F n (ϕ , ϕ ) and where F n = K∈Th FK K h h n+1 n in the second moment we will give an estimate for the term kBh − Yh kL2 (Ω) . Using (44) and the Green formula (34) we get Z Z Z 1 X 1 F n (ϕh , ϕh ) = − |Ane,K |[ϕh ]e · [ϕh ]e dΓ − C n ϕh · ϕh dxdt − ∇ · un kϕh k22 dx (75) 2 2 e Ω Ω e∈Eh

From the condition (3), and the equations (75) and (74) we get   ∆t ∆t n+1 kBh kL2 (Ω) ≤ 1 − α kBnh kL2 (Ω) − α kYhn kL2 (Ω) + kBn+1 − Yhn kL2 (Ω) . h 2 2

(76)

It remains to estimate the term kBn+1 − Yhn kL2 (Ω) . To do this we remark that after summing over h all elements K of the partition Th the sum (37)−(36)/2 gives Z  ∆t (Bn+1 − Yhn ) · ϕh dx = F n+1 (Yhn , ϕh ) − F n (Bnh , ϕh ) . (77) h 2 Ω Let us estimate the right hand side of (77). We have 2 Z X X   ∆t n+1 n n n ∆t An+1 Yhn − Ani Bnh · ∂i ϕh dx F (Yh , ϕh ) − F (Bh , ϕh ) = i 2 K∈Th i=1 K Z X X X Z   n+1 n n n − ∆t g (Yh )e − g (Bh )e · ϕK dx − ∆t C n+1 Yhn − C n Bnh · ∂i ϕh dx. K∈Th e∈∂K

e

K∈Th

K

(78) By seeing that there exists a time t\ and t† such that \ Ane,K = An+1 e,K − ∆t∂t Ae,K

and 18

† |Ane,K | = |An+1 e,K | − ∆t∂t |Ae,K |

where t\ = tn + θ\ ∆t and t† = tn + θ† ∆t, with 0 ≤ θ\ ≤ 1 and 0 ≤ θ† ≤ 1, then we get g

n+1

(Yhn )e

−g

n

(Bnh )e

n An+1 e,K Y h

n Ane,K Bh

[Yhn ]e

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

=



|An+1 e,K |

where n

g ] (Bnh )e = ∂t A\e,K Bh −

∂t |A†e,K |

+

|Ane,K | 2

[Bnh ]e

[Bnh ]e .

2

Then (78) can be rewritten as  ∆t ∆t F n+1 (Yhn , ϕh ) − F n (Bnh , ϕh ) = − F n+1 (Bnh − Yhn , ϕh ) + G(Bnh , ϕh ) 2 2 where G(Bnh , ϕh ) = ∆t2

2 Z X i=1



X X Z

∂t A‡i i Bnh · ∂i ϕh dx − ∆t2

K∈Th e∈∂K

g ] (Bnh )e · ϕK dΓ − ∆t2

e

Z

∂t C [ Bnh · ϕh dx



(79) with t‡i = tn + θ‡i ∆t and t[ = tn + θ[ ∆t such that Ani = An+1 − ∆t∂t A‡i i i

C n = C n+1 − ∆t∂t C [ .

and

First we estimate the term G(Bnh , ϕh ). From the Young and Cauchy-Schwarz inequalities and the inverse properties (11) we get Z

X Z 2 X

‡i n ‡i n 2 ∆t ∂ A B ∂ A B · ∂ ϕ dx ≤ ∆t

k∂i ϕh k2 dx t t i h h h i i 2 K∈Th K K∈Th K   Z

X ∆t4 ε 2

‡i 2 n 2 kBh k2 + h k∂i ϕh k2 dx C(ε) 2 ∂t Ai ∞ ≤ h 3 L ([0,T ]×R2 ) K K∈Th

≤ C(ε)

Z 2 X [ n ∆t ∂t C Bh · ϕh dx ≤ K∈Th K

ε ∆t4 kBnh k2L2 (Ω) + kϕh k2L2 (Ω) , 2 h 3

X

∆t2

Z

[ ∂ C

t K

K∈Th

L∞ ([0,T ]×

≤ C(ε)∆t4 kBnh k2L2 (Ω) +

19

R2 )

(80)

kBnh k2 kϕh k2 dx

ε kϕh k2L2 (Ω) , 3

(81)

and ! Z Z † X X X X 2 2 ∂t |Ae,K | n n \ ] n ∆t g (Bh )e · ϕK dΓ ≤ ∆t [Bh ]e · ϕK dΓ ∂t Ae,K Bh − 2 K∈Th e∈∂K e K∈Th e∈∂K e ! Z 2 X X ∂t A\e,K − ∂t |A†e,K | n ∂t A\e,K + ∂t |A†e,K | n BKe + BK · ϕK dΓ ≤ ∆t 2 2 K∈Th e∈∂K e Z   X X

n 

BK kϕK k + kBnK k kϕK k dΓ ≤ C k∂t Ae,K kL∞ ([0,T ]×R2 ) ∆t2 2 2 2 e 2 K∈Th e∈∂K

≤C

X X Z  K∈Th e∈∂K e ∆t4 X

≤ C(ε) ≤ C(ε)

h ∆t4 h2

e

  ε

∆t4  2 n 2 n 2

BKe 2 + kBK k2 + h kϕK k2 dΓ C(ε) h 3

kBnh k2L2 (∂K) + h

K∈Th

ε X kϕh k2L2 (∂K) 3 K∈Th

kBnh k2L2 (Ω) +

ε kϕh k2L2 (Ω) . 3

(82)

 ∆t4 4 + ∆t kBnh k2L2 (Ω) + ε kϕh k2L2 (Ω) h2

(83)

From (80)-(82), we get that |G(Bnh , ϕh )|

 ≤ C(ε)

Now we will estimate the term ∆tF n+1 (Vh , ϕh ). From the Young and Cauchy-Schwarz inequalities and the inverse properties (11) we get Z X X Z

n+1 ∆t

An+1 Vh k∂i ϕh k dx Ai Vh · ∂i ϕh dx ≤ ∆t 2 i 2 K∈Th K K∈Th K  X Z  ∆t2 ε ≤ C(ε) 2 kAi kL∞ ([0,T ]×R2 ) kVh k22 + h2 k∂i ϕh k22 dx h 3 K K∈Th

≤ C(ε)

X Z n+1 ∆t ≤ C V · ϕ dx h h K∈Th K

∆t2 ε kVh k2L2 (Ω) + kϕh k2L2 (Ω) , h2 3

X

Z ∆t K

K∈Th

n+1

C

∞ kVh k2 kϕh k2 dx L ([0,T ]×R2 )

≤ C(ε)∆t2 kVh k2L2 (Ω) +

20

(84)

ε kϕh k2L2 (Ω) , 3

(85)

and ! Z Z n+1 X X X X |Ae,K | n+1 n+1 ∆t [Vh ]e · ϕK dΓ g (Vh )e · ϕK dΓ ≤ ∆t Ae,K Vh − 2 K∈Th e∈∂K e K∈Th e∈∂K e ! X X Z An+1 − |An+1 | An+1 + |An+1 | e,K e,K e,K e,K ≤ ∆t V Ke + VK · ϕK dΓ 2 2 K∈Th e∈∂K e Z   X X ≤ C kAe,K kL∞ ([0,T ]×R2 ) ∆t (kVKe k2 kϕK k2 + kVK k2 kϕK k2 ) dΓ K∈Th e∈∂K

≤C

X X Z  K∈Th e∈∂K e ∆t2 X

≤ C(ε) ≤ C(ε)

h ∆t2 h2

e

  ε ∆t2  2 2 2 C(ε) kVKe k2 + kVK k2 + h kϕK k2 dΓ h 3

kVh k2L2 (∂K) + h

K∈Th

ε X kϕh k2L2 (∂K) 3 K∈Th

kVh k2L2 (Ω) +

ε kϕh k2L2 (Ω) . 3

(86)

From (84)-(86), we get that |∆tF(Vh , ϕh )| ≤ C(ε)

∆t2 kVh k2L2 (Ω) + ε kϕh k2L2 (Ω) h2

(87)

By taking ϕh = Bn+1 − Yhn in (77), and from (83) and (87) we obtain h

n+1

2

B − Yhn L2 (Ω) ≤ ∆tF(Bnh − Yhn , Bn+1 − Yhn ) + G(Bnh , Bn+1 − Yhn ) h h h

∆t2 n 2 + C(ε) ≤ 2ε Bn+1 − Y kBnh − Yhn k2L2 (Ω) h L2 (Ω) h 2 h  4  ∆t 4 + C(ε) + ∆t kBnh k2L2 (Ω) . h2

(88)

Let estimate kBnh − Yhn k2L2 (Ω) . By taking ϕh = Bnh − Yhn in (36) and from the estimate (87) we get kBnh − Yhn k2L2 (Ω) = ∆tF n (Bnh , Bnh − Yhn ) ≤ ε kBnh − Yhn k2L2 (Ω) + C(ε)

∆t2 kBnh k2L2 (Ω) . h2 (89)

If ε is small enough, from (89) we deduce ∆t2 kBnh k2L2 (Ω) . h2 Therefore, if ε is small enough from (88) and (90) we get  

n+1

∆t4 ∆t4 n 2 4

B − Yh L2 (Ω) ≤ C(ε) ∆t + 4 + 2 kBnh k2L2 (Ω) . h h h kBnh − Yhn k2L2 (Ω) ≤ C(ε)

(90)

(91)

From (91) and (76) we get   ∆t ∆t4 ∆t4 ∆t n+1 4 1−α kBh kL2 (Ω) ≤ + C(ε)∆t + C(ε) 4 + C(ε) 2 kBnh kL2 (Ω) − α kYhn kL2 (Ω) 2 h h 2   4 4 ∆t ∆t ∆t ≤ 1−α + C(ε)∆t4 + C(ε) 4 + C(ε) 2 kBnh kL2 (Ω) . (92) 2 h h 21

If we suppose ∆t ≤ C(ε)h4/3 then (92) leads to n kBn+1 h kL2 (Ω) ≤ (1 + C∆t)kBh kL2 (Ω) .

which implies that CT kBn+1 kBh (0)kL2 (Ω) h kL2 (Ω) ≤ e

and finishs the proof of the L2 -stability. 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 Y(t, x) = B(t, x) + ∂t B(t, x)∆t (93) 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

(94)

Using Taylor expansion and from the equation (2) and (93) we obtain ∂t B(t + ∆t, x) = − = −

2 X i=1 2 X

∂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 X

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

i=1

(95) From (93) and (94) we get 1 1 ∆t B(t + ∆t, x) = B(t, x) + Y(t, x) + ∂t B(t + ∆t, x) + O(∆t3 ) 2 2 2

(96)

From (95) and (96) we get 2

B(t + ∆t, x) =

1 ∆t X 1 B(t, x) + Y(t, x) − ∂xi (Ai (t + ∆t, x)Y(t, x)) 2 2 2 i=1



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

(97)

Finally from (93) and (97) we obtain Y(tn , x) = B(tn , x) − ∆t

2 X

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

(98)

i=1

B(t

n+1

, x) =

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

(99) 22

If we take the scalar product of (98) and (99) with a test function ϕh , integrate over a cell K and use the Green formula (34) then we get Z Z n n Y · ϕh dx = Bn · ϕh dx + ∆tFK (Bn , ϕh ) (100) K K  Z Z Z  ∆t n+1 n 1 n 1 n n+1 B + Y · ϕh dx + FK (Y , ϕh ) + (tn , x) · ϕh dx(101) B · ϕh dx = 2 2 2 K K K where k(tn , x)kl∞ (0,T ;L2 (Ω)) = O(∆t3 ). We set the following notations n = π Yn − Yn , δY h h

n = π Bn − B n , δB h h

enY = Yn − Yhn ,

n = π Y n − Y n , I n = π Bn − Bn , e n = Bn − B n . IY h h B B h

If we substract (36) to (100) and (37) to (101) we obtain Z Z n n n δY · ϕh dx = δB · ϕh dx + HK (ϕh ) K K  Z Z  1 n 1 1 n n+1 δY + δB · ϕh dx + JKn+1 (ϕh ) δB · ϕh dx = 2 2 2 K K where n HK (ϕh )

Z

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

=

(102) (103)

(104)

K

and JKn+1 (ϕh )

Z = K

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

We set Hn =

X

n HK 1K ,

and

J n+1 =

K∈Th

X

JKn+1 1K .

(106)

K∈Th

n in (102) and ϕ = δ n in (103), after summing over all the element K of the If we take ϕh = δB h Y partition Th the sum (102)/2+(103) gives

n+1

n+1

n n n n

δ

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

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

n+1

n n n

I − I + kIY − IB kL2 (Ω) ≤ C∆thm+1 |∂t B|L∞ (0,T ;Hm+1 (Ω)) (108) B B L2 (Ω) n ) the Using (108), the Cauchy-Schwarz and Young inequalities we have for the first term of Hn (δB estimate Z X Z n n n n n n ) · δB dx ≤ kIY − IB k2 kδB k2 (IY − IB K∈Th

K

Ω n 2 ≤ C∆th2m+2 + ∆t kδB kL2 (Ω) .

23

(109)

Now we estimate the term ∆tF n (Bn − πh Bn , ϕh ). X Z X Z ∆t Ani (Bn − πh Bn ) · ∂i ϕh dx =∆t (Ph Ani − Ani )(πh Bn − Bn ) · ∂i ϕh dx K

K∈Th

K

K∈Th

− ∆t

X Z K

K∈Th

= ∆t

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

X Z K∈Th

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

(110)

K

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 (110). Thanks the approximation properties (9) and (55), the inverse properties (11), the Cauchy-Schwarz and Young inequalities we obtain X Z ∆t (Ph Ani − Ani )(πh Bn − Bn ) · ∂i ϕh dx K∈Th

K

≤ ∆th

X Z K∈T

h X Z 



K∈Th K 2m+2

≤ Ch

kAi kL∞ (0,T ;W 1,∞ (Ω)) kπh Bn − Bn k2 k∂i ϕh k2 dx

K

 ∆tkAi k2L∞ (0,T ;W 1,∞ (Ω)) kπh Bn − Bn k22 + ∆th2 k∂i ϕh k22 dx

∆t|B|2L∞ (0,T ;Hm+1 (Ω)) + ∆tkϕh k2L2 (Ω)

≤ Ch2m+2 ∆t + ∆tkϕh k2L2 (Ω) .

(111)

Next we have X Z X Z n n n kCkL∞ (0,T ;L∞ (Ω)) kπh Bn − Bn k2 kϕh k2 dx ∆t C (πh B − B ) · ϕh dx ≤ ∆t K

K∈Th

K∈Th

X Z 



K∈Th K 2m+2

≤ Ch

K

 ∆tkCk2L∞ (0,T ;L∞ (Ω)) kπh Bn − Bn k22 + ∆tkϕh k22 dx

∆t + ∆tkϕh k2L2 (Ω) .

(112)

As it has been done for the continuous case (see estimates (59)-(60)) we obtain Z XZ ∆t X n n n 2m+1 |Ane |[ϕh ]e · [ϕh ]e dΓ. ∆t g (πh B − B )e · [ϕh ]e dΓ ≤ C∆th + 4 e e e∈Eh

(113)

e∈Eh

From (111)-(113) we get n

n

n

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

2m+1

+

3∆tkϕh k2L2 (Ω)

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

(114)

e∈Eh

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

(115)

From (114) we obtain n

n

∆tF (B − πh B

n

n , δB )

≤ C∆th

2m+1

+

n 2 3∆tkδB kL2 (Ω)

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

24

(116)

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

n

n n (δB , δB )



n 2 −α∆tkδB kL2 (Ω)

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

(117)

e∈Eh

From (104), (106), (109), (115)-(117) H

n

n (δB )

≤ C∆t(h

2m+1

+h

2m+2

) + (4 −

n 2 kL2 (Ω) α)∆tkδB

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

(118)

e∈Eh

n) Estimate of J n+1 (δY n) Using (108), the Cauchy-Schwarz and Young inequalities we have for the first term of J n+1 (δY the estimate Z X Z n+1 n+1 n n n n n n − IY − IB + (tn x)) · δY dx ≤ k2IB − IY − IB + (tn x))k2 kδY k2 dx (2IB K∈Th

K



Z ≤ Ω

 n n+1 n n n k2 + k(tn x))k2 kδY k2 dx k2(IB − IB )k2 + kIB − IY

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

(119)

n ) − ∆tF n+1 (Y n , δ n ). First we notice that Now we estimate the term ∆tF n+1 (Yn , δY h Y 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 ).

(120)

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

n+1

n

n

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

2m+1

+ ε∆tkϕh k2L2 (Ω)

Z ∆t X + |Ane |[ϕh ]e · [ϕh ]e dΓ. (121) 4 e e∈Eh

From (121) we obtain n

n

∆tF (Y − πh Y

n

n , δY )

≤ C(ε)∆th

2m+1

+

n 2 ε∆tkδY kL2 (Ω)

Z ∆t X n n + |Ane |[δY ]e · [δY ]e dΓ. (122) 4 e e∈Eh

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

n

n n (δY , δY )



n 2 kL2 (Ω) −α∆tkδY

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

(123)

e∈Eh

From (105), (106), (119), (120), (122) and (123) Z  ∆t X n n n n 2 |Ane |[δY ]e ·[δY ]e dΓ. J n+1 (δY ) ≤ C(ε)∆t (∆t2 + hm+1 )2 + h2m+1 ) +(4ε−α)∆tkδY kL2 (Ω) − 4 e e∈Eh

(124) n+1 kδB

nk δY L2 (Ω) .

Estimate of − First we notice that after summing over all elements K of the partition Th the sum (103)−(102)/2 gives Z  1 n+1 n (δB − δY ) · ϕh dx = J n+1 (ϕh ) − Hn (ϕh ) . (125) 2 Ω

25

where  1 1 J n+1 (ϕh ) − Hn (ϕh ) = 2 2

Z Ω

1 n+1 n n (2IB −IY −IB +(tn x))·ϕh dx− 2

Z

n n (IY −IB )·ϕh dx+Ln (ϕh ).



(126) with  Ln (ϕh ) = ∆t F n+1 (Yn , ϕh ) − F n+1 (Yhn , ϕh ) − F n (Bn , ϕh ) + F n (Bnh , ϕh ) .

(127)

From (108) and (109) we have n n |(IY − IB , ϕh )Ω | ≤ C∆th2m+2 + ∆tkϕh k2L2 (Ω)  n n 2I n+1 − IY − IB + (tn x), ϕh Ω ≤ C∆t(∆t2 + hm+1 )2 + ∆tkϕh k2L2 (Ω) . B

(128) (129)

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 )) .

(130)

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

ε kϕ k2 2 + C(ε)∆t2 h2m kBk2L∞ (0,T ;Hm+1 (Ω)) , 4 h L (Ω) ε kϕ k2 2 + C(ε)∆t2 h2m kYk2L∞ (0,T ;Hm+1 (Ω)) . 4 h L (Ω)

(131) (132)

n , ϕ )−∆tF n+1 (δ n , ϕ ). Now let us estimate the remaining term of Ln (ϕh ), that is to say ∆tF n+1 (δY h h B First we 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 )

(133)

where G(·, ϕh ) is defined by (79). From (87) we get ε ∆t2 n n n n 2 ∆tF n+1 (δB − δY , ϕh ) ≤ kϕh k2L2 (Ω) + C(ε) 2 kδB − δY kL2 (Ω) 4 h

(134)

and from (83) we get n |G(δB , ϕh )|

 ≤ C(ε)

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

(135)

n − δ n k2 n n Now let us estimate the term kδB Y L2 (Ω) which appears in (134). If we take ϕh = δB − δY in (102) we get n n 2 n n kδB − δY kL2 (Ω) = Hn (δY − δB ) (136)

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 ). (137) n − δ n ) the inequality The estimate of the three terms in the right hand side of (137) gives for Hn (δY B n n n n 2 Hn (δY − δB ) ≤ ∆tkδY − δB kL2 (Ω) + C∆th2m+2 n n 2 + C(ε)∆t2 h2m + εkδY − δB kL2 (Ω) n n 2 + εkδY − δB kL2 (Ω) + C(ε)

26

∆t2 n 2 kδ k 2 . h2 B L (Ω)

(138)

If ε and ∆t is small enough, (136) and (138) leads to   ∆t2 n 2 n n 2 2 2m 2m+2 kδY − δB kL2 (Ω) ≤ C(ε) ∆t h + ∆th + 2 kδB kL2 (Ω) . h

(139)

n+1 n in (125) and thanks to relation (126)-(135) and (139) we get By taking ϕh = δB − δY

 n+1 n 2 kδB − δY kL2 (Ω) ≤ C(ε)∆t h2m+2 + (∆t2 + hm+1 )2 + ∆th2m + ∆t3 h2m−2 + ∆t2 h2m   ∆t4 ∆t4 n+1 n 2 n 2 + 2(∆t + ε)kδB − δY kL2 (Ω) + C(ε) ∆t4 + 2 + 4 kδB kL2 (Ω)(140) . h h If ∆t and ε is small enough and ∆t ≤ C(ε)h4/3 then (140) leads to   n+1 n 2 kδB − δY kL2 (Ω) ≤ C(ε)∆t h2m+2 + (∆t2 + hm+1 )2 + h2m+4/3 + h2m+8/3   n 2 + C(ε) ∆t4 + ∆t5/2 + ∆t kδB kL2 (Ω)   n 2 kL2 (Ω) . ≤ C(ε)∆t (∆t2 + hm+1 )2 + h2m+4/3 + C(ε)∆tkδB

(141)

Finally the estimates (107), (118) , (124) and (141) leads to Z Z ∆t X ∆t X n+1 2 n n n n n kδB kL2 (Ω) ≤ − |Ae |[δB ]e · [δB ]e dΓ − |Ane |[δY ]e · [δY ]e dΓ 4 4 e e e∈Eh

+ (1 +

n 2 C(α, ε)∆t) kδB kL2 (Ω)

+ (4ε −

e∈Eh n 2 α)kδY kL2 (Ω)



 + C(ε)∆t (∆t2 + hm+1 )2 + h2m+1 + h2m+4/3 .

(142)

If we suppose ε small enough so that (4ε − α) ≤ 0 then (142) leads to n+1 2 n 2 kL2 (Ω) + C(ε)∆t (∆t2 + hm+1 )2 + h2m+1 kδB kL2 (Ω) ≤ (1 + C(α, ε)∆t) kδB



(143)

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

(144)

0 k2 0 0 2 2m+2 ). Finally, using (144) we get because kδB L2 (Ω) = kB − πh B kL2 (Ω) = O(h n 2 kBnh (x) − B(tn , x)k2L2 (Ω) ≤ kπh Bn (x) − B(tn , x)k2L2 (Ω) + kδB kL2 (Ω)  2 ≤ Ch2m+2 + C ∆t2 + hm+1/2 2  ≤ C ∆t2 + hm+1/2

and therefore

  kBnh (x) − B(tn , x)kL2 (Ω) ≤ C ∆t2 + hm+1/2 .

To end the proof we notice that (142) implies that Z  2 ∆t X n n n 2 |Ane |[δB ]e · [δB ]e dΓ ≤ (1 + C(α, ε)∆t) kδB kL2 (Ω) + C∆t ∆t2 + hm+1/2 4 e e∈Eh

27

(145)

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

NT X Z X

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

e

n=0 e∈Eh

n ), J n+1 (δ n ) and Now we suppose that Bh ∈ Wh1 . Then we see how change the estimates for Hn (δB Y n+1 nk kδB − δY L2 (Ω) . If we replace πh by Πh then the estimate (108) becomes

n+1

n n n

I − IB + kIY − IB kL2 (Ω) ≤ C∆th|∂t B|L∞ (0,T ;H1 (Ω)) B L2 (Ω) and therefore (109) becomes X Z K∈Th

K

(146)

n n n n 2 (IY − IB ) · δB dx ≤ C∆th2 + ∆t kδB kL2 (Ω)

(147)

Now we give a new estimate for the term ∆tF n (Bn − πh Bn , ϕh ). Using the Green formula (34) we obtain − ∆t

2 Z X i=1

Ani (Πh Bn

Z

n

− B ) · ∂i ϕh dx + ∆t



= ∆t



X Z K∈Th

− ∆t

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

K

2 X

∂i (Ani (Πh Bn − Bn )) + ∆tC n (Πh Bn − Bn )

· ϕh dx

(148)

i=1

X X Z K∈Th e∈∂K

!

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

(149)

e

First we look at the term (149). Using the approximation properties (18)-(20), inverse properties (11), the fact that div Πh Bn = div Bn = 0 on K, the Young and the Cauchy-Schwarz inequalities then we get X X Z ∆t Ane,K (Πh Bn − Bn ) · ϕh dΓ K∈Th e∈∂K



X Z K∈Th

∂K

e

2kukL∞ (0,T ;L∞ (Ω)) ∆t kΠh Bn − Bn k2 kϕh k2 dΓ

o X n ≤ 2kuk2L∞ (0,T ;L∞ (Ω)) ∆th−1 kΠh Bn − Bn k2L2 (∂K) + ∆th kϕh k2L2 (∂K) K∈Th

o X n 2 2 −1 n n 2 ≤ 2kukL∞ (0,T ;L∞ (Ω)) ∆th kΠh B − B kW(K) + ∆t kϕh kL2 (K) K∈Th



o   X n 2kuk2L∞ (0,T ;L∞ (Ω)) ∆th−1 kΠh Bn − Bn k2H(curl;K) + k(Πh Bn − Bn ) × νk2L2 (∂K) + ∆t kϕh k2L2 (K) K∈Th

 ≤ C T, kukL∞ (0,T ;L∞ (Ω)) , kBkL∞ (0,T ;H2 (Ω)) ∆th + ∆t kϕh k2L2 (Ω) .

28

(150)

As div Πh Bn = div Bn = 0 on K the term (148) can be recast in ! 2 X Z X ∆t ∂i (Ani (Πh Bn − Bn )) + C n (Πh Bn − Bn ) · ϕh dx K∈Th

K

i=1

X Z

= ∆t

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

K

K∈Th

X Z

= ∆t

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

K

K∈Th

(151) Now we have to estimate the three terms of (151). For the two first term of (151) we get X Z ∆t ((Πh Bn − Bn ) · ∇) un · ϕh dx K∈Th

K

X Z

≤ ∆t

K

K∈Th

2kukL∞ (0,T ;W 1,∞ (Ω)) kΠh Bn − Bn k2 kϕh k2 dx

o X n ≤ 4kuk2L∞ (0,T ;W 1,∞ (Ω)) ∆t kΠh Bn − Bn k2L2 (K) + ∆t kϕh k2L2 (K) K∈Th

 ≤ C T, kukL∞ (0,T ;W 1,∞ (Ω)) , |B|L∞ (0,T ;H1 (Ω)) h2 + ∆t kϕh k2L2 (Ω) .

(152)

and ∆t

X Z K

K∈Th

≤ ∆t

X Z K∈Th

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

K

kdivukL∞ (0,T ;L∞ (Ω)) kΠh Bn − Bn k2 kϕh k2 dx

o X n ≤ kdivuk2L∞ (0,T ;L∞ (Ω)) ∆t kΠh Bn − Bn k2L2 (K) + ∆t kϕh k2L2 (K) K∈Th

 ≤ C T, kdivukL∞ (0,T ;L∞ (Ω)) , |B|L∞ (0,T ;H1 (Ω)) h2 + ∆t kϕh k2L2 (Ω) .

(153)

Using the approximation properties (18)-(20), the continuous imbedding (28), the fact that div Πh B = div B = 0 on K and the Cauchy-Schwarz inequality then we get X Z ∆t (un · ∇)(Πh Bn − Bn ) · ϕh dxdt K∈Th

K

o X n ≤ 4kuk2L∞ (0,T ;L∞ (Ω)) ∆t kΠh Bn − Bn k2H1 (K) + ∆t kϕh k2L2 (K) K∈Th



C∆t kΠh Bn − Bn k2H(curl;K) +

K∈Th

)

!

( X

X

k(Πh Bn − Bn ) × ν e k2H1/2 (e)

+ ∆t kϕh k2L2 (K)

e∈∂K

 ≤ C T, kukL∞ (0,T ;L∞ (Ω)) , kBkL∞ (0,T ;H2 (Ω)) h∆t + ∆t kϕh k2L2 (Ω) . As it has been done for the continuous case (see estimates (59) and (64)) we obtain Z XZ ∆t X n n n g (Πh B − B )e · [ϕh ]e dΓ ≤ C∆th + ∆t |Ane |[ϕh ]e · [ϕh ]e dΓ. 4 e e e∈Eh

e∈Eh

29

(154)

(155)

From (148)-(155) we get ∆tF n (Bn − Πh Bn , ϕh ) ≤ C∆th + 3∆tkϕh k2L2 (Ω) +

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

(156)

e∈Eh

n) and finally, following what have been done when B ∈ Vhm , the new estimate for the term Hn (δB becomes Z ∆t X n n 2 n 2 n n H (δB ) ≤ C∆t(h + h ) + (4 − α)∆tkδB kL2 (Ω) − |Ane |[δB ]e · [δB ]e dΓ 4 e e∈Eh

n ). The same proof for the estimate (156) leads to Now we treat the term J n+1 (δY Z ∆t X n+1 n n 2 ∆tF (Y − Πh Y , ϕh ) ≤ C(ε)∆th + ε∆tkϕh kL2 (Ω) + |Ane |[ϕh ]e · [ϕh ]e dΓ. 4 e e∈Eh

n ) when B ∈ V m we obtain Following what have been done to obtain an estimate for J n+1 (δY h Z  ∆t X n n 2 n n J n+1 (δY ) ≤ C(ε)∆t (∆t2 + h)2 + h) + (4ε − α)∆tkδY kL2 (Ω) − |Ane |[δY ]e · [δY ]e dΓ. 4 e e∈Eh

n+1 nk Now we look at the term kδB − δY L2 (Ω) . From the approximation properties of Πh (see (18)) the estimates (128) and (129) become n n |(IY − IB , ϕh )Ω | ≤ C∆th2 + ∆tkϕh k2L2 (Ω)  n n n 2I n+1 − IY ≤ C∆t(∆t2 + h)2 + ∆tkϕh k2 2 , − I + (t x), ϕ h B L (Ω) B Ω

and the estimates (131) and (132) change into |∆tF n (πh Bn − Bn , ϕh )| ≤ ∆tF n+1 (Yn − πh Yn , ϕh ) ≤

ε kϕ k2 2 + C(ε)∆t2 kBk2L∞ (0,T ;H1 (Ω)) , 4 h L (Ω) ε kϕ k2 2 + C(ε)∆t2 kYk2L∞ (0,T ;H1 (Ω)) . 4 h L (Ω)

n+1 nk m Following what have been done to obtain an estimate for the term kδB − δY L2 (Ω) when B ∈ Vh we obtain   ∆t2 n 2 n n 2 2 2 kδY − δB kL2 (Ω) ≤ C(ε) ∆t + ∆th + 2 kδB kL2 (Ω) , h

  n+1 n 2 n 2 kδB − δY kL2 (Ω) ≤ C(ε)∆t (∆t2 + h)2 + h4/3 + C(ε)∆tkδB kL2 (Ω) , and n+1 2 kδB kL2 (Ω)

Z Z ∆t X ∆t X n n n n n ≤ − |Ae |[δB ]e · [δB ]e dΓ − |Ane |[δY ]e · [δY ]e dΓ 4 4 e e e∈Eh

n 2 C(α, ε)∆t) kδB kL2 (Ω)

+ (4ε −  + C(ε)∆t (∆t2 + h)2 + h + h4/3 . + (1 +

e∈Eh n 2 α)kδY kL2 (Ω)



If we suppose ε small enough so that (4ε − α) ≤ 0 then (157) leads to  n+1 2 n 2 kδB kL2 (Ω) ≤ (1 + C(α, ε)∆t) kδB kL2 (Ω) + C(ε)∆t (∆t2 + h)2 + h . 30

(157)

Then the end of the proof still remains the same and leads to the estimates   kBnh (x) − B(tn , x)kL2 (Ω) ≤ C ∆t2 + h1/2 . and ∆t

NT X Z X n=0 e∈Eh

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

e

which complete the proof.



References [1] S. Alinhac, P. G´erard, Op´erateurs pseudo-diff´erentiels et th´eor`eme de Nash-Moser, CNRS ´edition, 1991. [2] G. A. Baker, W. N. Jureidini, O. A. Karakashian, A piecewise solenoidal vector fields and the stokes problem, SIAM J. Numer. Anal., Vol 27 (1990), pp. 1466-1485. [3] D. S. Balsara, D. S. Spicer, A piecewise solenoidal vector fields and the stokes problem, J. Comput. Phys., 149 (1999), pp. 270-292. [4] J. U. Brackbill, D. C. Barnes The effect of nonzero ∇ · B on the numerical solution of the magnetohydrodynamic equations, J. Comput. Phys. 35, 426 (1980). [5] P. G. Ciarlet, Basic error estimates for elliptic problems, in Handbook of numerical analysis, 17-351, P.G. Ciarlet and J.L. Lions editors, North-Holland, 1991. [6] B. Cockburn, Discontinuous Galerkin methods for convection dominated problems, in Highorder methods for computational physics, 69-224, Lect. Notes Comput. Sci. Eng., 9, Springer, Berlin, 1999. [7] B. Cockburn, F. Li, C.-W. Shu, Locally divergence-free discontinuous Galerkin-methods for the Maxwell equations, J. of Comput. Phys. 194, 588-610, (2004). [8] M. Costabel, 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, t. 327, S´erie I, p. 849-854, (1998). [9] M. Costabel, A remark on the regularity of solutions of Maxwell’s equations on Lipschitz domains, Math. Methods Appl. Sci. 12 (4) (1990), 365-368. [10] W. Dai, P. R. Woodward, On the divergence-free condition and conservation laws in numerical simulations for supersonic magnetohydrodynamic flow, Astrophys. J. 494, 317 (1998). [11] A. Dedner, F. Kemm, D. Kr¨oner, C.-D. Munz, T. Schnitzer, M. Wesenberg, Hyperbolic Divergence cleaning for the MHD equations, J. Comput. Phys. 175, 645-673 (2002). [12] C. R. Evans, J. F. Hawley, Simulation of magnetohydrodynamic flows, a constrained transport method, Astrophys. J. 332, 659 (1988). [13] K. O. Friedrichs, Symmetric positive linear differential equations, Comm. on pure and appl. math., VOL. XI, (1958), 333-418.

31

[14] V. Gilrault, P.-A. Raviart, Finite element methods for the Navier-Stokes equatons, Theory and algorithms, Springer Series in Comput. Math. 5, Springer-Verlag, Berlin, 1986. [15] P. Grisvard, Elliptic problems in nonsmooth domains, Monographs and studies in Mathematics, 24 Pitman, London, 1985. [16] O. A. Karakashian, W. N. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations, SIAM J. Numer. Anal., Vol 35, No 1 (1998), pp. 93-120. [17] F. Li, C.-W. Shu, Locally divergence-free discontinuous Galerkin methods for MHD equations, SIAM J. Sci. Comput., to appear. [18] J. C. N´ed´elec, Mixed finite element in

R3 , Numer. Math. 35, 315-341 (1980)

[19] K. G. Powell, An approximate Riemann solver for magnetohydrodynamics (that works in more than one dimension), ICASE-Report 94-24 (NASA CR-194902) (NASA Langley Research Center, Hampton, VA 23681-0001, 8 april 1994). [20] P.-A. Raviart, J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, In A. Dold, B. Eckmann (eds). Mathematical aspects of finite elements methods. Proceedings of the conference held in Rome, 10-12 Dec, 1975. Springer, Berlin Heidelberg New York (Lectures Notes in Mathematics vol 606) [21] G. T´oth, The ∇·B = 0 constraint in shock-capturing magnetohydrodynamics codes, J. Comput. Phys. 161, 605 (2000). [22] J.-P. Vila, P. Villedieu, Convergence of an explicit finite volume scheme for first order symmetric systems, Numer. Math. (2003) 94, pp. 573-602. [23] Q. Zhang, C.-W. Shu, Error Estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws , SIAM J. Numer. Anal., to appear.

32