Philippe Chartier
email: philippe.chartier@inria.fr

Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

Florian Méhats
email: florian.mehats@univ-rennes1.fr

Gilles Vilmart

Christian Lubich

A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations

Keywords: highly-oscillatory problems, stroboscopic averaging, standard averaging, micro-macro method, pullback method, uniform accuracy AMS subject classification (2010): 65L20, 74Q10, 35K15

We introduce a new methodology to design uniformly accurate methods for oscillatory evolution equations. The targeted models are envisaged in a wide spectrum of regimes, from non-stiff to highly-oscillatory. Thanks to an averaging transformation, the stiffness of the problem is softened, allowing for standard schemes to retain their usual orders of convergence. Overall, high-order numerical approximations are obtained with errors and at a cost independent of the regime.

Introduction

In this article, we consider time-dependent differential equations with (possibly high-) oscillations of the form uε (t) = f t/ε (u ε (t)) ∈ X, u ε (0) = u 0 ∈ X, t ∈ [0, 1], (1.1) where the function (θ, u) ∈ T × X → f θ (u) is assumed to be 1-periodic w.r.t. θ on the torus T = R/Z ≡ [0, 1], and where X denotes a Banach space endowed with a norm • X . The set X may be simply R d , in which case (1.1) is a simple ordinary differential equation in finite dimension, or it may be a functional space, say for instance the space of square-integrable functions L 2 (R), in applications to partial differential equations such as the Schrödinger or the Klein-Gordon equations. Here, ε is a parameter whose value can vary in an interval of the form ]0, 1]. We emphasize that we do not restrict (1.1) to its asymptotic regime where ε tends to zero, but also allow ε to have a non-vanishing value, say close to 1.

Remark 1.1. The autonomous problem

vε = 1 ε Av ε + f (v ε ) ∈ X, v ε (0) = v 0 ∈ X, t ∈ [0, 1], (1.2)
where the map θ → e θA is assumed to be 1-periodic, can be straightforwardly rewritten within the format (1.1). As a matter of fact, the change of variable v ε = e -t ε A u ε transforms (1.2) into (1.1) with f θ (u) = e -θA f e θA u .

Methods from the (recent and not so recent) literature are suitable either for the highlyoscillatory regime where ε is small, or for the non-stiff one where ε is close to 1. However, they usually fail for certain values of ε in the range ]0, 1], indeed becoming inaccurate or inefficient. For instance, the authors of this article have introduced several schemes [START_REF] Calvo | Numerical stroboscopic averaging for ODEs and DAEs[END_REF][START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF][START_REF] Chartier | Higher-order averaging, formal series and numerical integration I: B-series[END_REF][START_REF] Chartier | Higher-order averaging, formal series and numerical integration II: The quasi-periodic case[END_REF][START_REF] Chartier | Higher-order averaging, formal series and numerical integration III: error bounds[END_REF], inspired by the theory of averaging or the theory of modified equations, which happen to be efficient for small values of ε and quite inaccurate for larger ones. This degraded behaviour appears as a consequence of the use of truncated series with no consideration of induced error terms. Recently, a new class of uniformly accurate methods, namely two-scale methods [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Chartier | Highly-oscillatory evolution equations with multiple frequencies: averaging and numerics[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF], has emerged (see also the methods developed in [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF][START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF] using a different approach) that are capable of solving (1.1) with an accuracy and at a cost both independent of the value of ε ∈]0, 1]. Such methods rely on an explicit separation of scales u ε (t) = U ε (t, θ)| θ=t/ε with θ ∈ T and the resolution of the transport equation

∂ t U ε (t, θ) + 1 ε ∂ θ U ε (t, θ) = f θ (U ε (t, θ)) (1.3)
with a well-chosen initial condition U ε (0, θ) = U ε 0 (θ). A numerical approximation u ε k of the solution u ε (t k ) of (1.1) at time t k is then obtained from a p th -order approximation of U ε (t k , θ k ) by taking θ k = t k /ε (see [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF] for the details of the technique). At a given computational cost χ ≈ 1 ∆t , this yields approximations on an uniform grid t k = k∆t, k = 0, . . . , K, that satisfy the following error estimate max

0<ε≤1 max 0≤k∆t≤1 u ε k -u ε (k∆t) X ≤ C∆t p , (1.4)
where the constant C may depend on p but is independent of ∆t and ε. Although the error behaviour (1.4) is unarguably favourable, the resolution of the transport equation (1.3) as introduced in [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Chartier | Highly-oscillatory evolution equations with multiple frequencies: averaging and numerics[END_REF] requires to design specific schemes, a task that is not straightforward for high orders. This aspect of two-scale methods has left them open to criticisms and has motivated the search for a technique dealing with (1.1) in a more direct and more affordable way.

In this work, we adopt the following strategy: through a suitable time-dependent (and periodic) change of variable, we transform equation (1.1) into a smooth differential equation which is then amenable to a standard numerical treatment. By standard treatment, we mean here that any non-stiff numerical method (i.e. any integrator which could be applied to eq. (1.1) with ε = 1) can be used for all values of ε ∈]0, 1] and will exhibit its usual order of convergence when applied to the transformed equation. Unsurprisingly, the aforementioned change of variable is a key-ingredient of the theory of averaging. However, and this constitutes the main advantage of this strategy over the methods discussed above, no abrupt truncation is applied in the process, so that all the information originally contained in (1.1) is retained. In this way, uniformly accurate schemes are elaborated simply from standard numerical schemes by adjoining a discretization of the change of variable.

In Section 2, we discuss the main ingredients of averaging theory. Firstly, we introduce in Subsection 2.1 the periodic change of variable and the averaged vector field for both stroboscopic and standard averaging. In particular, we write down, (i) the equation (2.3) that they formally satisfy, (ii) the inferred iteration (2.4) from which they can be recursively built and, (iii) the corresponding defect (2.6). Subsection 2.2 completes this expository section with the regularity estimates obtained so far in previous works.

Section 3 is the core of our paper. Subsection 3.1 presents the two methods that are later on analysed, namely the micro-macro method and the pullback method, together with possible variants. In Subsection 3.2, the defect and its time-derivatives after n iterations are upper bounded by powers of ε. Finally, we show in Subsection 3.3 that the transformed equations obtained both by the micro-macro and pullback methods are smooth and indeed can be solved by non-stiff standard methods.

Section 4 describes the actual implementation of the two most interesting methods, that is to say on the one hand, the micro-macro method with standard averaging, and on the other hand, the pullback method with stroboscopic averaging. Both methods are uniformly accurate. The first one appears to be easier to implement and cheaper. The second one has the advantage of being geometric, allowing in particular for the conservation of first integrals.

Finally, Section 5 presents numerical experiments with the two techniques of Section 4 used in combination with standard methods of orders 2, 3 and 4. The Hénon-Heiles problem and the Klein-Gordon equation are used as test-cases. As expected, the stroboscopic pullback technique exhibits a clear-cut preservation of invariants. In practice, all methods are evidently uniformly accurate. Up to our knowledge, this section demonstrates the first example of geometric, high-order and uniformly accurate methods available in the literature.

Numerical methods pertaining to averaging theory

In this section, we recall the main result of averaging. Since our technique is intented for a large range of ε-values, the following assumption on (1.1) appears as a natural prerequisite.

Assumption 2.1.

There exists an open set K ⊂ X with bounded closure K and a open set U 0 ⊂ K, such that for all initial condition u 0 in U 0 and all ε ∈]0, 1], the solution u ε (t) of (1.1) exists and remains in K for all times t ∈ [0, 1].

Averaging in a nutschell

We have chosen here to follow closely the presentation in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF] as it appears to be more convenient for the forthcoming developments. Other presentations can be found in classical references, e.g. [START_REF] Perko | Higher order averaging and related methods for perturbed periodic and quasi-periodic systems[END_REF] and [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF], or in recent papers, with a more algebraic point of view [START_REF] Chartier | Higher-order averaging, formal series and numerical integration I: B-series[END_REF][START_REF] Chartier | Higher-order averaging, formal series and numerical integration II: The quasi-periodic case[END_REF][START_REF] Chartier | Higher-order averaging, formal series and numerical integration III: error bounds[END_REF]. Generally speaking, the purpose of averaging methods is to write the exact solution of (1.1) as the composition

u ε (t) = Φ ε t/ε • Ψ ε t • (Φ ε 0 ) -1 (u 0 ) (2.1) of a periodic change of variable (θ, u) ∈ T × K → Φ ε θ (u) and of the flow map (t, u) ∈ [0, 1] × K → Ψ ε t (u) of an autonomous differential equation of the form Ψε t (u) = F ε (Ψ ε t (u)), Ψ ε 0 (u) = u, (2.2)
where the map u ∈ K → F ε (u) is a smooth vector field. The rationale behind averaging is apparent: the dynamics is sought after as the superposition of a slow drift, captured by Ψ ε t , and of high oscillations, captured by Φ ε t/ε . Various design choices are possible for the map Φ ε and we refer to [START_REF] Chartier | Higher-order averaging, formal series and numerical integration II: The quasi-periodic case[END_REF] for a thorough investigation of their implications. In this paper, we shall impose that either Φ ε 0 or the average Φ ε coincides with the identity map id X , where

g (u) = 1 0 g θ (u)dθ
stands for the average on T of any continuous function g : T × X → X, two procedures referred to respectively as stroboscopic averaging and standard averaging:

• Stroboscopic averaging Φ ε 0 = id X .
Owing to periodicity, Φ ε t/ε = id X for times t ∈ εN that are integer multiples of the period, hence the name "stroboscopic". When the problem (1.1) under consideration embeds a structure (e.g. Hamiltonian structure) or possesses invariants (e.g. quadratic first integrals), stroboscopic averaging, being geometric, is to be preferred.

• Standard averaging Φ ε = id X . This choice is more widely used and leads to a more simple implementation. In the context of kinetic theory, it is related to Chapman-Enskog method (see for instance in [START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF]).

For the decomposition in (2.1) to exist, the maps Φ ε and F ε would be required to satisfy the equation

1 ε ∂ θ Φ ε θ (u) + ∂ u Φ ε θ (u) F ε (u) = f θ • Φ ε θ (u) with F ε = ∂ u Φ ε -1 f • Φ ε , (2.3)
obtained by separation of scales. Unfortunately, it is a known fact that equation (2.3) has no rigorous solution: only a formal solution under the form of a divergent1 series in ε can be exhibited. Truncation thereof generates errors which can be made polynomially or exponentially small w.r.t. ε, depending on the smoothness of the vector field f in (1.1) w.r.t. u. In averaging theory, as presented in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF], approximations Φ [n] and F [n] of Φ ε and F ε are thus constructed from equation (2.3), starting from Φ [0] = id X and

F [0] = f , as follows Φ [k+1] θ = id X -εG [k+1] + ε θ 0 f τ • Φ [k] τ -∂ u Φ [k] τ F [k] dτ, F [k] = ∂ u Φ [k] -1 f • Φ [k] , (2.4)
for k = 0, . . . , n -1. The functions G [k+1] could in principle be chosen arbitrarily: they vanish for stroboscopic averaging, or are determined to ensure that Φ [k] = id X , i.e.

G [k+1] = θ 0 f τ • Φ [k] τ -(∂ u Φ [k] τ ) f • Φ [k] dτ , (2.5)
for standard averaging. Since equation (2.3) is impossible to satisfy exactly, a natural measure of success of the averaging idea is given by the size of the defect, defined by

