Error Estimates of Runge-Kutta Discontinuous Galerkin Methods for ...

Error Estimates of Runge-Kutta Discontinuous Galerkin Methods for the Vlasov-Maxwell System He Yang and Fengyan Li

∗

June 14, 2014

Abstract In this paper, error analysis is established for Runge-Kutta discontinuous Galerkin (RKDG) methods to solve the Vlasov-Maxwell system. This nonlinear hyperbolic system describes the time evolution of collisionless plasma particles of a single species under the self-consistent electromagnetic field, and it models many phenomena in both laboratory and astrophysical plasmas. The methods involve a third order TVD Runge-Kutta discretization in time and upwind discontinuous Galerkin discretizations of arbitrary order in phase domain. With the assumption that the exact solutions have sufficient regularity, the L2 errors of the particle number density function as well as electric and magnetic fields at any given 1 time T are bounded by Chk+ 2 + Cτ 3 under a CFL condition τ /h ≤ γ. Here k is the polynomial degree used in phase space discretization, satisfying k > dx2+1 (with dx being the dimension of spatial domain), τ is the time step, and h is the maximum mesh size in phase space. Both C and γ are positive constants independent of h and τ , and they may depend on the polynomial degree k, time T , the size of the phase domain, certain mesh parameters, and some Sobolev norms of the exact solution. The analysis can be extended to RKDG methods with other numerical fluxes and to RKDG methods solving relativistic Vlasov-Maxwell equations.

Keywords: Vlasov-Maxwell system, Runge-Kutta discontinuous Galerkin methods, error estimates AMS(MOS) subject classification: 65M15, 65M60 , 65M06, 35Q83, 35L50

1

Introduction

In this paper, we will establish error estimates of the Runge-Kutta discontinuous Galerkin (RKDG) methods for solving the dimensionless Vlasov-Maxwell (VM) equations    ∂t f + v · ∇x f + (E + v × B) · ∇v f = 0, ∂t E = ∇ × B − J, ∂t B = −∇ × E, (1.1)   ∇ · E = ρ − ρi , ∇ · B = 0, with

Z ρ(x, t) =

Z f (x, v, t)dv,

J(x, t) =

Ωv

f (x, v, t)vdv.

(1.2)

Ωv

This system describes the time evolution of collisionless plasma particles of a single species, such as electrons or ions, under the self-consistent electromagnetic field. Here f (x, v, t) ≥ 0 is the particle number density function in the phase space with (x, v) at time t, E(x, t) is the electric field, B(x, t) is the magnetic field, J(x, t) is the current density, ρ(x, t) is the charge density, and ρi is the charge density of the background particles. The system (1.1) is defined on the phase domain Ω = Ωx × Ωv , where Ωx = [Lx,1 , Lx,2 ]dx is the spatial domain and Ωv = [Lv,1 , Lv,2 ]dv is the velocity domain (dx , dv = 1, 2 or 3), with periodic boundary ∗ This research is partially supported by NSF CAREER grant DMS-0847241 and NSF grant DMS-1318409. Address: Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, United States. Email: [email protected], [email protected].

1

conditions in x. In v Rdirection, f is assumed to have compact support. We further assume the VM system is globally neutral, i.e. Ωx (ρ − ρi )dx = 0. Note that this is compatible with the periodic boundary conditions in x. The VM equations model many phenomena in both laboratory and astrophysical plasmas, and accurate and reliable numerical simulation of this system has fundamental importance. Particle methods [3, 10, 15] have been widely used since 60’s because of their low computational cost especially when the dimension of the phase space is high, yet with their numerical noise, it is hard for the methods to produce very accurate approximations. In recent years, many high order Eulerian methods have been developed in the context of Vlasov-Poisson equations or the VM equations. Some examples include semi-Lagrangian methods [6, 20, 17, 18], continuous finite element methods [24, 25], and methods based on Fourier transform [16, 12, 13]. In [5], one of the authors and her collaborators proposed and analyzed semi-discrete discontinuous Galerkin (DG) methods for the VM system. The methods were further combined with Runge-Kutta time discretizations, resulting in Runge-Kutta discontinuous Galerkin (RKDG) methods, and their performance in accuracy, stability, and conservation was demonstrated numerically. Note that DG methods were previously proposed and studied in [14, 1, 2, 4] for the Vlasov-Poisson system. DG discretizations are chosen for the phase domain in [5] due to their high accuracy, compactness, high efficiency in parallel implementations, flexibility with complicated geometry as well as boundary conditions and adaptive simulations. All aforementioned properties make the methods a competitive candidate to simulate the VM system accurately with reasonable computational cost especially for lower dimensional cases (e.g. the 1D2V or 2D2V system). This is even so if one further makes good use of the modern computer architectures. More importantly, with properly designed numerical fluxes and up to some boundary effect, semi-discrete DG methods have provable conservation property of both mass and total energy, and this is shared by very few high order methods. On the other hand, though RKDG methods have been widely used in many applications since they were introduced [9, 8], theoretical analysis for such fully discrete methods is relatively little. Error estimations based on Fourier analysis and some symbolic computations were carried out in [29] and [23] when RKDG methods are applied to linear problems on uniform meshes. Important developments were made by Zhang and Shu in [26, 27, 28], where error estimates to the smooth solutions on uniform or non-uniform meshes were developed for RKDG methods with the second order RK time discretization for scaler [26] and symmetrizable conservation laws [27], and with the third order RK method for scaler conservations laws in [28]. In [28], L2 -norm stability was also obtained for linear conservation laws, and such analysis so far is unavailable for nonlinear cases. In this paper, we use the idea in [28] to obtain the error estimates of the fully discrete RKDG methods for the VM system when the exact solutions have sufficient regularity. In particular, a third order TVD Runge-Kutta method [19], a commonly used Runge-Kutta method in RKDG schemes, is considered as time discretization together with DG discretizations of arbitrary accuracy in phase domain. The second order Runge-Kutta time discretization analyzed in [26] is not examined here partially due to their restrictive condition on timestep when the spatial accuracy is higher than second order [26, 27]. Our analysis is based on Taylor expansion, energy analysis, some techniques for analyzing the semi-discrete DG methods of the VM system [5] and for analyzing the third order Runge-Kutta time discretization within the method of lines framework. To treat the nonlinearity due to the nonlinear coupling of the Vlasov and Maxwell parts, the polynomial degree k is required to satisfy k > dx2+1 . In addition, a priori assumption is made for the L∞ error of the electric and magnetic fields. This assumption will be proved later by mathematical induction. Even though the Gauss’s laws are not considered in formulating the numerical methods, they are shown to be approximated accurately when the exact solutions have sufficient regularty. The remaining of this paper is organized as follows. In section 2, the formulations of RKDG methods are presented for the VM system. We also introduce notations and review some standard approximation results and inverse inequalities in finite element methods. In section 3, error estimates are established for the RKDG methods. Here we start with the error equations and energy equations. Based on these equations, the errors from the Vlasov and Maxwell parts are estimated and then combined. To better present the analysis, the proofs of some lemmas are given in section 4. Finally, we summarize and generalize our work in section 5.

2

2

Formulation of Runge-Kutta Discontinuous Galerkin Methods

In this section, we will introduce notations, review some standard approximation results and inverse inequalities, and present the formulation of Runge-Kutta discontinuous Galerkin methods for the Vlasov-Maxwell system (1.1).

2.1

Some preliminaries

Throughout this paper, standard notations are used for Sobolev spaces and norms: for a bounded domain D and any nonnegative integer m, we denote the L2 -Sobolev space of order m by H m (D) equipped with Sobolev norm || · ||m,D , and the L∞ -Sobolev space of order m by W m,∞ (D) with the Sobolev norm || · ||m,∞,D . When m = 0, L2 (D) is used instead of H 0 (D), so is L∞ (D) instead of W 0,∞ (D). For the brevity of notation, we use ? = x or v in this subsection. For the computational domain Ω = Ωx × Ωv , assume Th? = {K? } is a partition of Ω? , with K? being a (rotated) Cartesian element or a simplex, then Th = {K : K = Kx × Kv , ∀Kx ∈ Thx , ∀Kv ∈ Thv } defines a partition of Ω. Let E? be the set of the edges of Th? , then the edges of Th will be E = {Thx × Ev } ∪ {Ex × Thv }. In addition, let Ev = Evi ∪ Evb with Evi (resp. Evb ) consisting of all interior (resp. boundary) edges of Thv . The mesh size of Th is denoted as h = max(hx , hv ) = maxK∈Th hK , where h? = maxK? ∈Th? hK? with hK? = diam(K? ), and hK = max(hKx , hKv ) with K = Kx × Kv . When the mesh hv hx := minK h∈Tx v hK and hx,min := minK h∈Tv x hK are uniformly bounded above is refined, we assume both hv,min v

h

v

x

h

x

by a positive constant σ0 . Therefore in our analysis we do not always distinguish h, hx , hv , hKx and hKv . It is further assumed that {Th? }h is shape-regular. That is, if ρK? denotes the diameter of the largest sphere hK? ≤ σ? , ∀K? ∈ Th? , for a positive constant σ? independent of h? . included in K? , there is ρK ? Next we introduce two finite dimensional discrete spaces Ghk

=

{g ∈ L2 (Ω) : g|K ∈ P k (K), ∀K ∈ Th },

Uhk

=

{U ∈ [L2 (Ωx )]3 : U |Kx ∈ [P k (Kx )]3 , ∀Kx ∈ Thx },

where P k (D) is the set of polynomials of the total degree at most k on D, with k being any nonnegative integer. Note that functions in Ghk (resp. Uhk ) are piecewise defined with respect to Th (resp. Thx ). For such function, we would need the notations of jumps and averages. Given an edge e = (Kx+ ∩ Kx− ) ∈ Ex with n± x as the outward unit normal vector of Kx± , for any g ∈ Ghk and U ∈ Uhk with g ± = g|Kx± and U ± = U |Kx± , the averages of g and U across e are {g}x =

1 + (g + g − ), 2

{U }x =

1 + (U + U − ), 2

and the jumps are − − [g]x = g + n+ x + g nx ,

− − [U ]x = U + · n+ x + U · nx ,

− − [U ]tan = U + × n+ x + U × nx .

The averages and jumps across any interior edge e = (Kv+ ∩ Kv− ) ∈ Evi can be defined similarly. For a boundary edge in Evb with nv being the outward unit normal vector, we set [g]v = gnv and {g}v = 12 g. This is consistent with the exact solution f being compactly supported in Ωv . Below are some equalities which will be frequently used in the analysis and can be easily verified, 1 2 [g ]? = {g}? [g]? , 2 [g1 g2 ]x − {g1 }x [g2 ]x − {g2 }x [g1 ]x = 0,

(2.3a) (2.3b)

[U × V ]x + {V }x · [U ]tan − {U }x · [V ]tan = 0,

(2.3c) R K∈Th K ,

where g, g1 , g2 ∈ Ghk , and U, V ∈ Uhk . We also introduce some shorthand notations, Ω = Th = R R R R R P P 1 2 2 2 with = T? = K? ∈Th? K? , E? = e∈E? e . Additionally, ||g||0,E = (||g||0,Ex ×Thv + ||g||0,Thx ×Ev ) Ω? h R R  21 R R  12 R 1 ||g||0,Ex ×Thv = Ex T v g 2 dvdsx , ||g||0,Thx ×Ev = T x Ev g 2 dsv dx , and ||g||0,Ex = ( Ex g 2 dsx ) 2 . R

h

h

3

R

P

Let Πk denote the L2 projection onto Ghk , and Πkx be the L2 projection onto Uhk . In this paper, the following approximation properties in (2.4) and inverse inequalities in (2.5) will be used: there exists a constant C > 0, such that ∀g ∈ H k+1 (Ω), ∀U ∈ [H k+1 (Ωx )]3 ,  1 k+1 k k 2   ||g − Π g||0,K + hK ||g − Π g||0,∂K ≤ ChK ||g||k+1,K , 

∀K ∈ Th ,

1 2

||U − Πkx U ||0,Kx + hKx ||U − Πkx U ||1,Kx + hKx ||U − Πkx U ||0,∂Kx ≤ Chk+1 Kx ||U ||k+1,Kx ,    k+1 k x ||U − Πx U ||0,∞,Kx ≤ ChKx ||U ||k+1,∞,Kx , ∀Kx ∈ Th .

∀Kx ∈ Thx , (2.4)

In addition, there exists a constant C > 0, such that ∀g ∈ P k (K), ∀U ∈ [P k (Kx )]3 , ||∇x g||0,K ≤ Ch−1 Kx ||g||0,K ,

