Ludovic Moya

Jan Verwer

Temporal convergence of a locally implicit discontinuous Galerkin method for Maxwell's equations

Keywords: temporal convergence, discontinuous Galerkin method, time-domain Maxwell equations, component splitting, order reduction convergence temporelle, méthode Galerkin discontinue, équations de Maxwell, décomposition des éléments, réduction d'ordre

In this paper we study the temporal convergence of a locally implicit discontinuous Galerkin method for the time-domain Maxwell's equations modeling electromagnetic waves propagation. Particularly, we wonder whether the method retains its second-order ordinary differential equation (ODE) convergence under stable simultaneous space-time grid refinement towards the true partial differential equation (PDE) solution. This is not a priori clear due to the component splitting which can introduce order reduction.

Introduction

Nowadays, many different types of methods exist for the numerical resolution of time-domain Maxwell's equations modeling electromagnetic wave propagation. The most prominent method among physicists and engineers is still the finite difference time-domain (FDTD) method based on Yee's scheme [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF]. This popularity is mainly due to its simplicity and efficiency in discretising simple domain problems. However, its inability to effectively handle complex geometries has prompted to search for alternatives methods. Also one of the main features of numerical methods based on finite element meshes like finite element time-domain (FETD) [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF], finite volume time-domain (FVTD) [START_REF] Remaki | A new finite volume scheme for solving Maxwell's system[END_REF] or discontinuous Galerkin time-domain (DGTD) [START_REF] Cohen | A spatial high order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time-domain[END_REF][START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF][START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF] methods is the possibility of using locally refined and non-conformal space grids to easily deal with complex geometries. In recent years there has been an increasing interest in the DGTD method. The latter is particularly well suited to the design of hp-adaptive strategies (i.e. where the characteristic mesh size h and the interpolation degree p change locally wherever it is needed) [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]. Thus the DGTD method can achieve a high order of accuracy and is used in many applications [START_REF]Discontinuous Galerkin methods. Theory, computation and applications[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF].

In the same time the choice of the temporal integration method is a crucial step for the global efficiency of the numerical method. We distinguish two major families for the temporal integration: implicit and explicit methods. An implicit integration method to numerically solve a time-dependent PDE leads in general to unconditional stability. Then the time step can be chosen arbitrarily large. However, implicit methods require the solution of large linear systems resulting in a high computational effort. An explicit integration method results in less computational effort per time step, but readily leads to unduly step size restrictions caused by the smallest grid elements. For examples of implicit and explicit integration methods for the semi-discrete Maxwell's equations we refer to [START_REF] Catella | An Unconditionally Stable Discontinuous Galerkin Method for Solving the 2-D Time-Domain Maxwell Equations on Unstructured Triangular Meshes[END_REF][START_REF] Verwer | Unconditionaly stable integration of Maxwell's equations[END_REF] and [START_REF] Botchev | Numerical Integration of Damped Maxwell Equations[END_REF][START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF], respectively. A possible alternative to overcome step size limitation, induced by local mesh refinements in explicit methods, is to use smaller time steps precisely where the smallest elements are located, given by a local stability criterion. These local time-stepping techniques have been recently studied for second-order wave equations discretized in space by continuous or discontinuous finite element methods in [START_REF] Diaz | Energy conserving explicit local time stepping for second.order wave equations[END_REF] and more specifically for the time-domain Maxwell's equations discretized in space by a discontinuous Galerkin method in [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF][START_REF] Montseny | Dissipative terms and local time-stepping improvements in a spatial high order Discontinuous Galerkin scheme for the time-domain Maxwell's equations[END_REF][START_REF] Taube | A high order discontinuous Galerkin method with local time stepping for the Maxwell equations[END_REF]. Considering the strengths of implicit and explicit methods, the authors of [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF] and [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] have proposed another alternative: locally implicit time integration methods. More precisely, the smallest grid elements are treated implicitly and the remaining elements explicitly by a technique we call component splitting. If the ratio of fine to coarse elements is small, the most severe step size restrictions are overcome. The counterpart of this approach is having to solve per time step a linear system. But due to the assumed small fine to coarse cell size ratio, the overhead will also be small while the solution can be advanced in time with step sizes determined by the coarse elements. Consequently, these methods are particularly well suited if the local refinement is strongly localized. Note that the method from [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF] has been especially designed for a discontinuous Galerkin discretization, while that from [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] covers the common spatial discretizations like finite-difference and various finite-element discretizations.

In this paper we study the temporal convergence of the locally implicit DGTD method for Maxwell's equations initially proposed by Piperno in [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF]. In particular we examine whether the method retains its second-order ODE convergence towards the true PDE solution under stable simultaneous space-time grid refinement. This question is legitimate because component splitting can cause order reduction.

The paper is organized as follows. Section 1 presents the problem and the notations. In Section 2 we introduce the component splitting and the implicit-explicit method for time-domain Maxwell's equations spatially discretized with a discontinuous Galerkin method. This implicit-explicit integration method is a blend of the explicit one-step second-order Leap-Frog scheme (LF2) and the implicit one-step secondorder Crank-Nicolson scheme (CN2). Next we study the temporal convergence order. The different steps of this convergence analysis are the same as in [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] where the author proposes another implicit-explicit integration method for the semi-discrete Maxwell's equations. The latter method is also a blend of LF2 and CN2 and retains its second-order ODE convergence towards the true PDE solution under stable si-multaneous space-time grid refinement [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF]. In Section 3 we give some numerical results for 2D Maxwell problems to illustrate the previous convergence analysis. Section 4 concludes the paper with final remarks and future plans.

Problem statement

We consider time-domain Maxwell's equations

   ε ∂ t E = curl H -σ E -J E , µ ∂ t H = -curl E, (1)
where E and H denote the electric and magnetic field, respectively. J E is the given source current and ε, µ and σ are coefficients representing dielectric permittivity, magnetic permeability and conductivity, respectively. After discretization in space by a DG method we obtain the semi-discrete Maxwell system

   M ε ∂ t E = SH -DE + M ε f E , M µ ∂ t H = -S T E + M µ f H , (2)
where, for convenience, we use the same notation for the electric and magnetic fields E and H as in the space-continuous case. For more details on DG spatial discretization we refer to [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF]. The matrices M ε , M µ are the DG mass matrices which contain the values of the dielectric permittivity and magnetic permeability coefficient. The matrix S emanates from the discretization of the curl operator. The matrix D is associated with the dissipative conduction term -σ E. Throughout D may be assumed symmetric positive semi-definite. The functions f E and f H are associated with source terms. More precisely f E represents the given source current -J E , but f E and f H may also contain Dirichlet boundary data.

We can give an equivalent formulation of (2) without mass matrix. As in [START_REF] Botchev | Numerical Integration of Damped Maxwell Equations[END_REF] we introduce the Cholesky factorizations

M ε = L M ε L T M ε and M µ = L M µ L T M µ , (3)
where L M ε and L M µ are triangular matrices. Then with (2) we have

 

Thus we can write the semi-discrete Maxwell system equivalent to (2) as

   ∂ t Ẽ = S H -D Ẽ + f E , ∂ t H = -ST Ẽ + f H . (8)
For convenience of notation and presentation we use the same notation in (2) and ( 8) i.e.

   ∂ t E = SH -DE + f E , ∂ t H = -S T E + f H . (9)
We will proceed with ( 9), the meaning of E, H, S, D, f E and f H will always be clear from the context or will be precised. In particular results obtained for (9) apply to (2) and vice versa. Note that S emanates from an appropriate DG discretization for the Maxwell problem under consideration and further on we will prove that

S ∼ 1 h , for h → 0, ( 10
)
where the parameter h denotes the maximum diameter of the (non-uniform) grid elements. Throughout the remainder we will assume initial values at time t = 0 and the source functions

f E (t), f H (t) ∈ C 2 [0, T ] for an interval [0, T ], so that E (t), H (t) ∈ C 3 [0, T ].

The implicit-explicit DGTD method

The implicit-explicit integration considered in this note is issued from [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF]. As we have previously mentioned it is a blend of the second order LF2 scheme that we write in the three-stage form, emanating from Verlet's method, see [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF]

                   H n+ 1 2 -H n ∆t/2 = -S T E n + f H (t n ) , E n+1 -E n ∆t = SH n+ 1 2 - 1 2 D E n+1 + E n + 1 2 f E (t n+1 ) + f E (t n ) , H n+1 -H n+ 1 2 ∆t/2 = -S T E n+1 + f H (t n+1 ) , (11)
and the second order, unconditionally stable CN2 scheme that we also write in the three-stage form

                   H n+ 1 2 -H n ∆t/2 = -S T E n + f H (t n ) , E n+1 -E n ∆t = 1 2 S H n+1 + H n - 1 2 D E n+1 + E n + 1 2 f E (t n+1 ) + f E (t n ) , H n+1 -H n+ 1 2 ∆t/2 = -S T E n+1 + f H (t n+1 ) . (12)
which only differ in the middle stage in the time level for H (see [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF]). Herein ∆t = t n+1t n denotes the step size and upper indices time levels, as usual.

Component splitting

The set of grid elements is assumed to be partitioned into two subsets: one made of the smallest elements that will be treated implicitly using the CN2 method and the other one of the remaining elements that will be treated explicitly with the LF2 method. In line with this splitting the problem unknowns are reordered as

E =   E e E i   and H =   H e H i   , (13)
where the indices i and e are associated to the elements of the subset treated implicitly and explicitly, respectively. Likewise the semi-discrete curl operator S is split into the block form

S =   S e -A ei -A ie S i   , (14)
for the specific meaning of the block-entries of S we refer to [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF]. D is supposed to be split accordingly into

D =   D e 0 0 D i   , (15)
and we write

f E =   f E e f E i   , f H =   f H e f H i   . (16)
Inserting this splitting into (9) we obtain the system of ODEs

                   ∂ t E e = S e H e -A ei H i -D e E e + f E e (t) , ∂ t E i = S i H i -A ie H e -D i E i + f E i (t) , ∂ t H e = -S T e E e + A T ie E i + f H e (t) , ∂ t H i = -S T i E i + A T ei E e + f H i (t) . (17)
Inria

The implicit-explicit time integration method

The implicit-explicit method proposed in [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF] is a blend of LF2 and CN2 applied to [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]. It reads

         H n+1/
                     E n+1 i -E n i ∆t = S i H n+1 i + H n i 2 -A ie H n+1/2 e -D i E n+1 i + E n i 2 + f E i (t n+1 ) + f E i (t n ) 2 , H n+1 i -H n i ∆t = -S T i E n+1 i + E n i 2 + A T ei E n+1/2 e + f H i (t n+1 ) + f H i (t n ) 2 ,          E n+1 e -E n+1
For the stability analysis of this method we refer to [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF]. The proof is based on the conservation of a discrete electromagnetic energy. It is proved that this energy is a positive quadratic form of the numerical unknowns E n e , E n i , H n e and H n i under a condition on the time step size. Consequently the non-dissipative nature of the method yields to the stability of the method.

Matrix behavior for h → 0

In this subsection we are interested in the behavior of the matrices in [START_REF] Yu | Local theory of extrapolation methods[END_REF] for h → 0. This is an essential point for convergence analysis because some of these matrices can lie at the origin of order reduction. Let us consider the general case of dimension d (d = 1, 2 or 3). First we investigate the behavior of the matrices in the formulation with the mass matrices. Thereafter we will be able to deduce the behavior of the matrices in [START_REF] Yu | Local theory of extrapolation methods[END_REF].

