Philippe Chartier
email: philippe.chartier@irisa.fr

Ander Murua
email: ander.murua@ehu.es

Jesus Maria Sanz-Serna
email: sanzsern@mac.uva.es

Ernst Hairer

Arieh P Iserles

Chartier

Konputazio Zientziak

A A Saila

Higher-Order Averaging, Formal Series and Numerical Integration I: B-series

Keywords: Averaging, High-order stroboscopic averaging, Highly oscillatory problems, Hamiltonian problems, Multiscale numerical methods, Numerical integrators, Formal series, B-series, Trees, Fermi-Pasta-Ulam problem, Adiabatic invariants, Inverted Kapitsa's pendulum Mathematics Subject Classification (2000) 34C29, 65L06

We show how B-series may be used to derive in a systematic way the analytical expressions of the high-order stroboscopic averaged equations that approximate the slow dynamics of highly oscillatory systems. For first-order systems we give explicitly the form of the averaged systems with O(ǫ j ) errors, j = 1, 2, 3( 2 πǫ denotes the period of the fast oscillations). For second-order systems with large O(ǫ -1 ) forces, we give the explicit form of the averaged systems with O(ǫ j ) errors, j = 1, 2. A variant of the Fermi-Pasta-Ulam model and the inverted Kapitsa pendulum are used as illustrations. For the former it is shown that our approach establishes the adiabatic invariance of the oscillatory energy. Finally we use B-series to analyze multiscale numerical integrators that implement the method of averaging. We construct

Introduction

The aim of this series of papers is to relate the method of averaging to the formal series expansions that are nowadays used as a powerful tool in the analysis of numerical integrators of time-dependent problems. The present work is restricted to B-series and systems with a single fast frequency; Part II will deal with other types of formal series and with quasiperiodic problems.

The method of averaging [START_REF] Fedorenko | Derivation and justification of equations in slow time (the stroboscopic method)[END_REF][START_REF] Perko | Higher order averaging and related methods for perturbed periodic and quasi-periodic systems[END_REF][START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF] has a long history that goes back to the work in celestial mechanics of Gauss and Laplace [START_REF] Arnold | Geometrical Methods in the Theory of Ordinary Differential Equations[END_REF]. The aim is to study the longtime behavior of highly oscillatory systems by constructing an averaged system that approximately captures that behavior and ignores the details of the oscillations with small periods.

In 1974 Hairer and Wanner [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF] introduced the concept of B-series. They associated with each numerical integrator (within a very broad class) its B-series, a formal series in powers of the step-length, and showed that the properties of the integrator are easily studied by manipulating the corresponding B-series. The importance of the notion of B-series and other similar formal series [START_REF] Murua | Formal series and numerical integrators. Part I: systems of ODEs and symplectic integrators[END_REF] has grown steadily in recent years [START_REF] Hairer | Geometric Numerical Integration[END_REF], in particular since the discovery of their relevance in connection to symplectic integration [START_REF] Calvo | Canonical B-series[END_REF] and modified equations [START_REF] Hairer | Backward error analysis of numerical integrators and symplectic methods[END_REF] (see also the recent contribution [START_REF] Chartier | Algebraic structures of B-series[END_REF]). Even though B-series are formal series, they may lead to rigorous error bounds when suitably truncated [START_REF] Hairer | Geometric Numerical Integration[END_REF].

In Sects. 2-4 of this paper we show how B-series may be used to derive in a systematic way the analytical expressions of high-order averaged equations. This is probably the first example of the application of B-series outside the field of numerical ordinary differential equations. Section 2 describes the general idea and Sects. 3 and 4 provide the details for large classes of first-order and second-order differential systems respectively. For first-order systems we give explicitly the form of the averaged systems with O(ǫ j ) errors, j = 1, 2, 3(2 πǫ denotes the period of the fast oscillations). When the original oscillatory problem is Hamiltonian, so are the averaged systems and explicit expressions for the corresponding Hamiltonian functions are also provided. For second-order systems with large, O(ǫ -1 ) forces, we give the explicit form of the averaged systems with O(ǫ j ) errors, j = 1, 2. A family of Hamiltonian problems that includes a variant of the Fermi-Pasta-Ulam model is employed to illustrate the material in Sect. [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF]. We show that our methodology, in addition to determining averaged equations of high order, yields as a byproduct the adiabatic invariance of the oscillatory energy. The material in Sect. 4 is exemplified in the case of Kapitsa's inverted pendulum. The final Sect. 5 deals with the analysis of numerical multiscale methods such as those introduced by E and Engquist [START_REF]Analysis of the heterogeneous multiscale method for ordinary differential equations[END_REF][START_REF] Engquist | The heterogeneous multiscale methods[END_REF]. While the scope of application of those methods widely exceeds the oscillatory problems considered in the present article, they are of particular relevance to our research because, when applied to the classes of systems studied here, they may be considered as numerical implementations of the idea of averaging. We shall construct multiscale methods that are able to approximate not only the simplest, lowest-order averaged equation but also its high-order counterparts.

We end this introduction by pointing out that there are of course analogies between the B-series approach considered here and both the modulated Fourier expansion methodology pioneered by Hairer and Lubich [START_REF] Hairer | Long-time energy conservation of numerical methods for oscillatory differential equations[END_REF] and the WKB and Magnus expansion techniques [START_REF] Iserles | On the global error of discretization methods for highly-oscillatory ordinary differential equations[END_REF].

A Modified Equation Approach to Averaging