(

− dx

||U ||0,∞,Kx ≤ ChKx2 ||U ||0,Kx ,

||∇v g||0,K ≤ Ch−1 Kv ||g||0,K , −1

(2.5)

||U ||0,∂Kx ≤ ChKx2 ||U ||0,Kx .

Each positive constant C in (2.4) and (2.5) is independent of the mesh sizes hKx and hKv , and it depends on k and the shape regularity parameters σx and (or) σv of the mesh. One can refer to [7] for more details of such standard results. Throughout the paper, τ is used to denote the time step and tn = nτ . Without loss of generality, we assume the time steps are uniform and τ , h ≤ 1. The analysis in this paper also holds for non-uniform time steps. Even though the numerical methods and error analysis will be presented when both the (exact and approximated) electric and magnetic fields have three components, they can be easily adapted to reduced VM equations, such as the one to study Weibel instability in [5] when dx = 1 and dv = 2, where some components of E and B vanish and need not be approximated numerically.

2.2

Runge-Kutta discontinuous Galerkin methods

Now we are ready to present the RKDG methods for the VM system, where upwind DG methods of arbitrary accuracy are used as the spatial and phase discretization and a third order TVD Runge-Kutta method [19] is used as the time discretization. Note that on the PDE level, the two Gauss’s laws in (1.1) involving the divergence of the magnetic and electric fields can be derived from the remaining equations of the VM system as long as they are satisfied by the initial data, these equations will not be discretized numerically just as in [5]. For sufficiently smooth solutions as considered in this work, the divergence equations can be approximated accurately by the proposed methods (see Theorem 3.1). One should be aware that for general cases, imposing divergence equations in numerical simulations can be important. To initialize the simulation, let fh0 = Πk f0 , Eh0 = Πkx E0 and Bh0 = Πkx B0 , where f0 , E0 and B0 are the initial data of the VM system. Then for n ≥ 0, the approximate solutions at time tn+1 = (n + 1)τ are defined as follows. We look for fhn,1 , fhn,2 , fhn+1 ∈ Ghk , and Ehn,1 , Ehn,2 , Ehn+1 , Bhn,1 , Bhn,2 , Bhn+1 ∈ Uhk satisfying (fhn,1 , g)Ω = (fhn , g)Ω + τ ah (fhn , Ehn , Bhn ; g), (Ehn,1 , U )Ωx

(Bhn,1 , V )Ωx

(Ehn , U )Ωx

(2.6a)

(Bhn , V )Ωx

τ bh (Ehn , Bhn , fhn ; U, V

+ = + + ), 1 τ 3 (2.6b) (fhn,2 , g)Ω = ( fhn + fhn,1 , g)Ω + ah (fhn,1 , Ehn,1 , Bhn,1 ; g), 4 4 4 3 1 3 1 τ (Ehn,2 , U )Ωx + (Bhn,2 , V )Ωx = ( Ehn + Ehn,1 , U )Ωx + ( Bhn + Bhn,1 , V )Ωx + bh (Ehn,1 , Bhn,1 , fhn,1 ; U, V ), 4 4 4 4 4 1 n 2 n,2 2τ n,2 n,2 n,2 n+1 (fh , g)Ω = ( fh + fh , g)Ω + ah (fh , Eh , Bh ; g), (2.6c) 3 3 3 1 2 1 2 2τ (Ehn+1 , U )Ωx + (Bhn+1 , V )Ωx = ( Ehn + Ehn,2 , U )Ωx + ( Bhn + Bhn,2 , V )Ωx + bh (Ehn,2 , Bhn,2 , fhn,2 ; U, V ), 3 3 3 3 3

for any g ∈ Ghk and U, V ∈ Uhk , where

4

Z ah (fh , Eh , Bh ; g) =

fh v · ∇x g + fh (Eh + v × Bh ) · ∇v gdxdv Z Z Z X − fh\ v · nx gdsx dv + Ω

K=Kx ×Kv ∈Th

Kv

∂Kx

Kx

Z Bh · ∇ × U − Eh · ∇ × V dx +

bh (Eh , Bh , fh ; U, V ) = Ωx

X Z Kx ∈Thx

Z −

Z



 fh (Eh +\ v × Bh ) · nv gdsv dx ,

∂Kv

 nx\ × Bh · U − nx\ × Eh · V dsx

∂Kx

Z Jh · U dx,

Ωx

Jh (x, t) =

fh (x, v, t)vdv. Ωv

Here nx and nv are outward unit normal vectors of ∂Kx and ∂Kv , respectively. All the hat functions are upwinding numerical fluxes defined as   |v · nx | [fh ]x · nx , fh\ v · nx = {fh v}x + 2   |(Eh + v × Bh ) · nv | \ [fh ]v · nv , fh (Eh + v × Bh ) · nv = {fh (Eh + v × Bh )}v + 2     1 1 \ \ nx × Eh = nx × {Eh }x + [Bh ]tan , nx × Bh = nx × {Bh }x − [Eh ]tan , 2 2 and they further specify ah (fh , Eh , Bh ; g) = ah,1 (fh ; g) + ah,2 (fh , Eh , Bh ; g) with  Z Z Z  |v · nx | [fh ]x · [g]x dsx dv, ah,1 (fh ; g) = fh v · ∇x gdxdv − {fh v}x + 2 Ω Thv Ex Z ah,2 (fh , Eh , Bh ; g) = fh (Eh + v × Bh ) · ∇v gdxdv Ω  Z Z  |(Eh + v × Bh ) · nv | {fh (Eh + v × Bh )}v + − [fh ]v · [g]v dsv dx, 2 Thx Ev and Z bh (Eh , Bh , fh ; U, V ) =

Z (Bh · ∇ × U − Eh · ∇ × V ) dx − Jh U dx Ωx Ωx    Z  1 1 + {Bh }x − [Eh ]tan · [U ]tan − {Eh }x + [Bh ]tan · [V ]tan dsx . 2 2 Ex

Note that both ah,1 and bh are linear with respect to each argument, yet ah,2 (fh , Eh , Bh ; g) is linear with respect to fh and g only. The overall RKDG methods are consistent, and this will be used to derive the error equations in next section.

3

Error Estimates

This section is devoted to the main result of the paper, which is given in Theorem 3.1. More specifically, we will establish error estimates at any given time T > 0 for the fully discrete RKDG methods in section 2.2 when they are used to solve sufficiently smooth solutions. The analysis of the present work is based on the following assumptions on the regularity of the exact solutions. Regularity Assumption (a.1) f (·, ·, t), E(·, t), B(·, t) ∈ C 4 ([0, T ]); (a.2) ∂t4 f (·, ·, t) ∈ L2 (Ω), ∂t4 E(·, t), ∂t4 B(·, t) ∈ L2 (Ωx ) uniformly in t ∈ [0, T ]; 5

(a.3) (∂ts E+v×∂ts B)∇v ((∂t E+v×∂t B)·∇v (∂t f )), ∇x ((∂t E+v×∂t B)·∇v (∂t f )), ∂t ((∂t E+v×∂t B)∂t (∇v f )), (∂t2 E + v × ∂t2 B)∂t2 (∇v f ) ∈ L2 (Ω), ∂ts E, ∂ts B ∈ W 1,∞ (Ωx ), ∇v (∂ts f ), ∇v ((∂t E + v × ∂t B) · ∇v (∂t f )) ∈ L∞ (Ω) uniformly in t for s = 0, 1, 2; (a.4) ∂ts E(·, t), ∂ts B(·, t) ∈ H k+1 (Ωx ), ∂ts f (·, ·, t) ∈ H k+1 (Ω), (∂t E + v × ∂t B) · ∇v (∂t f ) ∈ H k+1 (Ω) uniformly in t ∈ [0, T ], s = 0, 1, 2. These, together with equations (3.8), will ensure f n,s ∈ H k+1 (Ω), E n,s , B n,s ∈ H k+1 (Ωx ) uniformly in n for s = 0, 1, 2. The assumption (a.3) can be replaced by the following, (a.30 ) All the mixed partial derivatives (in t, x, and v) of E and B up to the second order, in addition to ∂t3 E and ∂t3 B, are in L∞ (Ωx ) uniformly in t ∈ [0, T ]; ∇v (∂ts f ), s = 0, 1, 2 and ∇v (∇v ∂t f ) are in L∞ (Ω) uniformly in t ∈ [0, T ]; All the mixed partial derivatives (in t, x, and v) of f up to the third order are in L2 (Ω) uniformly in t ∈ [0, T ]. The assumption on L2 norms of functions together with (a.1) is for the third order accuracy in time discretization, while (a.4) is for the accuracy in the spatial and/or phase domain. The assumption on L∞ norms of functions is used to deal with the nonlinear coupling in the Vlasov part. Unless otherwise specified, C is used to denote a generic positive constant, and it can take different values at different occurrences. This constant is independent of n, h, τ , and may depend on polynomial degree k, mesh parameter σ0 , σx , σv , domain parameters Lx,i , Lv,i , i = 1, 2, and the time T . It may also depend on the exact solution in the form of its certain Sobolev norms or semi-norms. The constant γ1 in Theorem 3.10, γ2 in Theorem 3.16, and γ in Theorem 3.1 have similar dependence as the generic constant C. For convenience, we do not distinguish the upper indices n, 0 and n. For instance, we regard g n,0 = g n for any function g. With the analysis being very technical, to make it easier to follow, the proofs of some lemmas are given later in section 4. Theorem 3.1. Let (f, E, B) be a sufficiently smooth exact solution to the VM system (1.1). Let (fhn , Ehn , Bhn ) ∈ Ghk × Uhk × Uhk be the solution to the scheme (2.6) at time tn with k > dx2+1 . Under a CFL condition τ ≤ γh, with h < h0 for some h0 > 0, we have ||f (·, ·, tn ) − fhn ||20,Ω + ||E(·, tn ) − Ehn ||20,Ωx + ||B(·, tn ) − Bhn ||20,Ωx ≤ Ch2k+1 + Cτ 6 for any n ≤ T /τ . In addition, ||E(·, tm ) − Ehm ||0,∞,Ωx ≤ Ch,

||B(·, tm ) − Bhm ||0,∞,Ωx ≤ Ch,

∀m + 1 ≤ T /τ,

and the Gauss’s laws in (1.1) are approximated as below for any n ≤ T /τ , ||∇ · Ehn − ρnh + ρni ||0,Ωx ≤ C Here ρnh (·) =

3.1

R Ωv

1 k+ 1 (h 2 + τ 3 ), hx

||∇ · Bhn ||0,Ωx ≤ C

1 k+ 1 (h 2 + τ 3 ). hx

fhn (·, v)dv, and ρni (·) = ρi (·, tn ).

Error equations

Let (f (x, v, tn ), E(x, tn ), B(x, tn )) be the exact solution to the system (1.1)Rand (1.2) at time tn = nτ . We denote f n (x, v) = f (x, v, tn ), E n (x) = E(x, tn ), B n (x) = B(x, tn ), J n (x) = Ωv f n (x, v)vdv, and define  n,1 f = f n − τ (v · ∇x f n + (E n + v × B n ) · ∇v f n ) ,      E n,1 = E n + τ (∇ × B n − J n ), B n,1 = B n − τ (∇ × E n ),     
  f n,2 = 3 f n + 1 f n,1 − τ v · ∇x f n,1 + (E n,1 + v × B n,1 ) · ∇v f n,1 , 4 4 Z4 Z   n,1 n,1 n,2  J = f vdv, J = f n,2 vdv,    Ωv Ωv       E n,2 = 3 E n + 1 E n,1 + τ (∇ × B n,1 − J n,1 ), B n,2 = 3 B n + 1 B n,1 − τ (∇ × E n,1 ). 4 4 4 4 4 4 6

(3.7)

Equations in (3.7) are obtained by applying one step of the third order Runge-Kutta time discretization in [19] to the VM system (without the divergence conditions) from t = tn with the exact solution as the initial data at tn . From (1.1) and (3.7), one can further represent f n,] , E n,] and B n,] (] = 1, 2) in terms of f n , E n , B n and their derivatives as below,  n,1 f = f n + τ ∂t f n , E n,1 = E n + τ ∂t E n , B n,1 = B n + τ ∂t B n ,      n,2 τ2 τ3 τ f = f n + ∂t f n + ∂t2 f n − (∂t E n + v × ∂t B n ) · ∇v (∂t f n ), (3.8) 2 4 4   2 2    E n,2 = E n + τ ∂ E n + τ ∂ 2 E n , B n,2 = B n + τ ∂ B n + τ ∂ 2 B n . t t 2 4 t 2 4 t In the next lemma, local truncation errors from each step of the Runge-Kutta time discretization are given. The results can be verified straightforwardly based on (3.8) and Taylor expansion, and the proof is omitted. Lemma 3.2. If we define f (x, v, tn+1 )

=

E(x, tn+1 )

=

B(x, tn+1 )

=


1 n 2 n,2 2τ f + f − v · ∇x f n,2 + (E n,2 + v × B n,2 ) · ∇v f n,2 + Tfn (x, v), 3 3 3 1 n 2 n,2 2τ E + E + (∇ × B n,2 − J n,2 ) + TEn (x), 3 3 3 1 n 2 n,2 2τ B + B − (∇ × E n,2 ) + TBn (x), 3 3 3

(3.9)

where Tfn (x, v), TBn (x), TEn (x) are the local truncation errors in the n-th time step ∀n : (n + 1)τ ≤ T , then ||Tfn ||0,Ω , ||TEn ||0,Ωx , ||TBn ||0,Ωx ≤ Cτ 4 . n,] For each stage of the Runge-Kutta method, we denote the error by en,] − fhn,] = ξfn,] − ηfn,] , where f =f

ξfn,] = Πk f n,] − fhn,] and ηfn,] = Πk f n,] − f n,] , for ] = 0, 1, 2. Given that ηfn,] can be estimated in a standard way from the approximation results in (2.4) of the discrete spaces, our error analysis will focus on the term ξfn,] , which is also called the projected error due to ξfn,] = Πk en,] f . Similar comments and conventions on n,] n,] n,] n,] n,] k notation go to en,] E , eB , ξE , ξB , ηE and ηB . Next we multiply an arbitrary test function g ∈ Gh (resp. k U, V ∈ Uh ) on both sides of the equations corresponding to the Vlasov equation (resp. Maxwell equations) in (3.7) and (3.9), integrate over each mesh element K (resp. Kx ), take integration by parts, and sum up with respect to all K ∈ Th (or all Kx ∈ Thx ). We then subtract (2.6a) - (2.6c) from the resulting equations, and reach the error equations,  n,1 (ξf , g)Ω = (ξfn , g)Ω + τ J (g),      n,2 1 τ 3 (ξf , g)Ω = ( ξfn + ξfn,1 , g)Ω + K(g), (3.10) 4 4 4    1 2 2τ   (ξ n+1 , g)Ω = ( ξ n + ξ n,2 , g)Ω + L(g), f 3 f 3 f 3

 n,1 n,1 n n (ξE , U )Ωx + (ξB , V )Ωx = (ξE , U )Ωx + (ξB , V )Ωx + τ Q(U, V ),      n,2 3 n 1 n,1 3 n 1 n,1 τ n,2 (ξE , U )Ωx + (ξB , V )Ωx = ( ξE + ξE , U )Ωx + ( ξB + ξB , V )Ωx + R(U, V ), 4 4 4 4 4    1 2 1 2  n,2 n,2  (ξ n+1 , U )Ω + (ξ n+1 , V )Ω = ( ξ n + ξ , U )Ω + ( ξ n + ξ , V )Ω + 2τ S(U, V ). x x x x E B 3 E 3 E 3 B 3 B 3

7

(3.11)

Here J (g)

K(g)

L(g)

ηfn,1 − ηfn

=

τ

! + ah (f n , E n , B n ; g) − ah (fhn , Ehn , Bhn ; g),

,g Ω

4ηfn,2 − 3ηfn − ηfn,1

=

τ

! + ah (f n,1 , E n,1 , B n,1 ; g) − ah (fhn,1 , Ehn,1 , Bhn,1 ; g),

,g Ω

3ηfn+1 − ηfn − 2ηfn,2 + 3Tfn (x, v)

=

2τ

Q(U, V )

=

R(U, V )

=

S(U, V )

=

! + ah (f n,2 , E n,2 , B n,2 ; g) − ah (fhn,2 , Ehn,2 , Bhn,2 ; g),

