Simon Rohou

Abderahmane Bedouhene

Gilles Chabert

Alexandre Goldsztejn

Luc Jaulin

Bertrand Neveu

Victor Reyes

Gilles Trombettoni

Towards a Generic Interval Solver for Differential-Algebraic CSP

come   L'archive ouverte pluridisciplinaire

Introduction

Differential-Algebraic CSP (DACSP) systems include real variables, functional variables (also called trajectories or dynamical variables) describing the dynamics of the system, and differential and/or "static" numerical constraints among those variables. Because they are at the core of so many applications, such as biological systems, mechanism dynamics, astronomy, robotics, control, a lot of work has been dedicated to solving specific subclasses of the DACSP. Most of the approaches follow a probabilistic approach and are limited to linear constraints with Gaussian errors [START_REF] Kalman | Contributions to the Theory of Optimal Control[END_REF][START_REF] Thrun | Probabilistic Robotics[END_REF].

There are a number of advantages to using interval methods for handling dynamical and/or static systems. They can manage nonlinear constraints and approximate the solutions rigorously, whatever the uncertainties on the parameters (e.g., uncertainties related to measurements or to inaccurate physical models). Bounded intervals are used to characterize these uncertainties, as well as the errors due to operations over floating-point numbers.

With Jaulin et al. works [START_REF] Jaulin | Nonlinear Bounded-error State Estimation of Continuous-Time Systems[END_REF][START_REF] Jaulin | Range-Only SLAM with Indistinguishable Landmarks: A Constraint Programming Approach[END_REF], a significant step has been made towards a declarative constraint programming approach for dynamical systems. Contrary to dominant constraint approaches dealing with differential equations such as [START_REF] Deville | Consistency Techniques in Ordinary Differential Equations[END_REF][START_REF] Revol | Taylor models and floating-point arithmetic: proof that arithmetic operations are validated in COSY[END_REF][START_REF] Tucker | A Rigorous ODE Solver and Smale's 14th Problem[END_REF], the trajectories are viewed as variables of the CSP. They consider trajectories as variables, differential equations as constraints and tubes as domains. The solution and the domain are given by a tube representing a set of possible trajectories (see Fig. 1). This increases the level of abstraction and simplifies the formalization of the problem. For estimating a trajectory, a set of contractors, similar to propagators in CP, are applied to filter (contract) the bounds of the domain/tube. A benefit of the contractor programming approach [START_REF] Chabert | Contractor Programming[END_REF] is the variety of dynamical systems that can be handled. More recently, the Tubex library [START_REF] Rohou | The Tubex library -Constraint-programming for robotics[END_REF] has provided data structures for representing tubes and a catalogue of contractors, similarly to the catalogue of propagators available in CP solvers [START_REF] Prud'homme | Choco documentation[END_REF]. The user defines a sequence of contractors for modeling his problem. The set of contractors is applied iteratively until a quasi fixpoint in terms of contraction is reached. For instance, the resulting framework has been recently applied to actual data for autonomous robot localization. [START_REF] Rohou | Set-membership State Estimation by Solving Data Association[END_REF] describes the problem of localizing an underwater robot. Its evolution is depicted by an Ordinary Differential Equation (ODE) and bounded measurements, but the initial condition (position) of the system is unknown. During the mission, its sonar detects indistinguishable rocks on the seabed, that all look alike but are known to belong to a discrete point map embedded beforehand. The map of rocks positions is also modeled as a constraint. The constraint propagation approach is able to merge all data coming from the robot evolution and rocks observations. The identity of the rocks is finally associated to items in the map, and the trajectory of the robot is accurately estimated. This application highlights how a complex problem involving discrete, continuous and differential constraints can be solved easily following a CP/contractor approach.

Contractor programming is relevant when the problem is defined by heterogeneous constraints, provided they are redundant and numerous enough to enable the contraction phase alone to solve the problem.

Contributions.

In this paper, we introduce a generic solver using the Tubex library that can handle a DACSP instance made of differential constraints and/or static continuous constraints. The numerical constraints relate real variables that either represent states of the trajectories at given instants or are independent from the dynamics (e.g., the rocks positions in the previous example). Compared to the Tubex approach alone, the solver is endowed with an operator that can perform a choice point by splitting (bisecting) a tube in two at a chosen time. Thus, a tree search can accurately estimate distinct trajectories (for problems with several solutions) and can better handle several hard problems. To our knowledge, no other tool for solving dynamical systems performs a tree search.

Our solver can manage most of the contractors available in the literature. On the one hand, several contractors coming from CP, e.g., 3B [START_REF] Lhomme | Consistency Techniques for Numeric CSPs[END_REF] and CID [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF], are included naturally to improve the pruning of domains of functional variables, i.e., tubes. On the other hand, our generic framework enables to wrap existing solvers dedicated to specific types of (differential) constraints into contractors. These contractors efficiently reduce domains by taking into account space and time dependencies. In particular, we show in the experiments the benefits of using the contractor CtcVnode built upon the state-of-the-art VNODE guaranteed integration IVP solver [START_REF] Nedialkov | Validated Solutions of Initial Value Problem for Ordinary Differential Equations[END_REF][START_REF] Nedialkov | A Validated Solver for Initial Value Problems in Ordinary Differential Equations[END_REF].

This generic CP framework allows separating the problem description, which includes here all combinations of differential and algebraic problems, from its resolution. We also wanted to solve difficult standard problems coming from numerical analysis, e.g., BVP for integro-differential equations, that are out of the scope of previous CP/ODE frameworks and remain today difficult to solve in the numerical analysis framework.

Related Work. VNODE [START_REF] Nedialkov | A Validated Solver for Initial Value Problems in Ordinary Differential Equations[END_REF], CAPD [START_REF] Kapela | CAPD -a rigorous toolbox for Computer Assisted Proofs in Dynamics[END_REF], COSY [START_REF] Revol | Taylor models and floating-point arithmetic: proof that arithmetic operations are validated in COSY[END_REF] and DynIbex [START_REF] Chapoutot | Dynibex[END_REF] are state of the art interval analysis solvers dedicated to IVPs. 6 They are fully relevant for determining the guaranteed solution of a system at a final time, that has crucial applications such as the position determination of celestial bodies in astronomy [START_REF] Wilczak | Heteroclinic connections between periodic orbits in planar restricted circular three-body problem-a computer assisted proof[END_REF] or to characterize chaotic attractors [START_REF] Tucker | A Rigorous ODE Solver and Smale's 14th Problem[END_REF]. They use different algorithms to reliably simulate the initial information over time. In particular, the VNODE tool used in our solver combines a high-order interval Taylor form to integrate the state from an instant to a next one, and a step limiting the wrapping effect implied by interval calculation: it encloses the solution at the discrete times by an envelope sharper than a box, such as rotated boxes, zonotopes or polygons [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF].

A BVP is generally defined by an ODE, but the trajectory is not entirely determined at the initial or final times, which prevents a solving algorithm from propagating (integrating) the information from the known state to the rest of the trajectory. Instead, static constraints are defined on specific states (at specified instants). To deal with BVPs, the shooting method [START_REF] Goodman | The Numerical Solution of Two-Point Boundary Value Problems[END_REF] is a sampling-andoptimization algorithm that runs several integration processes from the initial state while trying to minimize the distance to a solution satisfying the ODE. This makes the approach incomplete and unable to determine several solutions, if any. Instead, our rigorous and deterministic solver can explore the whole search space and isolate distinct solutions.

The constraint programming community has contributed to dynamical systems. Janssen et al. propose a strategy [START_REF] Deville | Consistency Techniques in Ordinary Differential Equations[END_REF] dedicated to interval IVPs that has similarities with VNODE or CAPD. The integration step between two consecutive states (instants) used a box consistency algorithm [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF] to better contract the output state. A first attempt to create a constraint language extended to ODEs was made a few years later [START_REF] Cruz | Constraint Satisfaction Differential Problems[END_REF][START_REF] Cruz | Constraint Reasoning with Differential Equations[END_REF]. The language considers the entire ODE as a whole, but the unknown function has a special status that prevents direct manipulations: any constraint involving the unknown trajectory is managed by a specific operator. The fact that the function is bounded in an interval requires the introduction of an ad-hoc constraint called minimum/maximum restriction.

The link with real variables is also achieved via an ad-hoc constraint called value restriction. This language does not have the level of genericity targeted by our solver, where the concept of trajectories appears at the same level as other variables and ODEs are syntactic constructions among others. In our solver, we will rather use the C eval contractor [START_REF] Rohou | Reliable Nonlinear State Estimation Involving Time Uncertainties[END_REF] to handle value restriction constraints.