δ [k] θ (u) = 1 ε ∂ θ Φ [k] θ (u) + ∂ u Φ [k] θ (u) F [k] (u) -f θ • Φ [k] θ (u), k = 0, . . . , n. (2.6) Remark 2.2. Starting from the stroboscopic decomposition u ε (t) = Φ ε t/ε • Ψ ε t (u 0 ) with Φ ε 0 = id X ,
we may write (once again formally at this stage)

u ε (t) = Φ ε t/ε • Φ ε -1 Φε t/ε • Φ ε • Ψ ε t • Φ ε -1 Ψε t • Φ ε ( Φε 0 ) -1 (u 0 ).
It can indeed be checked, on the one hand, that

Φε = Φ ε • Φ ε -1 = Φ ε • Φ ε -1 = id X and Φε 0 = Φ ε 0 • Φ ε -1 = id X • Φ ε -1 = Φ ε -1 ,
and on the other hand, that

d dt Ψε t (u 0 ) = (∂ u Φ ε F ) • Φ ε -1 • Ψε t (u 0 ) = f • Φε • Ψε t (u 0 ) = F ε • Ψε t (u 0 ).
This mean that the maps Φε , Ψε and F ε are those corresponding to standard averaging and that we have the following relation on the defects,

δε θ = δ ε θ • Φ ε -1 .
Conversely, we may relate Φ ε to Φε through the formula

Φ ε θ = Φε θ • ( Φε 0 ) -1 .

Standard error estimates in averaging theory

The main result of [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF] now ensures that the maps Φ [n] and F [n] given by iteration (2.4) are well-defined and satisfy equation (2.3) up to small error terms of the form O(ε n ). It is derived under additional analyticity assumptions on f . As in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF], we thus introduce the complexification X C of X (which is just X C = C d in the case where X = R d ), defined as

X C = {U := u r + iu i , (u r , u i ) ∈ X 2 },
where u r = (U ) ∈ X and u i = (U ) ∈ X are the real and imaginary parts of U . The space X C is still a Banach space when endowed with the norm (which we will simply denote • in the sequel)

U X C := sup λ∈C * (λU ) X |λ| ,
and we can consider for ρ ≥ 0 the closed set

K ρ = {u + ũ : u ∈ K, ũ ∈ X C , ũ ≤ ρ}.
Now, given a bounded map (θ, u) ∈ T × K ρ → ϕ θ (u), we finally denote

ϕ ρ := sup θ∈T, u∈Kρ ϕ θ (u)
and if ϕ has p bounded derivatives w.r.t. θ ϕ ρ,p := max ν=0,...,p

∂ ν θ ϕ θ ρ .
The main hypothesis that was required in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF] is the following.

Assumption 2.3. The function (θ, u) → f θ (u) is continuous w.r.t. θ and there exists R > 0 such that for all θ ∈ T, the function f θ (•) can be extended to an analytic map from K 2R to X C . In addition, there exists M > 0 such that

f 2R ≤ M. (2.7)
An immediate consequence is that for all u 0 + ũ0 ∈ U 0 i where

U 0 i := u 0 + ũ0 : u 0 ∈ U 0 , ũ0 ∈ X C , ũ0 ≤ Re -M R ,
for all ε ∈]0, 1] and all t ∈ [0, 1], the solution u ε (t) of

uε (t) = f t/ε (u ε (t)), u ε (0) = u 0 + ũ0 , (2.8)
remains in K R . By Assumption 2.3, f θ is indeed Lipschitz on K R with Lipschitz constant M/R, so that by virtue of the Gronwall lemma and Assumption 2.1, u ε (t) remains in K R for ũ0 e M/R ≤ R. Now, the following result has been proved in the case of stroboscopic averaging.

Theorem 2.4 (See Castella, Chartier, Murua and Méhats [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF]). Suppose that Assumptions 2.1 and 2.3 are satisfied and let ε 0 := R/(8M ). For every positive integer n and positive ε such that (n+1)ε ≤ ε 0 , the maps Φ [n] and F [n] given by iteration (2.4) with all G [k] ≡ 0, are well-defined, Φ [n] is continuously differentiable w.r.t. θ and both Φ [n] and its average Φ [n] are invertible with analytic inverse on K R . Moreover, the following estimates hold for all 0 < ε < ε 0 ,

(i) Φ [n] -id R ≤ R 2(n + 1) (ii) F [n] R ≤ 2M (iii) δ [n] R ≤ 2M 2(n + 1) ε ε 0 n .
It is known that by choosing (n + 1) appropriately depending on ε, the upper-bound of δ [n] can be made exponentially small in ε, but this is not our main interest here.

The uniform accuracy paradigm for averaging methods

The vector field F [n] appearing in Theorem 2.4 has a Taylor expansion with respect to ε → 0 which coincides to (within error terms of size ε n ) with a formal series f + εF 2 + ε 2 F 3 + . . ., whose successive terms can be determined analytically [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF][START_REF] Chartier | Higher-order averaging, formal series and numerical integration II: The quasi-periodic case[END_REF]. Hence, it makes perfect sense for small values of ε to simply ignore the discrepancy in equation (2.1) and to approximate (1.1) from the differential equation u = F [n] (u) (complemented with the computation of Φ [n] ), yielding an error of size O(ε n ). However, if the value of ε is not so small -close to 1 for instance -, this procedure becomes completely inaccurate. The aim of the forthcoming techniques is precisely to resolve this difficulty.

Uniformly accurate averaging: micro-macro and pullback methods

A uniformly accurate method should thus satisfy two requirements: on the one hand, it should filter high-oscillations in one way or another -and thus rely on the principle of averaging -in order to be efficient for small values of ε, and on the other hand, it should incorporate the relevant information for large values of ε so as to remain accurate. This leads to the two following alternative methods considered in this paper.

Micro-macro method. If we insist on solving an autonomous equation of the form (2.2), i.e.

Ψ[n] t (u) = F [n] Ψ [n] t (u) , Ψ [n] 0 (u) = u,
then it is of paramount importance to retain the information contained in the remainder

G [n] t (u 0 ) = u ε (t) -Φ [n] t/ε • Ψ [n] t • (Φ [n] 0 ) -1 (u 0 ),
which obeys the following differential equation

d dt G [n] t (u 0 ) = uε (t) - 1 ε ∂ θ Φ [n] t/ε + ∂ u Φ [n] t/ε F [n] (Ψ [n] t (ǔ 0 )) = f t/ε Φ [n] t/ε • Ψ [n] t (ǔ 0 ) + G [n] t (u 0 ) -f t/ε Φ [n] t/ε • Ψ [n] t (ǔ 0 ) -δ [n] t/ε (Ψ [n] t (ǔ 0 )) with ǔ0 = (Φ [n] 0 ) -1 (u 0 )
. Upon using the Gronwall lemma, it can be shown (see Section 3.3) that G

[n] t (u 0 ) and its n first time-derivatives are uniformly bounded w.r.t. ε provided δ [n] (see relation (2.6)) and its first n derivatives with respect to θ are of size ε n . This leads to the following micro-macro

system v = F [n] (v), v(0) = (Φ [n] 0 ) -1 (u 0 ) ẇ = f t/ε Φ [n] t/ε (v) + w -1 ε ∂ θ Φ [n] t/ε + ∂ u Φ [n] t/ε F [n] (v), w(0) = 0 (3.1)
from which the solution of (1.1) can be recovered as u ε (t) = Φ

[n]

t/ε (v(t)) + w(t).

Pullback method. If we insist on satisfying an equation of the form (2.1) exactly, i.e.

u ε (t) = Φ [n] t/ε • Ψ [n] t • (Φ [n] 0 ) -1 (u 0 ), (3.2)
then an autonomous equation of the form (2.2), i.e.

Ψ[n] t (ǔ 0 ) = F • Ψ [n] t (ǔ 0 ), Ψ [n] 0 (u 0 ) = ǔ0 ,
where ǔ0 = (Φ

[n] 0 ) -1 (u 0 ), can not hold and should instead be amended to

Ψ[n] t (ǔ 0 ) = F [n] (Ψ [n] t (ǔ 0 )) + R [n] t/ε (Ψ [n] t (ǔ 0 )), Ψ [n] 0 (ǔ 0 ) = ǔ0 . (3.3)
As a matter of fact, by differentiating (3.2) as follows (with ǔ0 = (Φ

[n] 0 ) -1 (u 0 )) 0 = uε (t) - 1 ε ∂ θ Φ [n] t/ε (Ψ [n] t (ǔ 0 )) -∂ u Φ [n] t/ε (Ψ [n] t (ǔ 0 )) Ψ[n] t (ǔ 0 ) = ∂ u Φ [n] t/ε (Ψ [n] t (ǔ 0 )) (F [n] (Ψ [n] t (ǔ 0 )) - Ψ[n] t (ǔ 0 )) -δ [n] t/ε (Ψ [n] t (ǔ 0 )), we obtain (3.3) provided that R [n] θ = -∂ u Φ [n] θ -1 δ [n] θ . (3.4)
This leads to the following pulled back equation

v = ∂ u Φ [n] t/ε -1 f t/ε • Φ [n] t/ε (v) - 1 ε ∂ θ Φ [n] t/ε (v) , v(0) = (Φ [n] 0 ) -1 (u 0 ), (3.5)
from which the solution of (1.1) can be recovered as

u ε (t) = Φ [n] t/ε (v(t)).
Both the micro-macro (3.1) and the pullback (3.5) formulations will later prove to be useful for the design of uniformly accurate numerical schemes. They can be considered for stroboscopic and standard averagings, yielding four different approaches (see Table 1).

Stroboscopic averaging

Standard averaging Pullback method Φ

[n] 0 = id X , G [n] ≡ 0 Φ [n] = id X , G [n] ≡ 0 Micro-macro method Φ [n] 0 = id X , R [n] ≡ 0 Φ [n] = id X , R [n] ≡ 0
Table 1: Designing choices for uniformly accurate averaging.

The numerical resolution of either system (3.1) or system (3.5) by a standard non-stiff method now requires estimates of time-derivatives of their right-hand side. As discussed above, this boils down to getting estimates of ∂ k t δ

[n]

t/ε for k = 0, . . ., which is the objective of next subsection.

Bounds on the defect δ ε

We introduce the two non-linear operators, which, to a map (θ, u)

∈ T × K ρ → ϕ θ (u) with invertible partial derivative ∂ u ϕ θ associate the mappings (θ, u) ∈ T × K ρ → Λ(ϕ) θ (u) and (θ, u) ∈ T × K ρ → Γ ε (ϕ) θ (u) defined as Λ(ϕ) θ (u) = f θ • ϕ θ (u) - ∂ϕ θ ∂u (u) ∂ ϕ ∂u (u) -1 f • ϕ (u), (3.6)
Γ ε (ϕ) θ (u) = u + ε θ 0 Λ(ϕ) ξ (u) dξ -γ θ 0 Λ(ϕ) ξ (u) dξ , (3.7)
with γ = 0 for stroboscopic averaging and γ = 1 for standard averaging. In the sequel, we shall constantly need the following assumption (more stringent than our previous Assumption 2.3).

Assumption 3.1. The function (θ, u) → f θ (u) is p-times (p ≥ 1)
continuously differentiable w.r.t. θ. Moreover, there exists R > 0 such that for all θ ∈ T, the function u → f θ (u) and all its derivatives u → ∂ ν θ f θ (u) up to order ν = p can be extended to analytic maps from K 2R to X C that are uniformly bounded w.r.t. (θ, u) ∈ T × K 2R . More precisely, there exists a constant M > 0 such that for all 0 < σ ≤ 3, max ν=0,...,p