We reintroduce the notation with a tilde for the elements involved in the formulation without mass matrix (see [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]) in order to avoid confusion. The specific meaning of the block-entries of the different matrices involved in the formulation with mass matrices can be found in [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF]. First we observe that the mass matrices are only composed of volumic terms, hence we have

M ε , M µ h d , for h → 0. ( 19
)
Thus L M ε , L T M ε , L M µ and L T M µ h d 2 , for h → 0. ( 20
)
The matrices S e and S i are composed of volumic and surfacic terms, hence

S e and S i h d-1 , for h → 0. ( 21
)
The matrices A ei and A ie represent surfacic terms (interface matrices). Hence A ei and A ie h d-1 , for h → 0 [START_REF] Nédélec | A new dfamily of mixed finite elements in R 3[END_REF] and from the block form of the matrix S (see ( 14)) we deduce with ( 21) and ( 22) that

S h d-1 , for h → 0. ( 23
)
From [START_REF] Diaz | Energy conserving explicit local time stepping for second.order wave equations[END_REF] we get

S = L -1 M ε S L -1 M µ T , (24)
then with the behaviors above we deduce that for h → 0

Se , Si = O 1 h , Ãei , Ãie = O 1 h , (25)
and we have the expected behavior (10) for S.

Temporal convergence

In this section we are interested in the PDE convergence of method [START_REF] Yu | Local theory of extrapolation methods[END_REF]. More precisely, we will examine whether the method retains its second-order ODE convergence under stable simultaneous space-time grid refinement ∆t h, h → 0 towards the true PDE solution. This is not a priori clear due to the component splitting which can introduce order reduction through error constants which grow with h -1 , for h → 0. This section is organized in four subsections. In Subsection 2.4.1 we will eliminate the intermediate values of [START_REF] Yu | Local theory of extrapolation methods[END_REF] to get an equivalent one step formula from t n to t n+1 that we will use for our convergence analysis. In Subsection 2.4.2 we will introduce the perturbed method obtained by substituting the true PDE solution restricted to the assumed space grid into [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF], and defects (space-time truncation errors) obtained by substituting this true PDE solution into the equivalent one step formula of our method. In Subsection 2.4.3 we will define the common one-step recurrence relation for the global error. In Subsection 2.4.4 we will point out the order reduction mentioned above. Finally in Subsection 2.4.5 we will see that this order reduction, affecting the local error, may (partly) cancel in the transition from the local to the global error.

Elimination of intermediates values

First we treat H e . From the first and last equations of [START_REF] Yu | Local theory of extrapolation methods[END_REF] we get

H n+ 1 2 e = H n e - ∆t 2 S T e E n e + ∆t 2 A T ie E n i + ∆t 2 f H e (t n ) , H n+ 1 2 e = H n+1 e + ∆t 2 S T e E n+1 e - ∆t 2 A T ie E n+1 i + ∆t 2 f H e (t n+1 ) . (26)
Inserting the first equation of ( 26) into the last equation of ( 18) yields

H n+1 e = H n e - ∆t 2 S T e E n e + E n+1 e + ∆t 2 A T ie E n i + E n+1 i + ∆t 2 f H e (t n ) + f H e (t n+1 ) . (27)
Next we treat E e . From the second and fifth equations of [START_REF] Yu | Local theory of extrapolation methods[END_REF] we get

E n+ 1 2 e = E n e + ∆t 2 S e H n+ 1 2 e - ∆t 2 A ei H n i - ∆t 2 D e E n e + ∆t 2 f E e (t n ) , E n+ 1 2 e = E n+1 e - ∆t 2 S e H n+ 1 2 e + ∆t 2 A ei H n+1 i + ∆t 2 D e E n+1 e - ∆t 2 f E e (t n+1 ) . (28)

Inria

Inserting the first equation of ( 28) and half of each expression of ( 26) for

H n+ 1 2 e
into the fifth equation of ( 18) yields

E n+1 e = E n e + ∆t 2 S e H n e + H n+1 e - ∆t 2 A ei H n i + H n+1 i - ∆t 2 D e E n e + E n+1 e + ∆t 2 f E e (t n ) + f E e (t n+1 ) + ∆t 2 4 S e -S T e E n e + A T ie E n i + ∆t 2 4 S e S T e E n+1 e -A T ie E n+1 i + ∆t 2 4 S e f H e (t n ) + f H e (t n+1 ) . (29)
Now we consider H i . Inserting half of each expression of ( 28) for

E n+ 1 2 e
in the fourth equation of ( 18) gives

H n+1 i = H n i - ∆t 2 S T i E n i + E n+1 i + ∆t 2 A T ei E n e + E n+1 e + ∆t 2 f H i (t n ) + f H i (t n+1 ) + ∆t 2 4 A T ei A ei H n+1 i -H n i + ∆t 2 4 A T ei D e E n+1 e -E n e + ∆t 2 4 A T ei f E e (t n ) + f E e (t n+1 ) . (30)
Finally we treat E i . Inserting half of each expression of ( 26) for

H n+ 1 2 e
in the third equation of ( 18) yields

E n+1 i = E n i + ∆t 2 S i H n i + H n+1 i - ∆t 2 A ie H n e + H n+1 e - ∆t 2
D i E n i + E n+1 i + ∆t 2 f E i (t n ) + f E i (t n+1 ) - ∆t 2 4 A ie -S T e E n e + A T ie E n i - ∆t 2 4 A ie S T e E n+1 e -A T ie E n+1 i - ∆t 2 4 A ie f H e (t n ) + f H e (t n+1 ) . (31)
The equivalent method of [START_REF] Yu | Local theory of extrapolation methods[END_REF]

E n+1 i = E n i + ∆t 2 S i H n i + H n+1 i - ∆t 2 A ie H n e + H n+1 e - ∆t 2
D i E n i + E n+1 i + ∆t 2 f E i (t n ) + f E i (t n+1 ) - ∆t 2 4 A ie -S T e E n e + A T ie E n i -S T e E n+1 e -A T ie E n+1 i - ∆t 2 4 A ie f H e (t n ) + f H e (t n+1 ) , (32)
H n+1 e = H n e - ∆t 2 S T e E n e + E n+1 e + ∆t 2 A T ie E n i + E n+1 i + ∆t 2 f H e (t n ) + f H e (t n+1
) ,

H n+1 i = H n i - ∆t 2 S T i E n i + E n+1 i + ∆t 2 A T ei E n e + E n+1 e + ∆t 2 f H i (t n ) + f H i (t n+1 ) + ∆t 2 4 A T ei A ei H n+1 i -H n i + ∆t 2 4 A T ei D e E n+1 e -E n e + f E e (t n ) + f E e (t n+1 ) .

The perturbed method and defects for the PDE solution

Let E h e (t) denote at time t the true solution of the PDE problem restricted to the assumed space grid that we have approximated with the semi-discrete system [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]

                         d dt E h e (t) = S e H h e (t) -A ei H h i (t) -D e E h e (t) + f E e (t) + σ E e (t) , d dt E h i (t) = S i H h i (t) -A ie H h e (t) -D i E h i (t) + f E i (t) + σ E i (t) , d dt H h e (t) = -S T e E h e (t) + A T ie E h i (t) + f H e (t) + σ H e (t) , d dt H h i (t) = -S T i E h i + A T ei E h e (t) + f H i (t) + σ H i (t) , (33)
where

∆t 2 A ei H h i (t n ) + H h i (t n+1 ) - ∆t 2 D e E h e (t n ) + E h e (t n+1 ) + ∆t 2 f E e (t n ) + f E e (t n+1 ) + ∆t 2 4 S e -S T e E h e (t n ) + A T ie E h i (t n ) -S T e E h e (t n+1 ) -A T ie E h i (t n+1 ) + ∆t 2 4 S e f H e (t n ) + f H e (t n+1 ) + ∆t δ E e,n , E h i (t n+1 ) = E h i (t n ) + ∆t 2 S i H h i (t n ) + H h i (t n+1 ) - ∆t 2 A ie H h e (t n ) + H h e (t n+1 ) - ∆t 2 D i E h i (t n ) + E h i (t n+1 ) + ∆t 2 f E i (t n ) + f E i (t n+1 ) - ∆t 2 4 A ie -S T e E h e (t n ) + A T ie E h i (t n ) -S T e E h e (t n+1 ) -A T ie E h i (t n+1 ) - ∆t 2 4 A ie f H e (t n ) + f H e (t n+1 ) + ∆t δ E i,n , H h e (t n+1 ) = H h e (t n ) - ∆t 2 S T e E h e (t n ) + E h e (t n+1 ) + ∆t 2 A T ie E h i (t n ) + E h i (t n+1 ) + ∆t 2 f H e (t n ) + f H e (t n+1 ) + ∆t δ H e,n ,
H h i (t n+1 ) = H h i (t n ) - ∆t 2 S T i E h i (t n ) + E h i (t n+1 ) + ∆t 2 A T ei E h e (t n ) + E h e (t n+1 ) + ∆t 2 f H i (t n ) + f H i (t n+1 ) + ∆t 2 4 A T ei A ei H h i (t n+1 ) -H h i (t n ) + ∆t 2 4 A T ei D e E h e (t n+1 ) -E h e (t n ) + f E e (t n ) + f E e (t n+1 ) + ∆t δ H i,n , (34)
where

δ E i,n = δ E h i - ∆t 4 A ie d dt H h e (t n+1 ) -H h e (t n ) + 1 2 σ E i (t n ) + σ E i (t n+1 ) - ∆t 4 A ie σ H e (t n ) + σ H e (t n+1 ) , δ H e,n = δ H h e + 1 2 σ H e (t n ) + σ H e (t n+1 ) (35)
δ H i,n = δ H h i - ∆t 4 A T ei S e H h e (t n+1 ) -H h e (t n ) + ∆t 4 A T ei d dt E h e (t n+1 ) -E h e (t n ) + 1 2 σ H i (t n ) + σ H i (t n+1 ) + ∆t 4 A T ei σ E e (t n ) + σ E e (t n+1 ) ,
where δ E h e denotes the implicit trapezoidal (CN) rule defect (see [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF]) 34) we obtain the error scheme

for
ε E e,n+1 = ε E e,n + ∆t 2 S e ε H e,n + ε H e,n+1 - ∆t 2 A ei ε H i,n + ε H i,n+1 - ∆t 2 D e ε E e,n + ε E e,n+1 + ∆t 2 4 S e -S T e ε E e,n + A T ie ε E i,n -S T e ε E e,n+1 -A T ie ε E i,n+1 + ∆t δ E e,n , ε E i,n+1 = ε E i,n + ∆t 2 S i ε H i,n + ε H i,n+1 - ∆t 2 A ie ε H e,n + ε H e,n+1 - ∆t 2 D i ε E i,n + ε E i,n+1 - ∆t 2 4 A ie -S T e ε E e,n + A T ie ε E i,n -S T e ε E e,n+1 -A T ie ε E i,n+1 + ∆t δ E i,n , ε H e,n+1 = ε H e,n - ∆t 2 S T e ε E e,n + ε E e,n+1 + ∆t 2 A T ie ε E i,n + ε E i,n+1 + ∆t δ H e,n , ε H i,n+1 = ε H i,n - ∆t 2 S T i ε E i,n + ε E i,n+1 + ∆t 2 A T ei ε E e,n + ε E e,n+1 + ∆t 2 4 A T ei A ei ε H i,n+1 -ε H i,n + ∆t 2 4 A T ei D e ε E e,n+1 -ε E e,n + ∆t δ H i,n . (37)
Let