,g Ω

! n,1 n ηB − ηB + ,V + bh (enE , enB , enf ; U, V ), τ Ωx Ωx ! ! n,2 n,1 n,2 n,1 n n 4ηB − 3ηB − ηB 4ηE − 3ηE − ηE n,1 n,1 + + bh (en,1 ,U ,V E , eB , ef ; U, V ), τ τ Ωx Ωx ! ! n,2 n,2 n+1 n+1 n n n 3ηE − ηE − 2ηE + 3TE (x) 3ηB − ηB − 2ηB + 3TBn (x) ,U ,V + 2τ 2τ n,1 n ηE − ηE ,U τ

!

Ωx

n,2 n,2 +bh (en,2 E , eB , ef ; U, V

Ωx

),

with any test functions g ∈ Ghk and U, V ∈ Uhk . For the functional J (·), we denote J1 (g) =



ηfn,1 −ηfn ,g τ

 Ω

and J2 (g) = ah (f n , E n , B n ; g) − ah (fhn , Ehn , Bhn ; g). Similarly, one can define K] , L] , Q] , R] , S] , ] = 1, 2. We now take the test function g = ξfn , 4ξfn,1 and 6ξfn,2 in each equation of (3.10), respectively, sum them up and obtain the energy equations 3||ξfn+1 ||20,Ω − 3||ξfn ||20,Ω

= τ [J (ξfn ) + K(ξfn,1 ) + 4L(ξfn,2 )] +

||2ξfn,2

−

ξfn,1

−

ξfn ||20,Ω

+

3(ξfn+1

(3.12) −

ξfn , ξfn+1

−

2ξfn,2

+

ξfn )Ω .

Similarly, the following equation holds n+1 2 n+1 2 n 2 n 2 3(||ξE ||0,Ωx + ||ξB ||0,Ωx ) − 3(||ξE ||0,Ωx + ||ξB ||0,Ωx )

=

n,1 n,1 n,2 n,2 n,2 n,1 n,2 n n n 2 τ [Q(ξE , ξB ) + R(ξE , ξB ) + 4S(ξE , ξB )] + ||2ξE − ξE − ξE ||0,Ωx + ||2ξB n,2 n,2 n+1 n+1 n+1 n+1 n n n n +3(ξE − ξE , ξE − 2ξE + ξE )Ωx + 3(ξB − ξB , ξB − 2ξB + ξB )Ωx .

(3.13) −

n,1 ξB

−

n 2 ξB ||0,Ωx

The main error estimate will be established based on equations (3.12) - (3.13) which describe how the L2 norms of the projected errors are accumulated in one time step. In particular, in the next two subsections, we will estimate the errors from the Vlasov and Maxwell solvers, respectively, and the results will be combined in section 3.4 to get the main result of this paper. Before continuing, we will make a priori assumption for the L∞ error of the magnetic and electric fields, n,] L∞ -Assumption: For any integer n + 1 ≤ T /τ , ||en,] E ||0,∞,Ωx , ||eB ||0,∞,Ωx ≤ Ch, with ] = 0, 1, 2, hold for small enough h, where C is a fixed positive constant depending only on the initial conditions.

This assumption will be used in section 3.2 and Lemma 3.3-(2) to estimate terms with ah,2 as this is where the nonlinear coupling of the Vlasov and Maxwell parts lies, and this assumption will eventually be established rigorously by mathematical induction in section 3.4. Our presentation will also benefit from the following shorthand notations, Z Z Z Z STAB?f = |v · nx ||[ξf? ]x |2 dsx dv + |(Eh? + v × Bh? ) · nv ||[ξf? ]v |2 dsv dx, (3.14a) Thv

STAB?EB =

Z Ex

Ex

? ? (|[ξE ]tan |2 + |[ξB ]tan |2 )dsx ,

Thx

Ev

? 2 ? 2 N? = ||ξf? ||20,Ω + ||ξE ||0,Ωx + ||ξB ||0,Ωx .

8

(3.14b)

where ? = n or n, ]. The terms STAB with different subscripts or superscripts provide stability mechanism due to the upwind phase and spatial discretizations. Later on another type of stability mechanism will emerge which is due to the temporal discretization. In our analysis, we will frequently encounter certain linear combinations of ηn , ηn,1 , ηn,2 and ηn+1 ,  = f, E, B. In the next lemma, the estimates for such terms are summarized, with their proofs given in section 4. Lemma 3.3. Let dn = d0 ηn + d1 ηn,1 + d2 ηn,2 + d3 ηn+1 ,  = f, E, B, where d0 , d1 , d2 , d3 are four constants satisfying d0 + d1 + d2 + d3 = 0 and independent of n, h, τ . Then for any g ∈ Ghk and U, V ∈ Uhk , we have 1

1

(1) ||dnf ||0,Ω + h 2 ||dnf ||0,E ≤ Cτ hk+1 , ||dn? ||0,Ωx + hx2 ||dn? ||0,Ex ≤ Cτ hk+1 , ? = E, B, x τ 2k+2 τ k+1 n,s n,s 2 n ||g||0,Ω ≤ C (h + ||g||0,Ω ), s = 0, 1, 2, (2) |ah (df , Eh , Bh ; g)| ≤ C h h h (3) |bh (dnE , dnB , dnf ; U, V )| ≤ Cτ hk (||U ||20,Ωx + ||V ||20,Ωx ).

3.2

(3.15a) (3.15b) (3.15c)

The Vlasov equation part

We start with a key decomposition of the error change in one time step [28], namely, ξfn+1 −ξfn = Gn1 +Gn2 +Gn3 , where Gn1 = ξfn,1 − ξfn , Gn2 = 2ξfn,2 − ξfn,1 − ξfn and Gn3 = ξfn+1 − 2ξfn,2 + ξfn . It is obvious that ||Gni ||20,Ω ≤ C

2 X

||ξfn,] ||20,Ω ,

i = 1, 2, 3.

(3.16)

]=0

From equations (3.10), one gets (Gn2 , g)Ω

=

(Gn3 , g)Ω

=

τ τ (K(g) − J (g)) ≡ KRK (g), 2 2 τ τ (2L(g) − K(g) − J (g)) ≡ LRK (g), 3 3

(3.17a) (3.17b)

for any g ∈ Ghk . In addition, one can verify based on (3.12) that 3||ξfn+1 ||20,Ω − 3||ξfn ||20,Ω = Ξ1 + Ξ2 + Ξ3 ,

(3.18)

where Ξi = τ [Ji (ξfn ) + Ki (ξfn,1 ) + 4Li (ξfn,2 )], i = 1, 2 and Ξ3 = (Gn2 , Gn2 )Ω + 3(Gn1 , Gn3 )Ω + 3(Gn2 , Gn3 )Ω + 3(Gn3 , Gn3 )Ω . In particular, Ξ2 relies on the phase space discretizations, for which some results were essentially established in the analysis of the semi-discrete DG methods for the VM system in [5]. On the other hand, Ξ1 and Ξ3 characterize more the contribution of the time discretization. One will see that there are two mechanisms contributing to numerical stability, one is STAB?f (? can be n or n, ]) which comes from the phase space discretization and is also used in analyzing the semi-discrete method in [5], the other one is ||Gn2 ||2 which comes from the third order Runge-Kutta time discretization. We first summarize in Lemma 3.4 some estimates, which are based on the phase space discretization and are essentially available in the analysis of the semi-discrete upwind DG method in [5]. For completeness, the proofs are given in section 4. Lemma 3.4. For ] = 0, 1, 2, we have 1 (1) ah (ξfn,] , Ehn,] , Bhn,] ; ξfn,] ) = − STABn,] f , 2

(3.19a)

1 dx STABn,] , for k ≥ f , 16 2 n,] (3) |ah (f n,] , E n,] , B n,] ; g) − ah (f n,] , Ehn,] , Bhn,] ; g)| ≤ C(||en,] E ||0,Ω + C||eB ||0,Ω )||g||0,Ω ,

(2) ah (ηfn,] , Ehn,] , Bhn,] ; ξfn,] ) ≤ Ch2k+1 + CNn,] +

moreover |ah (f n,] , E n,] , B n,] ; ξfn,] ) − ah (f n,] , Ehn,] , Bhn,] ; ξfn,] )| ≤ Ch2k+2 + CNn,] . 9

(3.19b) ∀g ∈ Ghk , (3.19c) (3.19d)

Proposition 3.5. The following estimates hold for Ξ1 and Ξ2 , (1) Ξ1 ≤ Cτ (h2k+2 + τ 6 ) + Cτ

2 X

||ξfn,] ||20,Ω ,

(3.20a)

]=0

(2) Ξ2 ≤ Cτ h2k+1 + Cτ

2 X ]=0

Nn,] −

 dx 7  n,2 τ STABnf + STABn,1 . , ∀k ≥ f + 4 STABf 16 2

(3.20b)

Proof. Recall Ξ1 = τ [J1 (ξfn ) + K1 (ξfn,1 ) + 4L1 (ξfn,2 )], with similarity, we only estimate L1 (ξfn,2 ). Applying Cauchy-Schwarz inequality, (3.15a) in Lemma 3.3, and truncation error estimate in Lemma 3.2, we get ! 3ηfn+1 − ηfn − 2ηfn,2 + 3Tfn (x, v) n,2 n,2 , ξf L1 (ξf ) = 2τ Ω  1  n,2 n+1 n n ≤ ||3ηf − ηf − 2ηf ||0,Ω + 3||Tf ||0,Ω ||ξfn,2 ||0,Ω 2τ ≤ C(hk+1 + τ 3 )||ξfn,2 ||0,Ω ≤ C(h2k+2 + τ 6 ) + C||ξfn,2 ||20,Ω . (3.21) To estimate Ξ2 , due to similarity, we will only estimate J2 (ξfn ). Using the results in Lemma 3.4, one has J2 (ξfn )

=

ah (f n , E n , B n ; ξfn ) − ah (fhn , Ehn , Bhn ; ξfn )

=

ah (ξfn , Ehn , Bhn ; ξfn ) + (ah (f n , E n , B n ; ξfn ) − ah (f n , Ehn , Bhn ; ξfn )) − ah (ηfn , Ehn , Bhn ; ξfn ) 7 Ch2k+1 + CNn − STABnf . 16

≤

2

Next we will estimate Ξ3 . One key is to use −||Gn2 ||20,Ω to control (C hτ + C hτ 2 )||Gn2 ||20,Ω under some condition on the time step τ . To make the details tractable, we first give some preparatory results. Lemma 3.6. For r = 0, 1, 2, s = 0, 1, and any g ∈ Ghk τ k+1 )h ||g||0,Ω , h τ |ah (ξfn,r , Ehn,s+1 , Bhn,s+1 ; g) − ah (ξfn,r , Ehn,s , Bhn,s ; g)| ≤ C(1 + )||ξfn,r ||0,Ω ||g||0,Ω . h |ah (ηfn,r , Ehn,s+1 , Bhn,s+1 ; g) − ah (ηfn,r , Ehn,s , Bhn,s ; g)| ≤ C(1 +

(1) (2)

(3.22a) (3.22b)

Lemma 3.7. For any g ∈ Ghk , we have   2 X τ n,] n,] n,] k+1 3 LRK (g) ≤ C (1 + ) (||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξf ||0,Ω + h ) + τ  ||g||0,Ω + ah (Gn2 , Ehn,1 , Bhn,1 ; g) h ]=0

(3.23a) 



2

||Gn2 ||0,Ω  τ X n,] n,] ) (||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξfn,] ||0,Ω + hk+1 ) + τ 3 + ||g||0,Ω . (3.23b) h h ]=0   2 τ X n,] n,] n,] KRK (g) ≤ C(1 + ) (||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξf ||0,Ω + hk+1 ) ||g||0,Ω + ah (Gn1 , Ehn,1 , Bhn,1 ; g). h ≤ C (1 +

]=0

(3.23c) Lemma 3.8.

≤

|ah (Gn1 , Ehn,1 , Bhn,1 ; Gn2 ) + ah (Gn2 , Ehn,1 , Bhn,1 ; Gn1 )|  C 1 τ n 2 C 1+ ||ξf ||0,Ω + ||Gn2 ||20,Ω + (STABnf + STABn,1 f ). h h 16 10

(3.24)

With all the preparation in Lemmas 3.6 - 3.8, we are now ready to estimate Ξ3 . Proposition 3.9. Ξ3

2  τ τ  X n,] ≤ C τ (1 + ) + τ 2 (1 + )2 ( N + h2k+2 ) + Cτ 7 h h ]=0   2 τ τ τ + −1 + C + C 2 ||Gn2 ||20,Ω + (STABnf + STABn,1 f ). h h 16

(3.25)

Proof. First note that Ξ3

=

−||Gn2 ||20,Ω + 2(Gn2 , Gn2 )Ω + 3(Gn1 , Gn3 )Ω + 3(Gn2 , Gn3 )Ω + 3(Gn3 , Gn3 )Ω

= −||Gn2 ||20,Ω + τ (KRK (Gn2 ) + LRK (Gn1 ) + LRK (Gn2 )) + 3(Gn3 , Gn3 )Ω .

(3.26)

Based on (3.23a) and (3.23c), KRK (Gn2 ) + LRK (Gn1 ) ≤

ah (Gn2 , Ehn,1 , Bhn,1 ; Gn1 ) + ah (Gn1 , Ehn,1 , Bhn,1 ; Gn2 )   2 X τ n,] n,] n,] (||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξf ||0,Ω + hk+1 ) + τ 3  (||Gn1 ||0,Ω + ||Gn2 ||0,Ω ). +C (1 + ) h ]=0

We now apply Lemma 3.8 to estimate the first two terms on the right, and apply Cauchy-Schwartz inequality and (3.16) to estimate the last term, and get

≤

KRK (Gn2 ) + LRK (Gn1 )   2  τ  X n,] C 1 C 1+ N + h2k+2  + Cτ 6 + ||Gn2 ||20,Ω + (STABnf + STABn,1 f ). h h 16

(3.27)

]=0

In order to estimate LRK (Gn2 ) in (3.26), we apply (3.23b) of Lemma 3.7. In particular, take g = Gn2 in (3.23b) and use (3.16), we have   2 n X τ ||G || n,] n,] 2 0,Ω  LRK (Gn2 ) ≤ C (1 + ) (||ξE ||0,Ωx + ||ξB ||Gn2 ||0,Ω ||0,Ωx + ||ξfn,] ||0,Ω + hk+1 ) + τ 3 + h h ]=0   2 τ X n,] C ≤ C(1 + ) N + h2k+2  + Cτ 6 + ||Gn2 ||20,Ω . (3.28) h h ]=0