Although [START_REF] Goldsztejn | Including Ordinary Differential Equations Based Constraints in the Standard CP Framework[END_REF] improves [START_REF] Cruz | Constraint Satisfaction Differential Problems[END_REF] in both modeling possibilities and solving efficiency, it is still restricted to ODE constraints relating solutions values at specified times. Finally, our model can accept various differential constraints (like an integro-differential equation) and static constraints such as distance relations between states at different instants, as long as the corresponding contractors exist.

Outline. Section 2 introduces the notations used in the paper and the background useful to understand the following sections. Sections 3 and 4 detail our DACSP solver and its different procedures and parameters: contractors, choice point heuristics, etc. Section 5 reports first experimental results obtained on interval IVPs, but also integro-differential equations and BVPs.

Background and Notations

We first provide some background about intervals, inclusion functions and contraction. We then present several differential constraints and interval techniques adapted to dynamical systems.

Intervals

Contrary to numerical analysis methods that work with single values, interval methods can manage sets of values enclosed in intervals. Interval methods are known to be particularly useful for handling nonlinear constraint systems.

Definition 1 (Interval, box, box width, box hull) An interval [x i ] = [x i , x i ] defines the set of reals x i such that x i ≤ x i ≤ x i . IR denotes the set of all intervals. A box [x] denotes a Cartesian product of intervals [x] = [x 1 ]×...×[x n ]. The size, width or diameter of a box [x] is given by w([x]) ≡ max i (w([x i ])) where w([x i ]) ≡ x i -x i .
The hull of boxes approximates the union operator. It returns the smallest box enclosing all the boxes hulled.

Interval arithmetic [START_REF] Moore | Interval Analysis[END_REF] has been defined to extend to IR the usual mathematical operators over R For instance, the interval sum is defined by

[x 1 ] + [x 2 ] = [x 1 +x 2 , x 1 +x 2 ].
When a function f is a composition of elementary functions, an inclusion function [f ] of f must be defined to ensure a conservative image computation. There are several inclusion functions. The natural inclusion function of a real function f corresponds to the mapping of f to intervals using interval arithmetic. For instance, the natural inclusion function

[f ] N of f (x) = x(x + 1) in the domain [x] = [0, 1] computes [f ] N ([0, 1]) = [0, 1] • [1, 2] = [0, 2]
. Another inclusion function is based on an interval Taylor form [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF].

Interval arithmetics can be used for solving the numerical CSP (NCSP), i.e. finding solutions to an NCSP instance P = (x, [x], c), where x is an n-set of variables taking their real values in the domain [x] and c is an m-set of numerical constraints using operators like +, -, ×, a b , exp, log, sin, etc. NCSP solvers, like GlobSol [START_REF] Kearfott | GlobSol: History, Composition, and Advice on Use[END_REF], Gloptlab [START_REF] Domes | GLOPTLAB: A configurable framework for the rigorous global solution of quadratic constraint satisfaction problems[END_REF], RealPaver [START_REF] Granvilliers | RealPaver: An Interval Solver using Constraint Satisfaction Techniques[END_REF] or Ibex [START_REF] Chabert | IBEX -an Interval-Based EXplorer[END_REF] to name a few, follow a Branch and Contract method to solve an NCSP. The branching operation subdivides the search space by recursively bisecting variable intervals into two subintervals and exploring both sub-boxes independently. The combinatorial nature of this tree search is not always observed thanks to the contraction (filtering) operations applied at each node of the search tree. Informally, a contraction applied to an NCSP instance can reduce the domain without losing any solution. A contractor used in this paper is the well-known HC4-revise [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF][START_REF] Messine | Méthodes d'optimisation globale basées sur l'analyse d'intervalle pour la résolution des problèmes avec contraintes[END_REF], also called forward-backward. This contractor handles a single numerical constraint and obtains a (generally non optimal [START_REF] Collavizza | Comparing Partial Consistencies[END_REF]) contracted box including all the solutions of that constraint. To contract a box w.r.t. an NCSP instance, the HC4 algorithm performs a (generalized) AC3-like propagation loop applying iteratively the HC4-Revise procedure on each constraint individually until a quasi fixpoint is obtained in terms of contraction. 3B-consistency [START_REF] Lhomme | Consistency Techniques for Numeric CSPs[END_REF] and CID-consistency [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF] are two other stronger consistencies, enforced on an NCSP, that are exploited by our solver. The corresponding 3B and CID algorithms should call their Var3B or VarCID procedure on all the NCSP variables for enforcing the 3B or CID consistency. In practice however the algorithms implemented apply these procedures on a subset of the variables to get a better tradeoff between contraction and performance. VarCID splits a variable interval in k subintervals, and runs a contractor, such as HC4, on the corresponding sub-boxes. The k sub-boxes contracted are finally hulled. Var3B is somehow a dual operator that tries to remove subintervals at the bounds of a variable interval. If a contraction, like HC4, applied to a sub-box, where the interval is replaced by a subinterval at a bound, leads to an empty domain, then it proves that the subinterval contains no solution and can be removed safely from the variable domain.

Trajectories and Tubes

A trajectory, denoted x(•) = (x 1 (•), .., x n (•)), is a function from [t 0 , t f ] ⊂ R to R n . The input (argument) of x(•) is named time in this article (and denoted • or t) while the output (image) is called state. A tube [x](•)
is the interval counterpart of a trajectory and is defined as an envelope over the same temporal domain [t 0 , t f ]. The concept appeared in [START_REF] Filippova | Ellipsoidal State Estimation for Uncertain Dynamical Systems[END_REF][START_REF] Kurzhanski | On the Theory of Trajectory Tubes -A Mathematical Formalism for Uncertain Dynamics, Viability and Control[END_REF] in the context of ellipsoidal estimations. In our solver, it is used as a domain on which we apply operations of contractions and bisections. We employ them as intervals of trajectories, which is consistent with the aforementioned tools.

Hence, we will use the definition given in [START_REF] Bars | Set-membership state estimation with fleeting data[END_REF] where a tube

[x](•) : [t 0 , t f ] → IR n is an interval of two trajectories [x(•), x(•)] such that ∀t ∈ [t 0 , t f ], x(t)
x(t). We also consider empty tubes that depict an absence of solutions. A trajectory x(•) belongs to the tube 1 illustrates a one-dimensional tube enclosing a trajectory x(•). Our choice is to represent numerically a tube as a set of boxes corresponding to temporal slices. More precisely, an n-dimensional tube [x](•) with a sampling time δ > 0 is implemented as a box-valued function which is constant for all t inside intervals [kδ, kδ

[x] (•) if ∀t ∈ [t 0 , t f ], x (t) ∈ [x] (t). Fig.
δ • [x](•) t f t1 t3 t 0 x(•) [x](•) t 1 t 3 t 0 output gate of [[x]] (2) slice [[x]] (2)
+ δ], k ∈ N. The box [kδ, kδ + δ] × [x] (t k ), with t k ∈ [kδ, kδ + δ], is called the k th slice of the tube [x](•) and is denoted by [[x]] (k) .
This implementation takes rigorously into account floating-point precision when building a tube: computations involving [x](•) will be based on its slices, thus giving a reliable outer approximation of the solution set. The slices may be of same width as depicted in Fig. 1, but the tube can also be implemented with a customized temporal slicing. Finally, we endow the definition of a slice [[x]] (k) with the slice (box) envelope (blue painted in Fig. 1) and two input/output gates [x](t k ) and [x](t k+1 ) (black painted) that are intervals of IR n through which trajectories are entering/leaving the slice.

Once a tube is defined, it can be handled in the same way as an interval. We can for instance use arithmetic operations as well as function evaluations. If f is an elementary function such as sin, cos or exp, we define f ([x](•)) as the smallest tube containing all feasible values:

f ([x](•)) = {f (x(•)) | x(•) ∈ [x](•)} .

Dynamical Systems and Differential-Algebraic CSP

A differential constraint relates one or several trajectories and/or real variables. Numerous types of differential constraints can be considered in our approach, including:

1. ẋ(•) = f x(•) (ODE) 2. ẋ(t) = f x(t) + t t0 x(τ )dτ (integro-differential equation) 3. x(t k ) = y, ẋ(•) = v(•) (evaluation constraint) 4. ∀t ∈ [t], x(t) ∈ [y] 5. x(t) = y(t + δ) (delay constraint)
The first one is the most widespread differential constraint (see Def. 2). The second constraint is a little more complicated in that the state at a given time depends on the sum (integral) of the previous states. The third evaluation constraint y = x(t k ) states that the trajectory goes through an uncertain real value in [y] at an uncertain time in [t k ]. The fourth constraint is the complementary, although more complicated, constraint of the evaluation. The fifth one imposes a delay constraint of unknown real value δ between two trajectories and is particularly useful for clock calibration purposes in autonomous systems [START_REF] Voges | Timestamp offset calibration for an IMU-camera system under interval uncertainty[END_REF].