ε n =          ε E e,n ε E i,n ε H e,n ε H i,n          and δ n =          δ E e,n δ E i,n δ H e,n δ H i,n          , (38)
then from (37) we can write the global error in a more compact form (one-step recurrence relation)

ε n+1 = Rε n + ∆tρ n , R = R -1 L R R , ρ n = R -1 L δ n , (39)
where

R L =               I + ∆t 2 D e - ∆t 2 4 S e S T e ∆t 2 4 S e A T ie - ∆t 2 S e ∆t 2 A ei ∆t 2 4 A ie S T e I + ∆t 2 D i - ∆t 2 4 A ie A T ie ∆t 2 A ie - ∆t 2 S i ∆t 2 S T e - ∆t 2 A T ie I 0 - ∆t 2 A T ei - ∆t 2 4 A T ei D e ∆t 2 S T i 0 I - ∆t 2 4 A T ei A ei               , (40)
R R =               I - ∆t 2 D e - ∆t 2 4 S e S T e ∆t 2 4 S e A T ie ∆t 2 S e - ∆t 2 A ei ∆t 2 4 A ie S T e I - ∆t 2 D i - ∆t 2 4 A ie A T ie - ∆t 2 A ie ∆t 2 S i - ∆t 2 S T e ∆t 2 A T ie I 0 ∆t 2 A T ei - ∆t 2 4 A T ei D e - ∆t 2 S T i 0 I - ∆t 2 4 A T ei A ei               , (41)
and ε n , ∆tρ n and δ n are respectively the (space-time) global, local and truncation errors. Note that the recursion (39) has the standard form (see e.g. [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF]) for the convergence analysis of one-step integration methods. It transfers local errors to the global error, essentially by adding all local errors. Indeed, for a given time interval [0, T ] we make the usual stability hypothesis R n ≤ K for h → 0 and n ≥ 0, n∆t ≤ T.

(42)

On the other hand the elaboration of the error recursion (39) gives

ε n = R n ε 0 + R n-1 ∆tρ 0 + ... + R∆tρ n-2 + ∆tρ n-1 , (43)
which leads directly (with (42)) to

ε n ≤ K ε 0 + K∆t n-1 ∑ j=0 ρ j for n∆t ≤ T. (44)
Recall that ρ j = R -1 L δ j and because we assume stability we may consider R L inversely bounded for ∆t ∼ h, h → 0. Consequently if δ j = O ∆t k we have ρ j ≤ C∆t k . Assuming ε 0 = 0, we deduce from (44) that

ε n ≤ K∆t k for n∆t ≤ T, (45)
with constant K = KTC.

Error analysis

We assume that the true PDE solutions E h e , E h i , H h e and H h i are sufficiently differentiable. Then we can Taylor expand the trapezoidal rule defect (36) at the symmetry point t n+1/2 to get

δ E h e (t n ) = ∑ j=2 -j 2 j ( j + 1)! (∆t) j E h( j+1) e , (46)
Inria where j = 2 means even values for j only and E h( j) e denotes the j-th derivative of E h e (t) at time t = t n+1/2 . We obtain similar expressions for δ

E h i (t n ), δ H h e (t n ) and δ H h i (t n ).
Note that these defects start with ∆t 2 and the third solution derivative of the true PDE solution (in the first section we have assumed that the true PDE solution is three times differentiable). We conclude that the trapezoidal rule defects are O ∆t 2 for ∆t ∼ h and h → 0. Then, from (36), we write the truncation error as follows

δ n = b n + O ∆t 2 , ( 47
)
where

b n =          b E e,n b E i,n b H e,n b H i,n          =             ∆t 4 S e d dt H h e (t n+1 ) -H h e (t n ) - ∆t 4 A ie d dt H h e (t n+1 ) -H h e (t n ) 0 - ∆t 4 A T ei S e H h e (t n+1 ) -H h e (t n ) + ∆t 4 A T ei d dt E h e (t n+1 ) -E h e (t n )             , ( 48
)
and O ∆t 2 contains the trapezoidal rule defects. Note that we have voluntarily omitted the spatial error parts contained in the σ E e , σ E i , σ H e , σ H i contributions from (36) because our interest lies in temporal convergence order. Further carrying these spatial error contributions only complicates the formulas and will not lead to different conclusions for the temporal errors.

Next we Taylor expand the components of b n at the symmetry point

t n+1/2 b E e,n = ∆t 4 S e ∑ j=1 1 2 j-1 j! (∆t) j H h( j+1) e , b E i,n = - ∆t 4
A ie ∑ j=1 1 2 j-1 j! (∆t) j H h( j+1) e , b H i,n = - ∆t 4 A T ei S e ∑ j=1 1 2 j-1 j! (∆t) j H h( j) e + ∆t 4 A T ei ∑ j=1 1 2 j-1 j! (∆t) j E h( j+1) e , (49)
where j = 1 means odd values for j only. For a fixed spatial dimension we find the expected secondorder ODE convergence, since S e , A ie and A T ei are bounded for fixed dimension. Indeed, with (49) we have b n = O ∆t 2 and because we assume stability we may consider R L inversely bounded, consequently 2 and we conclude that we have the second-order convergence for a fixed dimension. Now we observe, with [START_REF] Suzuki | Fractal decomposition of exponential operators with applications to many-body theories and Monte-Carlo simulations[END_REF], that under stable simultaneous space-time refinement, ∆t ∼ h and h → 0, we might lose one unit of ∆t in b E e,n and b E i,n (due to S e and A ie , respectively) and two units of ∆t in b H i,n

ρ n = R -1 L δ n = O ∆t
(due to the product A T ei S e in the first terms). Then b n = O (1), for ∆t ∼ h and h → 0 and with (45) we should expect a severe order reduction. However, as mentioned in [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF], this result is based on standard local error analysis and in the transition from local to global errors it can happen that the order reduction for local errors is (partly) canceled. Often this cancellation can be shown to exist through a transformation of the global error recurrence to one by which we may gain one unit of ∆t in the transformed local error.

A transformed global error recursion

The transformation used in [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] emanates from [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF], Lemma II.2.3. We write the latter for our one-step global error recursion (39) and stability assumption (42) Lemma 2.1. Suppose the local error ∆ρ n can be written as

∆tρ n = (I -R) ξ n + η n ,
with ξ n ≤ C∆t k , η n ≤ C∆t k+1 and ξ n+1ξ n ≤ C∆t k+1 for all n. Then there is a constant C > 0, depending on C, K and T , such that ε n ≤ C ∆t k for n∆t ≤ T .

The proof can be found in [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF]. First we assume for ∆t ∼ h, h → 0 that

A T ei S e H h(1) e = O 1 ∆t . ( 50
)
Consequently we get b n = O (∆t) (in this case b H i,n = O (∆t), see (49)). With the above-mentioned Lemma we can assume that if the local error ∆tρ n allows a decomposition

∆tρ n = (I -R) ξ n + η n (51)
such that ξ n = O ∆t 2 , η n = O ∆t 3 for ∆t ∼ h, h → 0, then we have the desired second-order convergence for ε n . So we need to verify (51), or equivalently,

∆tδ n = (R L -R R ) ξ n + R L η n , (52)
such that ξ n = O ∆t 2 , η n = O ∆t 3 for ∆t ∼ h, h → 0. Now we deal with the condition η n = O ∆t 3 . Recall that δ n = b n + O ∆t 2 , then (52) can be written as ∆t b n + O ∆t 2 = (R L -R R ) ξ n + R L η n . (53)
Furthermore, R L is inversely bounded, then by assigning the O ∆t 2 terms present in (53) multiplied by ∆t R -1 L for η n (i.e. η n = ∆t R -1 L O ∆t 2 ), we have η n = O ∆t 3 . Consequently, we deduce from Lemma 2.1 that we have the desired second-order convergence if a vector ξ n exists such that ξ n = O ∆t 2 for ∆t ∼ h, h → 0 and

∆tb n = (R L -R R ) ξ n , i.e. ( 54
)          D e 0 -S e A ei 0 D i A ie -S i S T e -A T ie 0 0 -A T ei S T i 0 0                   ξ E e,n ξ E i,n ξ H e,n ξ H i,n          =          b E e,n b E i,n b H e,n b H i,n          . (55)
Equivalently, we have second-order convergence if a vector

ξ n = ξ E n T , ξ H n T T exists such that ξ n = O ∆t 2 for ∆t ∼ h, h → 0 and Dξ E n -Sξ H n = b E n , S T ξ E n = b H n , (56)
where

ξ E n = ξ E e,n T , ξ E i,n T T , ξ H n = ξ H e,n T , ξ H i,n T T , b E n = b E e,n T , b E i,n T T , b H n = b H e,n T , b H i,n T T
.

Inria

Now we will check the existence of a such vector ξ n . At this stage of the derivation we must be careful because the matrix S is not necessarily a square matrix (in 2D this is not the case) and consequently S may no be invertible. More precisely, if we denote n do f the number of degrees of freedom, and if we consider the two-dimensional transversal magnetic (TM) model, then the size of the matrix S is n do f × 2n do f . That is why we now use the notion of pseudo inverse.

Definition 2.1. Let A ∈ R m×n , b ∈ R m , x ∈ R n and A + the Moore-Penrose pseudo inverse of A which is a generalization of the inverse and exists for any m × n matrix. If A has full rank, then

A + = A -1 (m = n) A + = A T AA T -1 (m < n) , A + = A T A -1 A T (m > n) , (57)
and the solution of Ax = b is x = A + b.

Assume the size of the matrix S is m × n with m ≤ n (the case m > n can be treated similarly). With the second equation of (56) and the above definition we derive

ξ E n = S T + b H n . ( 58
)
With [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] we have

S T + = SS T -1 S ∼ h, for h → 0, ( 59
)
and recall that b H e,n = 0 and b H i,n = O (∆t) (due to the initial assumption (50)) for ∆t ∼ h, h → 0. Then we conclude with (58) and ( 59) that for ∆t ∼ h, h → 0

ξ E n = O ∆t 2 . ( 60
)
From the first equation of (56) we get

ξ H n = -S + b E n -Dξ E n , (61)
and with [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] we get

S + = S T SS T -1 ∼ h, for h → 0. ( 62
)
Recalling that b E e,n = O (∆t), b E i,n = O (∆t) and ξ E n = O ∆t 2 , we conclude from (61) and (62) that for

∆t ∼ h, h → 0 if D = 0 (no conduction term -σ E) or D = O h k with k ≥ -1, then ξ H n = O ∆t 2 . ( 63
)
With ( 60) and (63) we conclude that

ξ n = ξ E n T , ξ H n T T
= O ∆t 2 for ∆t ∼ h, h → 0 and consequently through the Lemma 2.1 we have the second-order convergence uniformly in h under the assumption (50). Note that if we do not assume (50) (i.e. b n = O (1)) a similar proof, based on the same Lemma, only guarantees the first-order convergence. We can now state the following theorem (t) = O ∆t -1 for h → 0.

Theorem 2.1. Let f H (t), f E (t) ∈ C 2 [0, T ]
To sum up, we can guarantee at least the first-order convergence of method [START_REF] Yu | Local theory of extrapolation methods[END_REF]. As might be feared, component splitting can be detrimental to the temporal convergence order (order reduction). We have also put forward a sufficient condition (50) on the the true solution of the PDE problem for second-order convergence. However it would have been better if this sufficient condition could be controlled through the source term, because in general the true solution is of course not (a priori) known.