For 3(Gn3 , Gn3 )Ω = τ LRK (Gn3 ), we take g = Gn3 in (3.23b) and obtain   2 n X || τ ||G n,] n,] 2 0,Ω  (||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξfn,] ||0,Ω + hk+1 ) + τ 3 + ||Gn3 ||0,Ω 3||Gn3 ||20,Ω ≤ Cτ (1 + ) h h ]=0   2 X τ τ2 ≤ Cτ 2 (1 + )2  Nn,] + h2k+2  + Cτ 8 + C 2 ||Gn2 ||20,Ω + ||Gn3 ||20,Ω . h h ]=0

Therefore, with a different constant C we have,   2 X τ2 τ Nn,] + h2k+2  + Cτ 8 + C 2 ||Gn2 ||20,Ω . 3(Gn3 , Gn3 )Ω ≤ Cτ 2 (1 + )2  h h ]=0

Finally, we complete the proof by combing (3.26)-(3.29). 11

(3.29)

Now we are ready to establish the main error estimate result for the Vlasov solver. Theorem 3.10. Let (f, E, B) be a sufficiently smooth exact solution to equations (1.1). Let (fhn , Ehn , Bhn ) ∈ Ghk × Uhk × Uhk be the solution to the scheme (2.6) at time tn with k ≥ d2x . Under the L∞ -Assumption, there exists a positive constant γ1 , such that for any hτ ≤ γ1 , the following estimate holds for ∀n : n + 1 ≤ T /τ 3||ξfn+1 ||20,Ω − 3||ξfn ||20,Ω ≤ Cτ h2k+1 + Cτ 7 + Cτ

2 X

Nn,] .

(3.30)

]=0

Proof. Since the constant C in the result (3.25) is independent of n, h, τ , there exists a positive constant γ1 2 independent of n, h, τ , such that −1 + C hτ + C hτ 2 ≤ − 12 as long as hτ ≤ γ1 . Under such condition, we combine (3.18) and the estimates in Proposition 3.5 and Proposition 3.9 and get 3||ξfn+1 ||20,Ω − 3||ξfn ||20,Ω = Ξ1 + Ξ2 + Ξ3 ≤

Cτ h2k+1 + Cτ 7 + Cτ

2 X ]=0

≤ Cτ h2k+1 + Cτ 7 + Cτ

2 X

 τ 1 n,2 3STABnf + 3STABn,1 + 14STAB Nn,] − ||Gn2 ||20,Ω − f f 2 8 Nn,] .

]=0

3.3

The Maxwell equations part

In this section, we estimate how error accumulates in the Maxwell solver. The procedure is parallel to the Vlasov part in section 3.2. We start with a decomposition of the error change in one time step n+1 n ξE − ξE = X1n + X2n + X3n ,

n+1 n ξB − ξB = Z1n + Z2n + Z3n ,

n,1 n,2 n,1 n,2 n,1 n,2 n,1 n+1 n n n n n −ξE , X3n = ξE where X1n = ξE −ξE , X2n = 2ξE −ξE −2ξE +ξE , Z1n = ξB −ξB , Z2n = 2ξB −ξB −ξB , n,2 n+1 n n and Z3 = ξB − 2ξB + ξB . It is obvious that

||Xin ||20,Ω ≤ C

2 X

n,] 2 ||ξE ||0,Ωx , ||Zin ||20,Ω ≤ C

]=0

2 X

n,] 2 ||ξB ||0,Ωx ,

i = 1, 2, 3.

(3.31)

]=0

Based on equations (3.11), the following hold for any U, V ∈ Uhk . (X2n , U )Ωx + (Z2n , V )Ωx

=

(X3n , U )Ωx + (Z3n , V )Ωx

=

τ τ (R(U, V ) − Q(U, V )) ≡ RRK (U, V ), 2 2 τ τ (2S(U, V ) − R(U, V ) − Q(U, V )) ≡ SRK (U, V ). 3 3

In addition, one can verify based on (3.13) that n+1 2 n+1 2 n 2 n 2 3(||ξE ||0,Ωx + ||ξB ||0,Ωx ) − 3(||ξE ||0,Ωx + ||ξB ||0,Ωx ) = Θ1 + Θ2 + Θ3 ,

(3.32)

where Θi =

  n,1 n,1 n,2 n,2 n n τ Qi (ξE , ξB ) + Ri (ξE , ξB ) + 4Si (ξE , ξB ) ,

Θ3 =

(X2n , X2n )Ωx + 3(X1n , X3n )Ωx + 3(X2n , X3n )Ωx + 3(X3n , X3n )Ωx

i = 1, 2

+(Z2n , Z2n )Ωx + 3(Z1n , Z3n )Ωx + 3(Z2n , Z3n )Ωx + 3(Z3n , Z3n )Ωx .

(3.33)

(3.34)

In particular, Θ2 depends on the spatial discretization, while Θ1 and Θ3 characterizes the contribution from the time discretization. Similarly as in the Vlasov part, there are two stability mechanisms, with one being 12

STAB?EB (? can be n or n, ]) from the spatial discretization, and the other is related to ||X2n ||2 , ||Z2n ||2 arising from the third order Runge-Kutta time discretization. Next we will estimate Θ1 and Θ2 . Some estimates for the spatial discretizations of the Maxwell part are summarized in Lemma 3.11. Lemma 3.11. For ] = 0, 1, 2, we have 1 n,] n,] n,] n,] n,] n,] 2 (i) bh (ξE , ξB , ξf ; ξE , ξB ) ≤ C(||ξfn,] ||20,Ω + ||ξE ||0,Ωx ) − STABn,] EB , 2 1 n,] n,] n,] n,] n,] n,] 2 (ii) |bh (ηE , ηB , ηf ; ξE , ξB )| ≤ Ch2k+1 + C||ξE ||0,Ex + STABn,] EB . 16 Proposition 3.12. The following estimates hold for Θ1 and Θ2 . Θ1 ≤ Cτ (h2k+2 + τ 6 ) + Cτ

2 X

Nn,]

(3.35a) (3.35b)

(3.36a)

]=0

Θ2 ≤ Cτ h2k+1 + Cτ

2 X

Nn,] −

]=0

7 n,2 τ (STABnEB + STABn,1 EB + 4STABEB ). 16

(3.36b)

  n,1 n,1 n,2 n,2 n n ) + R1 (ξE , ξB ) + 4S1 (ξE , ξB ) , with similarity, we only estimate , ξB Proof. Recall that Θ1 = τ Q1 (ξE n,2 n,2 S1 (ξE , ξB ). Applying Cauchy-Schwarz inequality, Lemma 3.2, Lemma 3.3 - (1), one gets ! ! n,2 n,2 n+1 n+1 n n − 2ηE + 3TEn (x) n,2 3ηE − ηE − 2ηB + 3TBn (x) n,2 3ηB − ηB n,2 n,2 S1 (ξE , ξB ) = + , ξE , ξB 2τ 2τ Ωx

Ωx

n,2 2 n,2 2 C(||ξE ||0,Ωx + ||ξB ||0,Ωx + h2k+2 + τ 6 ) ≤ C(Nn,2 + h2k+2 + τ 6 ).   n,1 n,1 n,2 n,2 n n ) + R2 (ξE , ξB ) + 4S2 (ξE , ξB ) , with similarity, we only need to look , ξB To bound Θ2 = τ Q2 (ξE

≤

n n ). By definition, and Lemma 3.11, , ξB at Q2 (ξE n n n n n n n n n n n n Q2 (ξE , ξB ) = bh (enE , enB , enf ; ξE , ξB ) = bh (ξE , ξB , ξfn ; ξE , ξB ) − bh (ηE , ηB , ηfn ; ξE , ξB ) 7 n 2 ≤ Ch2k+1 + C(||ξfn ||20,Ω + ||ξE ||0,Ωx ) − STABnEB 16 7 n 2k+1 n ≤ Ch + CN − STABEB . 16

Next we will estimate Θ3 . One key is to use −||X2n ||20,Ωx −||Z2n ||20,Ωx to control C hτ (||X2n ||20,Ωx +||Z2n ||20,Ωx ). To make the analysis easy to follow, we first present some preparatory results in two lemmas. Lemma 3.13. For any U, V ∈ Uhk , we have SRK (U, V ) ≤ C(hk+1 + τ hk + τ 3 )(||U ||0,Ωx + ||V ||0,Ωx ) + bh (X2n , Z2n , Gn2 ; U, V ) 1 ≤ C(hk+1 + τ hk + τ 3 + (||Z2n ||0,Ωx + ||X2n ||0,Ωx ))(||U ||0,Ωx + ||V ||0,Ωx ) h + C||Gn2 ||0,Ω ||U ||0,Ωx . k+1

RRK (U, V ) ≤ C(h

k

+ τ h )(||U ||0,Ωx + ||V ||0,Ωx ) +

bh (X1n , Z1n , Gn1 ; U, V

).

(3.37a)

(3.37b) (3.37c)

Lemma 3.14. bh (X1n , Z1n , Gn1 ; X2n , Z2n ) + bh (X2n , Z2n , Gn2 ; X1n , Z1n ) ≤

2 X C 1 (||X2n ||20,Ωx + ||Z2n ||20,Ωx ) + C (||Gnj ||20,Ω + ||Xjn ||20,Ωx ) + (STABnEB + STABn,1 EB ). h 16 j=1

13

Now we get ready to estimate Θ3 . Proposition 3.15. Θ3

≤

 τ τ2 Cτ (h + τ h + τ ) + −1 + C + C 2 (||X2n ||20,Ωx + ||Z2n ||20,Ωx ) h h 2  X τ  STABnEB + STABn,1 + Nn,] . EB + Cτ 16 2k+2

2 2k

6



(3.38)

]=0

Proof. First note that Θ3

= −||X2n ||20,Ωx − ||Z2n ||20,Ωx +τ (RRK (X2n , Z2n ) + SRK (X1n , Z1n ) + SRK (X2n , Z2n )) + 3||X3n ||20,Ωx + 3||Z3n ||20,Ωx .

(3.39)

Based on (3.37a), (3.37c), and Lemma 3.14, in addition to (3.16) and (3.31) we have RRK (X2n , Z2n ) + SRK (X1n , Z1n ) ≤ bh (X2n , Z2n , Gn2 ; X1n , Z1n ) + bh (X1n , Z1n , Gn1 ; X2n , Z2n ) + C(hk+1 + τ hk + τ 3 )(||X1n ||0,Ωx + ||Z1n ||0,Ωx ) + C(hk+1 + τ hk )(||X2n ||0,Ωx + ||Z2n ||0,Ωx ) 1 C ≤ (||X2n ||20,Ωx + ||Z2n ||20,Ωx ) + (STABnEB + STABn,1 EB ) h 16 2 X + C(h2k+2 + τ 2 h2k + τ 6 ) + C (||Gnj ||20,Ω + ||Xjn ||20,Ωx + ||Zjn ||20,Ωx ) j=1

≤

2 X C 1 2k+2 2 2k 6 Nn,] . (3.40) (||X2n ||20,Ωx + ||Z2n ||20,Ωx ) + (STABnEB + STABn,1 ) + C(h + τ h + τ ) + C EB h 16 j=1

Next we estimate SRK (X2n , Z2n ) by using (3.37b) in Lemma 3.13. 1 (||Z2n ||0,Ωx + ||X2n ||0,Ωx ))(||X2n ||0,Ωx + ||Z2n ||0,Ωx ) + C||Gn2 ||0,Ω ||X2n ||0,Ωx h 2 X C (3.41) ≤ C(h2k+2 + τ 2 h2k + τ 6 ) + C Nn,] + (||Z2n ||20,Ωx + ||X2n ||20,Ωx ). h j=1

SRK (X2n , Z2n ) ≤ C(hk+1 + τ hk + τ 3 +

Finally we turn to 3||X3n ||20,Ωx + 3||Z3n ||20,Ωx in (3.39). By applying (3.37b) in Lemma 3.13, 3||X3n ||20,Ωx + 3||Z3n ||20,Ωx = τ SRK (X3n , Z3n ) 1 ≤ Cτ (hk+1 + τ hk + τ 3 + (||Z2n ||0,Ωx + ||X2n ||0,Ωx ))(||X3n ||0,Ωx + ||Z3n ||0,Ωx ) + Cτ ||Gn2 ||0,Ω ||X3n ||0,Ωx h τ2 ≤ C(τ 2 h2k+2 + τ 4 h2k + τ 8 ) + C 2 (||Z2n ||20,Ωx + ||X2n ||20,Ωx ) + Cτ 2 ||Gn2 ||20,Ω + ||X3n ||20,Ωx + ||Z3n ||20,Ωx , h therefore, with a different value of C, we have τ2 (||Z2n ||20,Ωx + ||X2n ||20,Ωx ) + Cτ 2 ||Gn2 ||20,Ω . (3.42) h2 Now by combining the results in (3.39) - (3.42) and (3.16), we can conclude the estimate for Θ3 in (3.38). 3||X3n ||20,Ωx + 3||Z3n ||20,Ωx ≤ C(τ 2 h2k+2 + τ 4 h2k + τ 8 ) + C

The main error estimate result for the Maxwell solver is now established as following. Theorem 3.16. Let (f, E, B) be a sufficiently smooth exact solution to equations (1.1). Let (fhn , Ehn , Bhn ) ∈ Ghk × Uhk × Uhk be the solution to the scheme (2.6) at time tn . There exists a positive constant γ2 , such that for hτ ≤ γ2 , the following estimate holds for ∀n : n + 1 ≤ T /τ n+1 2 n+1 2 n 2 n 2 ||0,Ωx + ||ξB ||0,Ωx ) ≤ Cτ h2k+1 + Cτ 7 + Cτ 3(||ξE ||0,Ωx + ||ξB ||0,Ωx ) − 3(||ξE

2 X ]=0

14

Nn,] .

(3.43)