σ ν ν! ∂ ν θ f 2R ≤ M. (3.8)
Stating upper-bounds (polynomial and exponential in terms of ε) of G [n] or R [n] now boils down to getting upper-bounds, as well as regularity estimates, of the defect δ [n] . Proposition 3.2. Suppose that Assumptions 2.1 and 3.1 are satisfied and let ε 0 := R/(8M ) and denote r n = R n+1 and R k = 2R -kr n , k = 0, . . . , n + 1. For all positive integer n and positive ε such that (n + 1)ε ≤ ε 0 , the maps Φ [k] , k = 0, . . . , n, defined by iteration (2.4) are analytic on K R k and both Φ [n] and its average Φ [n] are invertible with analytic inverse on K R . Furthermore, all maps Φ [k] are (p + 1)-times continuously differentiable w.r.t. θ with derivatives analytic on K R k and upper-bounded:

∀k = 0, . . . , n, ∀1 ≤ ν ≤ p + 1, ∂ ν θ Φ [k] R k ≤ 8 ε M ν! Proof.
Owing to Remark 2.2, we can restrict ourselves to the case of stroboscopic averaging, for which the (rather technical) proof is given in Appendix section.

Theorem 3.3. Suppose that Assumptions 2.1 and 3.1 are satisfied and let ε 0 := R/(8M ). For all positive integer n and positive ε such that (n + 1)ε ≤ ε 0 , let the map Φ [n] and F [n] be defined by iteration (2.4). Then the defect δ [n] defined by (2.6) is p-times differentiable and satisfies the following estimate:

δ [n] R,p ≤ 2 M 2 Q p (n + 1) ε ε 0 n , where Q p ≤ 5 + 2 p+1 e 5/2 p!.
For the proof of Theorem 3.3, we need the following lemma with proof postponed in appendix.

Lemma 3.4. Let 0 < δ < ρ ≤ 2R and consider two periodic, near-identity mappings (θ, u) ∈ T × K ρ → ϕ θ (u) and (θ, u) ∈ T × K ρ → φθ (u), analytic on K ρ and satisfying

ϕ -Id ρ ≤ δ 2 and φ -Id ρ ≤ δ 2 .
(3.9)

Let ε 0 := R/(8M ) and suppose that they are p-times continuously differentiable w.r.t. θ and satisfy the estimates

∀ 0 < ε < ε 0 , ∀ 1 ≤ ν ≤ p, max ∂ ν θ ϕ ρ , ∂ ν θ φ ρ ≤ β ν := 8 M ε ν!
Then the following estimates hold true

∀ 0 < ε < ε 0 , ∀ 0 ≤ ν ≤ p, Λ(ϕ) -Λ( φ) ρ-δ,ν ≤ ξ Q ν (ξ) 8 ε ϕ -φ ρ,ν , Γ ε (ϕ) -Γ ε ( φ) ρ-δ,ν+1 ≤ ξ Q ν (ξ) 8 ϕ -φ ρ,ν , where Q ν (ξ) ≤ 1 + 2ξ + 2 ν+1 √ e e ξ ν! and ξ = 16M ε δ .
Proof of Theorem 3.3. Again, owing to Remark 2.2, we limit ourselves to the case of stroboscopic averaging. By virtue of Proposition 3.2, the iterates Φ [k] , k = 0, . . . , n are p + 1 times differentiable and satisfy the assumptions of Lemma 3.4. By definition, we have for n ≥ 1 δ

[n] θ = Λ(Φ [n-1] ) θ -Λ(Φ [n] ) θ . Hence, using Lemma 3.4 with δ = R n+1 and successively ρ = R n , ρ = R n-1 , ..., ρ = R 1 yields with ξ = 16M ε(n+1) R : δ [n] R,p ≤ 1 8ε ξQ p (ξ) Φ [n] -Φ [n-1] Rn,p = 1 8ε ξQ p (ξ) Γ ε Φ [n-1] -Γ ε Φ [n-2] Rn,p ≤ 1 ε ξQ p (ξ) 8 2 Γ ε Φ [n-2] -Γ ε Φ [n-3] R n-1 ,p ≤ . . . ≤ 1 ε ξQ p (ξ) 8 n Φ [1] -Φ [0] R 1 ,p .
In addition, a direct estimation gives Φ [1] -

Φ [0] R 1 ,p ≤ 2ε M, so that δ [n] R,p ≤ 2M ξQ p (ξ) 8 n .
Now, (n + 1)ε ≤ ε 0 , so that ξ ≤ 2 and we have

Q p (ξ) ≤ Q p := 1 + 4 + 2 p+1
p!e 5/2 ≤ 5 + 13 2 p+1 p!

Uniform accuracy of numerical schemes

We are now in position to prove the main result of Section 3 which permits to apply a non-stiff standard integrator to the pullback and micro-macro formulations. We consider the numerical solution by a standard scheme of the pullback problem (3.5) and of the micro-macro problem (3.1).

Theorem 3.5. Suppose that Assumptions 2.1 and 3.1 are satisfied for a given p ≥ 1 and let ε 0 := R/(8M ). For problems (3.5) and (3.1) with 0 < (n + 1)ε ≤ ε 0 , consider an approximation (v k ) 0≤k≤N of (3.5) or (v k , w k ) 0≤k≤N of (3.1) at times t 0 = 0 < t 1 < . . . < t N = 1 obtained by a standard stable method of nonstiff order q, i.e. a method which exhibits order q of convergence when applied to (1.1) with ε = 1. Then, these approximations yield a uniformly accurate approximation of order r = min(p, q, n) of the solution of (1.1).

Precisely, we have for the pullback method,

max 0≤k≤N Φ [n] (v k ) -u ε (t k ) X ≤ C∆t r ,
and for the micro-macro method,

max 0≤k≤N Φ [n] (v k ) + w k -u ε (t k ) X ≤ C∆t r ,
where the constant C is independent of ε and the time grid size ∆t = max k (t k+1 -t k ) assumed small enough.

Proof. We first note that according to Proposition 3.2, Φ [n] is analytic and bounded on K R , hence it is Lipschitz continuous, so that all we have to prove is that the approximations v k and (v k , w k ) to equations (3.5) and (3.1) are of order r. As far as the pullback method is concerned, equation (3.5) can be written as (see

(3.3)) v = F [n] (v) + R [n] t/ε (v), v(0) = (Φ [n] 0 ) -1 (u 0 ),
where F [n] is Lipschitz and where R

[n]

θ is defined in (3.4). The result thus follows from the existence of the derivatives of R

[n] θ w.r.t. θ up to order p and the estimates

max 1≤ν+ξ≤p ∂ ν t ∂ ξ u (R [n] t/ε ) R ≤ Cε n-ν ,
which are a consequence of Theorem 3.3 and of the Leibniz formula.

Consider now the micro-macro system (3.1). By definition of δ

[n]
θ , it may be expressed as

v = F [n] (v), v(0) = (Φ [n] 0 ) -1 (u 0 ), ẇ = L [n] (t/ε, t, w) w + b [n] (t/ε, t), w(0) = 0, (3.10)
with

L [n] (θ, t, w) = 1 0 ∂ u f θ Φ [n] θ (v(t)) + µw dµ , b [n] (θ, t) = -δ [n] θ (v(t)),
and where we have written in (3.1)

f t/ε Φ [n] t/ε (v(t)) + w -f t/ε Φ [n] t/ε (v(t)) = 1 0 ∂ u f t/ε Φ [n]
t/ε (v(t)) + µw dµ w.

Here, both the operator function L [n] (θ, t, w) and the source term b [n] (θ, t) = -δ

L [n] := max 0≤ν+µ+ξ≤p sup (θ,t)∈T×[0,T ] ∂ ν θ ∂ µ t ∂ ξ w L [n] (θ, t, •) R and B [n] := max 0≤ν+µ≤p sup (θ,t)∈T×[0,T ] ∂ ν θ ∂ µ t b [n] (θ, t) X ≤ b [n] ε n ,
where b [n] is a constant independent of ε. Exploiting the fact that b [n] (•, t) = 0 in the application of the Gronwall lemma allows to assert that there exists M 0 such that

∀t ∈ [0, 1], w(t) X ≤ M 0 b [n] exp(L [n] t) ε n+1 .
Then, by differentiating (3.10) up to r = min(p, n) times and using again the Gronwall lemma, it can be verified that there exists a constant M r ∈ R * such that

∀ν = 0, . . . , r + 1, ∀t ∈ [0, 1], ∂ ν t w(t) X ≤ M r b [n] exp(L [n] t) ε n+1-ν .
As for the numerical approximations w k of w(t k ), the boundedness of L and the stability of the scheme allow to write

w k+1 X ≤ (1 + αL [n] ∆t) w k X + βb [n] ε n ∆t, k = 0, . . . , N -1,
where α and β are two constants depending on the scheme, so that ∀k = 0, . . . , N, w k X ≤ β b [n] exp(αL [n] ) ε n

When the numerical scheme advances the solution from t k to t k+1 , it uses the right-hand side of (3.1) evaluated at w k . The fact that w k is of order ε n implies that all derivatives of this right-hand side up to order r + 1, when evaluated at w k , are bounded in ε. This allows to assert that the scheme retains its usual order q, provided q ≤ r.

Construction of uniformly accurate integrators

In this section, we explain how our framework allows for the derivation of effective uniformly accurate schemes.

Standard averaging micro-macro methods

We now detail the procedure for the construction of a scheme with accuracy of order 1 uniformly with respect to ε ∈]0, 1]. The simplest change of variable can be obtained after a single iteration in (2.4), i.e. for n = 1, which yields

Φ [1] θ (w) = w + ε(g θ (w) -g (w)) (4.1)
where we denote

g θ (w) = θ 0 f τ (w) -f (w) dτ, (4.2)
while the associated vector field F [1] is defined by

F [1] (v) = f (v).
We can then apply any standard non-stiff integrator to the micro-macro system (3.1), for instance the simplest Euler method. To derive a uniformly accurate method of order n, we follow the same methodology and consider the change of variable Φ [n] and the vector field F [n] defined by the iterations

F [n] = f • Φ [n] , Φ [n+1] θ = id + ε θ 0 f τ • Φ [n] τ -∂ u Φ [n] τ F [n] dτ -ε θ 0 f τ • Φ [n] τ -∂ u Φ [n] τ F [n] dτ
and then apply a standard non-stiff integrator of order n to the micro-macro system (3.1).

Stroboscopic averaging pullback methods

In this section, we construct two changes of variable for the pullback formulation (3.5) with order one and two respectively, for the construction of efficient geometric schemes that are uniformly accurate w.r.t. ε. The proposed changes of variables are perturbations of the maps Φ [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF] θ and Φ [START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF] θ defined in (2.4) and chosen in order to preserve quadratic first integrals of the original system (1.1). Their construction relies on the implicit midpoint rule, known to preserve quadratic first integrals.

A first order change of variable

We observe that the change of variable

Φ [1] θ (v) = v + εg θ (v)
where g θ is given by (4.2), is 1-periodic with respect to θ, and it may sound attractive because it is completely explicit. However, considering such an explicit change of variable is not desirable because it destroys the preservation of quadratic first integrals, and the Hamiltonian structure if f τ is assumed Hamiltonian, i.e. the transformed pullback system (3.5) looses these geometric properties in general. Alternatively, we shall consider the following change of variable based on the implicit midpoint rule, Φ [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF] θ