Numerical results

In this section we will conduct numerical tests in two dimensions by considering a wave propagation problem for which an analytical solution is available (see the subsection 3.1). The subsection 3.2 is devoted to the numerical convergence analysis of the locally implicit method [START_REF] Yu | Local theory of extrapolation methods[END_REF]. First in the ODE sense (i.e. on a fixed mesh) to illustrate that the method retains indeed its second-order convergence in time. Then in the PDE sense (i.e. for a simultaneous stable space-time grid refinement) to illustrate the reduction by one in the temporal convergence order, proved theoretically in the subsection 2.4. In the latter case we will also considered the second-order method [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF] to ensure that the reduction order is due to component splitting and not to the spatial discretization itself. Finally in the subsection 3.3 we will present an overview of the influence of the reduction order on the efficiency of the locally implicit method, when high order approximation polynomials are used within the DG method. For that purpose we will also considered the methods ( 18) -( 11).

Simulation setting

We solve the two-dimensional (2D) Transverse Magnetic (TM) model for the components E z (x, y,t), H x (x, y,t) and H y (x, y,t)

                 µ ∂ H x ∂t = - ∂ E z ∂ y , µ ∂ H y ∂t = ∂ E z ∂ x , ε ∂ E z ∂t = ∂ H y ∂ x - ∂ H x ∂ y -σ E z -J z E . (64)
In the following, we set ε = µ = 1 and σ = 0. Equations (64) are space discretized using a DG method formulated on quadrangular and triangular meshes, see Figure 1. In the preliminary implementation of this DG method, the approximation of the electromagnetic field components within a quadrangle c i or a triangle τ i relies on a nodal Q l or P l interpolation method, respectively. The a priori convergence analysis for this DGTD method based on a centered numerical flux and formulated on simplicial meshes shows that the convergence rate is O h l for a l-th interpolation order [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]. A quadrangle c i is characterized by the discretization parameters (∆x i , ∆y i ) and a triangle τ i by its height h i . The critical step size, denoted ∆t c , used in the numerical tests is given by

∆t c =              CFL 1 ∆x 2 k + 1 ∆y 2 k , on quadrangular meshes, CFL × h min k , on triangular meshes, (65)
where the parameters (∆x k , ∆y k ) characterize the smallest quadrangle c k treated explicitly and h min k the smallest height of a triangle τ k inside the region treated explicitly. The value of the CFL number corresponds to the numerical stability, i.e. the limit beyond which we observe a growth of the discrete energy. Finally for the component splitting we choose to treat implicitly the elements inside the red regions (see Figure 1).

We consider the propagation of an eigenmode in a unitary perfectly electrically conducting (PEC) Inria cavity. In this problem there is no source term i.e. J z E = 0 in (64) and the exact solution is given by

             H x (x, y,t) = - kπ ω sin (lπx) cos (kπy) sin (ωt) , H y (x, y,t) = lπ ω cos (lπx) sin (kπy) sin (ωt) , E z (x, y,t) = sin (lπx) sin (kπy) cos (ωt) , (66)
where the resonance frequencies are given by For numerical tests we put k = l = 1 and we initialize the electromagnetic field with the exact analyt-

ω = π k 2 + l 2 . ( 67
)
ical solution at t = 0 i.e.          H x (x, y,t = 0) = 0, H y (x, y,t = 0) = 0, E z (x, y,t = 0) = sin (πx) sin (πy) . (68)
For the boundary conditions, we consider a PEC cavity such that the tangential component of the electric field vanishes on the boundaries

n × E z = 0 on ∂ Ω, (69)
where the domain Ω = [0, 1] 2 (the unitary PEC cavity) and n denotes the unit outward normal to ∂ Ω.

Numerical convergence analysis

Throughout this subsection the total simulation time T is set to T = 3.34 × 10 -9 s, and we use the DGTD-Q 2 and DGTD-P 2 methods so that the spatial error is not detrimental to the temporal convergence orders.

For our integration method and the chosen spatial DG scheme the numerical CFL number, used in (65), is equal to 0.24 and 0.20 for DGTD-Q 2 and DGTD-P 2 methods, respectively.

Convergence in the ODE sense

To estimate the temporal convergence in the ODE sense we fix the discretization in space and we measure the electromagnetic field for different time steps (denotes by W ∆t ) at the final time T . We plot the sup-norm of W ∆t/n (T ) -W ∆t/2n (T ) as a function of ∆t/n, in logarithmic scale. The use of the logarithmic scale allows to visualize the convergence rates as the slopes of the curves with •-marker and -marker for uniform and non-uniform meshes, respectively. The dashed line has slope two for secondorder convergence. The results given in Figure 2 and Table 1 clearly confirm the expected behavior i.e. the second-order convergence of the method [START_REF] Yu | Local theory of extrapolation methods[END_REF] for ∆t → 0 on a fixed space grid. 18) based on the DGTD-Q 2 method.

Convergence in the PDE sense

First we conducted numerical investigations to check whether the condition (50) is satisfied. The following results have been obtained for ∆t ∼ h, h → 0

A T ei S e H h (1) e
=        O 1 ∆t 1.3 , for uniform meshes, O 1
∆t 1.2 , for non-uniform meshes. (70)
Consequently the condition (50) is not satisfied and we expect to observe a reduction PDE order by one. The possible loss in accuracy could originate from the spatial discretization itself, regardless of component splitting. To eliminate this latter possibility we also repeat the numerical tests using the fully explicit method [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF]. Furthermore in [START_REF] Buffa | Discontinuous Galerkin approximation of the Maxwell eigenproblem[END_REF] the authors have proven that DG methods on conformal tetrahedral/triangular meshes are spurious-free when the approximation spaces are made of elementwise polynomials of degree l in each variable or local Nédélec elements of the first type of degree l [START_REF] Nédélec | Mixed finite elements in R 3[END_REF]. They also remarked that this remains true for the local Nédélec elements of the first type of degree l on conformal hexahedral/quadrilateral meshes but not for elementwise polynomials of degree l in each variable or local Nédélec elements of the second type of degree l [START_REF] Nédélec | A new dfamily of mixed finite elements in R 3[END_REF]. In the latter case spurious modes can appear. Similarly the methods ( 18) -( 11) with a DG-Q 2 spatial discretization on quadrangular meshes can produce spurious modes. Therefore, we have also conducted numerical tests on triangular meshes (see Figure 1) with a DG-P 2 spatial approximation to ensure that the reduction of convergence is not due to the emergence of spurious modes, but indeed to component splitting.

To estimate the temporal convergence in the PDE sense we measure the maximal L 2 -norm of the error for different meshes of increased resolution. We plot this error as a function of the square root of the number of degrees of freedom (DOF), in logarithmic scale. The obtained results, given in Figure 3 and Table 2, clearly confirm the theoretical behavior i.e. the first-order convergence of the method ( 18) and the second-order of [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF], for ∆t ∼ h, h → 0; except for [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF] on non-uniform quadrangular meshes where we observe a super-convergence phenomenon (about 2.8 instead of 2.0).

High order approximation within the DG method

An attractive feature of the DG approach is the ability to easily increase the spatial convergence order. The reduction by one of the PDE order for the locally implicit method [START_REF] Yu | Local theory of extrapolation methods[END_REF] raises the question of the efficiency of the method when high order approximation polynomials are used within the DG method. We have conducted a numerical investigation for the propagation of an eigenmode in a unitary PEC cavity (66). The total simulation time T is set to T = 5 × 10 -8 s which corresponds to a propagation of the initial wave over 10 wavelengths. Note that the simulation time was chosen large enough in order that if we observe the evolution of the L 2 -norm error during time, the oscillations and the growth of the latter become very small from a time less than the final time. We plot the maximal L 2 -norm error as a function of the mesh spacing h and the total number of DOF for the DGTD-P We observe in Figure 4 that for a given error the locally implicit method needs a finer grid or significantly more DOF compared to the fully explicit case. In other words the method (18) needs more points per wavelength. We have also indicated in Tables 4 -5 the rates of decrease of the error relatively to the polynomial degree for both methods, for a given error tolerance (L 2 -norm error = 10 -3 ) or a given total number of DOF (= 10000). As expected these rates are lower for (18) than for [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF].

Because of its first order temporal convergence, an high order spatial discretization is less advantageous for method [START_REF] Yu | Local theory of extrapolation methods[END_REF] than for the fully explicit method [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF] which retains its second-order PDE convergence. Nevertheless increasing the polynomial order remains relatively useful for improving the accuracy as shown in Figure 4 and Tables 45. Furthermore, as previously mentioned in the introduction, another attractive feature of th DG approach is to handle geometrical details by using locally refined space grids. Consequently we can expect that the higher cost of the locally implicit method compared to the fully explicit method (in the sense of number of ppw for a given error) can be partly compensated because the unduly large step size restriction of ( 11) is overcome when a local refinement is necessary for a practical time dependent problem.

r (P k , P l ) |errorL 2 (P k ) -errorL 2 (P l ) | / |errorL 2 (P k ) | Method (18) Method (11)
r (P 2 , P 3 ) 0.14 0.34 r (P 2 , P 4 ) 0.28 0.52 r (P 3 , P 4 ) 0.17 0.28 In this paper, we have shown that component splitting can be detrimental to the temporal convergence uniformaly in h of a locally implicit DG method. We have presented a sufficient condition (see (50)) on the true PDE solution to recover the second-order convergence. Thus in the general case we can only guarantee first-order convergence. Numerical investigations have been conducted to illustrate this last theoretical result in 2D by considering a wave propagation problem for which an analytical solution is available. Note that the aim of these tests was not to illustrate the computational efficiency of the locally implicit DG method (for that purpose see [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF]) but the temporal convergence order of the method (18) in the ODE and PDE senses (i.e. for ∆t → 0 on a fixed space grid and ∆t ∼ h, h → 0, respectively).

In this paper we have mentioned another component splitting method from [START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] which cover the common spatial discretizations like finite difference and various finite element discretizations. We plan to compare the latter with the implicit-explicit integration method from [START_REF] Dolean | Locally implicit discontinuous Galerkin method for time domain electromagnetics[END_REF][START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF] for a DG discretization on unstructured meshes. The objective will be to illustrate the practical virtue of both component splitting methods using locally refined space grids, and to highlight the most efficient one which then can be advocated for future use.