Proof. Since the constant C in the estimate (3.38) is independent of n, h, τ , there exists a positive constant 2 γ2 independent of n, h, τ , such that −1 + C hτ + C hτ 2 ≤ − 21 as long as hτ ≤ γ2 . Under this condition on the time step τ , we also have τ 2 h2k ≤ γ22 h2k+2 . Now we combine the estimates for Θi , i = 1, 2, 3 in Proposition 3.12 and Proposition 3.15, and get n+1 2 n+1 2 n 2 n 2 3(||ξE ||0,Ωx + ||ξB ||0,Ωx ) − 3(||ξE ||0,Ωx + ||ξB ||0,Ωx ) = Θ1 + Θ2 + Θ3 2  X τ 1 n,2 3STABnEB + 3STABn,1 + 14STAB ≤ Cτ h2k+1 + Cτ 7 + Cτ Nn,] − (||X2n ||20,Ωx + ||Z2n ||20,Ωx ) − EB EB 2 8 ]=0

≤ Cτ h2k+1 + Cτ 7 + Cτ

2 X

Nn,] .

]=0

Remark 3.17. Unlike in Theorem 3.10, the a priori assumption about the L∞ norm of the error in the magnetic and electric fields, together with the condition of k ≥ d2x , are not needed in Theorem 3.16. This difference is due to the nonlinear coupling terms in the Vlasov equation.

3.4

Proof of the main result: Theorem 3.1

The following lemma is the final preparation. Lemma 3.18. Suppose hτ ≤ α, where α is a positive constant independent of n, h, τ . Then under the L∞ -Assumption, the following inequalities are satisfied ||ξfn,1 ||20,Ω

≤

n 2 n 2 Ch2k+2 + C(||ξfn ||20,Ω + τ 2 ||ξE ||0,Ωx + τ 2 ||ξB ||0,Ωx )

||ξfn,2 ||20,Ω

≤

n,1 2 n,1 2 Ch2k+2 + C(||ξfn ||20,Ω + ||ξfn,1 ||20,Ω + τ 2 ||ξE ||0,Ωx + τ 2 ||ξB ||0,Ωx )

n,1 2 ||ξE ||0,Ωx

≤

n 2 n 2 Ch2k+2 + C(τ 2 ||ξfn ||20,Ω + ||ξE ||0,Ωx + ||ξB ||0,Ωx )

n,2 2 ||ξE ||0,Ωx

≤

n,1 2 n,1 2 n 2 Ch2k+2 + C(τ 2 ||ξfn,1 ||20,Ωx + ||ξE ||0,Ωx + ||ξE ||0,Ωx + ||ξB ||0,Ωx )

n,1 2 ||ξB ||0,Ωx

≤

n 2 n 2 Ch2k+2 + C(||ξE ||0,Ωx + ||ξB ||0,Ωx )

n,2 2 ||ξB ||0,Ωx

≤

n,1 2 n,1 2 n 2 Ch2k+2 + C(||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξB ||0,Ωx ).

A direct consequence of these inequalities is 2 X

Nn,] ≤ CNn + Ch2k+2 .

(3.44)

]=0

Now we are ready for the proof of Theorem 3.1. Proof. First we make the a priori assumption for the L∞ error of the magnetic and electric fields as in L∞ -Assumption. Based on Theorem 3.10 and Theorem 3.16, for any hτ ≤ γ := min(γ1 , γ2 ), we get 3Nn+1 − 3Nn

≤

Cτ 7 + Cτ h2k+1 + Cτ

2 X


Nn,] ≤ Cˆ τ 7 + τ h2k+1 + τ Nn .

(3.45)

]=0

Here (3.44) is used to get the last inequality. Let Υn = Nn /(1 + ˆ 7 +τ h2k+1 ) C(τ . ˆ n 3(1+ C 3 τ)

ˆ n Cτ 3 ) ,

then (3.45) leads to Υn − Υn−1 ≤

We now sum up Υ] − Υ]−1 and use Υ0 = 0 to obtain Υn ≤

n ˆ 7 X C(τ + τ h2k+1 ) ]=1

3(1 +

ˆ C ] 3 τ)

15

≤ τ 6 + h2k+1 .

(3.46)

Therefore

ˆ ˆ ˆ Cnτ CT Cτ )n Υn ≤ e 3 (τ 6 + h2k+1 ) ≤ e 3 (τ 6 + h2k+1 ), 3 for any n ≤ T /τ . This can also be written as Nn = (1 +

(3.47)

n 2 n 2 ||ξfn ||20,Ω + ||ξE ||0,Ωx + ||ξB ||0,Ωx ≤ C(τ 6 + h2k+1 ).

(3.48)

∞

ξf0

0 ξE

0 Next we will establish the L -Assumption using mathematical induction. For n = 0, = = ξB =0 0 0 and the approximation property (2.4) implies that ||eE ||0,∞,Ωx , ||eB ||0,∞,Ωx ≤ Ch, where C depends only on 0,] the initial conditions and is independent of n, h, τ . Furthermore, Lemma 3.18 shows ||ξf0,] ||0,Ω , ||ξE ||0,Ωx , 0,] ||ξB ||0,Ωx ≤ Chk+1 for ] = 1, 2. Thus 0,] 0,] − ||e0,] E ||0,∞,Ωx ≤ ||ξE ||0,∞,Ωx + ||ηE ||0,∞,Ωx ≤ Ch

dx 2

0,] ||ξE ||0,Ωx + Chk+1 ≤ Chk+1− dx +1 2

dx 2

≤ Ch,

dx 2

for small enough h. Here the inverse inequality (2.5) and the condition k > > are used. Similarly 0,] n n ||eB ||0,∞,Ωx ≤ Ch. Suppose ||eB ||0,∞,Ωx , ||eE ||0,∞,Ωx ≤ Ch for all n ≤ m, therefore inequality (3.48) with m+1,] m+1,] n = m + 1 is satisfied. Combining this with Lemma 3.18 implies ||ξfm+1,] ||0,Ω , ||ξE ||0,Ωx , ||ξB ||0,Ωx ≤ 1 3 k+ 2 ) for ] = 0, 1, 2. Thus C(τ + h m+1,] ||eE ||0,∞,Ωx

≤ Ch−

dx 2

m+1,] ||ξE ||0,Ωx + Chk+1 ≤ Ch−

≤ Ch−

dx 2

τ 3 + Chk+ 2 −

1

dx 2

≤ Cγ 3 h3−

dx 2

dx +1 2

dx 2

1

(τ 3 + hk+ 2 ) + Chk+1 1

+ Chk+ 2 −

dx 2

≤ C0 h1+ ,

dx +1 2 ,2

(3.49)

dx 2 )

for k > (recall dx ≤ 3), under the condition hτ ≤ γ and with  = min(k − − > 0. The positive constant C0 depends on T in general and is independent of m, h, τ . One can see that there exists m+1,] h0 > 0, such that C0 h < C, hence ||eE ||0,∞,Ωx ≤ Ch for h < h0 . Up to now the a priori assumption is established, and this also completes the L2 error estimates of f , E and B. Finally we will show how the Gauss’s laws are approximated by the proposed methods. Based on the fact that the Gauss’s laws are satisfied by the exact solutions, using inverse inequalities (2.5) and approximating properties of the discrete spaces in (2.4), and applying the L2 error estimate of f and E, we get ||∇ · Ehn − ρnh + ρni ||0,Ωx = ||∇ · Ehn − ρnh + ρni − (∇ · E n − ρn + ρni )||0,Ωx =

||∇ · (Ehn − E n ) − (ρnh − ρn )||0,Ωx

n n ||∇ · ηE ||0,Ωx + ||∇ · ξE ||0,Ωx + C||enf ||0,Ω 1 C k+ 1 C n ||ξE ||0,Ωx + C(hk+ 2 + τ 3 ) ≤ (h 2 + τ 3 ). ≤ Chkx ||E n ||k+1,Ωx + hx hx

≤

Similarly one can get ||∇ · Bhn ||0,Ωx ≤

C k+ 12 hx (h

+ τ 3 ). This completes all the estimates in Theorem 3.1.

Remark 3.19. The condition k > dx2+1 is needed only to establish the L∞ -Assumption, while for all the other results, the requirement k ≥ d2x is enough.

4

Proofs of Some Lemmas

In this section, we will provide the proofs of some lemmas in section 3.

4.1

Proof of Lemma 3.3

To get (3.15a), we start with f n,1 and f n,2 given in (3.8), and get d0 f n + d1 f n,1 + d2 f n,2 + d3 f n+1   τ2 d2 τ3 ∂t f n + d2 ∂t2 f n − d2 (∂t E n + v × ∂t B n ) · ∇v (∂t f n ) =d3 (f n+1 − f n ) + τ d1 + 2 4 4   2 d2 τ τ3 =τ d3 ∂t f (x, v, t? ) + τ d1 + ∂t f n + d2 ∂t2 f n − d2 (∂t E n + v × ∂t B n ) · ∇v (∂t f n ), 2 4 4 16

(4.50)

for some t? ∈ [tn , tn+1 ]. Here we have used d0 + d1 + d2 + d3 = 0. Therefore, with I as the identity operator,     d2 τ2 τ3 dnf = (Πk − I) τ d3 ∂t f (x, v, t? ) + τ d1 + ∂t f n + d2 ∂t2 f n − d2 (∂t E n + v × ∂t B n ) · ∇v (∂t f n ) . 2 4 4 For sufficiently smooth solution, there is ||dnf ||0,Ω ≤ Cτ hk+1 max

∀t∈[0,T ]


||∂t f ||k+1,Ω + τ ||∂t2 f ||k+1,Ω + τ 2 ||(∂t E + v × ∂t B) · ∇v (∂t f )||k+1,Ω . 1

Recall that τ ≤ 1, we further have ||dnf ||0,Ω ≤ Cτ hk+1 . Similarly one can show ||dnf ||0,E ≤ Cτ hk+ 2 . These two estimates will lead to the upper bound in (3.15a). The proof of the results for E and B can be proceeded similarly. To get (3.15b), based on definition |ah (dnf , Ehn,s , Bhn,s ; g)| Z Z ≤ |dnf v · ∇x g|dxdv +

  Z |v · nx | n |{dnf v}x + [df ]x | |[g]x |dsx dv + |dnf (Ehn,s + v × Bhn,s )∇v g|dxdv v 2 Ω Th E x Ω   Z Z n,s n,s |(E + v × B ) · n | v h h [dnf ]v · [g]v |dsv dx | {dnf (Ehn,s + v × Bhn,s )}v + + 2 Thx Ev
≤ C ||dnf ||0,Ω ||∇x g||0,Ω + ||dnf ||0,Thv ×Ex ||g||0,Thv ×Ex
+C (||Ehn,s ||0,∞,Ωx + ||Bhn,s ||0,∞,Ωx ) ||dnf ||0,Ω ||∇v g||0,Ω + ||dnf ||0,Thx ×Ev ||g||0,Thx ×Ev   1 C ||g||0,Ω ||dnf ||0,Ω + h 2 ||dnf ||0,Thv ×Ex ≤ h   1 C + ||g||0,Ω (||Ehn,s ||0,∞,Ωx + ||Bhn,s ||0,∞,Ωx ) ||dnf ||0,Ω + h 2 ||dnf ||0,Thx ×Ev . h Z

By the a priori L∞ -Assumption, there is n,s ||Ehn,s ||0,∞,Ωx ≤ ||en,s ||0,∞,Ωx ≤ Ch + C ≤ C, E ||0,∞,Ωx + ||E

(4.51)

||Bhn,s ||0,∞,Ωx ≤ C.

(4.52)

τ τ |ah (dnf , Ehn,s , Bhn,s ; g)| ≤ C ||g||0,Ω hk+1 ≤ C (h2k+2 + ||g||20,Ω ). h h

(4.53)

and similarly, Finally we apply (3.15a) and conclude

To get (3.15c), bh (dnE , dnB , dnf ; U, V ) Z Z Z Z = dnB · (∇ × U )dx − dnE · (∇ × V )dx − ( dnf vdv) · U dx Ωx Ωx Ωx Ωv Z Z 1 1 n n n ({dE }x + [dnB ]tan ) · [V ]tan dsx + ({dB }x − [dE ]tan ) · [U ]tan dsx − 2 2 Ex ZEx Z Z Z 1 n 1 n n = − ( df vdv) · U dx + ({dB }x − [dE ]tan ) · [U ]tan dsx − ({dnE }x + [dnB ]tan ) · [V ]tan dsx , 2 2 Ωx Ωv Ex Ex R R where the second equality is due to ∇×U, ∇×V ∈ Uhk , which leads to Ωx dnB ·(∇×U )dx = Ωx dnE ·(∇×V )dx = 0. We now can apply inverse inequality (2.5) and (3.15a) to obtain (3.15c).

17

4.2

Proof of Lemma 3.4

To get (3.19a), based on definition, for any g ∈ Ghk , U, V ∈ Uhk  |v · nx | gv · ∇x gdxdv − [g]x · [g]x dsx dv {gv}x + 2 Ω Thv Ex  Z Z Z  X Z 1 |v · nx | v · nx g 2 dsx dv − {gv}x + [g]x · [g]x dsx dv 2 Thv 2 Thv Ex Kx ∈Thx ∂Kx    Z Z  1 |v · nx | v · [g 2 ]x − {gv}x + [g]x · [g]x dsx dv v 2 2 Th E x   Z Z  |v · nx | 1 2 [g ]x − {g}x [g]x · v − |[g]x |2 dsx dv 2 2 Thv Ex Z Z 1 − |v · nx ||[g]x |2 dsx dv. 2 Thv Ex Z

ah,1 (g; g)

= = = = =

Z



Z

(4.54)

Here the property ofRjump R and average in (2.3a) is used to get the last equality in (4.54). Similarly, there is ah,2 (g, U, V ; g) = − 21 T x Ev |(U + v × V ) · nv ||[g]v |2 dsv dx. Finally, taking g = ξfn,] , U = Ehn,] , V = Bhn,] h leads to (3.19a). To get (3.19b), we will proceed the proof in two steps. Step 1: To estimate ah,1 (ηfn,] ; ξfn,] ). By definition, ah,1 (ηfn,] ; ξfn,] )

Z = Ω

ηfn,] v

·

∇x ξfn,] dxdv

  |v · nx | n,] n,] [ηf ]x · [ξfn,] ]x dsx dv. {ηf v}x + 2 Ex

Z

Z − Thv

(4.55)

Let v0 be the L2 projection of the function v onto the piecewise constant space with respect to Thv , then Z Z Z ηfn,] v · ∇x ξfn,] dxdv = ηfn,] (v − v0 ) · ∇x ξfn,] dxdv + ηfn,] v0 · ∇x ξfn,] dxdv Ω Ω ZΩ n,] n,] = ηf (v − v0 ) · ∇x ξf dxdv. (4.56) Ω