This section relates well-known ideas from the theory of averaging [START_REF] Fedorenko | Derivation and justification of equations in slow time (the stroboscopic method)[END_REF][START_REF] Perko | Higher order averaging and related methods for perturbed periodic and quasi-periodic systems[END_REF][START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF]t o modified equations as used in the analysis of numerical integrators [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Sanz-Serna | Numerical Hamiltonian Problems[END_REF].

We are concerned with initial value problems for differential systems of the form

d dt y = f y, t ǫ ; ǫ , (1)
where y is a D-dimensional real vector, ǫ is a small parameter and the indefinitely differentiable function f is assumed to depend 2π -periodically on the variable t/ǫ. Our interest is in situations where, as ǫ → 0, the solutions or some of their derivatives with respect to t become unbounded. Such is the case, for instance, when f and its partial derivatives are O(1) and the partial derivative of f with respect to t/ǫ does not vanish (so that (1) is effectively non-autonomous); then differentiation in [START_REF] Ariel | A multiscale method for highly oscillatory ordinary differential equations with resonance[END_REF]sho ws that (d 2 /dt 2 )y = O(1/ǫ). But even in cases where (1) is autonomous, so that there is no effective dependence of f on the fast time t/ǫ, the derivatives of y will become unbounded as ǫ → 0iff itself behaves like a negative power of ǫ.

For the analysis it is useful to rewrite the system in terms of the scaled (nondimensional) time τ = t/ǫ:

d dτ y = ǫf (y , τ ; ǫ). (2)
If we denote by ϕ τ 0 ,τ ;ǫ : R D → R D the solution operator of (2), so that

y(τ) = ϕ τ 0 ,τ ;ǫ (y 0 )
is the solution that satisfies the initial condition y(τ 0 ) = y 0 , then a simple but important observation is that the map Ψ τ 0 ;ǫ = ϕ τ 0 ,τ 0 +2π;ǫ depends on τ 0 in a 2π -periodic manner; this follows from the fact that both ϕ τ 0 ,τ ;ǫ (y 0 ) and ϕ τ 0 +2π,τ+2π;ǫ (y 0 ) satisfy the initial value problem d dτ y(τ) = ǫf (y(τ ), τ ; ǫ) = ǫf (y(τ ), τ + 2π; ǫ), y(τ 0 ) = y 0 .

From this observation, it follows that, at the stroboscopic times τ n = τ 0 + 2πn, n = 0, ±1, ±2,...,

y(τ n ) = ϕ τ 0 ,τ n ;ǫ (y 0 ) = ϕ τ n-1 ,τ n ;ǫ ϕ τ 0 ,τ n-1 ;ǫ (y 0 ) = ϕ τ 0 ,τ 0 +2π;ǫ ϕ τ 0 ,τ n-1 ;ǫ (y 0 )
and, hence, we arrive at the fundamental formula:

y(τ n ) = (Ψ τ 0 ;ǫ ) n (y 0 ), n = 0, ±1, ±2,.... (3)
In Sects. 3 and 4 we shall describe general situations where the exact solution of (2) with initial value y(τ 0 ) = y 0 , sampled at τ n = τ 0 + 2πn, admits a formal expansion of the form

y(τ n ) = y 0 + ∞ j =1 ǫ j j l=1 (τ n -τ 0 ) l G j,l (y 0 ), (4)
with suitable indefinitely differentiable maps G j,l : R D → R D independent of ǫ.In particular, this implies that

Ψ τ 0 ;ǫ (y 0 ) = y 0 + ∞ j =1 ǫ j j l=1 (2π) l G j,l (y 0 ),
and thus Ψ τ 0 ;ǫ is a smooth near-to-identity map. Standard backward error analysis [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Sanz-Serna | Numerical Hamiltonian Problems[END_REF] then shows the existence of an autonomous system (the modified system of Ψ τ 0 ;ǫ )

d dτ Y = F(Y; ǫ) = ǫF 1 (Y ) + ǫ 2 F 2 (Y ) + ǫ 3 F 3 (Y ) +••• (5) or d dt Y = 1 ǫ F(Y; ǫ) = F 1 (Y ) + ǫF 2 (Y ) + ǫ 2 F 3 (Y ) +••• (6)
(F and the F j depend on τ 0 , but this has not been incorporated into the notation) whose (formal) solutions satisfy that Y(τ n ) = Ψ τ 0 ;ǫ (Y (τ n-1 )) for n = 0, ±1, ±2,... so that

Y(τ n ) = (Ψ τ 0 ;ǫ ) n (Y 0 )n = 0, ±1, ±2,.... (7)
We conclude from ( 3) and ( 7) that, if one chooses Y(τ 0 ) = y(τ 0 ), then Y(τ) exactly coincides with y(τ) at the stroboscopic times τ n = τ 0 + 2πn. In this way it is possible in principle to find y(τ n ) by solving the system (5)or [START_REF] Calvo | Carrying an inverted pendulum on a bumpy road[END_REF], where all t-derivatives of Y remain bounded as ǫ → 0. When τ does not coincide with one of the stroboscopic times, we note that obviously

y(τ) = (ϕ τ n ,τ ;ǫ • Φ τ n -τ ;ǫ )Y (τ ),
where τ n is the largest stroboscopic time ≤τ and Φ •;ǫ denotes the flow of [START_REF] Calvo | Instabilities and inaccuracies in the integration of highly oscillatory problems[END_REF]. In this way, y is 'enslaved' to Y through the mapping ϕ τ n ,τ ;ǫ • Φ τ n -τ ;ǫ whose dependence on τ is easily seen to be 2π -periodic.

It is well known that the series (5) does not converge in general, and in order to get rigorous results one has to consider truncated versions

d dτ Y = ǫF 1 (Y ) + ǫ 2 F 2 (Y ) + ǫ 3 F 3 (Y ) +•••+ǫ J F J (Y ), (8)
whose solutions satisfy that Y(τ n ) -Ψ τ 0 ;ǫ (Y (τ n-1 )) = O(ǫ J +1 ).IfY solves [START_REF] Calvo | A stroboscopic numerical method for highly oscillatory problems[END_REF] with Y(τ 0 ) = y(τ 0 ), then Y(τ n ) and y(τ n ) differ by an O(ǫ J ) amount, where the constant implied in the O notation is uniform as the stroboscopic time τ n ranges in an interval τ 0 ≤ τ n ≤ τ n + T/ǫ, with T = O(1) as ǫ → 0. 1 The process of obtaining the autonomous system (8) from the original system (2) is referred to in the averaging literature [START_REF] Fedorenko | Derivation and justification of equations in slow time (the stroboscopic method)[END_REF][START_REF] Perko | Higher order averaging and related methods for perturbed periodic and quasi-periodic systems[END_REF][START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF] as high-order stroboscopic averaging. In this paper we shall show how a number of techniques and results currently used in the analysis of numerical integrators may be applied to the task of constructing explicitly the functions F j that define the averaged equations [START_REF] Calvo | A stroboscopic numerical method for highly oscillatory problems[END_REF].

Let us briefly discuss how to find the modified system (5) corresponding to Ψ τ 0 ;ǫ . There are several techniques that may be used to accomplish this task; here we shall apply Theorem 1 in [START_REF] Murua | Formal series and numerical integrators. Part I: systems of ODEs and symplectic integrators[END_REF]. It follows from ( 4) that the formal solution Y(τ) of the autonomous system [START_REF] Calvo | Instabilities and inaccuracies in the integration of highly oscillatory problems[END_REF] with initial value Y(τ 0 ) = y 0 can be written as

Y(τ)= y 0 + ∞ j =1 ǫ j j l=1 (τ -τ 0 ) l G j,l (y 0 ). (9)
(This is proved as follows. The difference (τ ) between the formal solution Y and the right-hand side of (9) vanishes when τ coincides with one of the stroboscopic times τ n . Therefore the polynomial P(τ) with degree ≤N that interpolates (τ ) at τ 0 ,...,τ N , vanishes identically. On the other hand the interpolation error -P = is easily seen to be O(ǫ N +1 ) and, since N is arbitrary, ≡ 0.) The vector field F is then retrieved by differentiating the solution flow:

2 d dτ Y(τ) τ =τ 0 = F(y 0 ; ǫ).
The following proposition lists some useful properties.

1 T h es i z eo fY(τ)y(τ) at non-stroboscopic times depends of course on the behavior of ϕ τ n ,τ ;ǫ .I f accurate approximations of y(τ) are required, they may be obtained by first finding Y(τ n ),whereτ n is the largest stroboscopic time ≤τ , and then integrating the original oscillatory system from τ n to τ . 2 This procedure is of course classical: old differential equations textbooks used to show that any given family of sufficiently differentiable functions Q = Ξ(τ,C) with values in R D and depending on a vector parameter C in R D was the 'general solution' of the differential equation obtained by eliminating the parameters from the equations Q = Ξ ,d Q/dτ = ∂Ξ/∂τ. In general the resulting differential equation is non-autonomous. If it is known beforehand-as is the case here-that the differential equation will turn out to be autonomous, then it is sufficient to consider Q = Ξ ,dQ/dτ = ∂Ξ/∂τ at τ = 0. In Sect. 4 we shall use the corresponding idea for second-order systems; these have as 'general solution' families Proposition 1

Q = Ξ(τ,C 1 ,C 2 ) with C 1 , C 2 in R D ,
(1) If the real function I(y,τ)depends 2π -periodically on τ and is conserved by all solutions of the original system (2) (I (y(τ ), τ ) = I(y(τ 0 ), τ 0 )), then I(Y,τ 0 ) is a (formal) conserved quantity for the averaged system (5). ( 2) If (2) is Hamiltonian, then so is [START_REF] Calvo | Instabilities and inaccuracies in the integration of highly oscillatory problems[END_REF]. More precisely, there is a formal Hamiltonian

H(Y ; ǫ) = ǫH 1 (Y ) + ǫ 2 H 2 (Y ) + ǫ 3 H 3 (Y ) +••• such that F(Y; ǫ) = J -1 ∇H(Y ; ǫ) and F j (Y ) = J -1 ∇H j (Y ) for each j( J is the canonical symplectic matrix). (3) If (2) is autonomous, then (5) is independent of τ 0 . (4) If (2)
is autonomous and possesses a first integral I(y), then I(Y) is a formal first integral of (5). ( 5) If (2) is autonomous and Hamiltonian with Hamiltonian function H , then H is a formal first integral of [START_REF] Calvo | Instabilities and inaccuracies in the integration of highly oscillatory problems[END_REF]. Furthermore, the Hamiltonian H of (5) is a formal first integral of (2).

Proof For (1) (which obviously implies (4)) note that S(τ) = I(Y(τ),τ 0 ) -I(Y(τ 0 ), τ 0 ) vanishes when τ coincides with one of the stroboscopic times τ n . Then the interpolation argument used to prove (9) implies S ≡ 0. If the original system is Hamiltonian (possibly non-autonomous), then Ψ τ ;ǫ is a symplectic map and therefore its modified equation (5) will be Hamiltonian [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Sanz-Serna | Numerical Hamiltonian Problems[END_REF]. This proves [START_REF] Arnold | Geometrical Methods in the Theory of Ordinary Differential Equations[END_REF].

The property (3) is true because for an autonomous system ϕ τ 0 ,τ ;ǫ only depends on the difference ττ 0 and Ψ τ 0 ;ǫ = ϕ τ 0 ,τ 0 +2π;ǫ is independent of τ 0 .

Under the hypotheses of ( 5), H is a first integral of the original system and, by (4), also of its averaged counterpart. But then H and H are in involution (their Poisson bracket vanishes) [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF] and therefore H is a formal invariant of the original [START_REF] Arnold | Geometrical Methods in the Theory of Ordinary Differential Equations[END_REF].

Let us finish this section with a remark. If the initial value τ 0 of τ is changed to τ ′ 0 , then, as noted before, the averaged system (5) changes. The relation

ϕ τ ′ 0 ,τ ′ 0 +2π;ǫ = ϕ τ 0 +2π,τ ′ 0 +2π;ǫ • ϕ τ 0 ,τ 0 +2π;ǫ • ϕ τ ′ 0 ,τ 0 ;ǫ , i.e. Ψ τ ′ 0 ;ǫ = (ϕ τ ′ 0 ,τ 0 ;ǫ ) -1 • Ψ τ 0 ;ǫ • ϕ τ ′ 0 ,τ 0 ;ǫ ,
shows that the new Ψ τ ′ 0 ;ǫ is conjugate to the old Ψ τ 0 ;ǫ by means of the change of variables ϕ τ ′ 0 ,τ 0 ;ǫ . Therefore the different averaged equations arising from different values of τ 0 are also conjugate to each other.

First-Order Differential Equations

Framework

Throughout this section we assume that the function f in (1) possesses an expansion in powers of ǫ of the form

f(y,τ; ǫ) = 1 ǫ ∞ j =1 ǫ j f j (y, τ ). (10)
If

f j (y, τ ) = ∞ k=-∞ exp(ikτ )f j,k (y) (11)
is the Fourier expansion of f j (y, τ ), then the Fourier coefficients f j,k (y) are in general complex vectors, but, in order to have a real system, we assume that, for each j and k, f j,k ≡ f * j,-k ( * denotes complex conjugate). The system (1) is supplemented by the initial condition y(0) = y 0 , where y 0 is a given vector of O(1) magnitude. Note that there is no loss of generality in assuming that the initial value of t is 0; the general case may be reduced to this by shifting t and redefining f . In terms of the scaled time τ = t/ǫ, the differential system (2) being solved is then

d dτ y = ǫf (y , τ ; ǫ) = ∞ j =1 ǫ j f j (y, τ ). (12)

Trees

As is customary in the analysis of numerical methods for ordinary differential equations, the solution y of (12) will be expanded in an appropriate B-series: a formal series whose terms are indexed by (rooted) trees. In this subsection we describe the trees that will be required in the expansion of y which will be carried out in the next. It is well known, see e.g. [START_REF] Hairer | Solving Ordinary Differential Equations I, Nonstiff Problems[END_REF], that in the standard analysis of numerical integrators for systems dy/dτ = f(y), the vertices of the trees correspond to f and its derivatives. In view of the structure (10)- [START_REF] Engquist | The heterogeneous multiscale methods[END_REF] of the right-hand side of [START_REF] Engquist | Heterogeneous multiscale methods: a review[END_REF], the vertices of the trees to be used here correspond to the functions f j,k and their derivatives; to keep track of this correspondence, each vertex possesses here a label, i.e. a pair of indices (j, k), j = 1, 2,..., k = 0, ±1, ±2,....NowthesetU of all trees may be defined recursively by the following two rules: (i) for each label (j, k), the corresponding tree with one vertex, jk00, belongs to U , (ii) if u 1 ,..., u n ∈ U , then, the result

u =[u 1 ,...,u n ] j,k (13)
of grafting their roots to a new root with label (j, k) belongs to U . The expansion of y will be graded according to powers of ǫ (see ( 14) below) and in this connection we introduce the notion of weight. The weight of a vertex with label (j, k) is defined to be j and the weight |u| of u ∈ U is the sum of the weights of its vertices.

The Expansion of y

The solution y of ( 12) possesses an expansion

y(τ) = y 0 + u∈U ǫ |u| α u (τ ) σ u F u (y 0 ) = y 0 + j =1 ǫ j |u|=j α u (τ ) σ u F u (y 0 ). (14)
Here F u , σ u and α u are, respectively, the elementary differential,thesymmetry and the elementary coefficient of the tree u. These will be described presently.

The elementary differential F u is, for each u ∈ U , a vector-valued function of a vector argument y, constructed in terms of the Fourier coefficients f j,k . Recursively, F ( j,k00) = f j,k and, for u in [START_REF] Engquist | Heterogeneous multiscale methods for stiff ordinary differential equations[END_REF],

F u (y) = f (n) j,k (y) F u 1 (y), . . . , F u n (y) , where f (n) j,k is the nth-order Fréchet derivative of f j,k .
The integer σ u counts the number of symmetries of the tree u ∈ U and may be defined recursively as follows. For trees with one vertex, σ( j,k00) = 1 and, for u in [START_REF] Engquist | Heterogeneous multiscale methods for stiff ordinary differential equations[END_REF],

σ u = r 1 !•••r m !σ u 1 •••σ u m ,
where u μ , μ = 1,...,m, are the pairwise distinct u ν , ν = 1,...,n, and r μ counts the number of times that u μ features among the u ν .

For u ∈ U , the elementary coefficient α u is a complex-valued function of the real variable τ . 3 Rewriting the differential equation ( 12) as an integral equation

y(τ) = y 0 + ǫ τ 0 f y(τ ′ ), τ ′ ; ǫ dτ ′
and expressing in both sides y by means of its expansion [START_REF] Fedorenko | Derivation and justification of equations in slow time (the stroboscopic method)[END_REF], leads to a Picard-like procedure for the recursive computation of the elementary coefficients. For the tree j,k00,

α u (τ ) = τ 0 exp(ikτ ′ ) dτ ′
and, for u in [START_REF] Engquist | Heterogeneous multiscale methods for stiff ordinary differential equations[END_REF],

α u (τ ) = τ 0 exp(ikτ ′ )α u 1 (τ ′ ) •••α u n (τ ′ ) dτ ′ . ( 15
)
Two important remarks follow. First, note that the elementary coefficients in the expansion of y are universal in the sense that they are independent of the particular f j,k in the differential equation. This is one of the main advantages of the B-series approach in numerical differential equations: with B-series there is a clear separation between, on the one hand, quantities depending on the particular differential equation being integrated but not on the specific integration method4 and, on the other hand, method-dependent quantities independent of the differential equation. The second observation is that the elementary coefficients only depend on the second component (wave number) of the label (j, k).

Computing the Elementary Coefficients

It does not seem possible to obtain a closed-form expression for the coefficients α u recursively defined in [START_REF] Hairer | Backward error analysis of numerical integrators and symplectic methods[END_REF]. 5 Nevertheless the use of this recursion may be systematized by considering the functions φ r,k (τ ) = τ r exp(ikτ ), r = 0, 1, 2,..., k = 0, ±1, ±2,... and the linear space Φ spanned by them, i.e. the set of all complex-valued functions ψ of a real variable that may be written as linear combinations

ψ(τ) = r,k a r,k φ r,k (τ ),
where the a r,k 's are complex constants.

Clearly, for each ψ ∈ Φ, the integral

I(ψ)(τ ) = τ 0 ψ(τ ′ ) dτ ′
will be given by the corresponding linear combination of the functions κ r,k = I(φ r,k ). These are in turn members of the space Φ, because, obviously

κ r,0 (τ ) = τ r+1 r + 1 , (16)
and, for k = 0, integration by parts leads to the recursion

κ r,k = ir k κ r-1,k - i k φ r,k ,r = 1, 2,.... ( 17
)
For future reference we list here the first few κ r,k with k = 0

κ 0,k = i k - i k exp(ikτ ), (18)
κ 1,k =- 1 k 2 + 1 k 2 exp(ikτ ) - i k τ exp(ikτ ), (19)
κ 2,k =- 2i k 3 + 2i k 3 exp(ikτ ) + 2 k 2 τ exp(ikτ ) - i k τ 2 exp(ikτ ). ( 20
) Since φ r,k φ r ′ ,k ′ = φ r+r ′ ,k+k ′ , (21)
the space Φ is also closed under multiplication (i.e. ψ 1 ψ 2 ∈ Φ,i fψ 1 ∈ Φ and ψ 2 ∈ Φ) and, by induction, we conclude from (15) that, for each u ∈ U , the elementary coefficient α u belongs to Φ and accordingly there exist complex constants a u,r,k such that

α u (τ ) = r,k a u,r,k φ r,k (τ ) = r,k a u,r,k τ r exp(ikτ ). (22)
With the help of ( 16), ( 17) and ( 21), it is easy to implement a recursive procedure to compute the elementary coefficients α u . Once the constants a u ν ,r,k corresponding to the trees in the right-hand side of ( 15) are known, we write, with the help of ( 21), the integrand as a linear combination of the φ r,k ,

exp(ikτ ′ )α u 1 (τ ′ ) •••α u n (τ ′ ) = r,k b u,r,k φ r,k (τ ),
and obtain

α u (τ ) = r,k b u,r,k κ r,k (τ ). (23)
Then the constants a u , r,k are found by expressing the functions κ r,k as linear combinations of the φ r,k (cf. ( 18)-( 20)).

After computing the elementary coefficients of all trees of weight j = 1, 2, we find the first terms in [START_REF] Fedorenko | Derivation and justification of equations in slow time (the stroboscopic method)[END_REF]:

y(τ) = y 0 + ǫ k κ 0,k (τ )f 1,k (y 0 ) + ǫ 2 k κ 0,k (τ )f 2,k (y 0 ) + ǫ 2 k,ℓ c k,ℓ (τ )f ′ 1,k (y 0 )f 1,ℓ (y 0 ) + O ǫ 3 , (24)
where

c k,ℓ (τ ) = i ℓ κ 0,k (τ ) - i ℓ κ 0,k+ℓ (τ ), ℓ = 0, and
c k,0 (τ ) = κ 1,k (τ ),
so that, from ( 16) and ( 18)-( 19),

c 0,0 (τ ) = τ 2 2 , c 0,ℓ (τ ) = 1 ℓ 2 - 1 ℓ 2 exp(iℓτ ) + i ℓ τ, ℓ = 0, c k,0 (τ ) =- 1 k 2 + 1 k 2 exp(ikτ ) - i k τ exp(ikτ ), k = 0, c k,-k (τ ) = 1 k 2 - 1 k 2 exp(ikτ ) + i k τ, k = 0, c k,ℓ (τ ) =- 1 k(k + ℓ) + 1 kℓ exp(ikτ ) - 1 ℓ(k + ℓ) exp i(k + ℓ)τ , k,ℓ,k+ ℓ = 0.

The Expansion of the Averaged Solution

In view of ( 22), the functions (polynomials)

ᾱu (τ ) = r,k a u,r,k τ r (25)
interpolate the values of α u (τ ) when τ is an integer multiple of 2π . Therefore, the formula

Y(τ)= y 0 + j =1 ǫ j |u|=j ᾱu (τ ) σ u F u (y 0 ), (26)
provides the B-series expansion of a function Y(τ) that interpolates the solution y(τ) with expansion ( 14) and is smooth, in the sense that all derivatives (d k /dt k )Y with respect to the original time t remain bounded as ǫ → 0 (cf. the bounded derivative principle [START_REF] Kreiss | Manifolds of slow solutions for highly oscillatory problems[END_REF]).

The leading terms of the B-series for y computed in [START_REF] Li | Deciding the nature of the coarse equation through microscopic simulations: the baby-bathwater scheme[END_REF] thus yield

Y(τ) = y 0 + ǫτf 1,0 (y 0 ) + ǫ 2 τ 2 2 f ′ 1,0 (y 0 )f 1,0 (y 0 ) + ǫ 2 τf 2,0 (y 0 ) + ǫ 2 τ k =0 i k f ′ 1,0 (y 0 )f 1,k (y 0 ) -f ′ 1,k (y 0 )f 1,0 (y 0 ) + ǫ 2 τ k =0 i k f ′ 1,k (y 0 )f 1,-k (y 0 ) + O ǫ 3 . (27
Y τ =0 = ǫf 1,0 (y 0 ) + ǫ 2 f 2,0 (y 0 ) + ǫ 2 k =0 i k f ′ 1,0 (y 0 )f 1,k (y 0 ) -f ′ 1,k (y 0 )f 1,0 (y 0 ) + ǫ 2 k =0 i k f ′ 1,k (y 0 )f 1,-k (y 0 ) + O ǫ 3
and, since Y(0) = y 0 , we conclude that the lowest-order (J = 1) averaged equation is

d dτ Y = ǫf 1,0 (Y )
(a result certainly expected) and that at the next, J = 2, order

d dτ Y = ǫf 1,0 (Y ) + ǫ 2 f 2,0 (Y ) + ǫ 2 k =0 i k f ′ 1,0 (Y )f 1,k (Y ) -f ′ 1,k (Y )f 1,0 (Y ) + ǫ 2 k =0 i k f ′ 1,k (Y )f 1,-k (Y ). (28)
Averaged equations of higher orders may in principle be found following the same methodology used for the cases J = 1, 2. The only difficulty stems from the need to determine the elementary coefficients α u corresponding to trees of weights |u|=3, 4,... that are required to compute, via [START_REF] Murua | Formal series and numerical integrators. Part I: systems of ODEs and symplectic integrators[END_REF], the coefficients ᾱu in the expansion of Y . To conclude this subsection we find explicitly the averaged equation with remainder O(ǫ 4 ). This is probably as far as one can reasonably go for the general format [START_REF] Engquist | Heterogeneous multiscale methods: a review[END_REF], but, of course, in particular cases it is possible to explicitly construct averaged equations with remainders of higher order.

Our task may be simplified by observing that, since we aim at computing the terms with weight j = 1, 2, 3 in the expansion of dY/dτ at τ = 0, it is not necessary to gain full knowledge of the polynomials ᾱu , |u|=1, 2, 3, but only of the coefficient of the first power τ 1 in such polynomials. Now each α u , |u|=1, 2, 3, is of the form [START_REF] Levi | Geometry and physics of averaging with applications[END_REF], with the summation index r restricted to the values 0, 1, 2 and the formulas ( 16) and ( 18)- [START_REF] Iserles | On the global error of discretization methods for highly-oscillatory ordinary differential equations[END_REF] reveal that in the corresponding functions κ r,k (τ ), r = 0, 1, 2, k = 0, ±1, ±2,..., the power τ 1 is only present in κ 0,0 and κ 1,k , κ 2,k , k = 0. This observation lowers the number of coefficients b u,r,k to be determined.

In Fig. 1 we have depicted the eight 'families' of trees u with weight ≤3 (each family has infinitely many trees corresponding to different values of the wave numbers). Accordingly, we write the averaged equation with J = 3 in the form where F μ stands for the contribution of the μth family. Using the methodology outlined above we find that in (29) (all functions evaluated at Y ):

d dτ Y = VIII μ=I F μ (Y ), (29)
F I = ǫf 1,0 ,F II = ǫ 2 f 2,0 ,F IV = ǫ 3 f 3,0 , F III = ǫ 2 k =0 i k (f ′ 1,0 f 1,k -f ′ 1,k f 1,0 ) + k =0 i k f ′ 1,k f 1,-k , F V = ǫ 3 k =0 i k (f ′ 2,0 f 1,k -f ′ 2,k f 1,0 ) + k =0 i k f ′ 2,k f 1,-k , F VI = ǫ 3 k =0 i k (f ′ 1,0 f 2,k -f ′ 1,k f 2,0 ) + k =0 i k f ′ 1,k f 2,-k , F VII = ǫ 3 k =0 1 k 2 f ′ 1,k f ′ 1,0 f 1,0 - ℓ =0 2 ℓ 2 f ′ 1,0 f ′ 1,ℓ f 1,0 + m =0 1 m 2 f ′ 1,0 f ′ 1,0 f 1,m + ℓ =0 1 ℓ 2 f ′ 1,0 f ′ 1,-ℓ f 1,ℓ - k =0 2 k 2 f ′ 1,k f ′ 1,0 f 1,-k + k =0 1 k 2 f ′ 1,k f ′ 1,-k f 1,0 - k =0,m =0 1 km (f ′ 1,k f ′ 1,-k f 1,m + f ′ 1,k f ′ 1,-m f 1,m ) - k =0,ℓ =0 k =-ℓ 1 (k + ℓ)ℓ f ′ 1,k f ′ 1,ℓ f 1,0 - ℓ =0,m =0 ℓ =-m 1 ℓ(ℓ + m) f ′ 1,0 f ′ 1,ℓ f 1,m + k =0,m =0 k =-m 1 km f ′ 1,k f ′ 1,-k-m f 1,m + k =0,m =0 k =-m 1 km f ′ 1,k f ′ 1,0 f 1,m , F VIII = ǫ 3 k =0 1 k 2 f ′′ 1,k [f 1,0 f 1,0 ]- ℓ =0 1 ℓ 2 f ′′ 1,0 [f 1,ℓ f 1,0 ] - ℓ =0 1 ℓ 2 f ′′ 1,-ℓ [f 1,ℓ ,f 1,0 ]+ k =0 1 k 2 f ′′ 1,0 [f 1,ℓ ,f 1,-ℓ ] - k =0,m =0 1 km f ′′ 1,k [f 1,-k f 1,m ]+ k =0,ℓ =0 k =-ℓ 1 k(k + ℓ) f ′′ 1,k [f 1,ℓ f 1,0 ] - ℓ =0,m =0 ℓ =-m 1 2ℓm f ′′ 1,-ℓ-m [f 1,ℓ f 1,m ]- ℓ =0,m =0 ℓ =-m 1 2ℓm f ′′ 1,0 [f 1,ℓ ,f 1,m ] .
Note that a comparison of the expressions for F I -F II -F IV or F III -F V -F VI bears out the independence-mentioned before-of the elementary coefficients of the first index j in the labels (j, k).

The Hamiltonian Case

Assume that H(y,τ; ǫ) is a Hamiltonian function with d degrees of freedom, 2πperiodic with respect to τ and possessing an expansion (cf. ( 10))

H(y,τ; ǫ) = 1 ǫ ∞ j =1
ǫ j H j (y, τ )

(here y is a 2d-dimensional vector). If f = J -1 ∇H , then (1)
is a non-autonomous Hamiltonian system whose right-hand side is of the form [START_REF]Analysis of the heterogeneous multiscale method for ordinary differential equations[END_REF] with f j = J -1 ∇H j . Moreover, if (cf. ( 11))

H j (y, τ ) = ∞ k=-∞ exp(ikτ )H j,k (y)
is the Fourier expansion of H j , then the Fourier coefficients f j,k in (11) satisfy f j,k = J -1 ∇H j,k . From Sect. 2 we know that the averaged equations are also Hamiltonian and furthermore, in the present section, we have expressed the averaged vector fields as combinations of elementary differentials. The general theory of symplectic integrators [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Sanz-Serna | Numerical Hamiltonian Problems[END_REF] shows that the averaged Hamiltonians will be combinations of the corresponding elementary Hamiltonians. For instance the Hamiltonian for ( 29)is

ǫH 1,0 + ǫ 2 H 2,0 + ǫ 2 k =0 i k ∇H T 1,0 J -1 ∇H 1,k + ǫ 2 k =0 i 2k ∇H T 1,k J -1 ∇H 1,-k + ǫ 3 H 3,0 + ǫ 3 k =0 i k ∇H T 2,0 J -1 ∇H 1,k +∇H T 1,0 J -1 ∇H 2,k + ǫ 3 k =0 i k ∇H T 2,k J -1 ∇H 1,-k + ǫ 3 k =0 1 k 2 ∇ 2 H 1,k J -1 ∇H 1,0 , J -1 ∇H 1,0 - ℓ =0 1 ℓ 2 ∇ 2 H 1,0 J -1 ∇H 1,ℓ , J -1 ∇H 1,0 - ℓ =0 1 ℓ 2 ∇ 2 H 1,-ℓ J -1 ∇H 1,ℓ , J -1 ∇H 1,0 + k =0 1 k 2 ∇ 2 H 1,0 J -1 ∇H 1,ℓ , J -1 ∇H 1,-ℓ - k =0,m =0 1 km ∇ 2 H 1,k J -1 ∇H 1,-k , J -1 ∇H 1,m + k =0,ℓ =0 k =-ℓ 1 k(k + ℓ) ∇ 2 H 1,k J -1 ∇H 1,ℓ , J -1 ∇H 1,0 - ℓ =0,m =0 ℓ =-m 1 2ℓm ∇ 2 H 1,-ℓ-m J -1 ∇H 1,ℓ , J -1 ∇H 1,m - ℓ =0,m =0 ℓ =-m 1 2ℓm ∇ 2 H 1,0 J -1 ∇H 1,ℓ , J -1 ∇H 1,m .

A Class of Autonomous Highly Oscillatory Hamiltonian Systems

We end this section by studying autonomous Hamiltonians of the form

1 2 p T 1 p 1 + 1 2 p T 2 p 2 + 1 2 q T 2 Kq 2 + U(q 1 ,q 2 ; ǫ), (30)
where p 1 ,q 1 are d 1 -vectors, p 2 ,q 2 are d 2 -vectors, U is a real-valued potential that may be expanded in non-negative powers of ǫ and K is a d 1 × d 1 symmetric positive definite matrix with eigenvalues of the form k 2 /ǫ 2 (k an integer that may change from eigenvalue to eigenvalue). Observe that when U is independent of q 2 , the system can be decoupled into d 2 harmonic oscillators with large frequencies k 1 /ǫ,...,k d 2 /ǫ and a Hamiltonian system with d 1 degrees of freedom and Hamiltonian (1/2)p T 1 p 1 + U(q 1 ; ǫ). When τ = t/ǫ is used as new independent variable, the Hamiltonian becomes

H(p 1 ,p 2 ,q 1 ,q 2 ; ǫ) = ǫ 2 p T 1 p 1 + 1 2 ǫp T 2 p 2 + 1 ǫ q T 2 Ω 2 q 2 + ǫU(q 1 ,q 2 ; ǫ), ( 31
)
where Ω is a symmetric positive definite matrix with integer eigenvalues, and the equations of motion are then d dτ p 1 =-ǫ ∇ 1 U(q 1 ,q 2 ; ǫ),

d dτ p 2 =- 1 ǫ Ω 2 q 2 -ǫ ∇ 2 U(q 1 ,q 2 ; ǫ), d dτ q 1 = ǫp 1 , d dτ q 2 = ǫp 2 .
(32)

While this system does not fit into the general family [START_REF] Engquist | Heterogeneous multiscale methods: a review[END_REF] considered so far in this section, its solutions y(τ) = (p 1 (τ ), p 2 (τ ), q 1 (τ ), q 2 (τ )) possess at stroboscopic times an expansion of the form (4) and are therefore amenable to our approach. In fact, let us introduce the time-dependent change of variables

p1 = p 1 , p2 = cos(τ Ω)p 2 + ǫ -1 Ω sin(τ Ω)q 2 , q1 = q 1 , q2 =-ǫΩ -1 sin(τ Ω)p 2 + cos(τ Ω)q 2 ,
that transforms (32) into a non-autonomous Hamiltonian system with Hamiltonian function

Ĥ p1 , p2 , q1 , q2 ,τ; ǫ = ǫ 2 pT 1 p1 + ǫU q1 , cos(τ Ω) q2 + ǫΩ -1 sin(τ Ω) p2 ; ǫ . (33)
Since the eigenvalues of Ω are positive integers, the change of variables is 2πperiodic in τ , and thus, it reduces to the identity map at stroboscopic times τ n = 2πn. Hence, for any solution y(τ) = (p 1 (τ ), p 2 (τ ), q 1 (τ ), q 2 (τ )) of (32), there exists a solution ŷ(τ) = ( p1 (τ ), p2 (τ ), q1 (τ ), q2 (τ )) of the system corresponding to (33) such that y(τ n ) =ŷ(τ n ) for all τ n = 2πn. We conclude that when trying to approximate y(τ n ) by stroboscopic averaging we may pretend that the solutions being treated are those of the Hamiltonian system with Hamiltonian (33), a system that does fit into the general framework of this section. This implies the existence of an averaged system for (32) that, furthermore, according to Sect. 2, will be Hamiltonian for a Hamiltonian function H(Y ; ǫ).

From our earlier formulas applied to (33), we find that (capital P 's and Q's denote the components of the averaged Y ) H(Y ; ǫ) satisfies

H(Y ; ǫ) = ǫ 1 2 P T 1 P1 + 1 2π 2π 0 U Q 1 , cos(Ωτ )Q 2 ; 0 dτ + O ǫ 2 .
The proposition in Sect. 2 implies that H(y) is a conserved quantity of the original system (32) and therefore the function

ǫ -1 H(y; ǫ) -H(y; ǫ) = 1 2 p T 2 p 2 + 1 ǫ 2 q T 2 Ω 2 q 2 - 1 2π
2π 0 U(q 1 ,q 2 ; 0) -U q 1 , cos(Ωτ )q 2 ; 0 dτ

+ O(ǫ)
Fig. 2 FPU problem. Oscillatory energies in each of the three stiff springs and total oscillatory energy. Exact solution will also be conserved. We have established in this way the adiabatic invariance [START_REF] Hairer | Geometric Numerical Integration[END_REF] of the oscillatory energy

1 2 p T 2 p 2 + 1 ǫ 2 q T 2 Ω 2 q 2
of the harmonic oscillators for solutions satisfying q 2 = O(ǫ) (in particular for solutions where the total energy (30) remains bounded as ǫ → 0). A well-known example of the family of Hamiltonians considered in this subsection is provided by the variant of the Fermi-Pasta-Ulam problem studied in [START_REF] Hairer | Geometric Numerical Integration[END_REF] where q 1 and q 2 are m-dimensional and the potential is given by

U(q 1 ,q 2 ) = 1 4 (q 1,1 -q 2,1 ) 4 + (-q 1,m -q 2,m ) 4 + m-1 j =1
(q 1,j +1q 2,j +1q 1,jq 2,j ) 4 .

We have integrated on 0 ≤ t ≤ 200 the FPU problem with m = 3 (the initial conditions were taken from [START_REF] Hairer | Geometric Numerical Integration[END_REF] and we used the implicit midpoint rule with the-very small-time step 0.001). Figure 2 shows the exchange of energy among the three stiff springs, a phenomenon that manifests itself in scales of time t ∼ 1/ǫ. Figure 3 shows the results for the corresponding integration of the averaged (J = 3) system (29).

The averaged system, in spite of not following the oscillations with O(ǫ) period in the elongations of the stiff springs, reproduces rather well the exchange between the associated energies. In this section we study highly oscillatory second-order differential equations

d 2 dt 2 q = f q, t ǫ ; ǫ , (34)
where q ∈ R d and

f(q,τ; ǫ) = 1 ǫ ∞ j =0 ǫ j f j (q, τ ) (35)
with

f j (q, τ ) = ∞ k=-∞ exp(ikτ )f j,k (q).
For each j and k, f j,k ≡ f * j,-k . Note that here f is of size O(1/ǫ) due to the j = 0 term in the sum (35)( c f .( 10)). 6 We assume throughout that f 0,0 ≡ 0: in this way all the leading O(1/ǫ) components of the force f in (35) are oscillatory and average out to zero, making it possible for q to undergo variations of size O(1) (rather than O(1/ǫ)) on time-intervals 0 ≤ t ≤ T of length O(1). The system (34) is supplemented by initial conditions q(0) = q 0 , p(0) = p 0 , where p(t) = dq(t)/dt and q 0 , p 0 are given vectors of O(1) magnitude.

In terms of the scaled time τ = t/ǫ,wehave

d 2 dτ 2 q = ǫ 2 f(q,τ; ǫ) = ∞ j =0 ǫ j +1 f j (q, τ ); (36)
and, initially, dq(τ)/dτ has the value ǫp 0 .

The Expansion of q

The solution q(τ) of (36) possesses an expansion

q(τ) = q 0 + u∈UN ǫ |u| α u (τ ) σ u F u (p 0 ,q 0 ) = q 0 + j =1 ǫ j |u|=j α u (τ ) σ u F u (p 0 ,q 0 ). ( 37
)
Let us briefly point out the most significant differences with [START_REF] Fedorenko | Derivation and justification of equations in slow time (the stroboscopic method)[END_REF]. Due to the fact that we are now dealing with a second-order differential system (36) where the force ǫ 2 f does not depend on dq/dτ ,weusespecial Nyström trees [START_REF] Hairer | Solving Ordinary Differential Equations I, Nonstiff Problems[END_REF], i.e. trees whose vertices are of two types, meagre and fat, in such a way that the root is meagre, a meagre vertex has at most one son, which is fat, and fat vertices have only meagre sons. Only fat vertices are labeled and the set of labels (j, k) includes the possibility j = 0. The notation UN refers to the set of all special Nyström trees with labeled vertices. The weight |u| of a tree u ∈ UN is again the sum of the weights of its vertices; a meagre vertex has unit weight and a fat vertex with label (j, k) has weight j .Note that there exist vertices with zero weight and that, as a consequence, the number of vertices in a tree may be larger than its weight.

The elementary differential F u is, for each u ∈ UN, a mapping from the space R d × R d of the variables (p, q) into R d constructed in terms of the Fourier coefficients f j,k . A fat vertex with label (j, k) and r (meager) sons is associated with the Fréchet derivative of order r of the function f j,k and a terminal meagre son gives rise toatermp.

We rewrite (36) as an integral equation

q(τ) = q 0 + ǫp 0 + ǫ 2 τ 0 (τ -τ ′ )f q(τ ′ ), τ ′ ; ǫ dτ ′ ,
to obtain a procedure for the recursive computation of the elementary coefficients. The recursion starts from the tree consisting only of the (meagre) root, for which α u (τ ) = τ .Foratreeu ∈ UN with two or more vertices, let (j, k) be the label of the son of the root and denote by u ν ∈ UN the trees obtained by removing from u the root and its son; then

α u (τ ) = τ 0 (τ -τ ′ ) exp(ikτ ′ ) ν α u ν (τ ′ ) dτ ′ .
Let us now turn to the effective computation of the α u . The role played by the integral operator I in first-order differential systems is taken here by the operator I 2 that maps each smooth complex-valued function ψ of a real variable into the function

τ 0 (τ -τ ′ )ψ(τ ′ ) dτ ′ .
The elementary coefficients are now linear combinations of the functions (cf. ( 23)) χ r,k = I 2 (φ r,k ), r = 0, 1, 2,..., k = 0, ±1, ±2,.... (Note that (d 2 /dτ 2 )χ r,k (τ ) = φ r,k .) Clearly

χ r,0 (τ ) = τ r+2 (r + 2)(r + 1)
and furthermore, from χ r,k (τ ) = τκ r,k (τ )κ r+1,k (τ ) k = 0, and ( 18)-( 20), we find that, for k = 0,

χ 0,k = 1 k 2 + i k τ - 1 k 2 exp(ikτ ), (38)
χ 1,k = 2i k 3 - 1 k 2 τ - 2i k 3 exp(ikτ ) - 1 k 2 τ exp(ikτ ), (39)
χ 2,k =- 6 k 4 - 2i k 3 τ + 6 k 4 exp(ikτ ) - 4i k 3 τ exp(ikτ ) - 1 k 2 τ 2 exp(ikτ ). ( 40
)
We have calculated the elementary coefficients of all trees of weight j = 1, 2, to find the first terms in (37):

q(τ) = q 0 + ǫτp 0 + ǫ k =0 χ 0,k (τ )f 0,k (q 0 ) + ǫ 2 χ 0,0 (τ )f 1,0 (q 0 ) + ǫ 2 k =0 χ 0,k (τ )f 1,k (q 0 ) + ǫ 2 k =0 χ 1,k (τ )f ′ 0,k (q 0 )p 0 + ǫ 2 k =0,ℓ =0 c k,ℓ (τ )f ′ 0,k (q 0 )f 0,ℓ (q 0 ) + O ǫ 3 . ( 41
)
Here

c k,ℓ (τ ) = 1 ℓ 2 χ 0,k (τ ) + i ℓ χ 1,k (τ ) - 1 ℓ 2 χ 0,k+ℓ (
τ ) or, after replacing the χ 's by their values found above,

c k,ℓ (τ ) = 1 k 2 ℓ 2 - 2 k 3 ℓ - 1 ℓ 2 (k + ℓ) 2 + i 1 kℓ 2 - 1 k 2 ℓ - 1 ℓ 2 (k + ℓ) τ - 1 k 2 ℓ 2 - 2 k 3 ℓ - 1 ℓ 2 (k + ℓ) 2 exp(ikτ ) - i k 2 ℓ τ exp(ikτ ), for k = -ℓ, while c k,-k (τ ) = 3 k 4 + 2i k 3 τ - 3 k 4 exp(ikτ ) + i k 3 τ exp(ikτ ) - 1 2k 2 τ 2 .

The Expansion of the Averaged Solution

Just as in Sect. 3, the expansion of the averaged solution is obtained by interpolating the elementary coefficients. The leading terms of the B-series for q provided in (41) thus yield

Q(τ ) = q 0 + ǫτp 0 + ǫ k =0 i k τf 0,k (q 0 ) + ǫ 2 τ 2 2 f 1,0 (q 0 ) + ǫ 2 k =0 i k τf 1,k (q 0 ) -ǫ 2 k =0 2 k 2 τf ′ 0,k (q 0 )p 0 + ǫ 2 k =0,ℓ =0 ck,ℓ (τ )f ′ 0,k (q 0 )f 0,ℓ (q 0 ) + O ǫ 3 , (42)
where, for

k = -ℓ, ck,ℓ (τ ) = i 1 kℓ 2 - 2 k 2 ℓ - 1 ℓ 2 (k + ℓ) τ while ck,-k (τ ) = 3i k 3 τ - 1 2k 2 τ 2 .

Averaged Differential Equations

Differentiation with respect to τ of (42) leads to

d 2 dτ 2 Q τ =0 = ǫ 2 f 1,0 (q 0 ) - k =0 1 k 2 f ′ 0,k (q 0 )f 0,-k (q 0 ) + O ǫ 3 ,
and therefore the lowest-order (J = 1) averaged equation is

d 2 dτ 2 Q = ǫ 2 f 1,0 (Q) - k =0 1 k 2 f ′ 0,k (Q)f 0,-k (Q) (43)
an expression not easily guessed from the mere inspection of (36).

It is important to observe that, from (42), so that there is an O(ǫ) discrepancy between the initial value ǫp 0 (itself of size O(ǫ)) of dq/dτ and the initial value of dQ/dτ . As a consequence, the (smooth) time derivative dQ/dτ of the smooth interpolant of q is not an approximation to the smooth interpolant of p(t) = dq/dτ . If the force f in the original oscillatory problem (35) is the gradient of a scalar potential, then the averaged (43) possesses also the conservative format (d 2 /dτ 2 )Q = ǫ 2 ∇W , with the potential W given by (cf. [START_REF] Levi | Geometry and physics of averaging with applications[END_REF])

d dτ Q τ =0 = ǫp 0 + ǫ k =0 i k f 0,k (q 0 ) + O ǫ 2 , (44)
W(Q)= V 1,0 (Q) - 1 2 k =0 1 k 2 ∇V 0,k (Q) 2 .
Here V 1,0 and V 0,k are the potentials for f 1,0 and f 0,k respectively, and |•|represents the Euclidean norm. We conclude this subsection by computing, as in Sect. 3, an averaged equation of higher accuracy. Here we have to face the extra difficulty brought in by the fact that there are low weight trees with a large number of nodes. In Fig. 4 we have listed the thirteen 'families' of trees u in UN with more than one vertex and weight ≤3. Accordingly, we write

d 2 dτ 2 Q τ =0 = XIII μ=I F μ + O ǫ 4 ( 45
)
where F μ stands for the contribution of the μth family. By using a methodology similar to that employed in Sect. 3 to reduce the number of coefficients b u,r,k to be computed, we find:

F I = F III = F VI = 0, F II = ǫ 2 f 1,0 , F IV =-ǫ 2 k =0 1 k 2 f ′ 0,k f 0,-k , F V = ǫ 3 f 2,0 , F VII =-ǫ 3 k =0 1 k 2 f ′′ 0,k [p 0 ,p 0 ],
F VIII =-ǫ 3 k =0 1 k 2 f ′ 0,k f 1,-k + k =0 1 k 2 f ′ 0,k f 1,0 , F IX = ǫ 3 - k =0 1 k 2 f ′ 1,k f 0,-k + k =0 1 k 2 f ′ 1,0 f 0,k , F X = ǫ 3 k =0 2i k 3 f ′ 0,k f ′ 0,-k p 0 , F XI =-ǫ 3 k =0,ℓ =0 2i k 2 ℓ f ′′ 0,k [p 0 ,f 0,ℓ ], F XII = ǫ 3 - k =0,ℓ =0 1 k 2 ℓ 2 f ′′ 0,k [f 0,ℓ ,f 0,-k ] - k =0,ℓ =0m =0 1 k 2 ℓm f ′′ 0,k [f 0,ℓ ,f 0,m ] + k =0,ℓ =0 k =-ℓ 1 2k 2 ℓ 2 (k + ℓ) 2 f ′′ 0,k [f 0,ℓ ,f 0,-(k+ℓ) ] , F XIII = ǫ 3 - k =0,m =0 2 k 3 m f ′ 0,k f ′ 0,-k f 0,m - k =0,m =0 1 k 2 m 2 f ′ 0,k f ′ 0,-k f 0,m + k =0,ℓ =0 1 k 2 ℓ 2 f ′ 0,k f ′ 0,ℓ f 0,-ℓ + k =0,ℓ =0 k =-ℓ 1 k 2 (k + ℓ) 2 f ′ 0,k f ′ 0,ℓ f 0,- (k+ℓ)
.

In (45) it is understood that the functions f j,k and their Fréchet derivatives

f ′ j,k , f ′′ j,k
that feature in the F μ 's are evaluated at q 0 . The averaged differential equation we are seeking is now obtained by discarding the O(ǫ 4 ) remainder in (45) and eliminating q 0 and p 0 with the help of (44), i.e.

d 2 dτ 2 Q = XIII μ=I F μ (46)
where the functions F μ are to be evaluated now at Q and

P = ω d dτ Q - k =0 i k f 0,k (Q)
rather than at q 0 and p 0 . Note that the right-hand sides (forces) of the original oscillatory system (36) and the averaged system (43) are independent of the corresponding velocities dq/dτ ,dQ/dτ ; this is not the case for the system (46).

An Example: The Kapitsa Pendulum

As an illustration of the preceding material, we consider the differential equation

d 2 dt 2 q = g L + 1 ǫ v max L cos t ǫ + θ 0 sin q (47)
that describes the motion of a pendulum whose suspension point is subjected to a fast vertical vibration. In (47), q is the angle between the pendulum rod and the upward vertical axis, L the length of the rod, g the acceleration of gravity, 1/ǫ the (large) angular frequency of the vibration of the suspension point, v max > 0theO(1) maximum vertical velocity of the suspension point and θ 0 a parameter that governs the initial phase of the vibration. The mechanical system described by ( 47) is sometimes referred to by the name of Kapitsa's pendulum [START_REF] Kapitsa | Dynamical stability of a pendulum when its point of suspension vibrates[END_REF] and attracts much interest in physics as an example of the possibility of stabilization by vibration that led to Paul's 1989 Nobel Prize [START_REF] Calvo | Carrying an inverted pendulum on a bumpy road[END_REF][START_REF] Levi | Geometry and physics of averaging with applications[END_REF][START_REF] Sanz-Serna | Stabilizing with a hammer[END_REF][START_REF] Sharp | Multiple time scale numerical methods for the inverted pendulum problem[END_REF]. At the lowest order in ǫ, the averaged equation (43) turns out to be

d 2 dt 2 Q = g L - v 2 max 2L 2 cos Q sin Q; (48)
the term with -v 2 max /(2L 2 ) opposes the gravity term g/L and, for v max suitable large, turns into a stable equilibrium the usually unstable configuration Q = 0, where the pendulum rod is above the suspension point.

At the next order in ǫ the averaged equation in formula (46) reads

d 2 dt 2 Q = g L - v 2 max 2L 2 cos Q sin Q + ǫ v max L cos(θ 0 )P 2 sin Q - v 2 max L 2 sin(2θ 0 )P sin 2 Q - v 3 max 4L 3 cos(3θ 0 ) -3 cos(θ 0 ) sin 3 Q , (49)
where

P = d dt Q + v max L sin θ 0 sin Q.
The O(ǫ) terms in (49) stem from the families of trees VII, XI and XII respectively. (The contributions from the families VIII and IX cancel each other and those from families V and XIII vanish.) If θ 0 =±π/2 (modulo 2π ), then (49) reduces to (48), a fact that could have been anticipated by noting that, for those values of θ 0 ,t h e equation ( 47) is not altered when ǫ is changed into -ǫ.

Application to the Analysis of Multiscale Methods

Heterogeneous Multiscale Methods

For the sake of brevity, we shall only study here the case of second-order problems (34); all our considerations may easily be extended to first-order systems of the form (1). The numerical integration of the highly oscillatory system (34) may be a very difficult task: standard explicit algorithms suffer from stability restrictions that limit the step-length to size O(ǫ) and more sophisticated methods are often hindered by the phenomenon of order reduction [START_REF] Calvo | Instabilities and inaccuracies in the integration of highly oscillatory problems[END_REF]. Heterogeneous multiscale methods (HMM) [1, 10-13, 29, 31]( c f . [START_REF] Li | Deciding the nature of the coarse equation through microscopic simulations: the baby-bathwater scheme[END_REF]) avoid these difficulties by aiming at finding only the slowly varying components of the solution q without keeping track of the rapidly oscillatory components. A full discussion of the various HMMs is completely out of the scope of this paper and we shall limit ourselves to the asynchronous approach suggested in [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF].

From the preceding section we know that the solution q of (34) with initial conditions q(0) = q 0 , p(0) = p 0 differs in terms of size O(ǫ) from the solution of the averaged problem

d 2 dt 2 Q = F(Q), (50)
F(Q) = f 1,0 (Q) - k =0 1 k 2 f ′ 0,k (Q)f 0,-k (Q), (51)
with initial conditions (see (44))

Q(0) = q 0 , d dt Q t=0 = p 0 + k =0 i k f 0,k (q 0 ).
(Note that, if the f 0,k are real, then (d/dt)Q(0) = p 0 .) HMMs obtain approximations to q by numerically integrating (50) without using the explicit analytic expression given in (51) for the averaged force F . Each time that the numerical integrator (macro-integrator) used to integrate (50) requires the value of the function F at a known value Q * of its argument, F(Q * ) is approximated numerically by an average F of values of the force f of the originally given oscillatory problem (34). More precisely:

F Q * = 2 η η/2 -η/2 K 2t η f q * (t), t ǫ ; ǫ dt, ( 52
)
where η is a (small) scaling parameter, the kernel or weight function K is an even, K(ξ) = K(-ξ), real-valued function of the real variable ξ , -1 ≤ ξ ≤ 1, with unit mass,

1 -1 K(ξ)dξ = 1, (53)
and q * (t) is the solution of the original (34) with initial conditions

q * (0) = q 0 , d dt q * t=0 = 0. ( 54
)
In practice, the integral in (52) is approximated by a quadrature rule and the required values of q * are found by numerically integrating (34) with the initial data (54)inthe interval -η/2 ≤ t ≤ η/2. The numerical method used to perform this task is referred to as micro-integrator; the step-length δ used by the micro-integrator will typically suffer from a stability/accuracy restriction δ = O(ǫ), but this is acceptable because micro-integrations are only performed over small windows of length η rather than over the whole integration range 0 ≤ t ≤ T . (In fact, due to symmetry considerations, it is enough to micro-integrate over the interval 0 ≤ t ≤ η/2, see [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF].)

The error in a HMM consists [START_REF] Sanz-Serna | Modulated Fourier expansions and heterogeneous multiscale methods[END_REF]oftheaveraging error introduced by approximating the solution q by its averaged counterpart Q and the numerical error resulting from computing Q numerically. The numerical error may be analyzed by standard techniques after noting that it is nothing but the error of the algorithm used as macrointegrator when applied to (50) with approximate values of the force F . Therefore the numerical error has two origins: (a) the truncation error of the macro-integrator and (b) the use of inexact force values. In turn, the contribution (b) is the result of (b1) the use of the approximation (52) rather than the true averaged force in F in (51), (b2) the use of numerical quadrature to compute the integral in (52), (b3) the use of inexact values of q * obtained by numerical micro-integration. The analysis of (b2) and (b3) is standard, and accordingly we are only interested in the analysis of (b1), a task that is performed in the subsections that follow.

The Filtering Error

In view of (34), the formal series expansion for the function f(q * (t), t/ǫ; ǫ) that appears in (52) may be obtained by differentiating twice with respect to t the B-series for q * (t), which is given in (41) with q 0 = Q * , p 0 = 0(see(54)). After recalling that (d 2 /dτ 2 )χ r,k (τ ) = φ r,k (τ ), we find in this way

f q * (t), t ǫ ; ǫ = 1 ǫ k =0 φ 0,k t ǫ f 0,k Q * + φ 0,0 t ǫ f 1,0 Q * + k =0 φ 0,k t ǫ f 1,k Q * + k =0,ℓ =0 1 ℓ 2 φ 0,k t ǫ + i ℓ φ 1,k t ǫ - 1 ℓ 2 φ 0,k+ℓ t ǫ f ′ 0,k Q * f 0,ℓ Q * + R, (55)
where R stands for the series consisting of the higher-order terms. We may write R = R 1 (t) + R 2 (t), where R 1 and R 2 comprise respectively the non-oscillatory and the oscillatory terms. It is easily checked that R 2 = O(ǫ), while R 1 is a power series in t 2 , t 4 , . . . whose coefficients are O(1).

Next, substitution of (55)into(52) leads to

F Q * = 1 ǫ k =0 φ0,k f 0,k Q * + φ0,0 f 1,0 Q * + k =0 φ0,k f 1,k Q * + k =0,ℓ =0 1 ℓ 2 φ0,k + i ℓ φ1,k - 1 ℓ 2 φ0,k+ℓ f ′ 0,k Q * f 0,ℓ Q * + R1 + R2 , (56)
where φr,k denotes the result of filtering the function φ r,k ,

φr,k = 2 η η/2 -η/2 K 2t η φ r,k t ǫ dt,
and Rj stands for the result of filtering the remainders R j , j = 1, 2,

Rj = 2 η η/2 -η/2 K 2t η R j (t) dt.
Clearly R1 = O(η 2 ) and R2 = O(ǫ). This means that, for accuracy, the value η should be so small that the filtering procedure in (52) retrieves satisfactorily the nonoscillatory components of the micro-force f being filtered. The normalization (53) implies that φ0,0 = 1, and then a comparison of (56) with (51) reveals that for an 'ideal' filter K for which φr,k = 0,r = 0, 1,k = 0, (57)

we have F(Q * ) -F(Q * ) = R1 + R2 = O(η 2 + ǫ).
Reference [START_REF] Engquist | Heterogeneous multiscale methods for stiff ordinary differential equations[END_REF] considers the use of filter functions K(ξ) of class C ν that vanish with their derivatives at ξ =±1 (the exponential filter used in the experiments of [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF][START_REF] Sanz-Serna | Modulated Fourier expansions and heterogeneous multiscale methods[END_REF][START_REF] Sharp | Multiple time scale numerical methods for the inverted pendulum problem[END_REF] satisfies this requirement for any value of ν =, 1, 2,...). Such filters 7 only satisfy the target (57) in an approximate manner. In fact, since

φr,k = η 2ǫ r 1 -1 K(ξ)ξ r exp ikη 2ǫ ξ dξ,
ν integrations by parts reveal that

φr,k = O k -ν ǫ η ν-r ,k = 0. ( 58
)
It follows that

F Q * -F Q * = R1 + R2 = O η 2 + ǫ + 1 ǫ ǫ η ν
(an expression found in [START_REF] Sanz-Serna | Modulated Fourier expansions and heterogeneous multiscale methods[END_REF] by a different approach). This shows that, for filters in this family, the window length η should be chosen large with respect to the period 2πǫ of the fast oscillations. Since at the same time, as pointed out above, η has to be small with respect to the characteristic time of the smooth, non-oscillatory components of q, it follows that algorithms based on smooth filters cannot operate unless ǫ is much smaller than such a characteristic time, i.e. unless the problem has wellseparated time scales.

Simple Filtering

A simple alternative to the smooth filters we have just discussed was suggested in [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF]. It has K ≡ 1/2 and η = 2πǫ, i.e. it filters through the familiar mean value

F Q * = 1 2πǫ πǫ -πǫ f q * (t), t ǫ ; ǫ dt. (59)
Since here η = O(ǫ) while smooth filters require η ≫ ǫ, the micro-integrations are now much cheaper. However it should be emphasized that this alternative technique exploits the fact that for the problems under consideration the period of the fast oscillations is known beforehand; the scope of the use of smooth filters does not suffer from such a limitation. As is the case for the smooth filters, the simple filtering formula (59) does not satisfy the conditions (57). In fact, while it is clear that φ0,k = 0, k = 0, the quantities φ1,k are non-zero. However φ1,-k =-φ1,k and therefore under the additional hypothesis f 0,k ≡ f 0,-k it will be true that, after comparing (51) and (56),

F(Q * ) -F(Q * ) = R1 + R2 or F Q * -F Q * = O(ǫ).
Since f 0,k ≡ f * 0,-k , the hypothesis f 0,k ≡ f 0,-k is equivalent to the requirement that the Fourier coefficient f 0,k be real or, in other words, that f 0 (q, t/ǫ) be an even function of t. When f 0 (q, t/ǫ) is not even, the filtering error F(Q * ) -F(Q * ) is O(1) and simple filtering is inconsistent as already noted in [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF].

Higher-Order Averaged Equations

The HMMs considered so far in the literature have been based on the lowest-order averaged equation ( 50)-(51) with O(ǫ) error at stroboscopic times. 8 We now illustrate the possibility of devising multiscale numerical methods based on the averaged equation ( 46), with O(ǫ 2 ) error. We write (46) compactly as

d 2 dt 2 Q = F 2 Q, dQ dt . (60)
(F 2 depends also on ǫ, but this is not reflected in the notation.) To simplify matters we assume that (i) f(q,t/ǫ; ǫ) is an even function of t and (ii) p 0 = 0. Then the solution Q of (60) with initial condition Q(0) = q 0 , (d/dt)Q(0) = 0i sa nO(ǫ 2 ) approximation to the smooth interpolant of the solution q of (34) with q(0) = q 0 , p(0) = 0.9 Our aim is then to construct HMM algorithms that integrate (60) in such a way that the filtering errors perpetrated when computing F 2 arealsoofsizeO(ǫ 2 ).

When performing the macro-integration of (60) it is necessary to evaluate approximately the force F 2 at known values Q * , Q * of its arguments, a task that we suggest may be performed with the formula

F2 Q * , Q * = 2 η η/2 -η/2 K 2t η f q * (t), t ǫ ; ǫ dt, (61)
where, following the ideas in [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF], we take q * (t) to be the solution of (34) with initial conditions

q * (0) = q 0 , d dt q * t=0 = Q * .
In the terminology of [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF], this is an asynchronous approach because the microintegrations are always performed on the interval -η/2 ≤ t ≤ η/2, regardless of the current value of t in the macro-integration. The algorithm just suggested may be analyzed in a way similar to that employed above. Now it is necessary to consider the O(ǫ 3 ) terms in the series (41); the elementary coefficients of those terms are linear combinations of the functions χ 0,k , 8 At non-stroboscopic times Qq is still O(ǫ). χ 1,k (as are the O(ǫ 2 ) terms) and χ 2,k . A function χ 2,k (t/ǫ), k = 0 in the series for q(t) gi v esrisetoafunctionǫ -2 φ 2,k (t/ǫ) in the series for f = (d 2 /dt 2 )q and, after filtering, to a coefficient ǫ -2 φ2,k in the series for F2 . On the other hand, the contribution of χ 2,k (t/ǫ) to the exact averaged force is obtained, as discussed in Sect. 4, by first replacing in (40) the factor exp(ikt/ǫ) by 1 and then computing the second derivative with respect to t at t = 0. It is then obvious that such a contribution equals -ǫ -2 (2/k 2 ). Therefore, the conditions (57) for an 'ideal' filter have now to be supplemented by the additional requirement

φ2,k =- 2 k 2 ,k = 0. ( 62
)
The bound in (58) suggests that it is not possible for the smooth filters considered above to satisfy (62). For the simple filter K ≡ 1/2, a trite computation yields

φ2,k = 2 k 2 exp(ikπ), k = 0,
and (62) holds only for odd k. Hence, the simpler filter will not achieve the O(ǫ 2 ) error we aim at, unless f 0,k vanishes for even k, i.e. unless the leading component f 0 (q, t/ǫ) in the force (35) consists only of Fourier modes exp(ikt/ǫ)f 0,k (q) for which the filtering interval -πǫ ≤ t ≤ πǫ comprises an odd number of periods, a conclusion not easily guessed in advance! We have integrated numerically on the interval 0 ≤ t ≤ π/4 the inverted pendulum problem (47) for which the leading component f 0 of the force f consists only of the odd wave numbers k =±1. The values of the parameters were g = 9.8, L = 0.2, v max = 4, θ 0 = 0, q(0) = 0.5, p(0) = 0.5, with ǫ ranging from 10 -2 to 10 -5 .T h e 'classical' fourth-order Runge-Kutta method 10 was used both as a macro-and microintegrator and the macro-force was taken from (61) with the simple, constant filter function K ≡ 1/2. The macro-step-size was taken from the sequence (π/20)/2 ν , ν = 0, 1, 2,... and the micro-step was determined from the formula δ = (πǫ/4)/2 ν (see [START_REF] Calvo | Heterogeneous multiscale methods for mechanical systems with vibrations[END_REF]). Since each micro-integration covers the interval -πǫ ≤ t ≤ πǫ, the number of micro-steps in a single micro-integration equals 8 × 2 ν (independently of ǫ). There are four evaluations of F2 /micro-integrations per macro-step and therefore the total number of micro-steps in 0 ≤ t ≤ π/4 equals 4 × 5 × 8 × 4 ν . The simulations were performed only for those combinations of and ǫ for which the macro-step-length exceeds the vibration period 2πǫ;for ≤ 2πǫ, HMMs do not make much sense, and it is better to use a conventional numerical method.

In Table 1 we give the maximum over 0 ≤ t ≤ π/4 of the absolute value of the difference between the computed Q and the exact solution of (60). Consider first the right-most column of the table (ǫ = 10 -5 ). For the coarser values of , the error behaves like O( 4 ) as corresponds to the fourth-order method used to macro/microintegrate. For = π/160 the algorithm yields errors of roughly one part in ten thousand: a remarkable achievement since the macro-step-length is more than 300 times larger than the vibrational period 2π × 10 -5 . However, towards the bottom of the column the error saturates: there the discretization error in the macro/microintegration is dominated by the filtering error. A comparison of the saturated errors at the bottom of the columns bears out neatly the O(ǫ 2 ) behavior of the filtering error predicted by our analysis. Note also that in the top rows of the table, where the filtering error is negligible, the errors are independent of ǫ.

We integrated again the same problem with the same parameter values, but now regarding the function cos(t/ǫ) in (47)ashavingperiod4πǫ, so that micro-integrations are performed in the interval -2πǫ ≤ t ≤ 2πǫ. With respect to this new period the wave numbers involved are k =±2 and (62) is not satisfied. The numerical results are given in Table 2: now the saturated errors at the bottom of the table are O(ǫ) in agreement with our analysis. A comparison between both tables shows that the violation of the condition (62) results in a degraded overall performance of the algorithm.

Finally, we integrated once more the same problem but now regarding cos(t/ǫ) in (47)a sh a v i n gp e r i o d6 π/ǫ, so that the relevant wave numbers are ±3 and (62) holds once more. The results of the experiment (not included here) show that the filtering error is restored to being O(ǫ 2 ). This dissipates any possible suspicion that the degraded performance in Table 2 was due to the wider micro-integration window, rather than to the fact that averaging was performed over an even number of cycles of the pendulum vibration.

Discussion

While, as pointed out before, the material in this section is not meant to discuss exhaustively HMMs, we hope it has illustrated two points:

1. The B-series/modified equation approach suggested in Sects. 2-4 may be advantageously used to analyze the behavior of the error in a variety of HMMs. 2. It is possible to devise multiscale methods that attain, for small macro-steplengths, errors of size O(ǫ 2 ) rather than merely O(ǫ). In fact, we have suggested above an asynchronous algorithm (based on the simple filter function) that achieves O(ǫ 2 ) errors under the hypotheses that the force f is an even function of t and consists only of Fourier modes with odd wave numbers. Presented in [START_REF] Calvo | A stroboscopic numerical method for highly oscillatory problems[END_REF] (by M.P. Calvo and the present authors) is an algorithm based on finite-difference techniques that results in O(ǫ ν ) errors ν = 1, 2,... in a wide range of oscillatory problems with a single high frequency.

Fig. 1

 1 Fig. 1 Families of trees in U with weight ≤3

Fig. 3

 3 Fig. 3 FPU problem. Oscillatory energies in each of the three stiff springs and total oscillatory energy. Solution of averaged problem

Fig. 4

 4 Fig. 4 Families of trees in UN with weight ≤3

  and these are retrieved from Ξ by eliminating the parameters from the relations Q = Ξ , (d/dτ)Q= ∂Ξ/∂τ, (d 2 /dτ)Q= ∂ 2 Ξ/∂τ 2 .

Table 1

 1 Errors in Q for the inverted pendulum, simple filtering Micro-steps ǫ = 10 -2 ǫ = 10 -3 ǫ = 10 -4 ǫ = 10 -5

	π/20	160	1.94(-1)	2.39(-1)	2.46(-1)	2.47(-1)
	π/40	640	3.24(-2)	1.74(-2)	1.75(-2)	1.75(-2)
	π/80	2,560	***	1.28(-3)	1.02(-3)	1.01(-3)
	π/160	10,240	***	3.15(-4)	7.17(-5)	6.70(-5)
	π/320	40,960	***	2.34(-4)	7.31(-6)	4.47(-6)
	π/640	163,840	***	***	2.67(-6)	3.22(-7)
	π/1280	655,360	***	***	2.26(-6)	4.48(-8)
	π/2560	2,621,440	***	***	2.21(-6)	2.45(-8)

Table 2

 2 Errors in Q for the inverted pendulum, simple filtering on a window of length 4πǫ (two vibrational periods)

		Micro-steps	ǫ = 10 -2	ǫ = 10 -3	ǫ = 10 -4	ǫ = 10 -5
	π/20	2 × 160	4.92(-1)	2.52(-1)	2.47(-1)	2.47(-1)
	π/40	2 × 640	5.22(-1)	6.76(-2)	1.35(-2)	1.71(-2)
	π/80	2 × 2,560	***	7.57(-2)	7.01(-3)	5.28(-4)
	π/160	2 × 10,240	***	7.65(-2)	7.79(-3)	7.23(-4)
	π/320	2 × 40,960	***	7.69(-2)	7.84(-3)	7.82(-4)

Comparing the general term α u /σ u F u of the B-series used here with that of the standard B-series, see[START_REF] Hairer | Geometric Numerical Integration[END_REF], Chap. III, one sees that our elementary coefficients are the counterpart of the product elementary weight × τ m ,wherem is the number of vertices in the tree.

In the theory of B-series the exact solution counts as a particular instance of numerical solution.

The structure of the elementary coefficients will be analyzed further in part II of the present work.

Of course, it would have been possible to include in Sect. 3 first-order differential equations with O(1/ǫ) right-hand sides. However, we shall see in this section that the terms with j = 0 introduce a number of significant complications that would have hindered the presentation there.

These filters will be referred to as 'smooth' filters in view of the fact that, after defining K(ξ) = 0f o r |ξ | > 1, they become C ν functions in the whole real line.

The hypotheses (i) and (ii) are not essential, but without them the presentation is encumbered by the formula that relates (d/dt)Q(0) and p 0 through differentiation of (42), cf. (44).

More precisely, we used the Runge-Kutta-Nyström method for second-order differential equations implied by the classical Runge-Kutta formula, see[START_REF] Hairer | Solving Ordinary Differential Equations I, Nonstiff Problems[END_REF], Chap. II.14.

Acknowledgements This research has been supported by 'Acción Integrada entre España y Francia' HF2008-0105. A. Murua is also supported by projects MTM2007-61572 (Ministerio de Educación, España) and EHU08/43 (Universidad del Pais Vasco/Euskal Herriko Unibertsitatea). J.M. Sanz-Serna is also supported by project MTM2007-663257 (Ministerio de Educación, España).