(v) = v + εg θ v + Φ[1] θ (v) 2 . ( 4.3)
We note that Φ [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF] θ as well as differentials of Φ [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF] θ involved in (3.5) can be computed by fixed point iterations based on the following identities. The quantities involved in the right-hand side of (3.5) can be obtained as

∂ θ Φ[1] θ (v) = ε f θ -f v + Φ[1] θ (v) 2 + ε 2 ∂ u g θ v + Φ[1] θ (v) 2 ∂ θ Φ[1] θ (v), (4.4)
(∂ u Φ[1] θ (v)) -1 K = K - ε 2 ∂ u g θ v + Φ[1] θ (v) 2 (K + (∂ u Φ[1] θ (v)) -1 K) (4.5)
This permits us to propose an algorithm for computing the transformed vector field involved in (3.5) and defined as

Fε,θ (v) = ∂ u Φ[1] θ (v) -1 f θ • Φ[1] θ (v) - 1 ε ∂ θ Φ[1] θ (v) . (4.6)
The resulting algorithm for the computation of v = Fε,t/ε (v) at time t and point v as given by equation (3.1) is detailed below. In practice, it assumes that the function (θ, v) → f θ (v) and its directional derivative ∂ u f θ (v)P in the direction P are given and that their Fourier coefficients as periodic functions of θ can be computed cheaply, for instance by using a Fast Fourier Transform (FFT). We compute by fixed point iterations the approximations

P Φε θ (v), U v+ Φ[1] θ (v) 2 , Q 1 ε ∂ θ Φ[1] θ (v), S Fε,θ (v).
Algorithm 1 for computing the vector field Fε,θ (w) in (4.3),(4.6) for given input v, θ, ε.

Set

P [0] = v, S [0] = Q [0]
= 0, and for i = 0, 1, 2, . . . compute

U [i] = v + P [i] 2 P [i+1] = v + εg θ (U [i] ) Q [i+1] = f θ (U [i] ) -f (U [i] ) + ε 2 ∂ u g θ (U [i] )Q [i] K [i+1] = f θ (P [i+1] ) -Q [i+1] S [i+1] = K [i+1] - ε 2 ∂ u g θ (U [i] )(S [i] + K [i+1] )
The output is S [i] that converges to Fε,θ (v) as i grows to infinity.

A second order change of variable

Carrying one additional iteration of (2.4) we obtain

Φ [2] θ (v) = v + ε θ 0 f τ • Φ [1] τ (v) -(∂ u Φ [1] τ (v)) ∂ u Φ [1] (v) -1 f • Φ [1] (v) dτ.
Using the notation (4.2), we have Φ [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF] θ = id + εg θ and Φ [START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF] θ may be rewritten as

Φ [2] θ = id + εg θ + ε 2 θ 0 ∂ u f τ g τ -∂ u f g -(∂ u g τ -∂ u g ) f dτ + O(ε 3 ), (4.7)
where we have used

∂ u Φ [1] τ = id + ε τ 0 ∂ u (f s -f ) ds, ∂ u Φ [1] -1 = id -ε 1 0 τ 0 ∂ u (f s -f ) dsdτ + O(ε 2 ),
and

f • Φ [1] = f + ε 1 0 τ 0 ∂ u f τ (f s -f ) dsdτ + O(ε 2 ).
Following the idea introduced in Section 4.2.1 (in order to conserve exactly quadratic first integrals), we search for a change of variables satisfying

Φ[2] θ (v) = v + εh θ Φ[2] θ (v) + v 2 . (4.8) and satisfying Φ[2] θ (v) = Φ [2]
θ (v) + O(ε 3 ). Expansing (4.8) and in Taylor series with respect to ε,

Φ[2] θ (v) = v + εh θ (v) + ε 2 2 ∂ u h θ (v)h θ (v) + O(ε 3 ),
and comparing with (4.7), we deduce the vector field h θ given by

h θ = g θ + ε θ 0 ∂ u f τ g τ -∂ u f g -(∂ u g τ -∂ u g ) f dτ - ε 2 ∂ u g θ g θ . (4.9) Remark 4.1. Since Φ[2] θ = Φ [2]
θ + O(ε 3 ), the result of Theorem 3.5 with n = 2 remains true.

Theorem 4.2. Consider the pullback formulation (3.5) with change of variable Φ [START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF] θ (v) defined in (4.8). If the original system (1.1) has a quadratic first integral I(u(t)) = I(u(0)), then the pullback formulation (3.5) has the same quadratic first integral

I(v(t)) = I(v(0)). If f θ (u) = J -1 ∇ u H θ (u) in (1.
1) is a Hamiltonian vector field, then the pullback formulation (3.5) is again Hamiltonian.

Proof. The advantage of formulation (4.8)-(4.9) becomes apparent if one rewrites h θ as

h θ = g θ + ε 2 θ 0 [f τ + f , g τ ] -[f + f , g] dτ, (4.10)
where we expand

-∂ u g θ g θ = θ 0 (∂ u g τ f + ∂ u f g τ -∂ u g τ f τ -∂ u f τ g τ )dτ.
Here, [•, •] denotes the usual Lie-bracket of two vector fields, i.e.

[k

1 , k 2 ](v) = ∂ u k 1 (v)k 2 (v) -∂ u k 2 (v)k 1 (v).
As a matter of fact, the Lie-bracket operation is inherently geometric and F θ preserves first integrals as soon as f θ does, and it is Hamiltonian as soon as f θ is. Since the midpoint rule is known to preserve quadratic first integrals and to be a symplectic method for Hamiltonian vector fields, then Φ [START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF] θ preserves these properties.

We then have the following formulas for the computation of the right-hand side of equation (3.5):

Φε θ (v) = v + εh θ Φ[2] θ (v) + v 2 , (4.11) ∂ θ Φ[2] θ (v) = ε∂ θ h θ Φ[2] θ (v) + v 2 + ε 2 ∂ u h θ Φ[2] θ (v) + v 2 ∂ θ Φ[2] θ (v), (4.12) (∂ u Φ[2] θ (v)) -1 K = K -ε∂ u h θ Φ[2] θ (v) + v 2 K + (∂ u Φ[2] θ (v)) -1 K 2 . (4.13)
The resulting algorithm for the computation of ẇ at time t and point w as given by equation (3.1) is detailed below2 :

Algorithm 2 for computing the vector field Fε,θ (v) in (4.8), (4.6) for given input v, θ, ε.

Set

P [0] = v, S [0] = Q [0]
= 0, and for i = 0, 1, 2, . . . compute

U [i] = v + P [i] 2 P [i+1] = v + εh θ (U [i] ) Q [i+1] = ∂ θ h θ + ε 2 ∂ u h θ (U [i] )Q [i] K [i+1] = h θ (P [i+1] ) -Q [i+1] S [i+1] = K [i+1] - ε 2 ∂ u h θ (U [i] )(S [i] + K [i+1] ).
The output is S [i] that converges to Fε,θ (v) as i grows to infinity.

Remark 4.3. The above methodology is not restricted to the preservation of quadratic first integrals. Consider the special case where the vector field in (1.1) has the form

f θ (u) = S θ ∇I(u)
where S θ = -S T θ is a skew-symmetric matrix of size d × d depending periodically on the variable θ ∈ T. In this case, I : R d → R turns out to be a first integral and needs not to be quadratic. Then, in place of the implicit midpoint rule in (4.8), one can consider the averaged vector field method [START_REF] Quispel | A new class of energy-preserving numerical integration methods[END_REF] defined as

Φ[2] θ (v) = v + ε 1 0 h θ µ Φ[2] θ (v) + (1 -µ)v dµ, (4.14)
and the corresponding pullback formulation (3.5) has the same first integral (I(v(t)) = I(v(0))).

Numerical experiments

In this section, we present some numerical experiments that confirm our theoretical analysis and illustrates the efficiency of our strategy. We test our method on two different models, the Hénon-Heiles equations and the Klein-Gordon equation.

The numerical schemes

To present our numerical schemes, let us write the standard averaging micro-macro equations defined in Subsection 4.1 as well as the stroboscopic averaging pullback equation defined in Subsection 4.2 under the following generic form:

u = G [n] ε,t/ε (t, u), (5.1)
where the superscript n refers to the order of the change of variable Φ [n] used in the method.

Recall that the time derivatives up to order n of the vector field G

[n] ε,t/ε are uniformly bounded with respect to ε ∈]0, 1]. In our numerical experiments, we shall use the following numerical schemes that are of integral type. This enables to gain one order of regularity and to use the simpler change of variable Φ [n-1] for schemes of order n, instead of the change of variable Φ [n] normally required in Theorem 3.5. 3 -Order 2 Integral Runge-Kutta scheme (explicit):

u k+1 = S RK2 ∆t u k , with S RK2 ∆t u k = u k + t k+1 t k G [1] ε,t/ε t k+1/2 , u k + t k+1/2 t k G [1] ε,s/ε (t k , u k )ds dt,
where t k = k∆t. This scheme will be applied to the standard averaging micro-macro equations.

-Order 2 Integral midpoint scheme (implicit):

u k+1 = S midpoint ∆t u k , with S midpoint ∆t u k = u k + t k+1 t k G [1] ε,t/ε t k+1/2 , u k + u k+1 2 dt.
This implicit scheme can be simply implemented using fixed point iterations. It will be applied to the stroboscopic averaging pullback equation. Recall that G

[n] ε,θ is periodic w.r.t. θ, so in these two methods the integrals can be either analytically precomputed, or easily computed numerically by using FFT. Moreover, it can be shown that, thanks to the uniform boundedness of the first time derivative of G [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF] ε,t/ε , these two schemes are uniformly accurate (UA) of order 2.

-Order 3 extrapolation scheme : If S RK2 ∆t denotes the above integral RK2 scheme constructed with G

[2] ε,t/ε instead of G [1]
ε,t/ε , we construct a UA method S extrap3 ∆t of order 3 by implementing the Richardson extrapolation:

S extrap3 ∆t = 1 3 4S RK2 ∆t/2 • S RK2 ∆t/2 -S RK2 ∆t .
This scheme will be applied to the standard averaging micro-macro equations.

-Order 3 composition scheme : If S midpoint ∆t denotes the above integral midpoint scheme constructed with G

[2] ε,t/ε instead of G [1]
ε,t/ε , we construct a UA method S compos3 ∆t of order exactly 3 by composition, setting:

S compos3 ∆t = S midpoint α 3 ∆t • S midpoint α 2 ∆t • S midpoint α 1 ∆t
, where (for instance) α 2 = -1.8 and α 1 > α 3 are uniquely defined such that the order three conditions hold, α 1 + α 2 + α 3 = 1, α 3 1 + α 3 2 + α 3 3 = 0. We emphasise that this non-symmetric composition scheme is considered for illustrating the theory only, because efficient high-order symmetric composition methods could also be considered. This scheme will be applied to the stroboscopic averaging pullback equation.

-Order 4 extrapolation scheme : If S RK2 ∆t denotes the above integral RK2 scheme constructed with G

[3] ε,t/ε instead of G [1]
ε,t/ε , we construct a UA method S extrap4 ∆t of order 4 by implementing the Bulirsch and the Stoer extrapolations:

S extrap4 ∆t = 1 12 27S RK2 ∆t/3 • S RK2 ∆t/3 • S RK2 ∆t/3 -16S RK2 ∆t/2 • S RK2 ∆t/2 + S RK2 ∆t .
This scheme will be applied to the standard averaging micro-macro equations.