The last equality is satisfied because v0 · ∇x ξfn,] ∈ Ghk and the L2 projection of ηfn,] onto Ghk vanishes. We further have Z X n,] n,] η n,] v · ∇x ξ n,] dxdv ≤ ||v − v0 ||0,∞,Ωv (h−1 Kx ||ηf ||0,K )(hKx ||∇x ξf ||0,K ) f f Ω

K∈Th

≤ Chv ||v||1,∞,Ωv

X

−1 n,] hk+1 ||k+1,K ||ξfn,] ||0,K ≤ Chk+1 ||ξfn,] ||0,Ω . K hKx ||f

(4.57)

K∈Th

Here Cauchy-Schwarz inequality and approximate properties (2.4) are used for the first and second inequality

18

above, respectively. Applying the similar technique to the second term of ah,1 (ηfn,] ; ξfn,] ), we get Z Z   |v · nx | n,] n,] n,] {ηf v}x + [ηf ]x · [ξf ]x dsx dv T v Ex 2 h ! Z Z |[ηfn,] ]x | ≤ |v · nx |(|{ηfn,] }x | + ) |[ξfn,] ]x |dsx dv 2 Thv Ex ! ! 21 Z Z |[ηfn,] ]x | 2 n,] 2 ≤ 2 |{ηf }x | + ( ) |v · nx |dsx dv 2 Thv Ex ! 21 Z Z ≤ C||ηfn,] ||0,Thv ×Ex

Thv

Ex

|v · nx ||[ξfn,] ]x |2 dsx dv

Z

! 21

Z

nx ||[ξfn,] ]x |2 dsx dv

|v · Thv k+ 12

Ex

Z

Z

≤ Ch

Thv

Ex

! 21 |v · nx ||[ξfn,] ]x |2 dsx dv

. (4.58)

By (4.57), (4.58) and Young’s inequality, there is ah,1 (ηfn,] ; ξfn,] ) ≤ C||ξfn,] ||20,Ω + Ch2k+1 +

1 16

Z Thv

Z Ex

|v · nx ||[ξfn,] ]x |2 dsx dv.

(4.59)

Step 2: To estimate ah,2 (ηfn,] , Ehn,] , Bhn,] ; ξfn,] ). Let E0 = Π0x E n,] , B0 = Π0x B n,] be the L2 projection of E n,] , B n,] onto piecewise constant vector space with respect to Thx , then Z η n,] (E n,] + v × B n,] ) · ∇v ξ n,] dxdv f h h f Ω Z = ηfn,] (Ehn,] − E0 + v × (Bhn,] − B0 )) · ∇v ξfn,] dxdv Ω

≤ (||Ehn,] − E0 ||0,∞,Ωx + C||Bhn,] − B0 ||0,∞,Ωx )||ηfn,] ||0,Ω ||∇v ξfn,] ||0,Ω .

(4.60)

The first equality of (4.60) is satisfied because (E0 + v × B0 ) · ∇v ξfn,] ∈ Ghk , combining with the fact that the L2 projection of ηfn,] onto Ghk is equal to zero. Since the operator Πkx is bounded in any Lp (1 ≤ p ≤ ∞) norm [11], we can further estimate ||Ehn,] − E0 ||0,∞,Ωx ≤ ||Ehn,] − Πkx E n,] ||0,∞,Ωx + ||Πkx E n,] − E0 ||0,∞,Ωx , and ||Πkx E n,] − E0 ||0,∞,Ωx = ||Πkx (E n,] − E0 )||0,∞,Ωx ≤ C||E n,] − E0 ||0,∞,Ωx ≤ Chx ||E n,] ||1,∞,Ωx . Thus, n,] ||Ehn,] −E0 ||0,∞,Ωx ≤ ||ξE ||0,∞,Ωx +Chx . Similar treatment can be applied to ||Bhn,] −B0 ||0,∞,Ωx . Therefore, Z η n,] (E n,] + v × B n,] ) · ∇v ξ n,] dxdv f h h f Ω

n,] n,] ≤ (||ξE ||0,∞,Ωx + ||ξB ||0,∞,Ωx + Chx )||ηfn,] ||0,Ω ||∇v ξfn,] ||0,Ω n,] n,] ≤ (||ξE ||0,∞,Ωx + ||ξB ||0,∞,Ωx + Chx )Chk ||ξfn,] ||0,Ω

≤ Chk−

dx 2

n,] n,] (||ξE ||0,Ωx + ||ξB ||0,Ωx )||ξfn,] ||0,Ω + Chk+1 ||ξfn,] ||0,Ω .

(4.61)

Here we have applied the inverse inequality in (2.5) to the last inequality above. The second term in

19

ah,2 (ηfn,] , Ehn,] , Bhn,] ; ξfn,] ) can be estimated as Z Z T x Ev h Z Z ≤

! |(Ehn,] + v × Bhn,] ) · nv | n,] n,] +v× + [ηf ]v · [ξf ]v dsv dx 2 ! |[ηfn,] ]v | |(Ehn,] + v × Bhn,] ) · nv |(|{ηfn,] }v | + ) |[ξfn,] ]v |dsv dx 2 ! 12 Z Z

{ηfn,] (Ehn,]

Thx

Ev

Z

Z

≤ Thx

Ev

Bhn,] )}v

2|{(ηfn,] )2 }v ||(Ehn,]

+v×

1

≤

2 C||ηfn,] ||0,Thx ×Ev (||Ehn,] ||0,∞,Ω x

k+ 12

≤ Ch

+

· nv |dsv dx Thx

+

2 ||Bhn,] ||0,∞,Ω ) x

Z Thx

Ev

|(Ehn,]

+v×

Bhn,] )

! 12 ·

nv ||[ξfn,] ]v |2 dsv dx ! 21

Z

Thx

Ev

|(Ehn,]

+v×

Bhn,] )

·

nv ||[ξfn,] ]v |2 dsv dx ! 12

Z

1

1

Z

1

2 ||Bhn,] ||0,∞,Ω ) x

1

1

2 (||Ehn,] ||0,∞,Ω x

Bhn,] )

(4.62)

Ev

|(Ehn,]

+v×

Bhn,] )

·

nv ||[ξfn,] ]v |2 dsv dx

1

1

.

1

n,] 2 n,] 2 2 2 2 ||0,∞,Ωx + ||Πkx E n,] ||0,∞,Ω ||0,∞,Ωx + C||E n,] ||0,∞,Ω ≤ ||ξE ≤ ||ξE . From the inverse Note ||Ehn,] ||0,∞,Ω x x x n,] inequality in (2.5), there is ||ξE ||0,∞,Ωx ≤ Ch−

C(1 + h−

dx 4

dx 2

1

n,] 2 ||ξE ||0,Ωx . Since ||E n,] ||0,∞,Ωx is bounded, ||Ehn,] ||0,∞,Ω ≤ x

1

1

n,] 2 2 ||ξE ||0,Ωx ). A similar estimate also holds for ||Bhn,] ||0,∞,Ω . Combing the results, we have x

  1 1 1 1 n,] 2 n,] 2 n,] 2 − d4x 2 + ||B || ≤ C 1 + h + ||ξ || ) . ||Ehn,] ||0,∞,Ω (||ξ || 0,∞,Ωx 0,Ωx B 0,Ωx E h x

(4.63)

Equations (4.61) - (4.63) lead to ah,2 (ηfn,] , Ehn,] , Bhn,] ; ξfn,] ) ≤ Chk−

dx 2

k+ 12

+ Ch

n,] n,] (||ξE ||0,Ωx + ||ξB ||0,Ωx )||ξfn,] ||0,Ω + Chk+1 ||ξfn,] ||0,Ω



− d4x

1+h

1

n,] 2 (||ξE ||0,Ωx

1

+

 Z

n,] 2 ||ξB ||0,Ωx )

Z

Thx

Ev

! 21 |(Ehn,]

+v×

Bhn,] )

·

nv ||[ξfn,] ]v |2 dsv dx

n,] 2 n,] 2 ≤ Ch2k−dx (||ξE ||0,Ωx + ||ξB ||0,Ωx ) + C||ξfn,] ||20,Ω + Ch2k+2 Z Z   dx 1 n,] n,] ||0,Ωx ) + |(Ehn,] + v × Bhn,] ) · nv ||[ξfn,] ]v |2 dsv dx. ||0,Ωx + ||ξB + Ch2k+1 1 + h− 2 (||ξE 16 Thx Ev dx

n,] 2 n,] n,] n,] 2 n,] 2 Since h2k+1 h− 2 (||ξE ||0,Ωx + ||ξB ||0,Ωx ) ≤ C(h4k+2−dx + ||ξE ||0,Ωx + ||ξB ||0,Ωx ) ≤ C(h2k+2 + ||ξE ||0,Ωx + n,] 2 dx 2k−dx ||ξB ||0,Ωx ) and h ≤ 1 for k ≥ 2 , we further have

ah,2 (ηfn,] , Ehn,] , Bhn,] ; ξfn,] ) ≤ Ch2k+1 + CNn,] +

1 16

Z Thx

Z Ev

|(Ehn,] + v × Bhn,] ) · nv ||[ξfn,] ]v |2 dsv dx.

(4.64)

Now with the results in (4.59) and (4.64), we can conclude (3.19b). To get (3.19c) and (3.19d), note that with f being smooth with compact support in Ωv , one has [f n,] ]v = 0, {f n,] }v = f n,] for any e ∈ Ev . Thus, applying the divergence theorem gives ah,2 (f n,] , Ehn,] , Bhn,] ; g) Z Z Z n,] n,] n,] = f (Eh + v × Bh ) · ∇v gdxdv − f n,] (Ehn,] + v × Bhn,] ) · [g]v dsv dx Ω Ωx Ev Z = − ∇v f n,] · (Ehn,] + v × Bhn,] )gdxdv. Ω

20

Similarly, ah,2 (f n,] , E n,] , B n,] ; g) = −

R Ω

∇v f n,] · (E n,] + v × B n,] )gdxdv. Therefore,

|ah (f n,] , E n,] , B n,] ; g) − ah (f n,] , Ehn,] , Bhn,] ; g)| = |ah,2 (f n,] , E n,] , B n,] ; g) − ah,2 (f n,] , Ehn,] , Bhn,] ; g)| Z = | ∇v f n,] · (Ehn,] − E n,] + v × (Bhn,] − B n,] ))gdxdv| Ω

≤

n,] ||∇v f n,] ||0,∞,Ω (||en,] E ||0,Ω + C||eB ||0,Ω )||g||0,Ω

n,] ≤ C(||en,] E ||0,Ω + C||eB ||0,Ω )||g||0,Ω .

(4.65)

This gives (3.19c). Taking g = ξfn,] , and with the approximation property in (2.4), we further get (3.19d), |ah (f n,] , E n,] , B n,] ; ξfn,] ) − ah (f n,] , Ehn,] , Bhn,] ; ξfn,] )| n,] n,] n,] n,] ≤ C(||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ηE ||0,Ωx + ||ηB ||0,Ωx )||ξfn,] ||0,Ω ≤ Ch2k+2 + CNn,] .

4.3

Proof of Lemma 3.6

To get (3.22a), based on the definition of ah , approximation property in (2.4) and inverse inequality (2.5), |ah (ηfn,r , Ehn,s+1 , Bhn,s+1 ; g) − ah (ηfn,r , Ehn,s , Bhn,s ; g)| = |ah,2 (ηfn,r , Ehn,s+1 , Bhn,s+1 ; g) − ah,2 (ηfn,r , Ehn,s , Bhn,s ; g)| Z ≤ | ηfn,r (Ehn,s+1 − Ehn,s + v × (Bhn,s+1 − Bhn,s ))∇v gdxdv| Ω Z Z +| {ηfn,r }v (Ehn,s+1 − Ehn,s + v × (Bhn,s+1 − Bhn,s )) · [g]v dsv dx| Tx

Ev

Zh Z 1 |[ηfn,r ]v | |(Ehn,s+1 + v × Bhn,s+1 ) · nv | − |(Ehn,s + v × Bhn,s ) · nv | |[g]v |dsv dx + 2 Thx Ev   ≤ C ||Ehn,s+1 − Ehn,s ||0,∞,Ωx + ||Bhn,s+1 − Bhn,s ||0,∞,Ωx (||ηfn,r ||0,Ω ||∇v g||0,Ω + ||ηfn,r ||0,Thx ×Ev ||g||0,Thx ×Ev )   ≤ Chk ||g||0,Ω ||Ehn,s+1 − Ehn,s ||0,∞,Ωx + ||Bhn,s+1 − Bhn,s ||0,∞,Ωx . By the a priori assumption and equation (3.8), we have ||Ehn,s+1 − Ehn,s ||0,∞,Ωx = ||Ehn,s+1 − E n,s+1 + E n,s − Ehn,s + E n,s+1 − E n,s ||0,∞,Ωx ≤

n,s+1 ||en,s+1 ||0,∞,Ωx + ||en,s − E n,s ||0,∞,Ωx ≤ C(h + τ ). E E ||0,∞,Ωx + ||E

(4.66)

Similarly, ||Bhn,s+1 − Bhn,s ||0,∞,Ωx ≤ C(h + τ ). Therefore, |ah (ηfn,r , Ehn,s+1 , Bhn,s+1 ; g) − ah (ηfn,r , Ehn,s , Bhn,s ; g)| ≤ C(1 +

τ k+1 )h ||g||0,Ω . h

To get (3.22b), we follow the similar lines as to prove (3.22a) and obtain |ah (ξfn,r , Ehn,s+1 , Bhn,s+1 ; g) − ah (ξfn,r , Ehn,s , Bhn,s ; g)|   ≤ C ||Ehn,s+1 − Ehn,s ||0,∞,Ωx + ||Bhn,s+1 − Bhn,s ||0,∞,Ωx (||ξfn,r ||0,Ω ||∇v g||0,Ω + ||ξfn,r ||0,Thx ×Ev ||g||0,Thx ×Ev )   τ ≤ C(h + τ ) ||ξfn,r ||0,Ω ||∇v g||0,Ω + ||ξfn,r ||0,Thx ×Ev ||g||0,Thx ×Ev ≤ C(1 + )||ξfn,r ||0,Ω ||g||0,Ω h Here (4.66) is used for the second inequality, and the inverse inequality (2.5) is used for the third one. 21

4.4

Proof of Lemma 3.7

By definition, we have LRK (g) = 2L(g) − K(g) − J (g) Z 3ηfn+1 + 3ηfn − 6ηfn,2 + 3Tfn (x, v)