In Section 3.3 we have seen that the reduction by one in the PDE convergence order for the method [START_REF] Yu | Local theory of extrapolation methods[END_REF] does not allow to fully exploit the gain in accuracy and efficiency that one can expect when the interpolation order is increased. Thus a future objective is to examine higher order implicit-explicit methods. Several strategies based on one of the schemes from [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF][START_REF] Verwer | Component splitting for semi-discrete Maxwell equations[END_REF] can be envisaged. One possibility is to exploit composition methods [START_REF] Hairer | Geometric Numerical Integration[END_REF], known to be accurate with well-designed composition coefficients [START_REF] Mclachlan | On the numerical integration of ordinary differential equations by symmetric composition methods[END_REF][START_REF] Suzuki | Fractal decomposition of exponential operators with applications to many-body theories and Monte-Carlo simulations[END_REF][START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF] to minimize truncation errors. Note that these composition methods with orders beyond two are restricted to problems with small (non-stiff) dissipative terms [START_REF] Botchev | Numerical Integration of Damped Maxwell Equations[END_REF]. Furthermore in the presence of source functions, the convergence order may be lower than the chosen composition order [START_REF] Verwer | Composition methods, Maxwell's and source term[END_REF]. In the latter reference the author used the scheme [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF] as the basis method in composition and a source function perturbation to solve Maxwell's equations. With this perturbed scheme the author obtained at least a third-order scheme, and a fourth-order one with additional sufficient conditions. Another possibility is to rely on local or global Richardson extrapolations [START_REF] Botchev | Numerical Integration of Damped Maxwell Equations[END_REF][START_REF] Faragó | Richardson-extrapolated sequential splitting and its application[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II -Stiff and Differential-Algebraic problems[END_REF][START_REF] Yu | Local theory of extrapolation methods[END_REF][START_REF] Verwer | Global extrapolation of a first order splitting method[END_REF] which are easy to implement and straightforward to parallelize. These high order extensions are well-known in the numerical solution of ODEs and have been already considered with [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF] as the basis methods to solve damped Maxwell equations [START_REF] Botchev | Numerical Integration of Damped Maxwell Equations[END_REF]. With Dirichlet boundary conditions the authors of [START_REF] Botchev | Numerical Integration of Damped Maxwell Equations[END_REF] observed that the global approach does not suffer from order reduction; this is not the case of composition methods and local extrapolations. However when the reduction order is inevitable, they advocate the local approach that allows to eliminate error terms instantaneously. The fact that these different methods, with an implicit-explicit approach, retain their ODE convergence order for simultaneous space-time grid refinement towards the true PDE solution is doubtful. Moreover, due to their accuracy, a comparison with high order local time stepping approaches, as in [START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF][START_REF] Taube | A high order discontinuous Galerkin method with local time stepping for the Maxwell equations[END_REF], will certainly be very interesting.

Figure 1 :

 1 Figure 1: Examples of quadrangular/triangular (uniform/non-uniform) meshes used in numerical tests.

  . / unif Loc. imp. / non-unif Slope two

Figure 2 :

 2 Figure 2: Temporal convergence (in the ODE sense) of the method (18) based on DGTD-Q 2 method.

Figure 4 :

 4 Figure 4: Maximal L 2 -norm error as function of the mesh spacing h and the total number of DOF (leftright) for methods (18) -(11) (top -bottom) based on DGTD-P 1,2,3,4 spatial discretization.

Contents 1 statement 4 2 2 . 1 2 . 2 2 . 3 2 . 4 21 4

 142122232421 Problem The implicit-explicit DGTD method 5 Component splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 The implicit-explicit time integration method . . . . . . . . . . . . . . . . . . . . . . . 7 Matrix behavior for h → 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Temporal convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.1 Elimination of intermediates values . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.2 The perturbed method and defects for the PDE solution . . . . . . . . . . . . . . 9 2.4.3 The error scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.4 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.5 A transformed global error recursion . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Numerical results 16 3.1 Simulation setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Numerical convergence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1 Convergence in the ODE sense . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.2 Convergence in the PDE sense . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.3 High order approximation within the DG method . . . . . . . . . . . . . . . . . . . . . Conclusion and future work 22 Inria RESEARCH CENTRE SOPHIA ANTIPOLIS -MÉDITERRANÉE 2004 route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

  . E h e (t n ) thus represents the vector that is approximated by E n

e . Assume the same notation for E i , H e and H i . Substituting E h e (t), E h i (t), H h e (t) and H h i (t) into (17) reveals the spatial truncation errors

  and suppose a Lax-Richtmyer stable space-time grid refinement ∆t ∼ h, h → 0. On [0, T ] the approximations H n

	H h e (t), H h i (t), E h e (t) and E h i (t)	e , H n i , E n e and E n i of method (18) then converge to
	(i) at least at first order,
	(ii) at least at second order, if in addition A T ei S e H e h(1)

Table 1 :

 1 Temporal convergence (in the ODE sense), on uniform and non-uniform quadrangular meshes, of the method (

	Inria

Table 2 :

 2 Temporal convergence (in the PDE sense), on uniform and non-uniform meshes, of methods (18) and (11) based on DGTD-Q 2 (on left) and DGTD-P 2 (on right) methods. Temporal convergence (in the PDE sense), on uniform and non-uniform quadrangular/triangular meshes (top/bottom), of methods (18) and (11) based on DGTD-P 2 method.

	Max. error (L 2 -norm)
	Loc. imp. / unif
	Fully exp. / unif.
	Loc. imp. / non-unif.
	Fully exp. / non-unif.
	Slope one
	Slope two
	1/2

Table 3 :

 3 1,2,3,4 methods on uniform triangular meshes. The critical time step ∆t c is determined by (65) where the numerical CFL number is given in Table3. Numerical value of the CFL number in (65), for triangular meshes.

	Method	DGTD-P 4	DGTD-P 3	DGTD-P 2	DGTD-P 1
	Numerical CFL	0.09	0.13	0.20	0.30

Table 4 :

 4 Rate of decrease of the error relatively to the polynomial degree, for methods (18) -[START_REF] Grote | Explicit Local Time Stepping Methods for Maxwell's Equations[END_REF], for a fixed number of DOF to 10000 errorL 2 (P l ) = [1r (P k , P l )] × errorL 2 (P k ) . r (P k , P l ) |DOF (P k ) -DOF (P l ) | / |DOF (P k ) |

		Method (18)	Method (11)
	r (P 2 , P 3 )	0.11	0.34
	r (P 2 , P 4 )	0.33	0.52
	r (P 3 , P 4 )	0.26	0.28

Table 5 :

 5 Rate of decrease of the error relatively to the polynomial degree, for methods (18) -(11), for a fixed error tolerance to 10 -3 (DOF (P l ) = [1r (P k , P l )] × DOF (P k )).

RR n°7533

Acknowledgments

The author wishes to acknowledge the many and important contributions to this work by Prof. dr. Jan Verwer of Centrum Wiskunde en Informatica (CWI, 1090 GB Amsterdam, The Netherlands), who passed away on 16 February 2011 before completion of this paper. The author wishes also to thank the CWI for its hospitality during the month of December 2010, period during which much of the theory presented in this paper was conducted in collaboration with Jan Verwer. Stéphane Descombes and Stéphane Lanteri are acknowledged for valuable discussions and important assistance throughout this work.


Pierre Gabriel
email: gabriel@ann.jussieu.fr

Léon Matar
email: leonmatar@yahoo.fr

High-order WENO scheme for Polymerization-type equations

Keywords: Aggregation-fragmentation equations, polymerization process, size repartition, long-time asymptotic, mass conservation, WENO numerical scheme AMS Class. No. 35B40, 35F50, 35L65, 35M30, 35Q92, 35R09, 65-06, 65M06, 65R20

Polymerization of proteins is a biochemical process involved in different diseases. Mathematically, it is generally modeled by aggregation-fragmentation-type equations. In this paper we consider a general polymerization model and propose a high-order numerical scheme to investigate the behavior of the solution. An important property of the equation is the mass conservation. The fifth-order WENO scheme is built to preserve the total mass of proteins along time.

Introduction

The central mechanism of amyloid diseases is the polymerization of proteins : PrP in Prion diseases, APP in Alzheimer, Htt in Huntington. The abnormal form of these proteins is pathogenic and has the ability to polymerize into fibrils. In order to well understand this process, investigation of the size repartition of polymers is a crucial point. To this end, we discuss in this paper the mathematical modeling of these polymerization processes and we propose numerical methods to investigate the mathematical features of the models.

Mathematical models are already widely used to study the polymerization mechanism of Prion diseases [START_REF] Calvez | Size distribution dependence of prion aggregates infectivity[END_REF][START_REF] Doumic | Scaling limit of a discrete prion dynamics model[END_REF][START_REF] Engler | Analysis of a model for the dynamics of prions ii[END_REF][START_REF] Greer | A mathematical analysis of the dynamics of prion proliferation[END_REF][START_REF] Greer | Effects of general incidence and polymer joining on nucleated polymerization in a model of prion proliferation[END_REF][START_REF] Laurençot | Well-posedness for a model of prion proliferation dynamics[END_REF][START_REF] Masel | Efficient inhibition of prion replication by prp-fc2 suggests that the prion is a prpsc oligomer[END_REF][START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF][START_REF] Lenuzza | Modélisation de la réplication des Prions: implication de la dépendance en taille des agrégats de PrP et de l'hétérogénéité des populations cellulaires[END_REF][START_REF] Pruss | Analysis of a model for the dynamics of prion[END_REF], Alzheimer [START_REF] Craft | A mathematical model of the impact of novel treatments on the a[beta]burden in the alzheimer's brain, csf and plasma[END_REF][START_REF] Lomakin | Kinetic theory of fibrillogenesis of amyloid-protein[END_REF][START_REF] Pallitto | A mathematical model of the kinetics of [beta]-amyloid fibril growth from the denatured state[END_REF] or Huntington [START_REF] Cajavec | Aggregation in huntington's disease: insights through modelling[END_REF]. Such models are also used for other biological polymerization processes [START_REF] Banasiak | On a coagulation and fragmentation equation with mass loss[END_REF][START_REF] Biben | Dynamics of Bio-Polymeric Brushes Growing from a Cellular Membrane: Tentative Modelling of the Actin Turnover within an Adhesion Unit; the Podosome[END_REF] and even for cell division [START_REF] Brikci | Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle[END_REF][START_REF] Doumic | Analysis of a population model structured by the cells molecular content[END_REF][START_REF] Perthame | Transport equations in biology[END_REF] or in neurosciences [START_REF] Pakdaman | Dynamics of a structured neuron population[END_REF].

Another field where we find aggregation-fragmentation equations is the physics of aggregates (aerosol and raindrop formation, smoke, sprays...). Among these models (see [START_REF] Laurençot | On coalescence equations and related models[END_REF] for a review), one can mention the Smoluchowsky coagulation equation [START_REF] Escobedo | Dust and self-similarity for the Smoluchowski coagulation equation[END_REF][START_REF] Filbet | Mass-conserving solutions and non-conservative approximation to the Smoluchowski coagulation equation[END_REF][START_REF] Filbet | Numerical simulation of the Smoluchowski coagulation equation[END_REF][START_REF] Laurençot | Liapunov functionals for Smoluchowski's coagulation equation and convergence to self-similarity[END_REF][START_REF] Mischler | Existence globale pour l'équation de Smoluchowski continue non homogène et comportement asymptotique des solutions[END_REF] with fragmentation [START_REF] Escobedo | Gelation and mass conservation in coagulation-fragmentation models[END_REF][START_REF] Escobedo | Gelation in coagulation and fragmentation models[END_REF][START_REF] Escobedo | On self-similarity and stationary problem for fragmentation and coagulation models[END_REF][START_REF] Laurençot | From the discrete to the continuous coagulation-fragmentation equations[END_REF][START_REF] Laurençot | Convergence to equilibrium for the continuous coagulationfragmentation equation[END_REF][START_REF] Laurençot | The continuous coagulation-fragmentation equations with diffusion[END_REF] and the Lifshitz-Slyosov system [START_REF] Carrillo | A numerical study on large-time asymptotics of the Lifshitz-Slyozov system[END_REF][START_REF] Collet | On solutions of the Lifshitz-Slyozov model[END_REF][START_REF] Collet | Some remarks on large-time asymptotic of the Lifshitz-Slyozov equations[END_REF][START_REF] Filbet | Numerical approximation of the Lifshitz-Slyozov-Wagner equation[END_REF][START_REF] Niethammer | The LSW model for domain coarsening: asymptotic behavior for conserved total mass[END_REF][START_REF] Niethammer | Well-posedness for measure transport in a family of nonlocal domain coarsening models[END_REF][START_REF] Niethammer | Global well-posedness for an inhomogeneous LSW-model in unbounded domains[END_REF][START_REF] Niethammer | On the convergence to the smooth self-similar solution in the LSW model[END_REF]. In [START_REF] Collet | Lifshitz-Slyozov equations: the model with encounters[END_REF][START_REF] Herrmann | Self-similar solutions with fat tails for a coagulation equation with nonlocal drift[END_REF] a Smoluchowsky coagulation term is added to the Lifshitz-Slyosov equation.

In this paper we are interested in a model including polymerization, coagulation and fragmentation phenomena. We consider a medium where there are monomers (normal proteins for instance) characterized by the concentration V (t) at time t and polymers (aggregates of abnormal proteins) of size x with the concentration u(t, x). The dynamics of the density function u(t, x) is driven by the system

                   d dt V (t) = - ∞ 0 T V (t), x u(t, x) dx, ∂ ∂t u(t, x) = - ∂ ∂x T V (t), x u(t, x) + Q(u)(t, x),
u(t, 0) = 0, u(0, x) = u 0 (x) ≥ 0 and V (0) = V 0 ≥ 0.

(

) 1
The monomers are attached by polymers of size x with the polymerization rate k on (x). Depolymerization occurs when monomers detach from polymers with a rate k off (x). Hence the transport term writes

T (V, x) = V k on (x) -k off (x). (2)
The two functions k on and k off are piecewise derivable but can be discontinuous. They are positive except that k on can vanish at zero. In this case, or more generally when T (V (t), 0) ≤ 0, the boundary condition on u(t, 0) is not necessary since the characteristic curves outgo from the domain. The choice of the boundary condition u(t, 0) = 0 is justified later. The coalescence of two polymers and the fragmentation of a polymer into two smaller ones are taken into account by the operator Q. More precisely, denoting by A x an aggregate of size x we have

A x + A y kc(x,y) ----→ A x+y coagulation A x+y k f (x,y) ----→ A x + A y fragmentation.
Thus the coagulation-fragmentation operator is

Q = Q c -Q f with Q c (u)(x) = 1 2 x 0 k c (y, x -y) u(y)u(x -y) dy -u(x) ∞ 0 k c (x, y) u(y) dy, (3)
Q f (u)(x) = 1 2 u(x) x 0 k f (y, x -y) dy - ∞ 0 k f (x, y) u(x + y) dy. (4)
The coalescence of two polymers of size x and y occurs with the symmetric rate

k c (x, y) = k c (y, x).
This rate is a nonnegative function as the fragmentation symmetric rate k f (x, y) = k f (y, x) with which a polymer of size x + y produces two fragments of size x and y.

There is a difference between V (t) and u(t, x = 0). In biochemical polymerization processes, small polymers are very unstable and thus do not exist. When they appear by detachment from a longer polymer, they are immediately degraded into monomers. Thus, the quantity of small polymers vanishes while the quantity of monomers is very high. To reflect this in the mathematical model, a quantity V (t) of monomers is introduced, which is different from the quantity of small polymers u(t, x = 0). The evolution of the first one is given by an ODE while the second one is forced to be equal to zero through the boundary condition u(t, 0) = 0. A consequence of this distinction is that starting from u 0 (x) = 0 and V 0 > 0 there is no evolution : the concentration of monomers is constant in time, V (t) = V 0 , and the concentration of polymers remains null, u(t, x) = 0. This is a very intuitive and natural behaviour which is important to preserve for biological applications. In the modeling, the distinction between V and u(x = 0) induces a separation of the polymerizationdepolymerization process from the coagulation-fragmentation. Indeed the aggregation of a monomer to a polymer can be seen as a coagulation but the resulting polymer has same size x than the initial one, since a monomer is very small compared to the typical size of a polymer. So a transport term is more accurate to model this phenomenon than an integral term (see [START_REF] Doumic | Scaling limit of a discrete prion dynamics model[END_REF]). There is also the fact that when a small polymer is degraded into monomers, it increases the quantity of monomers. In a discrete model, this term appears in the equation on V (see n 0 in [START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF]). In the continuous model [START_REF] Banasiak | On a coagulation and fragmentation equation with mass loss[END_REF] this term can be neglected since the quantity of monomers produced by degradation of small polymers is very small compared to the total quantity of monomers.

Mass Conservation

The mechanism of polymerization is nothing but a rearrangement of the proteins, there is no creation and no disparition. So the total quantity of proteins has to be constant in time and this is the case in model [START_REF] Banasiak | On a coagulation and fragmentation equation with mass loss[END_REF]. We define the total mass of the system as

P (t) = V (t) + ∞ 0 xu(t, x) dx, (5)
since a polymer of size x "contains x monomers". Integrating the equation on u(t, x) multiplied by x and adding the equation on V we obtain

∀t > 0, dP dt (t) = 0, (6)
so the total mass is conserved along time. This is a very important property that we want to keep in the numerical scheme and for this we rewrite equation (1) under a conservative form.

Conservative formulation

The classical discretization methods for transport equations are mass preserving. So the idea is to write the coagulation-fragmentation operator Q, which preserves the mass, under a conservative form in order to use a transport scheme. For this we follow the paper [START_REF] Filbet | An asymptotically stable scheme for diffusive coagulation-fragmentation models[END_REF] where such a transformation is made :

       xQ c (u)(x) = - ∂C(u) ∂x (x), xQ f (u)(x) = - ∂F(u) ∂x (x),
where the operator C(u) is given by

C(u)(x) := x 0 ∞ x-y yk c (y, z)u(y)u(z) dzdy, (7)
and F(u) is

F(u)(x) := x 0 ∞ x-y yk f (y, z)u(y + z) dzdy. (8)
Under this form, the mass conservation is clearer and the use of conservative schemes possible.

A useful consequence of the property (6) for the numerical scheme is that the ODE on V can be replaced by a mass conservation equation (see [START_REF] Herrmann | Self-similar solutions with fat tails for a coagulation equation with nonlocal drift[END_REF])

∀t > 0, V (t) = V 0 + ∞ 0 x(u 0 (x) -u(t, x)) dx. (9)
Numerically, this equation is much easier to compute than the ODE to be solved. Moreover [START_REF] Collet | On solutions of the Lifshitz-Slyozov model[END_REF] provides an explicit expression for V as a function of u. So we set

G(u)(x) := V 0 + ∞ 0 y u 0 (y) -u(y) dy k on (x) -k off (x) (10)
and we obtain a new equation equivalent to (1)

     x ∂ ∂t u(t, x) + ∂ G(u)xu ∂x (t, x) + ∂C(u) ∂x (t, x) - ∂F(u) ∂x (t, x) = G(u)u(t, x), u(t, 0) = 0, u(0, x) = u 0 (x). (11)
In this equation [START_REF] Craft | A mathematical model of the impact of novel treatments on the a[beta]burden in the alzheimer's brain, csf and plasma[END_REF], we have written the transport term as

x ∂ G(u)u ∂x (t, x) = ∂ G(u)xu ∂x (t, x) -G(u)u(t, x). ( 12
)
This formulation enhances the relation

d dt ∞ 0 xu(t, x) dx = ∞ 0 G(u)u(t, x) dx = - d dt V (t) (13)
and allows to preserve this property numerically when using conservative transport schemes.

Domain truncation

Numerically, equation ( 11) is solved on a truncated domain [0, R] so the integration bounds have to be changed in order to keep the mass preservation. For the coagulation term, we introduce as in [START_REF] Filbet | An asymptotically stable scheme for diffusive coagulation-fragmentation models[END_REF] C R (u)(x) := With this truncation, we have

C R (u)(0) = C R (u)(R) = F R (u)(0) = F R (u)(R) = 0.
So the total mass does neither increase nor decrease with respect to time if we consider the coagulation and fragmentation processes. If we look at the effects of this truncation on the original coagulation and fragmentation operators we have

Q R c (u)(x) := - 1 x ∂ x C R (u)(x) = 1 2 x 0 k c (y, x -y)u(u)u(x -y) dy -u(x) R-x 0 k c (x, y)u(y) dy and Q R f (u)(x) := - 1 x ∂ x F R (u)(x) = 1 2 u(x) x 0 k f (y, x -y) dy - R x k f (x, y -x)u(y) dy.
In the coagulation term, the truncation corresponds to the assumption that a polymer of size x cannot coagulate with a polymer of size greater than R -x. Concerning the fragmentation term, it is nothing but the assumption that polymers of size greater than R cannot split. Biologically they are the natural assumptions to avoid the loss of mass.

Concerning the transport term, the only way to avoid the loss of mass is to set

G R (u)(R, t) = 0. (14)
The meaning we give to this relation in the numerical scheme is exposed in Section 2.2. It is useless to do such a truncation for x = 0 since xu(t, x) vanishes when x = 0.

Finally we obtain a conservative truncated equation for x ∈ (0, R)

     x ∂ ∂t u R (t, x) + ∂ G R (u R )xu R ∂x (t, x) + ∂C R (u R ) ∂x (t, x) - ∂F R (u R ) ∂x (t, x) = G(u R )u R (t, x), u R (t, 0) = 0, u R (0, x) = u 0 (x). ( 15
)
When there is no transport term, convergence of the solution of Equation ( 15) to the solution of Equation (1) when R → ∞ is proved in [START_REF] Dubovskiȋ | Existence, uniqueness and mass conservation for the coagulation-fragmentation equation[END_REF][START_REF] Laurençot | On coalescence equations and related models[END_REF][START_REF] Laurençot | From the discrete to the continuous coagulation-fragmentation equations[END_REF][START_REF] Laurençot | The continuous coagulation-fragmentation equations with diffusion[END_REF][START_REF] Stewart | A global existence theorem for the general coagulationfragmentation equation with unbounded kernels[END_REF] under growth conditions on k c and k f .

A High Order WENO Scheme

In order to obtain a mass preserving scheme, we consider equation ( 11) as a transport equation and for high order accuracy we choose a fifth-order WENO (Weighted Essentially Non Oscillatory) reconstruction for the fluxes. This high order scheme is comonly used [START_REF] Devys | A model describing the growth and the size distribution of multiple metastatic tumors[END_REF][START_REF] Shen | A high order WENO scheme for a hierarchical size-structured population model[END_REF] since it is not more complicated to implement than a third order WENO one for instance.

Numerical fluxes

Before using the WENO reconstruction we have to know if the fluxes are positive or negative in order to appropriately upwind the scheme. Concerning the coagulation and the fragmentation terms, we consider a positive upwinding as suggested in [START_REF] Filbet | An asymptotically stable scheme for diffusive coagulation-fragmentation models[END_REF]. For the transport term ∂ x G(u)xu we have to make a flux splitting because G has no sign. A natural splitting here is to put the terms of G that are preceded by a plus sign in the positive part and the terms preceded by a minus sign in the negative part, namely

G = G + 1 + G - 1 where    G + 1 (u)(x) = V 0 + ∞ 0 yu 0 (y)dy k on (x), G - 1 (u)(x) = - ∞ 0 yu(y)dy k on (x) -k off (x). ( 16
)
An other decomposition is the polymerization-depolymerization one

   G + 0 (u)(x) = V 0 + ∞ 0 y u 0 (y) -u(y) dy k on (x), G - 0 (u)(x) = -k off (x). ( 17
)
The term

G + 0 is necessarily positive because V 0 + ∞ 0 y u 0 (y) -u(t, y) dy = V (t) ≥ 0.
With these two flux splittings, we built others by convex combination. For any λ ∈ [0, 1] we set

G λ = λG 1 + (1 -λ)G 0 which gives    G + λ (u)(x) = V 0 + ∞ 0 y u 0 (y) -u(y) dy + λ ∞ 0 yu(y)dy k on (x), G - λ (u)(x) = -λ ∞ 0 yu(y)dy k on (x) -k off (x). ( 18
)
We also consider the Lax-Friedrichs scheme which corresponds to

G + LF (u) = 1 2 G(u) + m , G - LF (u) = 1 2 G(u) -m , (19)
with m = max x≥0 |G(u)|. This term has to be computed at each time step because G(u) depends on time.

Finally, the WENO reconstruction is done with the fluxes

   H + (u) = G + (u)xu + C(u) -F (u), H -(u) = G -(u)xu, (20)
and the choice among the different flux splittings is discussed in Section 3.2.

WENO reconstruction

The point of view adopted here is the finite difference one, as recommanded in [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF], because it is better than the finite volume in terms of operation counts. We assume the spatial domain [0, R] is divided into N uniform cells and we denote x i = i∆x for 0 ≤ i ≤ N with ∆x = R N . We use the WENO formulation of Jiang and Peng [START_REF] Jiang | Weighted eno schemes for hamilton-jacobi equations[END_REF] which consists in applying WENO to approach the spacial derivative directly on the nodes of the grid. The spatial derivative

∂ x (H + (v)+H -(v)) is approximated at the point x i by 1 ∆x H + i+ 1 2 + H - i+ 1 2 -H + i-1 2 -H - i-1 2
where the fifth order accurate numerical flux H + i+ 1 2 is given by the WENO reconstruction. For each node x i we denote by H + i the numerical approximation of H + (v(x i )). The stencil choice for each flux is specified in Figure 1, and the fluxes H ± i± 1 2 are expressed as convex combination of the H ± k of the stencil. Let us detail how we proceed :

-for

H - i-1 2 we set W 1 = H - k , W 2 = H - k+1 , W 3 = H - k+2 , W 4 = H - k-1 , W 5 = H - k-2 ,
-for H -

i+ 1 2 we set W 1 = H - k+1 , W 2 = H - k+2 , W 3 = H - k+3 , W 4 = H - k , W 5 = H - k-1 , -for H + i-1 2 we set W 1 = H + k-3 , W 2 = H + k-2 , W 3 = H + k-1 , W 4 = H + k , W 5 = H + k+1 , -for H + i+ 1 2 we set W 1 = H + k-2 , W 2 = H + k-1 , W 3 = H + k , W 4 = H + k+1 , W 5 = H + k+2 . H - i-1 2 H - i+ 1 2 H + i-1 2 H + i+ 1 2 i-3 i-3 i-2 i-2 i-1 i-1 i i i+1 i+1 i+2 i+2 i+3 i+3 Figure 1: stencil choice
For the regularity coefficients we define for each previous flux

S 1 = 13 12 (W 1 -2W 2 + W 3 ) 2 + 1 4 (W 1 -4W 2 + 3W 3 ) 2 , S 2 = 13 12 (W 2 -2W 3 + W 4 ) 2 + 1 4 (W 2 -W 4 ) 2 , S 3 = 13 12 (W 3 -2W 4 + W 5 ) 2 + 1 4 (3W 3 -4W 4 + W 5 ) 2 .
Then we take the weights

w r = a r 3 j=1 a j , with a r = d r (ǫ + S r ) 2 , d 1 = 3 10 , d 2 = 6 10 , d 3 = 1 10 r = 1, 2, 3
where ǫ is introduced to prevent the denominator from vanishing. Finally we take the different flux parts given by

           H - i± 1 2 = w 1 W 3 3 - 7W 2 6 + 11W 1 6 + w 2 -W 2 6 + 5W 1 6 + W 4 3 + w 3 W 1 3 + 5W 4 6 - W 5 6 , H + i± 1 2 = w 1 W 1 3 - 7W 2 6 + 11W 3 6 + w 2 -W 2 6 + 5W 3 6 + W 4 3 + w 3 W 3 3 + 5W 4 6 - W 5 6 .
Concerning the boundaries x = 0 and x = R, we compute the fluxes using the WENO reconstruction with ghost points x -3 , x -2 , x -1 , and x N +1 , x N +2 , x N +3 . In the first three points we use that for all time t ≥ 0, xu(t, x)

x=0 = C R (u)(x = 0, t) = F R (u)(x = 0, t) = 0 to set H ± -3 = H ± -2 = H ± -1 = 0.
For the last three ones we use the truncation

G R (u)(x = R, t) = C R (x = R, t) = F R (x = R, t) = 0 to put H ± N +1 = H ± N +2 = H ± N +3 = 0.

Integration method

For the integral terms, we use a fifth order composite rule introduced in [START_REF] Shen | A high order WENO scheme for a hierarchical size-structured population model[END_REF]. If f k denotes an approximation of f (x k ), the method can be written as

j∆x i∆x f (x) dx ≃ ∆x j k=i ′ f k where j k=i ′ f k = 251 720 f i + 299 240 f i+1 + 211 240 f i+2 + 739 720 f i+3 + 739 720 f j-3 + 211 240 f j-2 + 299 240 f j-1 + 251 720 f j + j-4 k=i+4 f k if j -i > 6.
This method is based on polynomial interpolations of the function f. On the first interval, we integrate without using the boundary value f 0 = 0 because the solution can be discontinuous at x = 0. So we use the fifth accurate approximation

1 k=0 ′ f k = 55 24 f 1 - 59 24 f 2 + 37 24 f 3 - 9 24 f 4 .
Finally for the intervals at the boundaries we have

2 0 ′ f k = 8 3 f 1 - 5 3 f 2 + 4 3 f 3 - 1 3 f 4 , 3 0 ′ f k = 21 8 f 1 - 9 8 f 2 + 15 8 f 3 - 3 8 f 4 , 4 0 ′ f k = 21 8 f 1 - 7 6 f 2 + 29 12 f 3 + 1 6 f 4 - 1 24 f 5 , 5 0 ′ f k = 21 8 f 1 - 7 6 f 2 + 19 8 f 3 + 17 24 f 4 + 1 2 f 5 - 1 24 f 6 , 6 0 ′ f k = 21 8 f 1 - 7 6 f 2 + 19 8 f 3 + 2 3 f 4 + 25 24 f 5 + 1 2 f 6 - 1 24 f 7 , 7 0 ′ f k = 21 8 f 1 - 7 6 f 2 + 19 8 f 3 + 2 3 f 4 + f 5 + 25 24 f 6 + 1 2 f 7 - 1 24 f 8 , N N -1 ′ f k = 9 4 f N + 19 24 f N -1 - 5 24 f N -2 + 1 24 f N -3 , N N -2 ′ f k = 1 3 f N + 4 3 f N -1 + 1 3 f N -2 , N N -3 ′ f k = 1 3 f N + 31 24 f N -1 + 7 8 f N -2 + 13 24 f N -3 - 1 24 f N -4 , N N -4 ′ f k = 1 3 f N + 31 24 f N -1 + 5 6 f N -2 + 13 12 f N -3 + 1 2 f N -4 - 1 24 f N -5 , N N -5 ′ f k = 1 3 f N + 31 24 f N -1 + 5 6 f N -2 + 25 24 f N -3 + 25 24 f N -4 + 1 2 f N -5 - 1 24 f N -6 , N N -6 ′ f k = 1 3 f N + 31 24 f N -1 + 5 6 f N -2 + 25 24 f N -3 + f N -4 + 25 24 f N -5 + 1 2 f N -6 - 1 24 f N -7 .
We use this quadrature method to discretize the operators C R and F R with

C R i -F R i = (∆x) 2 i j=0 ′ N l=i+1 ′ x l k c j,l-j u j u l-j -k f j,l-j u l (21)
as suggested in [START_REF] Filbet | An asymptotically stable scheme for diffusive coagulation-fragmentation models[END_REF]. Grouping the two terms in an unique summation is lighter regarding to operation counts.

Time discretization

The time step is denoted by ∆t and changes along time because of the CFL stability condition that is time dependent. For the time discretization we choose a third order Runge-Kutta method. To approach the time evolution of an equation ∂ t u = L(u), we compute at time n∆t

u (1) = u n + ∆t L(u n ) and u (2) = 3 4 u n + 1 4 u (1) + 1 4 ∆t L(u (1) ),
where u n is an approximation of u(n∆t). Then the approximation of v at time (n + 1)∆t is given by 2) ).