3 Indeed, it can be shown that the local error and hence the global error of order two of the integral schemes S RK2 ∆t and S midpoint ∆t applied to any system ẏ = f (t, y) involves only the norms of the second order derivative ÿ of the solution and the derivatives ∂tf, ∂yf, ∂yyf, ∂ytf of the vector field, and in contrast to a standard Runge-Kutta method, it does not depend on ... y , ∂ttf .

The Hénon-Heiles model -micro-macro method

In this subsection and the next subsection, we consider the Hénon-Heiles model [START_REF] Hénon | The applicability of the third integral of motion: Some numerical experiments[END_REF][START_REF] Hairer | Geometric numerical integration[END_REF], parametrized with ε. This is a Hamiltonian system with the following Hamiltonian energy:

H(q 1 , q 2 , p 1 , p 2 ) = p 2 1 2ε + p 2 2 2 + q 2 1 2ε + q 2 2 2 + q 2 1 q 2 - 1 2 q 3 2 .
When ε is small, the variables q 1 , p 1 are highly oscillatory. Let us write the associated filtered system, satisfied by the variable w(t) ∈ R4 defined by

               w 1 (t) = cos t ε q 1 (t) -sin t ε p 1 (t), w 2 (t) = q 2 (t), w 3 (t) = sin t ε q 1 (t) + cos t ε p 1 (t), w 4 (t) = p 2 (t).
This variable w(t) satisfies the system

dw dt (t) = f t ε , w(t) , (5.2)
with          f 1 (θ, w) = 2 sin θ (w 1 cos θ + w 3 sin θ) w 2 f 2 (θ, w) = w 4 f 3 (θ, w) = -2 cos θ (w 1 cos θ + w 3 sin θ) w 2 f 4 (θ, w) = -2 (w 1 cos θ + w 3 sin θ) 2 + w 2 2 -w 2 .
For the standard averaging micro-macro method, all the integrals in θ in our numerical schemes can be pre-computed analytically. We have implemented these computations with the software Maple.

Let us present our numerical results obtained by the standard averaging micro-macro method. The initial data is (0.12, 0.12, 0.12, 0.12). In Figure 1a, we have represented the error at T f inal = 1 as a function of ∆t, for different values of ε, for our UA scheme S RK2 ∆t of order 2, for our UA scheme S extrap3 ∆t of order 3 and for our UA scheme S extrap4 ∆t of order 4. The reference solution is computed with the Matlab ode45 function applied to the filtered system (5.2). The error is independent of ε as shown in Figure 1a where the curves for different values of ε are nearly identical. This confirms that the schemes are uniformly accurate respectively of order 2, 3 and 4, as predicted by Theorem 3.5.

The Hénon-Heiles model -pullback method

In this section, we present the results obtained with the pullback method with stroboscopic averaging for the Hénon-Heiles model. In Figure 1b, we have represented the error at time T f inal = 1 as a function of ∆t for different values of ε, for our UA scheme S midpoint ∆t of order 2 and for our UA scheme S compos3 ∆t of order 3. These figures confirm that our schemes are UA, as stated by Theorem 3.5. Note that for these schemes, the variable θ of the vector field F

[n] ε,θ is discretized with 32 points and FFT is used to compute all the integrals 4 . The initial data is (0.12, 0.12, 0.12, 0.12). Let us now check that our schemes constructed on the stroboscopic averaging pullback method have the expected geometric properties. In Figure 2a, we have represented the evolution of the error on the Hamiltonian H for ε = 1 over long times t ≤ 30000 for the scheme S midpoint ∆t with ∆t = 0.1 and ∆t = 0.2 and for the scheme S compos3 ∆t with ∆t = 0.2. In Figure 2b, we have represented the results of the same calculations with ε = 0.001. The initial data is (0, 0, p 1 , p 2 ) with p 1 = 2ε/12 sin( π 8 ), p 2 = 2/12 cos( π 8 ). We observe no drift for the Hamiltonian.

Let us now observe the so-called Poincaré cuts [START_REF] Hénon | The applicability of the third integral of motion: Some numerical experiments[END_REF][START_REF] Hairer | Geometric numerical integration[END_REF], i.e. the points (q 2 , p 2 ) reached when the solution crosses the hyperplane q 1 = 0 with p 1 > 0.

In Figure 3a, we have represented 5 orbits for ε = 1, 0 ≤ t ≤ 30000 and for the energy H = 1/12. The initial data are (0, 0, p 1 , p 2 ) with p 1 = √ 2εH sin(θ π 2 ), p 2 = √ 2H cos(θ π 2 ) with θ = 0.025, 0.25, 0.75, 0.95 and 1.25. Our UA scheme S midpoint ∆t is used with the timestep ∆t = 0.1 for t ≤ 30000.

In Figure 3b, we have represented 3 orbits for ε = 0.001, 0 ≤ t ≤ 30000 and for the energy H = 1/6. The initial data are (0, 0, p 1 , p 2 ) with p 1 = √ 2εH sin(θ π 2 ), p 2 = √ 2H cos(θ π 2 ) with θ = 0.1, 0.5 and 0.9. Our UA scheme S midpoint ∆t is used with the timestep ∆t = 0.1 for t ≤ 30000.

On these two figures, we observe that the Poincaré curves are well reproduced, which indicate that the "formal second invariant" is well preserved by our scheme.

5.4

The nonlinear Klein-Gordon equation in the nonrelativistic limit regime -micro-macro method

We consider the nonlinear Klein-Gordon (NKG) equation in the nonrelativistic limit regime, also studied numerically in [START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF][START_REF] Baumstark | Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting[END_REF][START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Faou | Asymptotic preserving schemes for the Klein-Gordon equation in the non-relativistic limit regime[END_REF]:

ε∂ tt u -∆u + 1 ε u + f (u) = 0, x ∈ R d , t > 0, (5.3)
with initial conditions given as

u(0, x) = φ(x), ∂ t u(0, x) = 1 ε γ(x), x ∈ R d .
(5.4)

The unknown is the function u(t, x) : R 1+d → C and the parameter ε > 0 is inversely proportional to the square of the speed of light. In our simulations, we take d = 1 and the nonlinearity is f (u) = 4|u| 2 u. The limit ε → 0 in equation (5.3), (5.4), referred to as the nonrelativistic limit, has been studied in [START_REF] Machihara | The nonrelativistic limit of the nonlinear Klein-Gordon equation[END_REF][START_REF] Machihara | Nonrelativistic limit in the energy space for nonlinear Klein-Gordon equations[END_REF][START_REF] Masmoudi | From nonlinear Klein-Gordon equation to a system of coupled nonlinear Schrödinger equations[END_REF].

It is convenient to rewrite (5.3) under the equivalent form of a first order system in time (see e.g. [START_REF] Masmoudi | From nonlinear Klein-Gordon equation to a system of coupled nonlinear Schrödinger equations[END_REF]). Setting

v + = u -iε(1 -ε∆) -1/2 ∂ t u, v -= u -iε(1 -ε∆) -1/2 ∂ t u, (5.5)
and denoting

f ± (v) = f 1 2 (v ± + v ∓ ) ,
we obtain that (5.3), (5.4) is equivalent to the following system on the unknown v = (v + , v -): Let us introduce the filtered unknowns

i∂ t v ± = - 1 ε (1 -ε∆) 1/2 v ± -(1 -ε∆) -1/2 f ± (v), (5.6) v(0, •) = (v + (0, •), v -(0, •)) = φ -i(1 -ε∆) -1/2 γ , φ -i(1 -ε∆) -1/2 γ . ( 5
v ± = e -i t ε √ 1-ε∆ v ± .
These quantities satisfy the following equation:

∂ t v ± = i(1 -ε∆) -1/2 e -i t ε √ 1-ε∆ f ± e i t ε √ 1-ε∆ v . (5.8)
Note that (5.8) is under the form

∂ t v ± (t) = f ± t ε , t, v(t) ,
if we set f ± (θ, t, u) = i(1 -ε∆) -1/2 e -iθ e -itAε f ± e iθ e itAε u (5.9)

with A ε = 1 ε √ 1 -ε∆ -1 .
The system has two natural invariants:

-the conserved charge (quadratic invariant)

Q = Im R d ε u u dx = 1 4 R d (1 -ε∆) 1/2 v + v + dx - 1 4 R d (1 -ε∆) 1/2 v -v -dx = 1 4 R d (1 -ε∆) 1/2 v + v + dx - 1 4 R d (1 -ε∆) 1/2 v -v -dx -the energy E = ε R d | u| 2 dx + R d |∇u| 2 dx + 1 ε R d |u| 2 dx + 2 R d |u| 4 dx.

The nonlinear Klein-Gordon equation in the nonrelativistic limit regime -pullback method

We now present the results of the pullback method with the schemes S midpoint ∆t and S compos3 ∆t . In Figure 4b, we have plotted the error as a function of ∆t for different values of ε. For this figure, we have discretized the θ variable with 64 points and the x variable with 128 points, the initial data is (5.10) and the final time of the simulation is T f inal = 0.25. These numerical experiments confirm that our schemes S midpoint ∆t and S compos3 ∆t are respectively UA of order 2 and 3.

In Figure 5a, we have represented the long time (t ≤ 1000) conservation of the quadratic invariant Q for our two schemes S midpoint ∆t and S compos3 ∆t and in Figure 5b we have plotted the evolution of the relative error on the energy. The maximal absolute value of the error on E is 10 -8 for the second order scheme and is 2 × 10 -10 for the third order scheme. The initial data, with a non trivial value of Q, is

φ(x) = (1 + i)(cos x + sin x), γ(x) = (1 -i/2) cos x + (1/2 + i) sin x.
The blue curves correspond to S midpoint ∆t with ∆t = 0.01 and the red curves correspond to S compos3 ∆t with ∆t = 0.01. In all the simulations, we have ε = 0.0001 and we have taken 64 points for the discretization of the x variable as well as for the θ variable. The quadratic invariant Q is an exact invariant (up to round-off errors) of the two schemes, and we observe no energy drift over long times.

Conclusion

We have shown that by transforming the original problem (1.1) through a change of variables inspired from the theory of averaging, it becomes possible to apply standard nonstiff numerical schemes (i.e. integrators which would properly solve equation (1.1) with ε = 1) for all values of ε ∈]0, 1] with a computational cost and accuracy both independent of ε ∈]0, 1]. This has been done in four different ways whose merits have been discussed in details. Roughly speaking, stroboscopic averaging used in conjunction with the pullback method leads to geometric uniformly accurate integrators while standard averaging associated with a micro-macro decomposition leads to comparatively cheap uniformly accurate methods (which are not geometric though). As far as efficiency is concerned, the most relevant question is to estimate the critical threshold value ε c below which our techniques become superior to standard schemes. Although such an ε c clearly always exists, it is delicate to estimate its value as it strongly depends on the problem itself and on the choices of implementation made in the various methods: for instance, we have exploited here the fact that both the change of variable and the averaged vector field can be computed at first and second orders exactly for the Henon-Heiles problem. Whenever this is possible, most of the additional formulae required by our techniques can be pre-computed, making the formers indeed very competitive even for large values of ε (i.e. ε c is relatively large). The situation is different when this option is precluded (e.g. for the Klein-Gordon equation) and in that case the value of ε c is undoubtedly smaller. However, we wish to emphasize that whatever its value might be, our technique remains of interest for variants of the model equation (1.1) where ε may undergo large variations.

θ ) iq -∂ ν θ ∂ u Φ [k] θ F [k]
where the inner sum runs over multi-indices i = (i 1 , . . . , i q ) ∈ N q such that i 1 + 2i 2 + . . . + qi q = q and where |i| = i 1 + . . . + i q and