gdxdv + 2ah (f n,2 , E n,2 , B n,2 ; g) − 2ah (fhn,2 , Ehn,2 , Bhn,2 ; g) Th   − ah (f n,1 , E n,1 , B n,1 ; g) − ah (fhn,1 , Ehn,1 , Bhn,1 ; g) − (ah (f n , E n , B n ; g) − ah (fhn , Ehn , Bhn ; g))

=

τ

=

Λ1 + Λ 2 + Λ 3 + Λ 4 ,

where Λ1 =

1 τ

Z Th



 3ηfn+1 + 3ηfn − 6ηfn,2 + 3Tfn (x, v) gdxdv,

2ah (f n,2 , E n,2 , B n,2 ; g) − 2ah (f n,2 , Ehn,2 , Bhn,2 ; g)   − ah (f n,1 , E n,1 , B n,1 ; g) − ah (f n,1 , Ehn,1 , Bhn,1 ; g) − (ah (f n , E n , B n ; g) − a(f n , Ehn , Bhn ; g)) ,

Λ2

=

Λ3

= −2ah (ηfn,2 , Ehn,2 , Bhn,2 ; g) + ah (ηfn,1 , Ehn,1 , Bhn,1 ; g) + ah (ηfn , Ehn , Bhn ; g),

Λ4

=

2ah (ξfn,2 , Ehn,2 , Bhn,2 ; g) − ah (ξfn,1 , Ehn,1 , Bhn,1 ; g) − ah (ξfn , Ehn , Bhn ; g).

The term Λ1 can be estimated by applying Lemmas 3.2 and 3.3,  1 |Λ1 | ≤ ||3ηfn+1 + 3ηfn − 6ηfn,2 ||0,Ω + 3||Tfn (x, v)||0,Ω ||g||0,Ω ≤ C(hk+1 + τ 3 )||g||0,Ω . τ To estimate Λ2 , we apply (3.19c) in Lemma 3.4, the approximation property (2.4), and obtain |Λ2 | ≤ C

2  X

2    X n,] n,] n,] ||en,] ||ξE ||0,Ωx + ||ξB ||0,Ωx + hk+1 ||g||0,Ω . E ||0,Ωx + ||eB ||0,Ωx ||g||0,Ω ≤ C

]=0

]=0

For Λ3 , we first rewrite it as below. Λ3

= −2ah (ηfn,2 − ηfn,1 , Ehn,2 , Bhn,2 ; g) − ah (ηfn,1 − ηfn , Ehn , Bhn ; g) −2ah (ηfn,1 , Ehn,2 , Bhn,2 ; g) + 2ah (ηfn,1 , Ehn,1 , Bhn,1 ; g) −ah (ηfn,1 , Ehn,1 , Bhn,1 ; g) + ah (ηfn,1 , Ehn , Bhn ; g).

Applying Lemma 3.3 to the first line, and Lemma 3.6 - (1) to the second and the third lines, one has Λ3 ≤ C(1 +

τ k+1 )h ||g||0,Ω . h

For Λ4 , recall Gn2 = 2ξfn,2 − ξfn,1 − ξfn , and using Lemma 3.6 - (2), we have   Λ4 = 2ah (ξfn,2 , Ehn,2 , Bhn,2 ; g) − 2ah (ξfn,2 , Ehn,1 , Bhn,1 ; g) + ah (Gn2 , Ehn,1 , Bhn,1 ; g)   + ah (ξfn , Ehn,1 , Bhn,1 ; g) − ah (ξfn , Ehn , Bhn ; g)  τ  ≤ C(1 + ) ||ξfn,2 ||0,Ω + ||ξfn ||0,Ω ||g||0,Ω + ah (Gn2 , Ehn,1 , Bhn,1 ; g). h Now we combine the estimates for Λi , i = 1, · · · 4, and can get (3.23a),   2 X τ n,] n,] LRK (g) ≤ C (1 + ) (||ξE ||0,Ωx + ||ξB ||0,Ωx + ||ξfn,] ||0,Ω + hk+1 ) + τ 3  ||g||0,Ω + ah (Gn2 , Ehn,1 , Bhn,1 ; g). h ]=0

22

To further bound the last term and hence obtain (3.23b), we apply the inverse inequality (2.5) and ||Ehn,1 ||0,∞,Ωx + ||Bhn,1 ||0,∞,Ωx being bounded by (4.51)-(4.52), and get n,1 n,1 ah (Gn2 , Eh , Bh ; g) ≤ C(||Gn2 ||0,Ω ||∇x g||0,Ω + ||Gn2 ||0,Thv ×Ex ||g||0,Thv ×Ex )  
+C ||Ehn,1 ||0,∞,Ωx + ||Bhn,1 ||0,∞,Ωx ||Gn2 ||0,Ω ||∇v g||0,Ω + ||Gn2 ||0,Thx ×Ev ||g||0,Thx ×Ev C ||Gn2 ||0,Ω ||g||0,Ω . h

≤

The estimate for KRK (g) in (3.23c) can be proved very similarly.

4.5

Proof of Lemma 3.8

First we consider ah,1 (Gn1 ; Gn2 ) + ah,1 (Gn2 ; Gn1 ).

=

ah,1 (Gn1 ; Gn2 ) + ah,1 (Gn2 ; Gn1 ) Z Z (Gn1 v · ∇x Gn2 + Gn2 v · ∇x Gn1 )dxdv −

Thv

Ω

Z

Z Thv

Z

Ex

Z v·

= Thv

Ex

Z

Z

− Thv

([Gn1 Gn2 ]x

(v{Gn1 }x +

Ex

|v · nx | n [G1 ]x ) · [Gn2 ]x dsx dv 2

|v · nx | n [G2 ]x ) · [Gn1 ]x dsx dv 2

(v{Gn2 }x +

−

=

Z

−

{Gn1 }x [Gn2 ]x

−

{Gn2 }x [Gn1 ]x ) dsx dv

Z

Z

− Thv

|v · nx |[Gn1 ]x · [Gn2 ]x dsx dv

Ex

|v · nx |[Gn1 ]x · [Gn2 ]x dsx dv.

(4.67)

Ex

The third equality above is due to (2.3b). Hence, |ah,1 (Gn1 ; Gn2 ) + ah,1 (Gn2 ; Gn1 )| ≤ ≤

1 16

Z Thv

Z



|v · nx | |[ξfn ]x |2 +

Ex



Z

Z

|v · nx | Thv

|[ξfn,1 ]x |2

Ex



dsx dv +

 1 |[Gn1 ]x |2 + 8|[Gn2 ]x |2 dsx dv 32

C ||Gn2 ||20,Ω . h

(4.68)

The last inequality in (4.68) is satisfied because of the inequality |[Gn1 ]x |2 = |[ξfn,1 − ξfn ]x |2 ≤ 2(|[ξfn ]x |2 + |[ξfn,1 ]x |2 ) and the inverse inequality in (2.5). Similarly, we can also show n,1 n,1 n,1 n,1 ah,2 (Gn1 , Eh , Bh ; Gn2 ) + ah,2 (Gn2 , Eh , Bh ; Gn1 ) Z Z = |(Ehn,1 + v × Bhn,1 ) · nv |[Gn1 ]v · [Gn2 ]v dsv dx Thx

Z

Ev

Z

|(Ehn,1 + v × Bhn,1 ) · nv ||[ξfn,1 ]v ||[Gn2 ]v |dsv dx

≤ Thx

Ev

Z

Z

+ Thx

Ev

|(Ehn,1 + v × Bhn,1 ) · nv ||[ξfn ]v ||[Gn2 ]v |dsv dx.

(4.69)

Denote each term on the right side of (4.69) as Λ1 and Λ2 , respectively, then   Z Z 1 n,1 2 n,1 n,1 n 2 |[ξ ]v | + 4|[G2 ]v | dsv dx Λ1 ≤ |(Eh + v × Bh ) · nv | 16 f Thx Ev Z Z ||Gn2 ||20,Ω 1 ≤ |(Ehn,1 + v × Bhn,1 ) · nv ||[ξfn,1 ]v |2 dsv dx + C||Ehn,1 + v × Bhn,1 ||0,∞,Ωx 16 Thx Ev h Z Z 1 C ≤ |(Ehn,1 + v × Bhn,1 ) · nv ||[ξfn,1 ]v |2 dsv dx + ||Gn2 ||20,Ω , (4.70) 16 Thx Ev h 23

here we have used the fact that ||Ehn,1 + v × Bhn,1 ||0,∞,Ωx ≤ C due to (4.51) - (4.52). For the estimate of Λ2 , Z Z n,1 n,1 Λ2 ≤ Eh − Ehn + v × (Bh − Bhn ) [ξfn ]v |[Gn2 ]v | dsv dx Thx

Ev

Z

Z

+ Tx

|(Ehn + v × Bhn ) · nv | [ξfn ]v |[Gn2 ]v | dsv dx

Ev

 h  ≤ C ||Ehn,1 − Ehn ||0,∞,Ωx + ||Bhn,1 − Bhn ||0,∞,Ωx ||ξfn ||0,Thx ×Ev ||Gn2 ||0,Thx ×Ev Z Z 2 1 C + |(Ehn + v × Bhn ) · nv | [ξfn ]v dsv dx + ||Gn2 ||20,Ω . 16 Thx Ev h

(4.71)

As implied in the proof of Lemma 3.6, ||Ehn,1 − Ehn ||0,∞,Ωx ≤ C(h + τ ). So we further have Z Z  2 C 1 τ n ||ξf ||0,Ω ||Gn2 ||0,Ω + ||Gn2 ||20,Ω + |(Ehn + v × Bhn ) · nv | [ξfn ]v dsv dx Λ2 ≤ C 1 + h h 16 Thx Ev   Z Z   2 τ τ 1 1 ≤ C 1+ |(Ehn + v × Bhn ) · nv | [ξfn ]v dsv dx. ||ξfn ||20,Ω + C 1 + + ||Gn2 ||20,Ω + h h h 16 Thx Ev (4.72) Note that 1, hτ ≤ lemma.

4.6

1 h,

hence 1 +

τ h

+

1 h

≤

C h.

We now combine (4.68), (4.70) and (4.72) and conclude this

Proof of Lemma 3.11

We only consider the case when ] = 0. The proof for other cases follows the same line. For part (i), with divergence theorem and equality (2.3c),

n n n n bh (ξE , ξB , ξfn ; ξE , ξB ) Z Z Z Z n n n n n = ξB · (∇ × ξE )dx − ξE · (∇ × ξB )dx − ( vξfn dv) · ξE dx Ωx Ωx Ωx Ωv   Z  Z  1 n 1 n n n n n + {ξB }x − [ξE ]tan · [ξE ]tan dsx − {ξE }x + [ξB ]tan · [ξB ]tan dsx 2 2 Ex Ex Z Z Z 1 n n n n n n n = ([ξE × ξB ]x + {ξB }x [ξE ]tan − {ξE }x [ξB ]tan ) dsx − ( vξfn dv) · ξE dx − STABnEB 2 Ex Ωx Ωv Z Z 1 1 n n n n 2 n 2 =− ( vξf dv) · ξE dx − STABEB ≤ C(||ξf ||0,Ω + ||ξE ||0,Ωx ) − STABnEB . 2 2 Ωx Ωv R R n n n n n For part (ii), since ∇ × ξE ∈ Uhk , there is Ωx ηB · (∇ × ξE )dx = 0. Similarly, Ωx ηE · (∇ × ξB )dx = 0. Then n n n n |bh (ηE , ηB , ηfn ; ξE , ξB )| Z Z Z n =| ( vηfn dv) · ξE dx + Ωx

Ωv

n ≤ C||ηfn ||0,Ω ||ξE ||0,Ωx

   Z  1 n 1 n n n n n {ηB }x − [ηE ]tan · [ξE ]tan dsx − {ηE }x + [ηB ]tan · [ξB ]tan dsx | 2 2 Ex Ex Z  12 n n n n + C(||ηE ||0,Ex + ||ηB ||0,Ex ) |[ξB ]tan |2 + |[ξE ]tan |2 dsx Ex

2k+1

≤ Ch

+

n 2 C||ξE ||0,Ex

1 + STABnEB . 16

24

4.7

Proof of Lemma 3.13

Based on definitions of SRK (U, V ), X2n , Z2n , and Gn2 , in addition to Lemma 3.2 and Lemma 3.3 - (1), we have

≤

SRK (U, V ) = 2S(U, V ) − R(U, V ) − Q(U, V ) ! ! n,2 n,2 n+1 n+1 n n 3ηE − 6ηE + 3ηE + 3TEn (x) 3ηB − 6ηB + 3ηB + 3TBn (x) + ,U ,V τ τ +bh (2en,2 E k+1

≤ C(h

−

en,1 E

−

enE , 2en,2 B

−

en,1 B

−

3

+ τ )(||U ||0,Ωx + ||V ||0,Ωx ) +

Ωx enB , 2en,2 f

n − en,1 f − ef ; U, V ) bh (X2n , Z2n , Gn2 ; U, V ) − bh (dnE , dnB , dnf ; U, V

Ωx

),

(4.73)

where dn? = 2η?n,2 − η?n,1 − η?n with ? = E, B, f . We further apply (3.15c) in Lemma 3.3, and this gives (3.37a), SRK (U, V ) ≤ C(hk+1 + τ hk + τ 3 )(||U ||0,Ωx + ||V ||0,Ωx ) + bh (X2n , Z2n , Gn2 ; U, V ).

(4.74)