The idea behind our approach is to decompose a differential-algebraic system into a set of such primitive constraints associated to contractors (similar to propagators in CSP solvers [START_REF] Prud'homme | Choco documentation[END_REF]) that belong, or can be added to, the Tubex library. We formally define below the following differential constraints used in the experimental part.

Definition 2 (ODE and integro

-differential equation) Consider x(•) : [t 0 , t f ] → R n , its derivative ẋ(•) : [t 0 , t f ] → R n ,
and an evolution function f : R n → R n , possibly non-linear. An ODE7 is defined by:

ẋ(•) = f x(•)
An integro-differential equation is defined by:

ẋ(t) = f x(t) + t t0
x(τ )dτ.

The ODEs considered are explicit, i.e. the evolution function f computes ẋ directly. These differential constraints can define dynamical systems.

Definition 3 (IVP, interval IVP, BVP)

The initial value problem (IVP) is defined by an ODE ẋ(•) = f x(•) and an initial condition x(t 0 ) = x 0 , where x 0 is a constant in R n . In an interval IVP, the initial condition is bounded by an interval, i.e. x(t 0 ) ∈ [x 0 ]. A boundary value problem (BVP) is defined by an ODE and a numerical constraint c x(t 1 )..x(t n ) = 0, where c : R n → R n and ∀i ∈ {1..n},

t i ∈ [t 0 , t f ].
A BVP generalizes an IVP in that no initial condition fully determines the trajectory at a unique instant. Instead, n algebraic constraints relate several states at times t 1 ..t n and enable the trajectory determination. Note that a condition known at any instant is equivalent to knowing the state at time t 0 . Indeed, numerical integration can propagate the information forward or backward in time indifferently. We are now in position to define the DACSP.

if MaxDiam([x](•)) ≤ maxT ubeDiam then return [x](•) L ← { ([x](•), null) } while L = ∅ /* Depth-first tree search */ do ([x](•), gate) ← Pop(L) ([x](•), [y]) ← Contraction(P , specialT imes, Integration, ε, true, gate) if MaxDiam([x](•)) ≤ maxT ubeDiam then solutionsList ← solutionsList ∪ {[x](•) } else ([x1](•), [x2](•), gate) ← Bisect ([x](•), bisectionP olicy) L ← {([x1](•), gate)} ∪ {([x2](•), gate)} ∪ L solutionsList ← TubeMerge(solutionsList, [x](•)) return solutionsList
Algorithm 1: The DACSP solver.

sibly different behaviors. This may impair the effectiveness of a unique generic algorithm. In practice however, the user may already have an intuition of some instants from which things should propagate. As a consequence, in addition to P , we allow the user to provide a set of special times, i.e. elements of the temporal domain that involve states of the trajectory x(•) and other static constraints. It allows this first solver to perform contraction more incrementally. The solver works in two main phases: a so-called slicing step splitting the temporal domain into time slices, followed by a tree search subdividing the vectorial tube [x](•). The last TubeMerge function compensates a potential clustering effect and merges together pairs of solution tubes that intersect along the temporal domain (on all the slices) in all the (x) dimensions. It is necessary when bisection is not used for identifying different solutions, but helps the solver to compute accurate trajectories.

The main precision parameter of the solver is maxTubeDiam, a size expressed as the maximum width over all the slices envelopes of the tube. The solver can indifferently compute "thin" trajectories of theoretical null volume (e.g., when dealing with pure IVPs) or "thick" trajectories (i.e., continua of trajectories, e.g. when dealing with interval IVPs). In the latter case, the user has to tune this precision parameter to get a good approximation of thick trajectories. A second user-defined #slicesM ax parameter is a maximum number of slices created during the solving, especially during the slicing phase. A large slice number leads to a better trajectory accuracy at the cost of worse performance. Note that the CPU time generally grows linearly in the number of slices.

The first slicing phase is performed by the first do..while loop interleaving contraction of P and time slicing ("discretization"). The latter splits several slices of [x](•) into two slices of equal temporal size. Three main slicing policies have been tested:

-(all) Split all the slices in two.

-(median) Compute for all the bounded slices a dx difference between the middle points of 2 consecutive gates, maximum over all the dimensions, i.e.

dx = max i |m([x k i ]) -m([x k-1 i ])|
, where m([x i ]) denotes the middle of [x i ]. Split half of the slices with the largest dx and all the unbounded slices.

-(average) Split the slices having a dx greater than the average value and all the unbounded slices.

If the loop terminates because the number of slices reaches #slicesM ax, the tree search will be in charge to get the maxT ubeDiam precision. This slicing phase seems to contradict the principles of most numerical algorithms that decide to subdivide a given time step adaptively. However, the Integration procedure called by the Contraction method carries out these adaptive time discretization steps, so that both mechanisms, i.e. integration and slicing phases, perform time discretization in a complementary manner.

The second phase performed by the second while loop is combinatorial. It implements a tree search branching on the domains of the trajectory variables, i.e. tubes. Although depth-first search is well-known in the CP community, to our knowledge, no prior work proposed to make choice points on tubes, defined as follows.

Definition 5 (Tube bisection) Let [x](•) be a tube of a trajectory x(•) defined over [t 0 , t f ]. Let t k be an instant in [t 0 , t f ], i a dimension in {1..n}, and

[x i ] the interval value of [x i ](•) at t k . Let mid(x i ) be xi+xi 2 . The tube bisection (t k , i) of [x](•) produces two tubes [x L ](•) and [x R ](•) equal to [x](•) except at time t k , where [x L i ] = [x i , mid(x i )] and [x R i ] = [mid(x i ), x i ].
In

Contractions in the Solver

The Contraction function consists of a simple propagation loop that calls the contractors corresponding to constraints in c until the relative gain in contraction volume is less than ε f p . Contractors can be of any type: HC4-Revise for a numerical constraint or the "map" contractor mentioned in introduction for the robotic application. The propagation loop is followed by a call to a contractor Dyn3B enforcing a strong consistency on the tube (see Algorithm 3).

Let us detail in Algorithm 2 an important contractor, called ExplicitDE, that carries out tube contractions based on ODEs or integro-differential equations. The procedure is mainly responsible for launching integration steps forward and backward in time through the tube. The actual integration method used is a parameter of Algorithm 2. Note that a guaranteed integration algorithm inferring a new information, like a value of the state known at a specific instant, is incremental in that it may contract only a subset of a tube if no more contraction is obtained at a given gate. The tube contraction is not incremental in the first slicing phase or at the top of the search tree (gateBis = null) because the Slicing procedure can subdivide numerous slices everywhere in the tube (parameter isIncremental set to false). Therefore integration is run from t 0 to t f (forward) and from t f to t 0 (backward). Conversely, during the tree search, integration is triggered by a tube bisection or a domain modification at a special time (whose state is related with a static variable). That is why incremental integrations start from each of these instants. 8 Our solver is endowed with two possible Integration procedures. The first one is an "internal" generic integration algorithm that will be incorporated in the Tubex library. Its signature is close to the procedure Integration shown in Algorithm 2. It can be triggered from any specified time in the tube, forward or backward, and with the

Algorithm ExplicitDE(f , [x](•), specialT imes, Integration, ε, isIncremental, gateBis) if gateBis = null or not isIncremental then [x](•) ← Integration(f , [x](•), t0, FORWARD, ε, f alse) [x](•) ← Integration(f , [x](•), t f , BACKWARD, ε, f alse) else gates ← Sort({gateBis} ∪ specialT imes)
forall gate ∈ gates do // forward and incremental simulation:

[x](•) ← Integration(f , [x](•), gate, FORWARD, ε, true)
forall gate ∈ gates, in reverse order do // backward and incremental simulation:

[x](•) ← Integration(f , [x](•), gate, BACKWARD, ε, true)
Algorithm 2: A generic contractor for ODE and integro-diff equation possibility of running the simulation incrementally, i.e. stopping it if no sufficient contraction volume gain has been obtained in a gate box. This procedure Integration is generic in that it can accept an evolution function f describing either an ODE or an integro-differential equation (see Definition 2). It can also be specialized by a "slice integration contractor" called at each time step of the simulation. Two slice contractors are highlighted in this paper. The first one, called DynBasic hereafter, wraps at the slice level two simple contractors available since the very first version of Tubex: CtcDeriv and an evaluation of the evolution function f called iteratively. CtcDeriv (denoted C d dt in the literature [START_REF] Rohou | Guaranteed computation of robot trajectories[END_REF]) is a tube contractor treating the constraint ẋ(•) = v(•), where x(•) and v(•) are two trajectories and v(•) is the derivative of x(•) over time. The fundamental theorem of calculus that relates differentiation and integration, is used by CtcDeriv for contracting the tube [x](•). The second slice contractor, called DynCIDGuess hereafter, generates for each integration step a "slice" contractor graph, where the variables correspond to the two gates and the slice envelope (see Sec. 2.1). Based on the input gate box and the envelope, DynCIDGuess can improve the output gate box using sophisticated singleton consistencies based on 3B and CID (see Section 2.1). This contractor will be detailed in another article. This generic integration contractor starts by calling a Picard operator that allows one to set non-infinite initial bounds on some tube [x](•), which is required for engaging contraction, and can create new slices adaptively [START_REF] Nedialkov | Validated Solutions of Initial Value Problem for Ordinary Differential Equations[END_REF].

A second Integration procedure wraps directly the state-of-the-art VN-ODE [START_REF] Nedialkov | A Validated Solver for Initial Value Problems in Ordinary Differential Equations[END_REF] guaranteed integration solver into a CtcVnode contractor. During the slicing phase, it calls VNODE simulations forward and backward from the smallest gate. After each bisection, it calls VNODE simulations forward and backward starting from the bisected gate. To make CtcVnode a contractor, the results obtained by the VNODE simulations are intersected with the current tube. VNODE performs its own slicing, especially in the first iterations, and the slices produced are added to those from the slicing phase. Finally, as we will see in the experiments, the contractor CtcVnode, and a slice integration contractor as DynBasic or DynCIDGuess, can be called successively inside the contraction loop performed by the Contraction procedure.

Another new and useful dynamic contractor is the Dyn3B contractor described in Algorithm 3. This is a dynamic adaption of the 3B algorithm described in Section 2.1. It selects iteratively the instant (gate) t k with the largest interval (the tube is thus not contracted at all instants) and applies a VarDyn3B shaving procedure to all the [x i ] intervals at t k . VarDyn3B is a straightforward adaptation of the standard Var3B shaving procedure (see Section 2.1) to tubes. Subintervals at the bounds of [x i ] can be safely eliminated if an integration starting from the corresponding sub-tube leads to an empty domain. This integration procedure can be achieved by DynBasic, or CtcVnode followed by DynBasic.

Algorithm Dyn3B (P , ε) do volumeSave ← volume([x](•)) t k ← SelectGate([x](•)) forall i ∈ {1..n} do [x](•) ← VarDyn3B(P , t k , i) while VolumeGain([x](•), volumeSave) > ε3B
Algorithm 3: The Dynamic 3B algorithm

Experiments

The goal of this section is to highlight that the DACSP model, the contractors available via Tubex and our DACSP solver can handle a large variety of systems that no competitor or a few ones can deal with. All the results have been obtained on a CPU computer using an x86-64 processor (1.6 GHz).

BVP for Integro-Differential Equation

Let us illustrate the versatility of our DACSP solver on the following problem.

It combines an integro-differential equation defined on the domain [0, 1] and a constraint between the initial and final values, as follows:

ẋ(t) = 1 -2x(t) -5 t 0 x(τ )dτ ; t ∈ [0, 1] x(0) 2 + x(1) 2 = 1 (1)
Our solver can find both solutions in 8.35 seconds and needs to resort to 66 bisections (and a search tree depth of 25) to isolate them at a good accuracy.

For both solutions, Table 1 reports some details. Note that only our generic Integration algorithm can be used in the solver for this particular problem since there is no ODE, contrarily to the following DACSP systems. It has been run with maxT ubeDiam = 0.02 and #slicesM ax = 400.

For the next two DACSP categories tested, we show the best combination of the CtcVnode (refered by vnode in the tables), DynBasic (basic), DynCIDGuess (CIDG) and Dyn3B (3Bvnode or 3Bbasic) contractors.

BVPs and Cruz & Barahona System

We have tested and reported in Table 2 five BVPs and the Cruz system close to a BVP because no state is fully determined at a given instant. However, note that this system has a thick tube solution.

       ẋ(•) = x(•) x(0) 2 + x(1) 2 = 1 [t0, t f ] = [0, 1] maxT ubeDiam = 0.0005 (2)        ẍ(•) = -x(•) x(0) = 0; x(π/2) = 2 [t0, t f ] = [0, π/2] maxT ubeDiam = 0.0005 (3)            ẍ(•) = 5 ẋ(•) x(0) = 1; x(1) = 0 ẋ(0) ∈ [-10, 10]; ẋ(1) ∈ [-10, 10] [t0, t f ] = [0, 1] maxT ubeDiam = 0.02 (4)            ẍ(•) = -10( ẋ(•) + x(•) 2 ) x(0) = 0; x(1) = 0.5 ẋ(0) ∈ [-20, 20]; ẋ(1) ∈ [-20, 20] [t0, t f ] = [0, 1] maxT ubeDiam = 0.05 (5)                ẍ(•) = -exp(x(•)) x(0) = 0; x(1) = 0 ẋ(0) ∈ [-20, 20] ẋ(1) ∈ [-20, 20] [t0, t f ] = [0, 1] maxT ubeDiam = 0.05 (6)                ẋ1(•) = -0.7x1(•) ẋ2(•) = 0.7x1(•) -(ln(2)/5)x2(•) x1(0) = 1.25 x2 ∈ [1.1, 1.3] during [1, 3] [t0, t f ] = [0, 6] maxT ubeDiam = 0.04 (7)
Fig. 2. Five BVPs and the Cruz system. (2) A one-dimensional problem with an algebraic constraint between the initial and final states; (3) Classical linear example cited in Wikipedia; (4) and (5) denote resp. Systems 2 and 23 in the BVPSolve benchmark [START_REF] Mazzia | Solving boundary value problems in the open source software R: Package bvpSolve[END_REF]; (6) the Bratu system, the only one with two solutions in the BVPSolve benchmark, (7) the Cruz system with a partial information in the middle of the temporal domain.

Interval IVPs

Although solving interval IVPs is not the primary purpose of the DACSP solver, we present results obtained on three interval IVPs (results in Table 3).

Table 2. Solutions obtained on BVP systems. For each system, strategy and solution (s1 and/or s2), we report the diameters of the two unknown states in t0 and t f (most of the systems tested are 2-dimensional, but 2 of the 4 bounds are provided as initial conditions), the volume of the solution tubes, the number of slices, the computational time and the number of choice points required (#bis.).

Sys

       ẋ = -x 2 x(0) ∈ [0.1, 0.4] [t0, t f ] = [0, 5] eps = 0.2 (8)
             ẋ1 = -x1 -2x2 ẋ2 = -3x1 -2x2 x1(0) ∈ [5.9, 6.1] x2(0) ∈ [3.9, 4.1] [t0, t f ] = [0, 1] eps = 0.5 (9)
             ẋ1 = -x2 + 0.1 x1 (1 -x 2 1 -x 2 2 ) ẋ2 = x1 + 0.1 x2 (1 -x 2 1 -x 2 2 ) x1(0) ∈ [0.7, 1.3] x2(0) = 0.0 [t0, t f ] = [0, 5] eps = 0.15 (10)
Fig. 3. Three interval IVPs tested. ( 8) and ( 9) were introduced in [START_REF] Deville | Consistency Techniques in Ordinary Differential Equations[END_REF]. (10) describes a limit cycle and is particularly sensitive to the wrapping effect caused by interval computation.

For the 2 examples from [START_REF] Deville | Consistency Techniques in Ordinary Differential Equations[END_REF], the exact solution is known, so we also report the relative error (column Gap in Table 3) on interval width of the final gates.

Discussion on Experiments

Note first that the VNODE solver alone (outside the DACSP solver) cannot cope with BVPs and is not efficient on the interval IVPs selected. CtcVnode (inside the DACSP solver) often provides a very good accuracy on the gates. The good performance is probably due to the high-order interval Taylor form used (order 11 has been set for the experiments). However, CtcVnode does generally not obtain good contraction on the whole tube (slice envelopes). Additional work is required to envisage obtaining a better tube volume accuracy using CtcVnode.

Table 3. Solutions obtained on interval IVPs systems. For each system, we report the best strategy, the diameters of the state at t f , the volume of the solution tube, the number of slices, the computational time and the number of bisections. DynCIDGuess alone is generally not efficient, except on Systems ( 7) and ( 8), because it requires too many slices to reach the precision (recall that the CPU time generally grows linearly in the number of slices). Finally, the best option for the Integration procedure is generally to call first CtcVnode and then DynBasic or DynCIDGuess. A final call to Dyn3B is useful for Systems (4) (using CtcVnode as subcontractor) and (5) (using DynBasic).

Sys

Overall, different solver strategies provide the best results on the different systems tested. All the devices offered by the solver can be useful on different instances: contractors, slicing, choice points. When the best strategy includes a number of bisections, this means that the #slicesM ax has been reached and the solver resorts to choice points to better approximate the solution tube. An issue for future work is to better study the interplay between slicing and bisection in order to obtain a more generic DACSP solver that can work without the #slicesM ax parameter.

Conclusion and Future Work

We have presented a new generic solver that can handle together differential and static numerical constraints. The originality of the approach lies both in the underlying model considering trajectories as variables and in a novel backtracking mechanism applicable to DACSP. Our DACSP solver is endowed with an exploration operator that enables to bisect a tube at a chosen time. This allows the DACSP solver to better handle hard DACSP systems and accurately estimate distinct trajectories of problems having several solutions. We have shown on first experimental results that our solver is versatile enough to solve DACSP instances for which no or a few algorithmic solutions currently exist. We have also demonstrated the benefits of wrapping the state-of-the-art VNODE in a CtcVnode contractor implemented in Tubex.

With regard to future work, we will first try to limit the number of userdefined parameters, in particular remove #slicesM ax. Also, we want to propose ideally only one combination of contractors in the DACSP solver for every DACSP subclass. Second, we will study a more general search tree branching static and dynamical variables domains indifferently, though we need to explore new ideas on large-scale problems. Finally, in the current solver, the propagation between functional and real variables domains is somewhat naive and is partly ensured by the "special times" specified by the user (see Algorithm 2). The quite recent Tubex 3.0 accepts bi-level slice/tube variables, which will enable a fully incremental contraction achieved by a propagation engine. Supplementary materials including the sources of the solver and the experiments are available on http://simon-rohou.fr/research/dacsp-solve/.

Fig. 1 .

 1 Fig. 1. A one-dimensional tube [x](•), interval of two functions [x(•), x(•)], enclosing a random trajectory x(•) depicted in orange. The tube is numerically represented by a set of δ-width slices illustrated by blue boxes.

Definition 4 (

 4 Differential-algebraic CSP -DACSP) A DACSP network is defined by a quintuplet (x(•), [x](•), y, [y], c), where x(•) is a trajectory variable of domain [x](•) (a tube R → IR n1 , as defined in the previous section), y ∈ R n2 denotes the static numerical variables with a domain/box [y] and c denotes the set of static or differential constraints. Solving a DACSP instance consists in finding the set of values in [x](•) and [y] satisfying c.

3 A

 3 Generic Solver for Differential-Algebraic CSP In Algorithm 1, we give a description of a first generic solver for DACSP. The solver works on a network P = (x(•), [x](•), y, [y], c) and returns a set of trajectories satisfying c. The input tube [x](•) is defined generally with one single slice [t 0 , t f ] × [-∞, ∞] n . We attempt to tackle a wide class of problems with pos-Algorithm DACSP-Solver (P = (x(•), [x](•), y, [y], c), specialT imes, maxT ubeDiam, #slicesM ax, Integration, ε = (ε f pt , εintegr, ε3B), slicingP olicy, bisectionP olicy) do /* Slicing loop: */ ([x](•), [y]) ← Contraction(P , specialT imes, Integration, ε, f alse) [x](•) ← Slicing([x](•), slicingP olicy) while MaxDiam([x](•)) > maxT ubeDiam and #Slices(tube) < #slicesM ax

  practice, a bisection (t k , i) is applied only to a gate of the tube. Two heuristics are proposed to the user for selecting the instant t k . The first one picks randomly one instant among the "special times" specified by the user. The second one selects the t k having the largest box [x](t k ). The dimension i ∈ {1..n}, on which the bisection is performed, is decided according to the largest component [x i ].Note that the DACSP solver is sound because no operator used in Algorithm 1 can eliminate a solution: Contraction, Slicing, Bisect, TubeMerge.

Table 1 .

 1 Solutions obtained on the integro-differential based system. The table reports the diameters of the initial and final gates, the tube volume and the slices number.

	Solution	Diam. of gate t0 = 0	Diam. of gate t f = 1	Tube volume	#slices
	1	0.015	0.030	0.018	400
	2	0.034	0.022	0.024	400

An Initial Value Problem is composed of an ODE and an initial condition. Numerical integration propagates the initial value through the whole trajectory by integrating the evolution function of the ODE.

Note that a high-order problem can be transformed automatically into a first-order ODE shown in Def. 2 by introducing auxiliary variables. Also note that non autonomous ODEs of the form ẋ(t) = f x(t), t can also be transformed into autonomous ODEs ẋ(t) = f x(t) whose derivative depends only on the state.

The actual code is a little bit more complicated. An instant is skipped if it is handled by the previous integration step.

Acknowledgements. This work was supported by the French Agence Nationale de la Recherche (ANR) [grant number ANR-16-CE33-0024].


G Nehmetallah

T Chaumont-Frelet

S Descombes

S Lanteri

A POSTPROCESSING TECHNIQUE FOR A DISCONTINUOUS GALERKIN DISCRETIZATION OF TIME-DEPENDENT MAXWELL'S EQUATIONS

Keywords: time-domain electromagnetics, Maxwell's equations, discontinuous Galerkin method, high-order method, postprocessing

We present a novel postprocessing technique for a discontinuous Galerkin (DG) discretization of time-dependent Maxwell's equations that we couple with an explicit Runge-Kutta time-marching scheme. The postprocessed electromagnetic field converges one order faster than the unprocessed solution in the H(curl)-norm. The proposed approach is local, in the sense that the enhanced solution is computed independently in each cell of the computational mesh, and at each time step of interest. As a result, it is inexpensive to compute, especially if the region of interest is localized, either in time or space. The key ideas behind this postprocessing technique stem from hybridizable discontinuous Galerkin (HDG) methods, which are equivalent to the analyzed DG scheme for specific choices of penalization parameters. We present several numerical experiments that highlight the superconvergence properties of the postprocessed electromagnetic field approximation.

Introduction

Maxwell's equations are the most general model of electrodynamic theory [START_REF] Griffiths | Introduction to Electrodynamics[END_REF]. As a result, they are employed in a variety of applications, ranging from telecommunication engineering [START_REF] Russer | Electromagnetics, microwave circuit and antenna design for communications engineering[END_REF] to nanophotonics [START_REF] Gaponenko | Introduction to nanophotonics[END_REF], to study the propagation of an electromagnetic field and its interaction with structures and matter.

Nowadays, numerical schemes are routinely employed to simulate the propagation of electromagnetic waves by computing approximate solutions to Maxwell's equations [START_REF] Bondeson | Computational Eelectromagnetics[END_REF]. While several approaches, such as finite difference methods [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF], are available, we focus here on discontinuous Galerkin methods [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 Maxwll's equations[END_REF][START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF], which have recently received a lot of attention, due to their great flexibility and ability to handle complex geometries.

Even if currently available computational power allows for useful and realistic simulations, modeling accurately the propagation of electromagnetic fields in complex geometries is still a challenging and very costly task. As a result, numerical schemes are expected to be accurate and robust, but also very efficient and adapted to modern computer architectures.

In the context of finite element methods, postprocessing techniques are an attractive way to improve the accuracy of an already computed discrete approximation. In many cases, these techniques can increase the order of convergence of the method at a very moderate cost. In addition, they often have a "local" nature, which allows for the design of embarrassingly parallel implementations. As a result, postprocessing techniques and superconvergence have attracted a considerable attention in the past decades [START_REF] Andreev | Superconvergence of the gradient for quadratic triangular finite elements[END_REF][START_REF] Babuska | Validation of recipes for the recovery of stresses and derivatives by a computer-based approach[END_REF][START_REF]Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of laplaces, poissons and the elasticity equations[END_REF][START_REF] Goodsell | A unified treatment of superconvergent recovered gradient functions for piecewise linear finite element approximations[END_REF].

In this work, we elaborate a novel postprocessing technique for time-dependent Maxwell's equations. Following [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF], Maxwell's equations are discretized with a first-order discontinuous Galerkin method coupled with an explicit Runge-Kutta time-integration scheme [START_REF] Carpenter | Fourth-order 2N-storage Runge-Kutta schemes[END_REF]. This postprocessing improves the convergence rate in the H(curl)-norm by one order. As with similar postprocessing techniques devised in the past, our proposed approach is local, in the sense that the enhanced solution is computed independently in each cell of the computational mesh, and at each time step of interest. This is a key property as (a) it enables the design of highly parallel numerical algorithms, and (b) when the targeted application only requires the knowledge of the electromagnetic field in a limited region of space and/or time, the amount of computations is greatly reduced. Our postprocessing technique is inspired by two recent works, namely, a postprocessing for an explicit HDG discretization of the 2D acoustic wave equation [START_REF] Stanglmeier | An explicit hybridizable discontinuous Galerkin method for the acoustic wave equation[END_REF], and a postprocessing for a HDG discretization of the 3D time-harmonic Maxwell's equations [START_REF] Abgrall | Handbook of numerical methods for hyperbolic problems[END_REF].

We do not carry out the mathematical analysis of the proposed postprocessing but instead, we present a number of numerical experiments highlighting its main features. As a result, our work is organized as follows: in Section 2, we recall the settings and key notations related to Maxwell's equations, discontinuous Galerkin methods, and Runge-Kutta schemes. We describe our postprocessing in Section 3, and Section 4 presents numerical illustrations of the resulting methodology.

Settings

2.1. Maxwell's equations. We consider Maxwell's equations set in a Lipschitz polyhedral domain Ω ⊂ R 3 and in a time interval (0, T ). Specifically, given J : (0,

T ) × Ω → R 3 , the electromagnetic field E, H : (0, T ) × Ω → R 3 satisfies (1a) ε∂ t E -∇ × H = J , µ∂ t H + ∇ × E = 0,
in (0, T ) × Ω, where the functions ε, µ : Ω → R respectively represent the electric permittivity and the magnetic permeability of the materials contained in Ω. We assume that 0 < c ≤ ε, µ ≤ C a.e. in Ω for fixed constants c and C. The boundary of Ω is split into two subdomains Γ A and Γ P , and we prescribe the boundary conditions

(1b)    E × n Ω + µ ε (H × n Ω ) × n Ω = G on (0, T ) × Γ A , E × n Ω = 0 on (0, T ) × Γ P ,
where n Ω denotes the unit vector normal to ∂Ω pointing outward Ω and G : (0, T ) × Γ A → R 3 is a tangential load term (i.e. G • n Ω = 0). The first relation of (1b) is a first-order absorbing boundary condition (ABC) known as the Silver-Muller ABC. It is the simplest form of ABC for Maxwell's equations, and one could alternatively consider higher order ABCs [START_REF] Hagstrom | Radiation boundary conditions for Maxwell's equations: A review of accurate time-domain formulations[END_REF] or perfectly matched layers [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF]. The second equation in (1b) models the boundary of a perfectly conducting material. Finally, initial conditions are imposed in Ω (1c)

E| t=0 = E 0 , H| t=0 = H 0 ,
where E 0 , H 0 : Ω → R 3 are given functions.

Classically [START_REF] Assous | Mathematical foundations of computational electromagnetism[END_REF], under the assumption that the data µ, ε, J , G, E 0 and H 0 are sufficiently smooth, there exists a unique pair of solution (E, H) to [START_REF] Abgrall | Handbook of numerical methods for hyperbolic problems[END_REF].

We finally mention that in many applications, G is defined in order to inject an "incident" field in the domain. In this case, we have

(2) G := E inc × n Ω + µ ε (H inc × n Ω ) × n Ω ,
where (E inc , H inc ) is a solution to Maxwell's equations in free space. An important example that we will consider in Section 4 is the case where the incident field is a plane wave.

2.2. Mesh and notations. The domain Ω is partitioned into a mesh T h . We assume that T h consists of straight tetrahedral elements K, but hexahedral and/or curved elements could be considered as well. We assume that ε and µ take constant values ε K and µ K in each element K ∈ T h . For the sake of simplicity, we restrict our attention to meshes that are conforming in the sense of [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. Specifically, the intersection K -∩K + of two distinct elements K ± ∈ T h is either a full face, a full edge, or a single vertex of both K -and K + . In particular, hanging nodes are not covered by the present analysis. This is not an intrinsic limitation of the method, but this assumption greatly simplifies the forthcoming presentation.

We denote by F h the faces of the partition. Recalling that T h is conforming, each face F ∈ F h is either the intersection ∂K -∩ ∂K -of two elements K ± ∈ T h , or is contained in the intersection ∂K ∩ ∂Ω of a single element K ∈ T h with the boundary of the domain. We respectively denote by F int h , F P h and F A h the set internal faces, and the sets of faces belonging to Γ P and Γ A .

We associate with each face

F ∈ F h a unit normal n F , with the convention that n F = n Ω if F ∈ F P h ∪ F A h . If F ∈ F int h
, the orientation of the normal is arbitrary, but fixed. If v : Ω → R 3 is a function admitting well-defined traces on F ∈ F h , the notations v F and {{v}} F denote the "jump" and the "mean" of v on F . If F ∈ F int h with F = ∂K -∩ K + , these quantities are defined by

v F := v + | F (n + • n F ) + v -| F (n -• n F ), {{v}} F := 1 2 (v + | F + v -| F ) ,
where v ± := v| K ± and n ± denotes the unit outward normal to K ± , while we simply set

v F := {{v}} F := v| F , if F ∈ F P h ∪ F A h .
In the remaining of this work, k is a fixed non-negative integer representing a polynomial degree. For every element K ∈ T h , P k (K) denotes the set of polynomials defined on K of degree less than or equal k, and P k (K) := (P k (K)) 3 denotes the space of vector-valued functions having polynomial components. We finally employ the notation

P k (T h ) := v : Ω → R 3 | v| K ∈ P k (K) ∀K ∈ T h ,
for the space of piecewise polynomial functions. We also employ the notation P t k (F ) for the set of vector-valued polynomial functions defined on F that are tangential to F . P t k (F h ) is then the set of tangential polynomial defined on the skeleton of the mesh that are piecewise in P t k (F ).

2.3. The discontinuous Galerkin scheme. We seek the discrete fields as piecewise polynomial functions, namely E h , H h ∈ P k (T h ). Following [START_REF] Arnold | Unified analysis of discontinuous Galerkin, methods for elliptic problems[END_REF], the first step is to multiply (1a) by two test functions v and w, and integrate by parts over each element K ∈ T h . We obtain

(3) (ε∂ t E h , v) T h -(H h , ∇ × v) T h + H t h , v × n F h = (J , v), (µ∂ t H h , w) T h + (E h , ∇ × w) T h -E t h , w × n F h = 0, where E t h , H t h ∈ P t
k (F h ) are face-based tangential fields called "numerical fluxes", and

M t h , u × n F h := F ∈F h F M t h • ( u F × n F ) , for M t h ∈ P t k (F h )
, and u h ∈ P k (T h ). We make use of numerical fluxes in the spirit of local DG methods that were originally introduced in [START_REF] Castillo | An a priori error analysis of the local discontinuous Galerkin method for elliptic problems[END_REF] for scalar elliptic equations, and later in [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwll's equations[END_REF] for Maxwell's equations. We follow [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] to define our numerical fluxes. Specifically, we set Z K := µ K /ε K and Y K := 1/Z K for each K ∈ T h , and we select

E t h | F := 1 {{Y }} {{Y E h }} t F + 1 2 H h F × n , H t h | F := 1 {{Z}} {{ZH h }} t F - 1 2 E h F × n , for all F = ∂K -∩ ∂K + ∈ F int h ,

and

E t h | F := 0 H t h | F := -Y E h × n + H t h , if F = ∂K ∩ Γ P ∈ F P h , and
E t h | F := 1 2 E t h + ZH h × n + G × n , H t h | F := Y 2 ZH t h -E h × n -G , when F = ∂K ∩ Γ A ∈ F A h .
2.4. Time discretization. We can rewrite problem (3) obtained after space discretization as

(4) M Uh (t) + KU h (t) = B(t), U h (0) = U h,0
where for each t ∈ [0, T ], the vector U h (t) contains the coefficients defining E h (t) and H h (t) in the nodal basis of P k (T h ), M and K are the usual mass and stiffness matrices associated with (3), and U h,0 is the interpolation of the initial conditions in the discretization space. Classically, the key asset of DG schemes is that the mass matrix is block-diagonal, and hence, easy to invert. Thus, we may safely rewrite (4) as ( 5)

Uh (t) = -GU h (t) + F (t), U h (0) = U h,0 ,
where G := M -1 K and F (t) := M -1 B(t). At this point, we recognize in (5) a system of ordinary differential equations that can be discretized with a time marching scheme.

Here, we focus on a low storage Runge-Kutta scheme, usually denoted by LSRK(5,4), presented in [START_REF] Carpenter | Fourth-order 2N-storage Runge-Kutta schemes[END_REF] [START_REF]Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of laplaces, poissons and the elasticity equations[END_REF].

U h (t n ), t n := n∆t. Specifically, we let U 0 h := U h,0 , and for n ≥ 0, U n+1 h is deduced from U n h through the following algorithm                  V 1 h = U n h V 2 h = a k V 2 h + ∆t GV 1 h + F (t n + c k ∆T ) V 1 h = V 1 h + b k V 2 h for k = 1, • • • , 5 U n+1 h = V 1 h
, where the coefficients a k , b k and c k are described in Table 1. Then, E h,n and H h,n are the element of P k (T h ) expended on the nodal basis with the coefficients stored in U n h . The above scheme is of particular interest as it is fourth-order accurate with respect to the time step ∆t while being memory efficient. Indeed, it only requires the storage of two coefficient vectors in memory.

Classically, as this time integration scheme is explicit, it is stable under a CFL condition linking together the mesh size h and the selected time step ∆t. Specifically, given a mesh T h , we fix the time step by [START_REF]Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of laplaces, poissons and the elasticity equations[END_REF] ∆t := α k min

K∈T h 1 c K V K A K
where, c K := 1/ √ ε K µ K is the wave speed in the element K, and V K and A K are respectively the volume and the area of K. The constant α k is selected according to the polynomial degree k. Here, we use the values listed in Table 2, that we obtained after testing the scheme on simple test-cases. Finally, to ease the discussions in numerical experiments below, we denote by N the number of time steps performed in each simulations.

A novel postprocessing

As discussed above, E h,n and H h,n are respectively meant to approximate E(t n ) and H(t n ). The purpose of this section is to introduce postprocessed solutions E h,n and H h,n that are more accurate representations of E(t n ) and H(t n ). This postprocessing is purely local in time, in the sense that the computation of E h,n and H h,n only involves E h,n and H h,n . It is also local in space as the computation are local to each element K ∈ T h . Actually,

E h,n | K (resp. H h,n | K ) only depends on E h,n | K (resp. H h,n | K )
, where K is the union of all elements K ∈ T h sharing (at least) one face with K.

Our approach closely follows previous works. Specifically, similar postprocessing strategies have been derived for the time-harmonic Maxwell's equations [START_REF] Abgrall | Handbook of numerical methods for hyperbolic problems[END_REF], as well as time-dependent acoustic wave equation [START_REF] Stanglmeier | An explicit hybridizable discontinuous Galerkin method for the acoustic wave equation[END_REF]. These works develop in the context of hybridizable discontinuous Galerkin (HDG) methods, but can be easily applied to the DG scheme under consideration, as we depict hereafter.

Our postprocessing hinges on element-wise finite element saddle-point problems. For each element K ∈ T h , there exists a unique pair (E h,n , p) ∈ P k+1 (K) × P k+2 (K)/R such that

     (∇ × E h,n , ∇ × w) K + (∇p, w) K = (∇ × E h,n , ∇ × w) K + E t h,n -E t h,n , n × ∇ × w ∂K , (E h,n , ∇v) K = (E h,n , ∇v) K ,
for all w ∈ P k+1 (K) and v ∈ P k+2 (K)/R. Similarly, for the magnetic field, there exists a unique pair (H h,n , q) ∈ P k+1 (K) × P k+2 (K)/R such that

     (∇ × H h,n , ∇ × w) K + (∇q, w) K = (∇ × H h,n , ∇ × w) K + H t h,n -H t h,n , n × ∇ × w ∂K , (H h,n , ∇v) K = (H h,n , ∇v) K ,
for all w ∈ P k+1 (K) and v ∈ P k+2 (K)/R. E h,n and H h,n are then our postprocessed approximations to E(t n ) and H(t n ).

The left-hand sides of the above definition lead to solve symmetric linear systems of small size. In addition, observing that the left-hand side is actually the same for the two postprocessing schemes, we deduce that only one matrix factorization is required per element.

The right-hand sides further show that for each K ∈ T h , the postprocessed field E h,n | K only depends on E h,n | K and the value at the flux E t h,n | F on each face F ∈ F K . In turn, since the flux is defined using the two elements sharing the face F , we see that E h,n | K depends on the values taken by E h,n on all the elements K sharing at least one face with K. A similar comment holds true for H h,n .

Numerical experiments

4.1. Standing wave in a cavity. We first consider a model problem given by the propagation of standing wave in unit cube Ω := (0, L) 3 , L := 1 m, with perfectly conducting walls (i.e. Γ P := ∂Ω and Γ A := ∅). Specifically, we consider Maxwell's equations (1) with right-hand sides J := 0, G := 0 and initial conditions where the angular frequency is given by ω := √ 3πc 0 /L, c 0 := 1/ √ ε 0 µ 0 being the speed of light.

E| t=0 :=   -cos(πx 1 ) sin(πx 2 ) sin(πx 3 ) 0 sin(πx 1 ) sin(πx 2 ) cos(πx 3 )   , 0 0.2 0.4 0.6 0.8 1 •10 -8 10 -1 10 -2 10 -3 Time (s) ∇ × (E(t n ) -E h,n ) Ω ∇ × (E(t n ) -E h,n ) Ω
We consider structured meshes T h that are obtained by first splitting Ω into n × n × n cubes (n := L/h), and then splitting each cube into 6 tetrahedra.

Figures 1 and2 show the behavior of the error for the original and postprocessed discrete solutions with respect to time on a fixed mesh built from a 8 × 8 × 8 Cartesian partition. The time step ∆t is selected following CFL condition [START_REF]Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of laplaces, poissons and the elasticity equations[END_REF]. Both the original and the postprocessed error exhibit an oscillatory behavior, which is typical of this particular test case. The postprocessed solution is about 10 times more accurate than the original one. Table 4.1 presents in more detail our results on a series of meshes and for different polynomial degrees. We see that in each case, the curl of the postprocessed solution converges with the expected order, namely k + 1.

4.2.

Plane wave in free space. We now consider the propagation of a plane wave in free space. Specifically, we consider Maxwell's equations [START_REF] Abgrall | Handbook of numerical methods for hyperbolic problems[END_REF] with Ω := (0, L) 3 , L := 1 m, Γ

P := ∅ 0 0.2 0.4 0.6 0.8 1 •10 -8 10 -1 10 -2 10 -3 Time (s) ∇ × (H(t n ) -H h,n ) Ω ∇ × (H(t n ) -H h,n )) Ω Figure 2.
Standing wave in a cubic cavity: time evolution of the error on the magnetic field.

h ∇ × (E(T ) -E h,N ) Ω ∇ × (E(T ) -E h,N ) Ω 1/4
Time (s) ∇ × (E(t n ) -E h,n ) Ω ∇ × (E(t n ) -E h,n ) Ω Figure 3.
Plane wave in free space: time evolution of the error on the electric field. and Γ A := ∂Ω. J := 0, and G is defined by ( 2) with

E inc (t, x) := p cos ω t - d • x c 0 , H inc (t, x) := ε 0 µ 0 d × E inc (t, x),
where p := (1, 0, 0) T is the polarization, d := (0, 0, 1) T is the direction of propagation and ω := 6πc 0 /L is the angular frequency. We impose the initial conditions (1c) with E 0 := E inc | t=0 and H 0 := H inc | t=0 . Then, since the medium under consideration is homogeneous, no reflection and/or diffraction occur, and the analytical solution is simply E = E inc and H = H inc . We select the simulation time T := 10 ns. As for the cubic cavity test, we consider structured meshes T h , that we obtain by first splitting Ω into n × n × n cubes (n := L/h), and then splitting each cube into 6 tetrahedra. As explained above, the time step is selected using [START_REF]Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of laplaces, poissons and the elasticity equations[END_REF]. Figures 3 and4 show the behaviour of the error for the original and postprocessed discrete solutions with respect to time on a fixed mesh based on a 12×12×12 Cartesian partition. The postprocessed solution is about 5 times more accurate than the original solution. Table 4.2 presents in more detail our results on a series of meshes and for different polynomial degrees. We see that in each cases, the curl of the postprocessed solution converges with the expected order, namely k + 1.

4.3.

Scattering of a plane wave by a dielectric sphere. We now consider a problem involving a dielectric sphere of radius 0.15 m with ε = 2ε 0 and µ = µ 0 . The computational domain is bounded by a cube of side 1 m on which the Silver-Muller absorbing condition is applied and the simulation time is T := 3 ns. We make use of an unstructured tetrahedral mesh, which consists of 32,602 elements with 565 elements in the sphere and ∆t is chosen via [START_REF]Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of laplaces, poissons and the elasticity equations[END_REF]. The right-hand sides J and G are the same than in Example 4.2, and the initial conditions are taken to be zero. We select P 2 elements, and denote by (E h , H h ) and (E h , H h ) the original and postprocessed solutions. As the analytical solution to the problem is unavailable, we compute a reference solution (E r , H r ) with P 4 elements on the same mesh and the time step is defined as ∆t r := ∆t/3. ∆t r is chosen as an integral division of ∆t to facilitate comparisons. We chose to divide ∆t by 3 since, following

Time (s) ∇ × (H(t n ) -H h,n ) Ω ∇ × (H(t n ) -H h,n )) Ω Figure 4.
Plane wave in free space: time evolution of the error on the magnetic field. for which CFL condition (6) holds true. We refer the reader to Figure 5 for a snapshot of the reference solution.

h ∇ × (E(T ) -E h,N ) Ω ∇ × (E(T ) -E h,N ) Ω 1/8
To assess the impact of the postprocessing, we consider a set of evaluation points A, and we compute relative errors

err(V ) 2 = N n=1 ||∇ × (V r (t n , A) -V h,n (A))|| 2 N n=1 ||∇ × (V r )(t n , A)|| 2 and err (V ) 2 = N n=1 ||∇ × (V r (t n , A) -V h,n (A))|| 2 N n=1 ||∇ × (V r )(t n , A)|| 2
with V := E or H. Table 4.3 shows that our postprocessing approach reduces the error by at least a factor of 2 for the 9 evaluation points that we have selected.

Conclusion

In this work we have presented a postprocessing approach for a discontinuous Galerkin discretization of the time-dependent Maxwell's equations in 3D. This postprocessing technique is inexpensive, and can be computed independently in each element of the mesh, and at every time step of interest. It is thus well adapted to parallel computer architectures. Moreover, it is particularly suited to applications requiring a higher accuracy in localized regions, either in time or space. We have presented numerical examples, both with analytical solution and in complicated geometries, that indicate that our postprocessing approach improves the convergence rate of the discrete solution in the H(curl)-norm by one order. Overall, this contribution is to be employed as an efficient way of reducing the H(curl)-norm error of discontinuous Galerkin discretizations.

Point

Field

Figure 1 .

 1 Figure 1. Standing wave in a cubic cavity: time evolution of the error on the electric field.

Figure 5 .

 5 Figure 5. Representation of |E r (T )| in the scattering example.

Table 1 .

 1 . After fixing a time-step ∆t, we iteratively construct approximations U n h of Values of the coefficients of the LSRK(5,4) scheme.

	Coeff Value	Coeff Value			Coeff Value
	a 1	0		b 1	1432997174477 9575080441755	c 1	0
	a 2	-	567301805773 1357537059087	b 2	5161836677717 1361206829357	c 2	1432997174477 9575080441755
	a 3	-	2404267990393 2016746695238	b 3	1720146321549 2090206949498	c 3	2526269341429 6820363962896
	a 4	-	3550918686646 2091501179385	b 4	3134564353537 4481467310338	c 4	2006345519317 3224310063776
	a 5	-	1275806237668 842570457699	b 5	2277821191437 14882151754819	c 5	2802321613138 2924317926251
			k	1	2	3	4
			α k 0.70 0.46 0.30 0.21

Table 2 .

 2 Values of α k in CFL condition

Table 3 .

 3 Standing wave in a cubic cavity: numerical convergence.

	7.99e-01	6.37e-01
	P 1 1/6 4.94e-01	(eoc 1.19) 2.69e-01	(eoc 2.13)
	1/8 3.65e-01	(eoc 1.05) 1.45e-01	(eoc 2.15)
	1/4 1.40e-01	3.80e-02
	P 2 1/6 6.55e-02	(eoc 1.87) 1.04e-02	(eoc 3.20)
	1/8 3.75e-02	(eoc 1.94) 4.24e-03	(eoc 3.12)
	1/4 2.05e-02	4.32e-03
	P 3 1/6 6.17e-03	(eoc 2.96) 9.29e-04	(eoc 3.74)
	1/8 2.62e-03	(eoc 2.98) 3.09e-04	(eoc 3.83)
	1/4 6.17e-01	4.18e-01
	P 1 1/6 3.76e-01	(eoc 1.22) 1.80e-01	(eoc 2.08)
	1/8 2.70e-01	(eoc 1.15) 9.71e-02	(eoc 2.15)
	1/4 9.94e-02	2.19e-02
	P 2 1/6 4.68e-02	(eoc 1.86) 6.00e-03	(eoc 3.19)
	1/8 2.71e-02	(eoc 1.90) 2.44e-03	(eoc 3.13)
	1/4 1.60e-02	2.46e-03
	P 3 1/6 4.83e-03	(eoc 2.95) 5.39e-04	(eoc 3.74)
	1/8 2.06e-03	(eoc 2.96) 1.82e-04	(eoc 3.77)

h ∇ × (H(T ) -H h,N ) Ω ∇ × (H(T ) -H h,N ) Ω

Table 2

 2

	, it is the smallest integer

Table 4 .

 4 Plane wave in free space: numerical convergence.

	5.37e-00	6.02e-00
	P 1 1/10 4.38e-00	(eoc 0.92) 3.99e-00	(eoc 1.84)
	1/12 3.75e-00	(eoc 0.86) 2.73e-00	(eoc 2.08)
	1/8 1.98e-00	7.92e-01
	P 2 1/10 1.36e-00	(eoc 1.70) 3.72e-01	(eoc 3.38)
	1/12 9.77e-01	(eoc 1.81) 2.08e-01	(eoc 3.18)
	1/8 4.63e-01	1.01e-01
	P 3 1/10 2.44e-01	(eoc 2.88) 4.25e-02	(eoc 3.87)
	1/12 1.43e-01	(eoc 2.93) 2.22e-02	(eoc 3.56)
	1/8 5.89e-00	6.01e-00
	P 1 1/10 4.68e-00	(eoc 1.03) 3.97e-00	(eoc 1.85)
	1/12 4.00e-00	(eoc 0.86) 2.75e-00	(eoc 2.03)
	1/8 2.16e-00	7.60e-01
	P 2 1/10 1.45e-00	(eoc 1.79) 3.71e-01	(eoc 3.21)
	1/12 1.03e-00	(eoc 1.89) 2.11e-01	(eoc 3.10)
	1/8 4.87e-01	1.01e-01
	P 3 1/10 2.54e-01	(eoc 2.93) 4.32e-02	(eoc 3.79)
	1/12 1.48e-01	(eoc 2.96) 2.29e-02	(eoc 3.48)

h ∇ × (H(T ) -H h,N ) Ω ∇ × (H(T ) -H h,N ) Ω

Table 5 .

 5 Scattering of a plane wave by a dielectric sphere: L 2 error between the reference solution and the solution with a P 2 interpolation with and without applying the postprocessing.

			err	err
		E	0.083 0.033
	A 1 (0, 0, 0.45)	H	0.103 0.048
		E	0.008 0.005
	A 2 (0.2, -0.3, 0.8)	H	0.008 0.006
		E	0.019 0.005
	A 3 (0.2, -0.3, 0.2)	H	0.020 0.006
		E	0.015 0.004
	A 4 (0.2, 0.3, 0.2)	H	0.017 0.005
		E	0.019 0.007
	A 5 (0.2, 0.3, 0.8)	H	0.027 0.007
		E	0.015 0.008
	A 6 (-0.2, -0.3, 0.8) H	0.014 0.008
		E	0.027 0.008
	A 7 (-0.2, -0.3, 0.2) H	0.028 0.008
		E	0.021 0.007
	A 8 (-0.2, 0.3, 0.2)	H	0.024 0.007
		E	0.010 0.005
	A 9 (-0.2, 0.3, 0.8)	H	0.011 0.005

† Laboratoire J.A. Dieudonné, Parc Valrose, 28 Avenue Valrose, 06108 Nice Cedex 02, 06000 Nice, France