C q i = q! i 1 !i 2 !2! i 2 . . . i q !q! iq . Owing to (ii), for u ∈ K R k+1 , Φ [k]
θ (u) takes its values in K R k+1 +rn/2 , a set on which ∂ j θ ∂ i u f θ is bounded (due to Assumption 3.1 and Cauchy estimates for analytic functions) as follows:

∂ j θ ∂ i u f R k+1 +rn/2 ≤ 1 (2R -(R k+1 + r n /2)) i M = M (k + 1 2 ) i r i n . Hence, if we denote α k,ν = ∂ ν θ Φ [k] R k /(εM ), for ν ≥ 1, then we have ∀k ≥ 1, ∀ν ≥ 1, α k+1,ν+1 ≤ 1 + 2εM α k,ν r n + ν q=1 ν q i C q i (εM ) |i| α i k r |i| n (k + 1 2 ) |i|
where we denote α i k = q j=1 α i j k,j and use that (owing to (ii) and a Cauchy estimate)

∂ u ∂ ν θ Φ [k] F [k] R k+1 ≤ 1 r n ∂ ν θ Φ [k] R k F [k] R k+1 ≤ 2M ∂ ν θ Φ [k] R k r n .
Finally, the assumption that (n + 1)ε ≤ ε 0 = R 8M leads to the inequality and an easy induction shows that the sequence (g k (1)) k≥1 is upper-bounded by 8. We conclude the proof by applying Cauchy estimates as follows:

∀1 ≤ k ≤ n + 1, ∀1 ≤ ν ≤ p + 1, α k,ν = d ν g k dξ ν (0) ≤ ν! sup |ξ|≤1 |g k (ξ)| 1 ν ≤ ν! g k (1) ≤ 8 ν!
Proof of Lemma 3.4. Under the assumptions (3.9), it has been shown in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF] that:

(a) ∂ u ϕ ρ-δ ≤ 3 2 , the mapping ∂ u ϕ -1 is analytic on K ρ-δ and obeys ∂ u ϕ -1 ρ-δ ≤ 2;

(b) the mappings Λ(ϕ) θ (u) and Γ ε (ϕ) θ (u) are analytic on K ρ-δ for all θ ∈ T, and

∀0 < ε < ε 0 , Λ(ϕ) ρ-δ ≤ 4 M and Γ ε (ϕ) -Id ρ-δ ≤ 4 M ε.
Now, on the one hand, for a multi-index i = (i 1 , . . . , i q ) ∈ N q , the algebraic identity

q j=1 X i j j - q j=1 Y i j j = q j=1 (X j -Y j ) j-1 l=1 Y i l l q l=j+1 X i l l i j r=1 X i j -r j Y r-1 j
leads, for any (∆ 1 , . . . , ∆ q ) ∈ X q C and any ( ∆1 , . . . , ∆q ) ∈ X q C , to

∂ |i| u f • φ (∆ i 1 1 , . . . , ∆ iq q ) -∂ |i| u f • φ ( ∆i 1 1 , . . . , ∆iq q ) ρ-δ ≤ ∂ |i| u f ρ-δ/2 q j=1 ∆ j -∆j j-1 l=1 ∆ l i l q l=j+1 ∆l i l i j r=1 ∆ j i j -r ∆j r-1 . (A.2)
On the other hand, the Faà di Bruno's formula with 1 ≤ ν < p gives

∂ ν θ (f θ • ϕ θ -f θ • φθ ) = ∂ ν θ f θ • ϕ θ -∂ ν θ f θ • φθ + ν q=1 ν q i C q i ∂ ν-q θ ∂ |i| u f θ • ϕ θ (∂ θ ϕ θ ) i 1 , . . . , (∂ q θ ϕ θ ) iq - ν q=1 ν q i C q i ∂ ν-q θ ∂ |i| u f θ • φθ (∂ θ φθ ) i 1 , . . . , (∂ q θ φθ ) iq
with the notations of Theorem 3.2. Hence, using the decomposition in (A.2) we have

∂ ν θ (f θ • ϕ θ -f θ • φθ ) ρ-δ ≤ ∂ u ∂ ν θ f θ ρ-δ/2 ϕ -φ ρ-δ + ν q=1 ν q i C q i β i ∂ ν-q θ ∂ |i|+1 u f ρ-δ/2 ϕ -φ ρ-δ + ν q=1 ν q i C q i β i ∂ ν-q θ ∂ |i| u f ρ-δ/2 q j=1 i j β j ∂ j θ (ϕ -φ) ρ-δ
Upon using Cauchy estimates, we obtain

∂ ν θ (f • ϕ -f • φ) ρ-δ ≤ 2M δ ϕ -φ ρ-δ + 2M δ ϕ -φ ρ-δ ν q=1 ν q i Cq i 16M ε δ |i| + 1 8ε ν q=1 ν q i Cq i 16M ε δ |i| q j=1 i j j! ∂ j θ (ϕ -φ) ρ-δ ≤ 2M δ P ν 16M ε δ ϕ -φ ρ,ν
where we denote Cq i = C q i q j=1 j! i

∂ ν θ (∂ u ϕ F -∂ u φ F ) ρ-δ ≤ ∂ ν θ ∂ u ϕ ρ-δ F -F ρ-δ + ∂ ν θ ∂ u ϕ -∂ ν θ ∂ u φ ρ-δ F ρ-δ ≤ 8M ε δ F -F ρ-δ + 2M δ ∂ ν θ (ϕ -φ) ρ
Denoting A θ (u) = ∂ u ϕ θ (u) and Ãθ (u) = ∂ u φθ (u), we have

F -F ρ-δ ≤ A -1 ρ-δ f • ϕ -f • φ ρ-δ + A -1 -Ã -1 ρ-δ f • φ ρ-δ ≤ 8M δ ϕ -φ ρ
where we have used the bound in (a) A -1 ρ-δ ≤ 2, Ã -1 ρ-δ ≤ 2 and the inequality

A -1 -Ã -1 ρ-δ ≤ A -1 ρ-δ Ã -1 ρ-δ A -Ã ρ-δ ≤ 4 δ ϕ -φ ρ .
Finally, we obtain

∂ ν θ (∂ u ϕ F -∂ u φ F ) ρ-δ ≤ 64M 2 ε δ 2 ϕ -φ ρ + 2M δ ∂ ν θ ϕ -∂ ν θ φ ρ . (A.3)
Collecting all terms, we finally have

∂ ν θ (Λ(ϕ) -Λ( φ)) ρ-δ ≤ 2M δ Q ν 16M ε δ ϕ -φ ρ,ν
with Q ν (ξ) = 1 + 2ξ + P ν (ξ). By a Cauchy estimate it holds that ∀ξ ≥ 0, P ν (ξ) ≤ 2 ν+1 √ e e ξ ν! and Q ν (0) ≤ Q ν (ξ) ≤ 2 ν+1 √ e e ξ ν! + 2ξ + 1.

The corresponding bound for Γ ε is then obtained straightforwardly.

θ

  (v(t)) are of class C p w.r.t. θ ∈ T, of class C ∞ w.r.t. t ∈ [0, T ] and analytic w.r.t. w on K R with bounds

Figure 1 :

 1 Figure 1: Hénon-Heiles model. Error as a function of ∆t for ε = 2 -k , k ∈ {0, 1, 2, • • • , 9}. Fig. 1a (micro-macro method): UA scheme S RK2 ∆t of order 2 (top curves), UA scheme S extrap3 ∆t

Figure 2 :

 2 Figure 2: Hénon-Heiles model, pullback method. Long time evolution of the error on the Hamiltonian. UA scheme S midpoint ∆t of order 2 (blue: ∆t = 0.2; red: ∆t = 0.1). UA scheme S compos3 ∆t

  ε = 0.001, H = 1/6.

Figure 3 :

 3 Figure 3: Poincaré cuts for the Hénon-Heiles model. Pullback method with S midpoint ∆t

Figure 5 :

 5 Figure 5: NKG model with ε = 10 -4 , pullback method. Long-time evolution of the relative error on the quadratic invariant and on the energy. Blue: UA scheme S midpoint ∆t with ∆t = 0.01. Red: UA scheme S compos3 ∆t

1 ) 0 G

 10 We now introduce the generating functions for ξ ∈ C,g k (ξ) = ν≥1 α k,ν ξ ν ν! , k = 1, . . .A direct estimation gives α 1,1 ≤ 2 and α 1,ν ≤ 1 for ν ≥ 2, so that 0 ≤ g 1 (ξ) ≤ ξ + e ξ -1 for all ξ ∈ R + and from (A.1), we obtain (owing to the Leibniz and the Faà di Bruno formulae)G k (ξ) = exp ξ + g k (ξ) 8k + 4so that, using α k+1,1 ≤ 4, we have for ξ ∈ R + g k+1 (ξ) k (s)ds.Noticing that sup 0≤ξ≤1 g k (ξ) = g k (1), the previous inequality leads to g k+1[START_REF] Bao | Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF]

.

  Considering now the second part of ∂ν θ (Λ(ϕ) θ -Λ( φ) θ ) and denoting F = ∂ u ϕ -1 f • ϕ and F = ∂ u φ -1 f • φ , we have for ν ≥ 1

Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland, Gilles.Vilmart@unige.ch

Note, however, that this series is convergent for a linear right-hand side in (1.1) with bounded operator and for all sufficiently small ε.

To facilitate the convergence of the fixed point iterations for large ε, for all ε one can multiply g θ in (4.3) or h θ in (4.11) by a damping term such as exp(-ε 2 ) = 1 + O(ε 2 ) which does not affect the uniform accuracy of order two.

This discretization turns out to be exact for the Hénon-Heiles problem, due to the low degree of the involved trigonometric polynomials.

Acknowledgements. This work was partially supported by the Swiss National Science Foundation, grants No: 200020 144313/1 and 200021 162404 and by the ANR project Moonrise ANR-14-CE23-0007-01.

In the simulations of this subsection, we take the following initial data:

(5.10)

For Figure 4a, the final time of the simulations is T f inal = 0.25. The computational domain in x is [0, 2π]. For the numerical evaluations of the function f in our scheme, a spectral method is used in the x variable and the fast Fourier transform is used in the practical implementation. For these tests, the reference solution is computed by our third order UA method S extrap3 ∆t applied to the standard averaging micro-macro method with 256 grid points in x, 128 grid points in θ and ∆t = T f inal /4096. We use the H 1 relative error of a given numerical scheme which we define as

where u num (t f inal , •) is the approximated solution obtained by the considered numerical scheme, at the final time T f inal of the simulation. Let us present the results for the standard averaging micro-macro method. In this case, the analytic computation of the integrals in θ involves a too heavy system of equations and it is more convenient to compute these integrals numerically at each iteration of the scheme, by using fast Fourier transform techniques. In Figure 4a, we have represented the error at T f inal = 0.25 as a function of ∆t, for different values of ε, for our UA scheme S RK2 ∆t of order 2, for our UA scheme S extrap3 ∆t of order 3 and for our UA scheme S extrap4 ∆t of order 4. These curves confirm that the schemes are UA respectively of order 2, 3 and 4. The numerical parameters are the following: we have taken 128 discretization points in the x variable and 64 discretization points in θ for the quadrature methods.

Appendix

This section is devoted to the proof of Proposition 3.2 and Theorem 3.3. We first recall the two following basic results from [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF]: Lemma A.1 (See Castella, Chartier, Murua and Méhats [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF]). Let 0 < δ < ρ ≤ 2R. Assume that the function (θ, u) ∈ T × K ρ → ϕ θ (u) ∈ X C is analytic, and that ϕ θ is a near-identity mapping, in the sense that