To further obtain (3.37b), we need to estimate bh (X2n , Z2n , Gn2 ; U, V ). bh (X2n , Z2n , Gn2 ; U, V ) Z Z Z Z Gn2 vdv) · U dx ( X2n · (∇ × V )dx − Z2n · (∇ × U )dx − = Ωv Ω Ωx Ωx   Z  Z x 1 n 1 n n + {Z2 }x − [X2 ]tan · [U ]tan dsx − {X2 }x + [Z2n ]tan · [V ]tan dsx 2 2 Ex Ex n n n ≤ ||Z2 ||0,Ωx ||∇ × U ||0,Ωx + ||X2 ||0,Ωx ||∇ × V ||0,Ωx + C||G2 ||0,Ω ||U ||0,Ωx +C(||Z2n ||0,Ex + ||X2n ||0,Ex )(||U ||0,Ex + ||V ||0,Ex ). Now we can apply inverse equalities in (2.5), get bh (X2n , Z2n , Gn2 ; U, V ) ≤

C (||Z2n ||0,Ωx + ||X2n ||0,Ωx )(||U ||0,Ωx + ||V ||0,Ωx ) + C||Gn2 ||0,Ω ||U ||0,Ωx , h

hence (3.37b). Similarly, one can get the estimate (3.37c) for RRK (U, V ).

4.8

Proof of Lemma 3.14

From the definition of bh , bh (X1n , Z1n , Gn1 ; X2n , Z2n ) + bh (X2n , Z2n , Gn2 ; X1n , Z1n ) Z = (Z1n · (∇ × X2n )dx − X1n · (∇ × Z2n ) + Z2n · (∇ × X1n )dx − X2n · (∇ × Z1n )) dx Ωx   Z  Z  1 n 1 n n n n {X1 }x + [Z1 ]tan · [Z2n ]tan dsx + {Z1 }x − [X1 ]tan · [X2 ]tan dsx − 2 2 Ex Ex   Z  Z  1 n 1 n n n n + {Z2 }x − [X2 ]tan · [X1 ]tan dsx − {X2 }x + [Z2 ]tan · [Z1n ]tan dsx 2 2 Ex Ex   Z Z Z Z − Gn1 vdv · X2n dx − Gn2 vdv · X1n dx. Ωx

Ωv

Ωx

Ωv

25

(4.75)

Using divergence theorem and equality (2.3c), in addition to inverse inequality (2.5) and Young’s inequality, we can simplify the inequality above as

=

bh (X1n , Z1n , Gn1 ; X2n , Z2n ) + bh (X2n , Z2n , Gn2 ; X1n , Z1n ) Z Z − ([X1n ]tan · [X2n ]tan + [Z1n ]tan · [Z2n ]tan ) dsx − (Gn1 v · X2n + Gn2 v · X1n ) dxdv Z

≤ Ex

≤

Ex

Ω



 2 X 1 1 |[X1n ]tan |2 + 8|[X2n ]tan |2 + |[Z1n ]tan |2 + 8|[Z2n ]tan |2 dsx + C (||Gnj ||20,Ω + ||Xjn ||20,Ωx ) 32 32 j=1

2  X C 1  STABnEB + STABn,1 (||X2n ||20,Ωx + ||Z2n ||20,Ωx ) + C (||Gnj ||20,Ω + ||Xjn ||20,Ωx ) + EB . h 16 j=1

(4.76)

The last inequality of (4.76) is due to that n,1 n,1 n n |[Z1n ]tan |2 = |[ξB − ξB ]tan |2 ≤ 2|[ξB ]tan |2 + 2|[ξB ]tan |2 ,

and a similar bound for |[X1n ]tan |2 .

4.9

Proof of Lemma 3.18

For any g ∈ Ghk , we first consider ah (f n,] , E n,] , B n,] ; g) − ah (fhn,] , Ehn,] , Bhn,] ; g), ] = 0, 1, 2, ah (f n,] , E n,] , B n,] ; g) − ah (fhn,] , Ehn,] , Bhn,] ; g) = −ah (ηfn,] , Ehn,] , Bhn,] ; g) + (ah (f n,] , E n,] , B n,] ; g) − ah (f n,] , Ehn,] , Bhn,] ; g)) + ah (ξfn,] , Ehn,] , Bhn,] ; g). With Cauchy-Schwarz inequality, approximation property in (2.4), inverse inequality (2.5), boundedness of ||Ehn,1 ||0,∞,Ωx + ||Bhn,1 ||0,∞,Ωx in (4.51) and (4.52), we have n,] n,] n,] n,] n,] ah (ηf , Eh , Bh ; g) ≤ C(||ηf ||0,Ω ||∇x g||0,Ω + ||ηf ||0,Thv ×Ex ||g||0,Thv ×Ex ) +C(||Ehn,] ||0,∞,Ωx + ||Bhn,] ||0,∞,Ωx )(||ηfn,] ||0,Ω ||∇v g||0,Ω + ||ηfn,] ||0,Thv ×Ex ||g||0,Thv ×Ex ) ≤ Chk ||g||0,Ω ,

(4.77)

n,] n,] n,] n,] n,] ah (ξf , Eh , Bh ; g) ≤ C(||ξf ||0,Ω ||∇x g||0,Ω + ||ξf ||0,Thv ×Ex ||g||0,Thv ×Ex ) +C(||Ehn,] ||0,∞,Ωx + ||Bhn,] ||0,∞,Ωx )(||ξfn,] ||0,Ω ||∇v g||0,Ω + ||ξfn,] ||0,Thv ×Ex ||g||0,Thv ×Ex ) ≤

C n,] ||ξ ||0,Ω ||g||0,Ω . h f

(4.78)

Following the derivation to get (4.65) in Lemma 3.4, we have ah (f n,] , E n,] , B n,] ; g) − ah (f n,] , Ehn,] , Bhn,] ; g) Z   = ∇v f n,] · Ehn,] − E n,] + v × (Bhn,] − B n,] ) gdxdv Ω

n,] n,] n,] k+1 ≤ C(||en,] )||g||0,Ω . E ||0,Ωx + ||eB ||0,Ωx )||g||0,Ω ≤ C(||ξE ||0,Ωx + ||ξB ||0,Ωx + h

(4.79)

Equations (4.77) - (4.79) lead to ah (f n,] , E n,] , B n,] ; g) − ah (fhn,] , Ehn,] , Bhn,] ; g) C n,] n,] ≤ C(||ξE ||0,Ωx + ||ξB ||0,Ωx + hk )||g||0,Ω + ||ξfn,] ||0,Ω ||g||0,Ω . h

26

(4.80)

Now based on (4.80) and Lemma 3.3, we can bound J (ξfn,1 ). ! ηfn,1 − ηfn n,1 n,1 + ah (f n , E n , B n ; ξfn,1 ) − ah (fhn, , Ehn , Bhn ; ξfn,1 ) J (ξf ) = , ξf τ Ω

n C(||ξE ||0,Ωx

≤

+

n ||ξB ||0,Ωx

+ hk )||ξfn,1 ||0,Ω +

C n ||ξ ||0,Ω ||ξfn,1 ||0,Ω . h f

(4.81)

Based on equations (4.81) and (3.10), we get ||ξfn,1 ||20,Ω

=

(ξfn,1 , ξfn )Ω + τ J (ξfn,1 )

≤

τ n n Cτ (||ξE ||0,Ωx + ||ξB ||0,Ωx + hk )||ξfn,1 ||0,Ω + (1 + C )||ξfn ||0,Ω ||ξfn,1 ||0,Ω , h

Cancel ||ξfn,1 ||0,Ω from both sides of the inequality, and take the square of both sides, we have ||ξfn,1 ||20,Ω

τ n 2 n 2 ≤ Cτ 2 (||ξE ||0,Ωx + ||ξB ||0,Ωx + h2k ) + (1 + C )2 ||ξfn ||20,Ω . h

(4.82)

By further using τ ≤ αh, we obtain the upper bound of ||ξfn,1 ||20,Ω . Similarly, ||ξfn,2 ||20,Ω can be estimated n,] n,] based on (4.80) and (3.10). Likewise, we can establish the estimates for ||ξB ||0,Ωx , ||ξE ||0,Ωx for ] = 1, 2.

5

Extension and Conclusion

In this paper we prove the error estimates of fully discrete methods, which involve a third order Runge-Kutta time discretization and upwind DG discretizations of arbitrary order of accuracy in phase domain, for solving the Vlasov-Maxwell system. When the exact solutions have enough regularity, we show that the L2 errors of 1 the numerical solutions by such methods are of O(hk+ 2 + τ 3 ) for k > dx2+1 . The third order Runge-Kutta 1 time integration contributes to the error O(τ 3 ), while the error from the DG approximation is O(hk+ 2 ), which is expected for hyperbolic systems with upwind numerical fluxes on general meshes. The techniques used in this paper can be applied to the RKDG methods which involve other numerical fluxes, such as central or alternating fluxes, (nx\ × Eh , nx\ × Bh ) = (nx × {Eh }, nx × {Bh }), (nx\ × Eh , nx\ × Bh ) = nx ×

(Eh− , Bh+ ), or

nx ×

(central)

(Eh+ , Bh− ),

(alternating),

in the Maxwell solver. It was shown that these fluxes will result better energy conservation in semi-discrete DG methods [5]. On the other hand, for RKDG methods with such fluxes, the stabilization mechanism STABn,] EB , ] = 0, 1, 2, in the form of the tangential jump Z n,] n,] |[ξE ]tan |20,Ωx + |[ξB ]tan |20,Ωx dsx Ex

is no longer available from the Maxwell solver (see part (i) of Lemma 3.11), and this in general will lead to a sub-optimal L2 -norm error estimate: Chk + Cτ 3 . (For alternating fluxes, better estimates can be obtained on Cartesian meshes and with the tensor structured polynomial discrete spaces.) With some insignificant modification to the details, almost the same error estimates can be established for the RKDG methods solving the smooth solutions of the relativistic Vlasov-Maxwell system of one species [21, 22],  v v  · ∇x f + (E + p × B) · ∇v f = 0, ∂t f + p   2  1 + |v| 1 + |v|2 ∂t E = ∇ × B − J, ∂t B = −∇ × E,     ∇ · E = ρ − ρi , ∇ · B = 0, with

Z ρ(x, t) =

Z f (x, v, t)dv,

J(x, t) =

Ωv

Ωv

27

v p

1 + |v|2

f (x, v, t)dv.

References [1] B. Ayuso, J.A. Carrillo and C.-W. Shu, Discontinuous Galerkin methods for the one-dimensional VlasovPoisson system, Kinetic and Related Models, 4 (2011), pp.955-989. [2] B. Ayuso, J.A. Carrillo and C.-W. Shu, Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem, Mathematical Models and Methods in Applied Sciences, 22 (2012), 1250042 (45 pages). [3] C. K. Birdsall and A. B. Langdon, Plasma Physics Via Computer Simulation, McGraw-Hill, New York, 1985. [4] Y. Cheng and I. M. Gamba, Numerical study of one-dimensional Vlasov-Poisson equations for infinite homogeneous stellar systems, Communications in Nonlinear Science and Numerical Simulation, 17 (2012), pp.2052-2061. [5] Y. Cheng, I. M. Gamba, F. Li and P. J. Morrison, Discontinuous Galerkin methods for Vlasov-Maxwell equations, SIAM Journal on Numerical Analysis, accepted, 2014. [6] C. Z. Cheng and G. Knorr, The integration of the Vlasov equation in configuration space, Journal of Computational Physics, 22 (1976), pp.330-351. [7] P. G. Ciarlet, Finite element method for elliptic problems, Noth-Holland, Amsterdam, 1978. [8] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, 52 (1989), pp.411435. [9] B. Cockburn and C.-W. Shu, The Runge-Kutta local projection p1-discontinuous Galerkin finite element method for scalar conservation laws, ESAIM Mathematical Modelling and Numerical Analysis, 25 (1991), pp.337-361. [10] G. H. Cottet and P. A. Raviart, On Particle in Cell methods for the Vlasov-Poisson equations, Transport Theory and Statistical Physics, 15 (1986), pp.1-31 [11] M. Crouzeix and V. Thom´ee, The stability in Lp and Wp1 of the L2 -projection onto finite element function spaces, Mathematics of Computation, 178 (1987), pp.521-532 [12] B. Eliasson, Numerical modelling of the two-dimensional Fourier transformed Vlasov-Maxwell system, Journal of Computational Physics, 190 (2003), pp.501-522 [13] B. Eliasson, Numerical simulations of the Fourier-transformed Vlasov-Maxwell system in higher dimensions-theory and applications, Transport Theory and Statistical Physics, 39 (2011), pp.387-465 [14] R. E. Heath, I. M. Gamba, P. J. Morrison and C. Michler, A discontinuous Galerkin method for the Vlasov-Poisson system, Journal of Computational Physics, 231 (2012), pp.1140-1174. [15] G. Jacobs and J.S. Hesthaven, Implicit-explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning, Computer Physics Communications, 180 (2009), pp.1760-1767 [16] A. J. Klimas and W. M. Farrell, A splitting algorithm for Vlasov simulation with filamentation filtration, Journal of Computational Physics, 110 (1994), pp.150-163. [17] M. C. Pinto and M. Mehrenberger, Convergence of an adaptive semi-Lagrangian scheme for the VlasovPoisson system, Numerische Mathematik, 108 (2008), pp.407-444.

28

[18] J.-M. Qiu and C.-W. Shu, Conservative high order semi-Lagrangian finite difference WENO methods for advection in incompressible flow, Journal of Computational Physics, 230 (2011), pp.863-889. [19] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, Journal of Computational Physics, 77 (1988), pp.439-471. [20] E. Sonnendr¨ ucker, J. Roche, P. Bertrand and A. Ghizzo, The semi-Lagrangian method for the numerical resolution of the Vlasov equation, Journal of Computational Physics, 149 (1999), pp.201-220. [21] A. Suzuki and T. Shigeyama, A conservative scheme for the relativistic Vlasov-Maxwell system, Journal of Computational Physics, 229 (2010), pp.1643-1660. [22] H. Yang and F. Li, Discontinuous Galerkin methods for relativitistic Vlasov-Maxwell equations, in preparation. [23] H. Yang, F. Li, and J. Qiu, Dispersion and dissipation errors of two fully discrete discontinuous Galerkin methods, Journal of Scientific Computing, 55 (2013), pp.552-574. [24] S. I. Zaki, L. R. T. Gardner and T. J. M. Boyd, A finite element code for the simulation of onedimensional Vlasov plasmas. i. theory, Journal of Computational Physics, 79 (1988), pp.184-199. [25] S. I. Zaki, T. J. M. Boyd and L. R. T. Gardner, A finite element code for the simulation of onedimensional Vlasov plasmas. ii. applications, Journal of Computational Physics, 79 (1988), pp.200-208. [26] Q. Zhang and C-W Shu, Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws, SIAM Journal on Numerical Analysis, 42 (2004), pp.641-666. [27] Q. Zhang and C-W Shu, Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for symmetrizable systems of conservation laws, SIAM Journal on Numerical Analysis, 44 (2006), pp.1703-1720. [28] Q. Zhang and C-W Shu, Stability analysis and a priori error estimates to the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws, SIAM Journal on Numerical Analysis, 48 (2010), pp.1038-1063. [29] X. Zhong and C.-W. Shu, Numerical resolution of discontinuous Galerkin methods for time dependent wave equations, Computer Methods in Applied Mechanics and Engineering, 200 (2011), pp. 2814-2827.

29