u n+1 = 1 3 u n + 2 3 u (2) + 2 3 ∆t L(u (
This method is an explicit one, so to ensure the stability we compute the time step ∆t at each iteration thanks to the CFL condition ∆t

≤ min (G + C + F ) -1 (22)
where

G = 1 ∆x sup i G + i -G - i , C = sup i N j=1 ′ k c i,j u j and F = sup i 1 2 i-1 j=1 ′ k f j,i-j .
For instance the Lax-Friedrichs decomposition leads to G LF = m/∆x.

Since we combine a fifth order WENO reconstruction and a third order time discretization, we predict that our scheme is convergent of third order. To validate it numerically, we compute the solution for different discretization grids with regular parameters and initial data. Comparing these solutions at time T = 20 in the L ∞ space norm (see Table 1 for the results), we obtain the numerical order 2.95 which validates the prediction. Table 1: Error between different discrete solutions and the reference computed with ∆x = 5/320.

Numerical Simulations

Parameters

Numerical values for the polymerization and fragmentation rates can be found in the biological literature (see [START_REF] Masel | Efficient inhibition of prion replication by prp-fc2 suggests that the prion is a prpsc oligomer[END_REF][START_REF] Rubenstein | Dynamics of the nucleated polymerization model of prion replication[END_REF][START_REF] Knowles | An Analytical Solution to the Kinetics of Breakable Filament Assembly[END_REF] for instance). It is of importance to note that the models considered in these papers are discrete ones, so some computations are necessary to deduce adimensional numerical values for the continuous model [START_REF] Banasiak | On a coagulation and fragmentation equation with mass loss[END_REF]. Another point is that the parameters of these models do not depend on the size of polymers, so we can only obtain mean values for the size-dependent parameters.

We choose to use the values of the recent paper [START_REF] Knowles | An Analytical Solution to the Kinetics of Breakable Filament Assembly[END_REF] to do numerical simulations. The mean length of polymers for the initial distribution is estimated to be 1380. With the continuous model we reduce this value to 0.2 by considering an initial profile equal to a positive constant on [0, 0.4] and null for x > 0.4 (see the first plot of Figure 5). Thus we define a parameter ε = 0.2/1380 ≈ 1, 4 × 10 -4 which allows to go from a discrete model to a continuous one (see [START_REF] Doumic | Scaling limit of a discrete prion dynamics model[END_REF] for more details). The values we find are for instance 2.9 × 10 -2 µM -1 s -1 for the polymerization rate and 2.1 × 10 -9 s -1 for the fragmentation (where M represents the concentration in mol and s the time in second). The polymerization rate appears in a derivative term, so the value of the discrete model has to be multiplied by ε to obtain the continuous accurate value 4 × 10 -6 µM -1 s -1 . Conversely, the fragmentation rate which appears in an integral term has to be divided by ε and we find 1.5 × 10 -5 s -1 . Concerning the depolymerization and coagulation, they are neglected in the models of [START_REF] Masel | Efficient inhibition of prion replication by prp-fc2 suggests that the prion is a prpsc oligomer[END_REF][START_REF] Rubenstein | Dynamics of the nucleated polymerization model of prion replication[END_REF][START_REF] Knowles | An Analytical Solution to the Kinetics of Breakable Filament Assembly[END_REF]]. So we consider numerical values that seem to be reasonable compared to the previous ones.

In the present study, the parameters are assumed to be size dependent as suggested in [START_REF] Calvez | Prion dynamic with size dependency -strain phenomena[END_REF][START_REF] Calvez | Size distribution dependence of prion aggregates infectivity[END_REF] and their choice is now presented and motivated. Concerning the numerical coefficients, they are chosen in order to have mean values of the same order than the values previously obtained from [START_REF] Knowles | An Analytical Solution to the Kinetics of Breakable Filament Assembly[END_REF]. For the polymerization we assume that small polymers have a different behavior compared to the big ones. We consider a critical size x c = 0.5 such that polymers of size x < x c convert monomers with the rate k (1) on (x) = (4x + 0.2) × 10 -6 µM -1 s -1 , and for x > x c with a constant rate k (2) on (x) ≡ 4 × 10 -6 µM -1 s -1 .

For the fragmentation kernel, we use the classical assumption that the fragmentation probability depends only on the size x + y of the polymer and we set

k f (x, y) = 80(x + y) 10 + (x + y) × 10 -5 s -1 .
The depolymerization is assumed to be constant and of the same order as the fragmentation. We discuss the dependence on this parameter by considering different intensities

k off (x) ≡ η × 10 -6 s -1 with 2 ≤ η ≤ 8. ( 23
)
Concerning the coagulation kernel, we do not use a classical one. Even if there is no space in model ( 1), we use a kernel which reflect some space effects. Indeed we consider that small polymers are very mobile and that the big ones, plaques, are very attractant. So the coagulation occurs preferentially between big and small aggregates. The kernel we choose is of the form

k c (x, y) = 4|x -y| 3 2 1 + (x + y) × 10 -6 µM -1 s -1 .
This kernel satisfies the growth assumption that we can find in [START_REF] Dubovskiȋ | Existence, uniqueness and mass conservation for the coagulation-fragmentation equation[END_REF] to ensure the convergence of the solution when R → ∞ if we consider the only coagulation-fragmentation process. In [START_REF] Knowles | An Analytical Solution to the Kinetics of Breakable Filament Assembly[END_REF] we also find numerical values for the initial data V 0 = 98 µM and

R 0 xu 0 (x) dx = 0.21µM.
This last value and the fact that the initial distribution of polymers is assumed to be under the form

u 0 = cst × 1l [0,0.4] lead to u 0 (x) = 2.6 if 0 ≤ x ≤ 0.4 0 if x > 0.4. (24)
For the following simulations, the discretization is made on a domain [0, 5] with a number of nodes N = 200, so the mesh size is ∆x = 0.025.

Choice among the different flux splittings

First we deal with the CFL condition. Thanks to the triangular inequality, we obtain that G LF ≤ G 0 . Moreover, there is an explicit expression for G λ

G n λ = 1 ∆x sup i V 0 + ∆x N j=1 ′ x j u 0 j + (2λ -1)∆x N j=1 ′ x j u n j k on i + k off i which shows that G λ increases with λ. So if 0 ≤ λ < Λ ≤ 1 then at each time step we have G n LF ≤ G n λ ≤ G n Λ .
Notice also that, with the numerical values we have chosen, the quantity of polymers ∆x N j=1 ′ x j u n j increases with n. Indeed we can see in Figure 8 that the quantity of monomers V n ≃ V (n∆t) defined by the mass conservation

V n + ∆x N j=1 ′ x j u n j = V 0 + ∆x N j=1 ′ x j u 0 j decreases. The consequence on the CFL condition is that G n λ increases with n if λ > 1 2 , decreases if λ < 1
2 and is time independent when λ = 1 2 . Thus, regarding to the numerical computation, the fastest scheme is the Lax-Friedrichs one and then the computation time increases significantly with λ.

Let us now turn to the effects of the flux splitting on the size distribution. First we consider a depolymerization corresponding to η = 8 in [START_REF] Filbet | An asymptotically stable scheme for diffusive coagulation-fragmentation models[END_REF] and investigate the differences between the solutions associated to the decompositions G LF , G 0 and G 1 . We can see in Figure 3 that the solutions for G 0 and G 1 are close together for small times and then the behavior of G 0 becomes closer to the Lax-Friedrichs one. The less oscillating scheme for t = 6h is the the Lax-Friedrichs one, but it is also the most oscillating at time t = 12h. Finally the solutions associated to the three flux decompositions are quite similar and they all present oscillations at some times, so we do not find with this simulation any reason to discard one of them. x x x If we change the depolymerization rate by considering k off = 6 × 10 -6 s -1 , we remark that the Lax-Friedrichs scheme is unstable (see Figure 4) while G 0 is stable. If we continue to decrease η, we find with k off = 2 × 10 -6 s -1 that the G 0 -scheme becomes unstable while G 0.2 is stable. Thus the Lax-Friedrichs scheme and the G λ -schemes with λ small has to be avoided to ensure stability when small depolymerization values are considered. x x Knowing that, we compare different stable schemes, namely G λ with 0.2 ≤ λ ≤ 1. We can see in Figure 5 that for large times (t = 20h), the three flux decompositions provide a good behavior where there are strong variations of the solution. These locations are x = 0 because of the boundary condition which enforces u(t, 0) to vanish, and x = 0.5 because the transport term k on is discontinuous at x = 0.5. If we look at smaller times (t = 6h for instance) we can see that the larger λ is, the less oscillating the curves are. But, as we already remarked, the quantity G λ is higher for λ close to 1 and it increases with time when λ > 1 2 . Thus it is penalizing for the computation time to use high values of λ. A good compromise could be to choose λ = 1 2 since G 1 2

t = 6 t = 6 t = 6 u(x, 6) u(x, 6) u(x, 6) x x x G LF G LF G LF G 0 G 0 G 1 40
G LF G LF G LF G 0 G 0 G 1
G 0 G 0 G 0.2
does not depend on time. The other solution is to adapt the λ when we change the parameters.

t = 6 t = 6 t = 6 u(x, 6) u(x, 6) u(x, 6) x x x G 0.2 G 0.2 G 0.2 G 0.5 G 0.5 G 1 0
x x x G 0.2 G 0.2 G 0.2 G 0.5 G 0.5 G 1

Interpretation of the numerical results

We have considered that the mean size of the polymers at the initial time t = 0 was 1380. This size can be multiplied by 5 along the polymerization process (see Figure 6 keeping in mind that the mean size is represented by 0.2 in this continuous model). So if we want to solve the discrete model, we have to consider a system of dimension close to 5000, and the computations are very heavy. If we limit this value to 200 keeping the discrete model, then we lose a lot of precision. It is the same for the continuous model if it is discretized with a first order scheme. That is why we use a high order discretization, and we can see the difference in Figures 6 and7 : the high order scheme is able to capture strong variations of amplitude while the first order flattens them. We also remark that the size repartition converges to a bimodal distribution, as observed by [START_REF] Silveira | The most infectious prion protein particles[END_REF]. This asymptotic profile is independent of the initial data (see the time t = 20 in Figures 6 and7) and can be seen as an eigenvector of the operator Q -∂ x T (see [START_REF] Michel | Existence of a solution to the cell division eigenproblem[END_REF][START_REF] Jauffret | Eigenelements of a general aggregation-fragmentation model[END_REF]). The evolution of the quantity of monomers V (t) is plotted in Figure 8 for different values of the depolymerization rate k off . This quantity decreases since the monomers aggregate to polymers. Thus the mass of polymers increases and the speed of this evolution is similar to those observed by [START_REF] Knowles | An Analytical Solution to the Kinetics of Breakable Filament Assembly[END_REF]. Concerning the dependence on k off , the difference between the three curves is more significant when the time increases. For small times, when V (t) is close to V 0 = 98µM, the depolymerization can be neglected since k off is small compared to the product V (t)k on (x). Conversely, the equilibrium is reached when d dt V (t) = 0, so when k off ≃ V k on (see Equation ( 1)). That is why variations of η influence essentially the ratio between the quantity of monomers and the mass of polymers at the equilibrium as we can see in Figure 8.

Conclusion and future work

We have written a high order conservative scheme for a polymerization-type equation. The choice of the flux splitting for the transport term has been discussed but the accurate decomposition remains unclear. It seems that unstabilities can be avoided by adapting the value of λ but the oscillations remain present for any choice of the flux decomposition, even for a regular initial size distribution. A possible explanation for these phenomena can be that the integration method is not "positive" for the intervals at the boundaries. These points remain to be investigated for a better understanding and improvement.

As we have remarked in Section 3.3, the size distribution converges toward an equilibrium which corresponds to an eigenvector. The high-order WENO scheme presented in this paper could be used to numerically compute such eigenvectors. Another application of the code is to solve inverse problems (see [START_REF] Doumic | Numerical solution of an inverse problem in sizestructured population dynamics[END_REF][START_REF] Perthame | On the inverse problem for a size-structured population model[END_REF]) in order to determine the size dependence of the different parameters.

  y, z -y)u(y)u(z -y) dzdy, yk f (y, z)u(y + z) dzdy = x 0 R x yk f (y, z -y)u(z) dzdy.