Then ∂ u ϕ ρ-δ ≤ 3 2 , the mapping ∂ u ϕ -1 is well-defined and analytic on K ρ-δ with ∂ u ϕ -1 ρ-δ ≤ 2, and the mappings (θ, u) ∈ T×K ρ-δ → Λ(ϕ) θ (u) and (θ, u) ∈ T×K ρ-δ → Γ ε (ϕ) θ (u) are well-defined and analytic. In addition, the following bounds hold for all ε ≥ 0,

Lemma A.1 shows that, starting from a function (θ, u) ∈ T×K 2R → ϕ θ (u) ∈ X C , we can consider iterates (Γ ε ) k (ϕ) θ at the cost of a gradual thinning of their domains of analyticity. The following contraction property holds.

Lemma A.2 (See Castella, Chartier, Murua and Méhats [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF]). Let 0 < δ < ρ ≤ 2R and consider two periodic, near-identity mappings (θ, u) ∈ T×K ρ → ϕ θ (u) and (θ, u) ∈ T×K ρ → φθ (u), analytic on K ρ and satisfying

Then the following estimates hold for all ε ≥ 0,

Proof of Proposition 3.2. We assume that all Φ [k] are defined by (2.4) in the stroboscopic case (all functions G [k] are taken null). The following properties were already proved in [START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF]:

(i) the maps Φ [k] and F [k] := ∂ u Φ [k] -1 f • Φ [k] are well-defined and analytic on K R k for all k = 0, . . . , n + 1;

(ii) the estimates


Natural Interaction for Object Hand-Over Mamoun Gharbi, Séverin Lemaignan, Jim Mainprice, Rachid Alami CNRS-LAAS, 7 av. du Colonel Roche, F-31077 Toulouse, France Université de Toulouse, UPS, INSA, INP, ISAE, LAAS, F-31077 Toulouse, France Email: surname.name@laas.fr

Abstract-The video presents in a didactic way several abilities and algorithms required to achieve interactive "pick and give" tasks in a human environment. Communication between the human and the robot relies on unconstrained verbal dialogue, the robot uses multi-modal perception to track the human and its environment, and implements real-time 3D motion planning algorithms to achieve collision-free and human-aware interactive manipulation.

We tackle here the issue of natural and legible behaviour for a robot intended to live in domestic environment, amongst human peers. Natural because we aim at unconstrained spoken English mixed with gestures and perspective taking to communicate with the robot ; legible because we apply human-aware motion planning algorithms that ensure a smooth human-robot interaction by accounting for implicit social rules which appear in placement or approach strategies.

Symbolic facts and beliefs management motion plan requests

Human-aware symbolic task planning shared plans

Geometric & Temporal Reasoning

A. Modelling of the environment

The video clip gives a overview of the techniques used online by the robot to build a model of its environment. We rely first on a geometric model, fed by Kinect-like sensors for human tracking and 2D barcodes for object identification and localization.

From this geometric model, we build in real-time a symbolic model, stored as an ontology [START_REF] Lemaignan | ORO, a knowledge management platform for cognitive architectures in robotics[END_REF]. This models holds spatial relations between objects and agents, along with a set of affordances (like reachability, visibility, etc.) [START_REF] Lemaignan | Grounding the interaction: Anchoring situated discourse in everyday human-robot interaction[END_REF].

B. Natural language processing

The video also introduces recent work on natural dialogue grounding [START_REF] Lemaignan | Grounding the interaction: Anchoring situated discourse in everyday human-robot interaction[END_REF]. The user verbal input is first parsed, then semantically grounded, in tight interaction with the symbolic model. It allows for multi-modal and perspective-aware communication: affordances (like visibility) and human gestures that are recognized (like pointing) and stored in the ontology are used during the grounding process.

C. Motion planning in human close vicinity

While planning the robot motions, it is important to account for interaction constraints such as those defined by the proxemics theory. Motion planning techniques such as the ones in [START_REF] Mainprice | Planning Human-aware motions using a sampling-based costmap planner[END_REF] enable to account for such constraints using a costbased representation where efficient algorithms incorporating features from stochastic optimization are applied to produce low-cost paths.

D. Sharing effort in handover tasks and proactive robot placement

The planner demonstrated in this video has been introduced in [START_REF] Mainprice | Sharing effort in planning human-robot handover tasks[END_REF], it also accounts for interaction constraints. It consists of sampling handover configurations encoding for the human and robot postures ensuring that they are accessible to them both. The configurations are evaluated regarding an interaction cost and the one that minimizes this cost is returned by the planner. This algorithm enables the robot to proactively decide where to give the object. It can produce plans either for a seated or standing human or through a window as shown in the video. We have introduced a notion of mobility to tune the amount of effort asked to the human.

Fig. 1 .

 1 Fig.1. Overview of the robot deliberative architecture[START_REF] Alami | When the robot considers the human[END_REF]. Coloured modules with plain borders are the specific focus of the video.

ACKNOWLEDGMENT

This work has been supported by EU FP7 "SAPHARI" under grant agreement no. ICT-287513.


Léopold Ouairy
email: leopold.ouairy@inria.fr

Hélène Le Bouder

Jean-Louis Lanet
email: jean-louis.lanet@inria.fr

Hélène Le-Bouder
email: helene.le-bouder@imt-atlantique.fr

Normalization of Java source codes

Keywords: applet security, identifier renaming, Chucky-ng, Java Card

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

An applet is a program embedded on Java Cards. Such applets contain secrets. Attackers can steal such secrets by exploiting wrong implementations of specifications. To do so, they try a fuzzing attack on the applet. This attack consists in sending many different messages to the card via the Application Protocol Data Unit (APDU) communication buffer. This APDU enables terminals to communicate with cards and vice versa. With this attack, attackers can detect and then exploit the wrong state machine of an applet. This exploit can lead the attackers to illegally obtain access to secret resources. This is possible if an applet omit to implement a check in a critical state, whereas the specification precises that such check shall be implemented. We want to protect applets by detecting those missing-checks in their source codes. Since applets are written in a subset of the Java language, we base our work on a version of Chucky-ng [START_REF] Maier | Assisted discovery of vulnerabilities in source code by analyzing program slices 1[END_REF][START_REF] Yamaguchi | Chucky: Exposing missing checks in source code for vulnerability discovery 1[END_REF] adapted for Java source files. This version is entitled ChuckyJava and it directly works on source codes without the need to modify or annotate them. Such tool aims at detecting the missing-checks in applets by outputing an anomaly score between 1.00 and -1.00. If this anomaly score is 1.00, the applet forgets to perform a test. On the contrary, if the anomaly score is -1.00, then the applet is the only one to perform a test. Since there are false positives, an analyst is required in order to interpret ChuckyJava's results.

We want to improve ChuckyJava by reducing the false positive rate and the entries of its report to analyze. Since it works with identifier names of constants, variables, parameters and methods, our idea is to normalize every identifier of an applet set. We suggest an algorithm to achieve this task and we need to develop the tool IdentNormalize. Such modifications can be adapted to other programming languages. However, our study focuses on Java Card applets. To benchmark our idea, we manually apply the modifications we recommend in the applet set. Then, we perform an analysis with ChuckyJava to verify if our idea improves the tool efficiency. Because our modifications are concluding, we suggest another idea to increase the ChuckyJava efficiency.

The context and state of the art are shown in section 2. Our contribution is presented in section 3. The results are explained in section 4. Our future work is shown in section 5. To finish, we conclude on our improvement in section 6.

Context and State of the art 2.1 How Chucky-ng operates

Chucky-ng's objective is to detect wrong implementations of source codes. To do so, it first parses source code within a project. Then the analyst chooses a function to process. Chucky-ng groups similar functions to the one selected. Finally, it compares this group by giving an anomaly score for every expression existing in the function set. This anomaly score's range is [-1.00, 1.00]. An anomaly score of -1.00 informs the analyst that the function executes an expression where none of the other do. On the opposite, a function with an anomaly score of 1.00 does not execute an expression whereas the others do. These expression are sensible to the name of identifiers. Two expressions not using the same identifier are flagged as not similar, even if the functions are semantically identical. ChuckyJava [START_REF] Ouairy | Protection des systemes face aux attaques par fuzzing 2 5. Pietig, A.: Functional Specification of the OpenPGP application on ISO Smart Card Operating Systems 6[END_REF] is based on Chucky-ng but parses Java files instead of C/C++. Moreover, we use a layer above ChuckyJava: FetchVuln. This last one automates ChuckyJava by requesting an analysis for every methods of the file. It reports the vulnerable methods, according to ChuckyJava in a single output file.

Tool requirements

This section highlights the necessity for us to design a normalization tool. It has to normalize the variable, constant, parameter and methods names.

identifier names A specification of a Java Card applet precises the commands to implement and the expected behavior of an applet once it receives the command. Since ChuckyJava works partly with identifiers, it would work optimized if every applet implementing a specification uses the exact same identifiers for constants, fields and parameter names. However, our work shows two applets which implement the OpenPGP 2.0.1 specification can use really different identifiers. For example, the command PUT-_DATA is named SET_DATA in another applet. Even if both names have the same meaning, they are completely different for ChuckyJava. It would warn the user for an anomaly wich is a false positive. Our tool has to rename such constants as either PUT_DATA or SET_DATA. method names Our tool has to normalize method names too. In ChuckyJava, there is a step where the tool gathers similar methods, using an unsupervised machine learning algorithm. The authors of Chucky-ng, the base of ChuckyJava, rely on a similarity of method names during its neighborhood discovery step. This step focuses on gathering the most similar methods (neighbors) within applet source codes, by using information retrieval techniques. Added to the cosine distance metric used, a distance reward bonus is added if the Jaccard distance between two method names are similar. The Jaccard distance measures the dissimilarity between two sets or two strings in our case. If the distance is closer to zero, then the strings are identical. On the opposite, if the distance is one, the strings have nothing in common. In our case, if they are not similar, a distance penalty is added to the method. As an example we have three methods putData. One of them forget to use a particular object of a meaningful type. ChuckyJava would not gather it as similar functions of putData, or even include a getData function. It would generate a lot of false positive. To prevent this, we want to force the gather of functions using the same names in order to compare methods implementing the same commands.

State of the art

There are different techniques or tools that aims at normalizing methods names in source codes. INFOX [START_REF] Zhou | Identifying features in forks 3[END_REF] is a tool that aims at tracking forks for a project on github. It then extracts features and clusters code snippets in order to highlight similar code.

Another approach is based on word study within methods [START_REF] Sridhara | Identifying word relations in software: A comparative study of semantic similarity tools 3[END_REF]. The authors have based their work on different methods to analyze the similarity of their semantics. It can be achieved by using a WordNet object structure to organize the words. This structure organizes the words in a hierarchical way. Based on this, one can use a similarity measure such as the Lowest Common Subsumer (LCS) to determine if two words are related. Another technique uses the definition of words to determine if a pair of word is similar. This last one is called glossed based technique. However, the authors conclude by precising that none of the tools seem to perform well and require improvements.

The paper [START_REF] Kuhn | Semantic clustering: Identifying topics in source code 3[END_REF] uses Latent Semantic Indexing (LSI). This technique groups term frequencies of documents in matrices. The tool works by gathering similar vocabulary words from comments and identifier names. The technique is to firstly retrieve information within source codes. Then, it is able to cluster the code snippets and label them. The analyst gets a comprehensive view of the code snippet topics, without needing to search in it a particular information. This plugin for the Microsoft's Phoenix framework [START_REF] Tairas | Phoenix-based clone detection using suffix trees 3[END_REF] enables an analyst to detect clone detection. It is based on techniques to detect biological identical ADN sequences to search for perfect matches. It works by using suffix trees (AST) and uses identifier names. It creates AST of functions and compares them. There are other plagiarism detection techniques available on the internet to discover clones of code.

The requirements of our tool The tool which clusters topics by extracting words from comments and identifier could fit to our prototype. In a future work, we have to adapt it for working with methods as the base unit. This could be possible since applets in general are production ones. In other words, their source code should be clean and commented in most of the cases. The Phoenix's plugin could match our requierements. It works only with perfect matches. However, the authors claim that for near exact matches, we should use the Smith-Waterman [START_REF] Greenan | Method-level code clone detection on transformed abstract syntax trees using sequence matching algorithms 4[END_REF] algorithm. It uses matrices to determine the similarity of code snippets.

Contribution

We have considered the Latent Semantic Analysis (LSA) approach for the method normalization step. However, such method is based on similarity and its results may not be always be true. For the rest, we have made the choice to rename identifiers by deducing the ones which fit almost every applet of the set. The tool works in four steps we develop:

1. remove of unused variables, 2. constant and field names normalization, 3. method names normalization, 4. parameter names normalization.

Remove unused variables

This step aims at removing unused variables. We work mainly on production applets. Some of them do not have any. However, it eliminates some useless output of ChuckyJava.

2. Constant and field names normalization IdentNormalizer tries to rename common constants within applet source codes. In most cases in specification implementation, the short constants defined at the beginning of the applets is the value of a specific command. This value may be unique. For example, the command CHANGE_REFERENCE _DATA has a value of 0x24 for its instruction byte. As it is declared once in an applet, we have a first step of comparison between constants of this value. We randomly choose the name for this constant by selecting a name in an applet and assigning it to the others. For constants using the same value at least twice in the source code should have different semantics, we prefer to not modify the names. We could replace all constants with their real values. However, Chucky-ng normalizes values as $NUM. For example, a case 0x20 and a case 0x30 would be transformed as case $NUM. By doing this, we would lose a lot of information. The tool normalizes the fields names too. For example, in the OpenPGP specification, there are three different passwords or card holder verifications which are OwnerPin objects. Some of them can be renamed as ch1 or pw1. The trick here is to rename the identifiers by deducing their relations based on the object type OwnerPin. This technique renames identifiers regardless of any specification. be flagged as similar. If two names are similar but semantically different such as setData and getData, then they could be flagged as similar, even if they do not have the same objective. We want to compare functions with the same objective together. One way to solve our problem is to use information retrieval techniques (IR). It is possible to rely on latent semantic indexing or analysis (LSI, LSA) on source code and comments to extract the main topic of methods. Kuhn et al. [START_REF] Kuhn | Semantic clustering: Identifying topics in source code 3[END_REF] propose a technique entitled Semantic Clustering. It is based on both LSI and clustering. The information retrieval step extracts identifier names and comments from the source code. Then, a clustering method gathers code snippets of similar topics. With this technique, it is possible to group similar functions together in sets of topics in order to rename those methods. However, since this technique is based on statistics, it could contain classification errors.

Method names normalization

Parameter names normalization

At this point, the method names are normalized. This is an essential condition since we want to normalize parameter names. Even if function can use a different number of parameters, we rename the one that are used by every similar method names together. This is based on the type name. However, two functions can use twice the same type of object as parameter. In this case, within the function, it is possible to gather in a set the methods used by the first parameter and in another set the ones using the second parameter. By comparing those sets, it should be possible to determine if two are similar and should be renamed the same way. If the sets are identical, we suggest to leave the parameters names unchanged.

Summary of the operation Our tool

IdentNormalize can be summed up as shown in Algorithm 1. It performs with a complexity of O(2n) with n, the total number of files for all applets. However, the step of method names normalization may be more complex than O(2n). The complexity of the tool depend on the LSI method's complexity. Appendix A shows two different code snippets as S 1 . Since they are semantically identical, our tool produces either two identical appendix A.1 or two appendix A.2 as S 1 .

Only the class name remains unchanged.

Results

Description of the experiments

We have divided the experiments in two different sets. The first is experimented in a controlled environment. It contains one-file applets which perform a match-on-card algorithm. They have a similar structure and nearly the same number of methods. On the contrary, the second set is composed of four applets coming from different github sources. All of them implement the OpenPGP v.2.0.1 specification [5]. However, their structure can be different and may interpret commands of the APDU in different ways.

For each experiment, we present the result of ChuckyJava before and after the normalization method. Moreover, we compare only the results for the anomaly scores of -1.00 and 1.00. The reason behind this decision is because those numbers are the priority to focus on for an analyst.

Controlled environment experimentation

This set gathers eleven applets which were designed by students. Every applet is written within a single file. All of these applets are constructed with globally the same program structure. It focuses on the importance of the identifier normalization. From this set, we have divided the experiment in two experiments. We first perform a ChuckyJava analysis on the applets without any modification. Then, we apply our renaming algorithm on the set. We operate a second analysis on this new set. We present the impact of such a renaming tool on ChuckyJava by comparing before and after our renaming step as shown in Table 1. This same table shows only the results we have obtained for the anomaly scores -1.00 and 1.00. However, Fig 1 . shows the values we get before and after renaming the identifiers for a few anomaly scores. We can see that the number of entries for the anomaly score -1.00 decreases on our custom set. However, it increases for the 1.00 anomaly score. After analysing at the results, it is because we have now the entries as anomalies for the tests/cases of switches. Such entries are lines which ChuckyJava reports to the analyst. Each line precises the anomaly score, the expression associated to this score and the location of the expression in its Java file. It adds a benefit for using a normalization tool.

Second experimentation

We have gathered four applets implementing the OpenPGP v2.0.1 specification. None of these applets come from the same author and have both different identifier names and program structures. Our applet set is composed of MyPGPid4 , OpenPGPApplet5 , JCOpenPGP6 , Gpg7 . Fig 2 . shows the result we have obtained from the first set to the last one which contains all the modifications. Table 2 summarizes the original number of results and with the improvements, for both -1.00 and 1.00 anomaly scores.

As for the controlled experiment, we can see that the number of entries for the anomaly score decreases for the anomaly score -1.00 but increase for the 1.00. To summarize number of entries for -1.00 number of entries for 1.00 original set 3769 255 improved set 3443 351

Table 2: Number of entries to analyse betwen the original set and the new one the result, we lose roughly 6% of entries to analyse. It corresponds to a decrease of 230 lines. This result may seem to be a bit low. However, this is mainly because Chucky-Java is really sensible to the program structure. Indeed, it compares identifiers between applets to establish the anomaly score. However, since our applets are made from different authors and from different quality, we have concluded that the gain of 6% reflects the lower limit of this approach. For example, industrials may have structure pattern for the applets they implement, reducing the impact of such structure on the results of ChuckyJava.

Fig. 2: The original set is on the left while the bar on the right represent the impact of IdentNormalizer, in a production applet set

Future work

This section presents the future work which completes the identifier normalization problem. We have noticed three different structural problems. The first one concerns an assignation problem. The second one focuses on the order of instruction within the source code. The last one is about the optional functionalities which an implementation offers.

Assignation type problem

We are aware of several assignation type problem. For example, the vast majority of the applets creates three different variables pw1,pw2 and pw3 for the three different Own-erPin(s). However, one applet has been designed to store them in an array of OwnerPin. When a method of OwnerPin is called, the identifier name of the array is used instead of one of the pw. Because of this, ChuckyJava generates false positives. It reports this applet with an anomaly score of 1.00 since it is the only one to use such an identifier. This creates one additional output every time the program uses this array.

Instructions call order

We are aware of a structural problem. Listing 1.1 shows a code snippet representing an applet using the APDU buffer as a field. Every time the applet receives a command, it first retrieves the buffer before calling the corresponding method. On the opposite, 1.2 shows another applet retrieving the APDU buffer within each command method. However, even if those applets are semantically identical, they can be reported as anomaly.

Because ChuckyJava analyzes methods as the base unit, it cannot see that the call to the getBuffer is made in another location in the source code. It adds additional false positive lines to the ChuckyJava's report.

Optional features

One last structure problem generating false positive is when a specification suggests optional features. As an example, the OpenPGP v2.0.1 specification allows developers to use optional data objects (DO). If one applet decides to not implement them while the others do it, then the applet is flagged with many additional entries with an anomaly score of 1.00. On the opposite, if one option is implemented in only one applet, Chuck-yJava generates entries with anomaly scores of -1.00. However, both the results are false positives since the implementation of the option is non mandatory.

Conclusion

We have developed IdentNormalizer which normalizes applet methods and identifiers names. It helps an analyst who uses ChuckyJava by reducing the number of entries between 6% to 15%. We have shown that the more the structure of an applet set is similar, the more IdentNormalizer is efficient. Our future work focuses mainly creating the tool that normalizes the identifiers. Moreover, we want to improve the efficiency of ChuckyJava by creating another tool which normalizes applet structures. As like IdentNormalizer, this new tool could be executed before performing an analysis with ChuckyJava.

Fig. 1 :

 1 Fig. 1: The original set is on the left while the bar on the right represent the impact of IdentNormalizer in a controlled applet set

Listing 1 . 1 :

 11 Snippet 1 b y t e [ ] b u f f e r = new b y t e [ APDU_MAX_SIZE ] ; p r o t e c t e d v o i d p r o c e s s (APDU apdu ) throws I S O E x c e p t i o n { b u f f e r = apdu . g e t B u f f e r ( ) ; s w i t c h ( b u f f e r [ ISO7816 . OFFSET_INS ] ) { c a s e VERIFY : v e r i f y ( apdu ) ; break ; [ . . . ] } p r i v a t e v o i d v e r i f y (APDU apdu ) { / / Some code , u s i n g w i t h o u t c a l l i n g t h e APDU } Listing 1.2: Snippet 2 p r o t e c t e d v o i d p r o c e s s (APDU apdu ) throws I S O E x c e p t i o n { b y t e [ ] b u f f e r = apdu . g e t B u f f e r ( ) ; s w i t c h ( b u f f e r [ ISO7816 . OFFSET_INS ] ) { c a s e VERIFY : v e r i f y ( apdu ) ; break ; [ . . . ] } p r i v a t e v o i d v e r i f y (APDU apdu ) { b y t e [ ] b u f f e r = apdu . g e t B u f f e r ( ) ; / / Some code , u s i n g t h e b u f f e r }

  It is not trivial to find which methods are similar. One way to gather similar methods is to analyze their name with the Jaccard distance. However, if two functions are semantically identical but use different names, they would not Algorithm 1: How IdentNormalizer operates

	Data: source codes of applets S1
	Result: The applet modified S 1
	forall s1 ∈ S1 do
	remove of unused variables;
	constant and field names normalization;
	end
	V1 := methods normalization names list;
	forall s1 ∈ S 1 do
	method names normalization based on V1;
	parameter names normalization;
	end

Table 1 :

 1 Number of entries to analyse betwen the original set and the new one

		number of entries for -1.00	number of entries for 1.00
	original set	909	121
	improved set	726	151

https://github.com/CSSHL/MyPGPid

https://github.com/Yubico/ykneo-openpgp

https://sourceforge.net/projects/jcopenpgp/

https://github.com/FluffyKaon/OpenPGP-Card