∆x 5 /

 5 [START_REF] Masel | Quantifying the kinetic parameters of prion replication[END_REF] 5/80 5/160 error 81.84 12.44 1.37

Figure 2 :

 2 Figure 2: Profiles of the coagulation and fragmentation kernels.

Figure 3 :

 3 Figure 3: Comparison between the flux splittings G LF , G 0 and G 1 for η = 8 × 10 -6 s -1 .

Figure 4 :

 4 Figure 4: Unstability of some schemes when η decreases. Left: G LF becomes unstable for η = 6. Right: G 0 becomes unstable for η = 2.

Figure 5 :

 5 Figure 5: Comparison of the behavior of the solution for different λ with η = 5.

Figure 6 :

 6 Figure 6: Comparison between the WENO scheme and a first order scheme for the initial size distribution[START_REF] Filbet | Numerical approximation of the Lifshitz-Slyozov-Wagner equation[END_REF], with the depolymerization value k off = 8 × 10 -6 s -1 and the flux splitting parameter λ = 0.5

Figure 7 :

 7 Figure 7: Comparison between the WENO scheme and a first order scheme for a regular initial size distribution, with the depolymerization value k off = 8 × 10 -6 s -1 and the flux splitting parameter λ = 0.5

8 Figure 8 :

 88 Figure 8: Evolution of the quantity of monomers for different depolymerization rates, with the scheme G 0.2 .

Acknowledgement

The authors are thankfull to Vincent Calvez, Marie Doumic Jauffret, Frédéric Lagoutière and Natacha Lenuzza for their help and support during the Cemracs'09 research session. This work has been done with the financial supports of -the ANR contract TOPPAZ, allocation grant No. 4243, http://www-roc.inria.fr/bang/TOPPAZ/ -the CEA-Institute of Emerging Diseases and Innovative Therapies, Route du Panorama, Bat.60, F-92265 Fontenay-aux-Roses.