Averaging of time-periodic systems without a small parameter

In this article, we present a new approach to averaging in non-Hamiltonian systems with periodic forcing. The results here do not depend on the existence of a small parameter. In fact, we show that our averaging method fits into an appropriate nonlinear equivalence problem, and that this problem can be solved formally by using the Lie transform framework to linearize it. According to this approach, we derive formal coordinate transformations associated with both first-order and higher-order averaging, which result in more manageable formulae than the classical ones. Using these transformations, it is possible to correct the solution of an averaged system by recovering the oscillatory components of the original non-averaged system. In this framework, the inverse transformations are also defined explicitly by formal series; they allow the estimation of appropriate initial data for each higher-order averaged system, respecting the equivalence relation. Finally, we show how these methods can be used for identifying and computing periodic solutions for a very large class of nonlinear systems with time-periodic forcing. We test the validity of our approach by analyzing both the first-order and the second-order averaged system for a problem in atmospheric chemistry.


Introduction
The theory of averaging falls within the more general theory of normal forms; it arises from perturbation analysis and is thus formulated usually for equations which contain a small parameter . A vast literature treats this case (e.g. [2,5,10,13,14,29,30,32,33], and references therein). In this paper we develop a time-averaging theory for the case of timeperiodic nonautonomous systems without an explicit dependence on a small parameter. The theory is presented here formally, subject to a conditional statement to be explained later. It is shown to apply numerically to our test problem, which does contain a "large-amplitude" periodic perturbation.
Among the methods of derivation of averaged systems, those based on Lie transforms offer an efficient computational framework (cf. [5,31], or more recently [36]). The rigorous use of such techniques, however, is confined so far to the construction of near-identity transformations expanded in a small parameter.
In order to overcome the main issues raised by the lack of -parametrization, we consider the problem in terms of "differentiable equivalence" between the time-dependent vector field and a corresponding averaged form, allowing us to put the problem in an appropriate Lie transform framework.
In Hamiltonian dynamics, G. Hori [21] was the first to use Lie series for coordinate transformations, while A. Deprit [7] introduced the slightly different formalism of Lie transforms; the two were shown to be equivalent by J. Henrard and J. Roels [19]. They developed a new, efficient way to construct canonical transformations for perturbation theory in celestial mechanics. The Poisson bracket in their formulae can be replaced by the Lie bracket to obtain a representation of general, non-canonical transformations [5,16,23]. Nowadays, Lie transforms are applied in many research fields, such as artificial satellites [34], particle accelerator design [9] and optical devices [12], to mention just a few.
The object of the present article is to perform averaging analysis for a T -periodic Ndimensional system (N ≥ 2, T > 0) of ODEs where Ω is an open subset of R N . The precise framework is given in §2.
In the present paper, we consider the averaging process via an appropriate notion of differentiable equivalence (Definition 2.1). This equivalence concept allows us to determine the underlying time-dependent family of diffeomorphisms, as the solution of the nonlinear functional equation (2.6). This equation plays a central role in our approach, as it allows one to generalize the Lie transform formalism. We put it into a pullback form that falls into this formalism via suspension (Proposition 2.4 and Appendix). This process allows us to obtain a computable formal solution of equation (2.6) (Theorem 2.9), based on the solvability of a recursive family of linear PDEs, typically known as Lie's equations [5,18,35]. It is shown in §2.2 and §3.2 that Lie's equations can be solved easily in our framework.
This approach is consistent with other classical equivalence or conjugacy problems, in which the so-called homological equations linearizing such problems are a cornerstone in proving Anosov's theorem or Hartman-Grobman theorem (cf. [2] for details). A recent proof of Kolmogorov's theorem on the conservation of invariant tori also used this classical framework [22]. The essential difference between our approach and the classical one is that (2.6) involves a "conjugacy" (see Definition 2.1 for the right notion) between a nonautonomous field in the original problem and an autonomous one in the transformed problem. The formal developments in §2 to §4 are contingent upon Lemma 6.1. This Lemma will be proven in a subsequent paper (in preparation), and permits a rigorous proof of the existence of a solution to equation (2.6). The numerical results in §5 and §6 demonstrate the usefulness of the present approach and the plausibility of the forthcoming rigorous results.
Our approach provides, furthermore, simpler formulae for higher-order averaged systems, even in the classical context of perturbations scaled by a small parameter . This is shown in Proposition 3.1 and is due to the suspension associated with solving (2.6); see §3.2, §3.3 and [10,32].
We show that the classical change of coordinates on the initial data respects the differentiable equivalence relation (2.6) at time t = 0. It arises here in a more natural way from the generator of the family of inverse diffeomorphims at time t = 0; see Proposition 4.2 and [17].
We demonstrate that the Lie transform framework permits gaining CPU time with respect to the standard numerical methods by, loosely speaking, integrating the averaged system with a larger step size; this is followed by performing efficient corrections with a step size adapted to the oscillations of the forcing terms. We test the validity of this numerical approach by computing the second-order terms in our atmospheric-chemistry problem. More precisely, we build a second-order averaged system in this test problem and the relevant correction in the same framework, after applying the averaging-correction method at the first order.
In §2 we describe how Lie transforms allow us to resolve formally a differentiable equivalence problem between a time-periodic system and an autonomous one. In §3, we give a formal algorithm of computing higher-order averaged systems based on Lie transforms and related corrections; furthermore we compare our results to the classical ones. In §4 we clarify how to compute initial data which satisfy the differentiable equivalence relation. Numerical results and advantages of this technique are discussed in §5 for a simple model of atmospheric chemistry. In §6, we comment on the use of these methods for constructing periodic orbits, in nonautonomous dissipative systems with periodic forcing. In the Appendix, the Lie transform framework is presented as a tool for solving pullback problems in general.
2 The Lie transform setting and the differentiable equivalence problem

Fitting the differentiable equivalence problem into the Lie transform formalism
For all the rest of the paper Ω will be an open subset of R N , N will be a positive integer, and T will be a positive real. For subsequent computations we shall consider Q k (Ω) as the class of continuous functions f : R + × Ω → R N which are T -periodic in t for each x ∈ Ω and which have k-continuous partial derivatives with respect to x for x ∈ Ω and t in R + . More precisely we will work with Q ∞ (Ω) = k≥0 Q k (Ω) when we deal with Lie transforms as often it is requires in this framework for dependence on the spatial variables (cf. Appendix).
There are known results on the boundedness and the global existence in time of solutions of ODEs. For instance, by an application of Gronwall's lemma it is easy to show that if there exists α, β ∈ C 0 (R + , , for all (t, x) ∈ R + × Ω, then every solution ofẋ = Y (t, x) is bounded. Thus, using the continuation theorem (cf. theorem 2.1 of [14] for instance), we can prove that the solutions of the above equations are global in time.
Taking into consideration this fact, we make the following assumption.
(λ) All solutions remain in Ω, for all considered non-autonomous and autonomous vector fields defined on Ω.
We describe some concepts related to the subject of this paper, dealing with differentiable equivalence between a nonautonomous vector field and an autonomous one.
Definition 2.1. Let r be a positive integer, let Y ∈ P r (Ω) and Z ∈ C r (Ω), then Y and Z are said to be P k dif f -equivalent (k ≤ r), if there exists a map Φ ∈ P k (Ω) such that for all t ∈ R + , Φ t := Φ(t, ·) : Ω → Ω is a C k -diffeomorphism, which carries the solutions x(t, x 0 ) ofẋ = Y (t, x), for x(0, x 0 ) = x 0 varying in Ω, into the solutions z(t, z 0 ) ofż = Z(z), in the sense that: Using the notations in this definition, we shall say that if (2.1) is satisfied for a pair (x, z) of solutions, then x and z are in P k dif f -correspondence by Φ, or, sometimes that they are in Furthermore, we will often use for the underlying transformations the functional space where the smoothness assumption on the map t → Φ −1 t , will be apparent from the proof of Lemma 2.3.
Remark: It is important to note here that the P k dif f -equivalence of Definition 2.1 does not conserve the period of the trajectories although it preserves parametrization by time and sense, and, in this meaning, realizes a compromise between the classical notions of conjugacy and equivalence (e.g. [13]). This will be essential in §6.
Let us recall the classical definition of pullback of autonomous vector fields.
The problem of P k dif f -equivalence takes the form of a family of "pulled back" problems via the ∈ Ω, and z an integral curve of Z. Then, x and z are in Proof: The proof of this lemma is an obvious application of the chain rule formula and Definition 2.1.
Let us now introduce the projection π : We define also a map, defined for each t ∈ R + , as which will be T -periodic in time t. According to Lemma 2.3, if there exists a map Φ ∈ P k d (Ω) which satisfies the following nonlinear functional equation, (called in this work equivalence problem) then every integral curve of Z is in P k dif f -correspondence by Φ to a unique solution of: x = Y (t, x), that is the vector fields Y and Z are P k dif f -equivalent. The problem is then to solve in Φ ∈ P k d (Ω), the equation (2.6) for a given pair of fixed vector fields Y ∈ P r (Ω) (non-autonomous) and Z ∈ C r (Ω) (autonomous).
In order to show how Lie transforms formalism (e.g., [5,7,19] and Appendix) can be used for solving formally equation (2.6) with r = k = ∞, we put this equation in an appropriate pullback form by the Proposition 2.4 below.
For that we introduce some notations. We denote by Y the so-called t-suspended vector field associated with (1.1) written in terms of the enlarged vector of dynamical variables, Note that by assumption (λ), every solution x of (2.7) is contained in Ω := {(t, x) ∈ R + ×Ω}, for the rest of the paper.
Remark: Note that for the rest of the paper the t-suspended form of a (time-dependent or not) vector field X on R N will be defined as X = [1 R , X] T , and the flat form of X will be defined as [0 R , X] T . The same notation X will be used for the latter, by abusing the notation. No confusion will be made with regard to the context. Proposition 2.4 Let Z ∈ C ∞ (Ω) be an autonomous vector field and let Y ∈ P ∞ (Ω), with Z and Y as their respective t-suspended forms. Let Θ ∈ P ∞ d ( Ω). If Θ satisfies Θ −1 t * Y = Z, for all t ≥ 0, then the family of C ∞ -diffeomorphisms (π • Θ t • I t ) t∈R + , acting on Ω, satisfies the equation (2.6). The converse is also true.
Proof: The proof relies on basic calculus and the definitions listed above.
We illustrate now the strategy used to solve equation (2.6) for Z = Y := 1, Y T with Y being the classical first averaged system, namely 1 T T 0 Y ds. The same reasoning will be used for higher-order averaged systems in §3. By Proposition 2.4, the problem of finding a solution of (2.6) falls into the Lie transform formalism if Y and Y are given by τ -series (cf. Appendix and Theorem 7.1). This point of view, according to the solvability of Lie's equations, will then allow us to solve formally (2.6) in §2.3.
Starting from that point, denote Y by Y ave . In order to write Y ave and Y as τ -series we introduce a natural auxiliary real parameter τ , by the following expressions with the obvious definitions for the vector fields Y Consider ϕ τ,t the semiflow at τ generated by a τ -suspended smooth field G t for each non-negative t, i.e. ϕ τ,t is the solution at τ of with a formal expansion of G(t, ξ) in powers of τ given by with G and G n smooth vector fields on R N .
Remark: The generator system would be denoted by another notation, the τ -suspension and the t-suspension being not equivalent generally. We slightly abuse the notation.
Applying the discussion of the Appendix with p = N + 1, A τ = Y τ and B τ = Y ave,τ for all τ , and utilizing the fact that Lie's equations are independent of τ , we obtain that if there exists ϕ τ,t such that (ϕ τ,t ) * Y τ = Y ave,τ for every non-negative real t, then by Corollary 7.2 with H τ,t = G τ,t , Lie's equations are solvable for the G n,t . Note that Y (resp. Y ave ) corresponds to τ = 1 in (2.9) (resp. (2.8)), so if (ϕ 1,t ) * Y = Y ave for every non-negative real t, then by Proposition 2.4, the one-parameter family (π • ϕ 1,t • I t ) t∈R + solves the equation (2.6). We solve Lie's equations (2.14) in the next section and we will show how this solvability will allow us to compute formally a solution of (2.6) in §2.3.

Solving
Lie's equations associated to the equivalence problem (2.6)-first-order averaged system We treat in this section the problem of solvability of Lie's equations associated with (2.6) for Z = Y ave . In that respect, let M be the map from Z * + × Z * + to the set of the vector fields of R N +1 defined by where the terms Y (j) k are calculated from those of (2.9) and (2.10) by the recursive formula (7.5).
We have the following result: Proposition 2.5 For every integer j greater than or equal to two and for every integer l belonging to {0, .., j − 2}, we have M (j, j − l) = 0.
We make the following hypothesis, P(j), until a j > 3, Then, as a consequence of this predicate, we get We show the heredity of the predicate. The formula (7.5) for n = l and i = j − (l + 1), takes the form so is re-written, for M , = 0 according to (2.14), so M (j, j − l) = 0 for every integer l between 0 and j − 3. Making l = j − 3 in (2.16) we get M (j, 2) = 0 and therefore P(j + 1).
We make a useful remark for the clarity of some coming computations.
Remark: Let Y be the vector field associated with the system (1.1), with an expanded t-suspended form given by (2.9) and (2.10). Let G be a vector field on R N with G given on the model of (2.12) and (2.13). Let Y (m) j be the terms calculated by the recursive formula This last proposition allows us to solve Lie's equations easily in our case, which we can summarize in the following Proposition 2.6 Suppose that ϕ 1,t generated by G t given by (2.12) entirely determined by the sequence (Ỹ (0) j ) j∈Z * + given by (2.10) and (2.9).
More exactly we have, for every positive integer j, modulo a constant vector, and, 0 ) is defined by the first Lie equation and (2.14).
Proof: First we prove that G 0 (and G 0 ) is determined as follows. By definition of the Lie derivative in R N , the first Lie equation (cf. Appendix) takes the form Y given by (2.9) and (2.14), which gives G 0 by integration. We obtain Y (1) 1 as follows. By (7.5) we get is given in (2.14) and Y (2) 0 is equal to zero. Thus we get Y 0 , which yields (2.18) for j = 1.
Computing G 1 by (7.5) gives: and, after easy manipulations as Y yielding (2.17) for j = 1 after integration.
We consider the following predicate, H(j), for j ≥ 1:

hal-00162628, version 1 -16 Jul 2007
We prove H(j) by recurrence. The preceding computations show H(1). Let j ≥ 2. Suppose H(j − 1). Again (7.5) gives, and, as previously, we have The terms Y (0) k are equal to zero for all k ≥ 2 (see (2.10)), so, using (2.12), (2.13), (2.23), and Y (0) 1 = Y , we get by integration the two first assertions of H(j), noting that the first component of the r.h.s of (2.23) is equal to zero (that gives α (1) j = 0). Finally, we obtain the last assertion of H(j) as follows. From (7.5), we have for every positive integer j: j−1 = 0 for every positive integer j, so we deduce the last assertion of H(j) from (2.25). Thus H(j) is a necessary condition of H(j − 1), that completes the demonstration.

Formal solution of the equivalence problem (2.6)
We describe in this section a procedure for obtaining formal series representation of a solution of (2.6), assuming that Lie's equations are solved. In other words, we assume, for this section, that Lie's equations associated with the problem of equivalence under consideration are solved.
Introduce for any sufficiently smooth vector fields F, H : (τ, ξ) ∈ R N +1 → R N , where D ξ F stands for the usual Jacobian of F in the local coordinates ξ ∈ R N . Naturally, Λ n H represents the iterated operation The following lemma is the first step toward efficiently describing P ∞ dif f -correspondence between a solution x of (1.1), and one z ofẏ = Z(y).
Proof: This lemma is a consequence of the Taylor formula. Let z be an integral curve ofẏ = Z(y), and x the solution of (1.1) being in P ∞ dif f -correspondence with z by (π • φ 1,t • I t ) t∈R + . We have formally for τ = 1: where ζ means the triplet (τ = 0, t, z(t)).
Consider the projection π on the N last components (projection which commutes with the operator d dτ ). If we compose (2.28) to the left by π we obtain: , it is therefore sufficient to prove for every integer n greater or equal than 1 that: We prove (2.30). Using the chain rule it is merely an exercise in basic calculus. By assumptions, Let us introduce the obvious notation φ τ,t = (Id R ; φ 1 τ,t ; ...; φ N τ,t ) T and, denoting the i th component of H t,τ by H i t,τ , the notation H t,τ = (1, H 1 t,τ , ..., H N t,τ ) T . Moreover, let us recall that the classical pullback by a diffeomorphism Φ of a map f : thus (2.31) can be viewed component-wise as: , by application of the chain rule, taking into account (2.32), we obtain finally the key formula: where ∇ stands for the gradient and ·, · for the usual Euclidean scalar product on R N . By introducing the following operator, which acts on families of maps from R N to R: we can write (2.34) in the form: and thus, by induction it is clear that, for every integer n, which, according to (2.33), leads to and taking this expression at τ = 0 we get, because φ 0,t = Id R N +1 : Let F be a map from R N +1 to R N . By definitions (2.32) and (2.26) we see that Consider now H t and F t time-dependent vector fields acting on R N , both admitting τseries representations. Omitting the dependence on t for the vector fields in the expansion of F t and H t , we consider the following recurrence formula: n = F n,t , F n,t being the n th term of the τ -series of F t ; and where H n−k is the (n − k) th term in the τ -series of H t .
With these notations, we formulate Definition 2.8. Let t ∈ R + fixed. The H t -transformation evaluated at τ , denoted by T Ht (τ ), of a time-dependent vector field F t acting on R N , by a time-dependent vector field H t acting on the same space, both admitting formal τ -series, is given by: 0,t is calculated using the recurrence formula (2.42), initialized by F n,t = F n,t ; F n,t being the n th term of the formal τ -series of F t . When F = H, the term F 0,t will be called the (n+1) th -corrector.
Remark: It is important to note here that T Ht (0) · F t = 0, allowing near-identity transformations for small τ (cf. Theorem 2.9).
We then have the following theorem: Theorem 2.9 If φ τ,t is generated by a time-dependent vector field of the form (1, H t ) T on a region R + × Ω, for which H t has a formal τ -series, n≥0 τ n n! H n,t , on this region for all τ ∈ I; I being an interval of reals such that 0 ∈ I, then: (2.44) Proof: First of all, note that the derivations in the proof of Lemma 2.7 show that, for all t ∈ R + and τ ∈ I, Omitting the time-dependence, simple calculus shows and defining the vector fields H [1] n as following a procedure similar to the one of [19], but for the operator Λ H , we get the recurrence formula: n is the n th term of the series n = H n for every non-negative integer n. We conclude therefore that: which proves the theorem after taking into account Definition 2.8.
In practice Theorem 2.9 and Corollary 7.2, permit to determine formally a solution of (2.6). For example, when the "target" autonomous field is Y , as shows Proposition 2.6, Lie's equations (2.14) are solvable, thus (Id R N + T Gt (1) · G t ) t∈R + can be computed using computer algebra implementation [3,25,6] and thereby gives a formal solution of (2.6) (with Z = Y ), where the generator G is determined by (2.17) in Proposition 2.6.

New higher-order averaged systems and related corrections
Higher-order averaged systems are linked to the way change of variables are performed. In the -dependent case, the usual formulae are obtained by using formulae for arbitrary-order derivatives of compositions of differentiable vector functions [11], on the one hand, and by using the implicit function theorem, on the other [10,30,32]. This form for higher-order averaging is not the only one, as it is pointed out in [29, pp. 36-37]. We describe here other higher-order averaged systems both for -dependent and independent cases, via Lie transforms formalism, but different from those of [31,36], and given by simple explicit formulae.

Formal algorithm of higher-order averaging based on Lie transforms
We describe here a way of computing formally the higher-order averaged systems according to the Lie transform of the vector field Y ∈ P ∞ ( Ω) in consideration.
In that respect we consider first the parametric standard formẋ = τ Y (t, x), where τ is a real parameter. Then we have the following proposition.
Proposition 3.1 Let n be an integer greater than or equal to two. Let Y (n) ∈ C ∞ (Ω) be the n th term of the finite sequence satisfying and W 0 = G 0 , with G 0 given from (2.22).
Proof: Let an integer n ≥ 2. If there exists φ τ,t generated by W t = [1 R , W t ] T having a τ -series representation on the model of (2.12) and (2.13), such that for every positive real t, the map π • φ τ,t • I t satisfies (3.1), then according to Proposition 2.4 the equivalence problem where Y (n) is given by assumption (a τ ).
If we applied (7.4) in this case, we get: 0,t ), and the t-suspended form for m = 0. By projecting onto the N -last components of R N +1 , by taking into account the remark before Proposition 2.6 we get for every m ∈ {2, · · · , n}, . By (7.5) we have for m ∈ {2, · · · , n}: and continuing the process to an arbitrary rank l ∈ {2, · · · , m − 1}: we get: is equal to Y (see (2.9)) and taking into account (2.10), we get and noting moreover that, as in (2.24), Integrating (3.11) with respect to time, from 0 to the minimal period T , we obtain according to the fact that Y Finally, noting that the first Lie equation is the same at each-order (cf. (3.18) and related discussion), we easily conclude that W 0 = G 0 , where G 0 is given by (2.22), which completes the proof of this proposition.
Note that when we talk about the n th averaged system, in any cases, we will consider naturally the following system of differential equations: where Y (n) will be given by (3.2) in Proposition 3.1.

Correction of higher-order averaged systems and solvability of related Lie's equations
We adopt here the procedure described in §2.3 for correcting higher-order averaged systems. We consider the system of differential equations (1.1), and for a given p ≥ 2 a higher- Consider the diffeomorphism φ 1,t generated by W t = [1, W t ] T in Proposition 3.1. Then, according to Theorem 2.9, the transformation π • φ 1,t • I t is given formally as the following formal series representation that allows computer algebra implementation [6], and where W t depends on the averaged system Y (p) (cf. again Proposition 3.1).
It is important to note here that for two given averaged systems, Y (p) and Y (q) , (p = q), the corresponding generators, say W and V for fixing the ideas, related to each system are not a priori equal. Therefore, considering Theorem 2.9 and Definition 2.8, the W t -and V t -transformation given by analogous formulae to (3.13), are also not equal. For instance for p = 1 (resp. p = 2), i.e. when W = G (resp. W = Γ), G (resp. Γ) being defined as the generator of time-dependent diffeomorphisms which transform Y into 0 ); the second corrector in T Gt (1) · G t and in T Γt (1) · Γ t are respectively G [1] 0,t = G 1,t + DG 0,t .G 0,t , (3.14) which is given by (2.42) with F = H = G, i = 0, n = 0 and, is given by (3.3) with m = 2. Same algebraic procedure for obtaining Γ 1,t and related computations, in the case p = 1, leads to

17)
modulo a constant. Therefore, we see that the G t -transformation and the Γ t -transformation are different. Indeed by substituting (3.17) and (3.16) respectively in (3.14) and (3.15), we note that the second correctors G [1] 0,t and Γ [1] 0,t are distinct. Same remarks hold for higher-order correctors related to two distinct higher-order averaged systems, allowing us to conclude the noninvariance of every higher-order corrector.
We consider now the invariance of the first corrector. Indeed, the relation (7.5) for i = 0, n = 0, A = Y and H 0 = W 0 or H 0 = V 0 yields for the same Lie's equation after projection. This demonstrates that the first corrector term W 0,t is independent from the averaged system considered.
We conclude this section concerning the solvability of Lie's equations related to higherorder averaged systems.
Fix p ≥ 2. For computer implementation of (3.13), we have to initialize (2.42) with H = W = F ; in other words we have to solve Lie's equations associated with (3.1). Those equations show a structure similar to that studied in §2.2. We outline how to solve them for p > 2, and we give the precise results for the case p = 2. Here, the key point of the Lie's equations solvability consists in extending Proposition 2.5 associated to the first-order equivalence problem (i.e. for n = 1 in (3.1)) to one associated to a higher-order equivalence problem (3.1). Indeed, based on (7.5), by introducing a map N on the model of (2.15) with elements of the diagonal given by (3.3), we can show that a proposition analogous to Proposition 2.5 holds, by shifting of (p − 1) ranges the first column of vanishing terms in the infinite matricial representation of N . This allows us to compute the W j (i.e. the j th -term of the τ -series representation of W in Proposition 3.1), by adapting computations giving the G j in Proposition 2.6. We illustrate this procedure in the case p = 2 for the convenience of the reader.
In this latter case, we can indeed show the two propositions below. First of all note that, by applying (3.3) in Proposition 3.1 with m = 2, we have where Γ 0 is given by integration of the r.h.s of (3.18) because of the invariance of the first corrector. The analogous of Proposition 2.5 associated to the second-order equivalence problem (i.e. n = 2 in (3.1)) is is given by (3.19), and N (j, j) = Y (j) 0 are null for every integer j greater than or equal to three. Then for every such integer j, and for every integer l belonging to {0, ..., j − 3}, we have N (j, j − l) = 0.
Following the model of the proof of Proposition 2.6 that was based on Proposition 2.5, we can show, taking into account Proposition 3.2, for all t ≥ 0. Then the sequence (Γ j , Y j ) j∈Z * + is entirely determined by the sequence (Ỹ (0) j ) j∈Z * + given by (2.10) and (2.9). More exactly we have for every positive integer j,

22)
with every first component of Y This concludes the solvability of Lie's equations related to the second-order averaged system. Such propositions can be extended for higher-order averaged systems, that rules, by induction, the solvability of related Lie's equations.

A heuristic comparison with the classical approach
We make in this section a few comments on previous works and ours as well, and discuss a rigorous setting concerning formulae derived in the preceding sections, in the case when τ is a small parameter, usually denoted by .
As explained in §3.2, Lie's equations are solvable to each order of the averaged systems, giving thereby the time-varying change of coordinates by Theorem 2.9. For instance, for n = 2 in Proposition 3.1, we can solve Lie's equations related to the second-averaged system (cf. Proposition 3.3). Consider the system (1.1), if we make the change of variables 0,t (y), (3.23) (that is a truncation up to order 2 of T Wt (1) · W t ), then following the procedure of [13, pp. 168-169] under sufficient smoothness assumptions, we obtain that y satisfieṡ This shows that our procedure of averaging-correction is coherent and rigorously justified up to order two, for sufficiently small. Extending this kind of reasoning to order greater than 2, following, for example, the procedure described in [30, §3], we can show that our procedure is rigorously justified to any order, for sufficiently small, making thereby near-identity transformations. This permits us to obtain convergence theorems on n th -order averaging, with suitable assumptions, such as those used in [32], for instance. These kind of results increase precision but not the timescale of validity. The latter can be improved if the system shows some attracting properties, as it is classically done [29,30,33]. ¿From a practical point of view, we obtain more manageable formulae than those usually obtained [10,32]. Indeed, our approach based on Proposition 3.1 gives an efficient procedure of construction of any higher-order averaged systems, a procedure which can be implemented on a computer using a symbolic computational software, as it is done for Lie transforms in general [3,25,6]. Actually, many authors [5,30,31] pointed out the fact that Lie transforms can be used to obtain higher-order averaged systems, but in our knowledge, a systematic explicit exposition like the one given here has not been published (although a related paper [36] presents some of the elements.) The fact that the classical formulae are not manageable may be more evident in [15], where the authors develop an algorithm of higher order averaging for linear equations with periodic coefficients, in order to simplify the one of [10].
Finally, it is interesting to note that usually the i th -correctors are thought to be the antiderivative in time of the oscillatory part associated to each order (cf. (6) via (7) and (5) of [32] for example). This structural form is the same for correctors W t∈R + to a unique solution x (n) of the n th averaged system given by (3.12) in Proposition 3.1. We recall that the latter means In a number of applications x 0 is given, and thus, in order to compute an approximation of x based on x (n) , we have to determine the initial condition appropriate for P ∞ dif fcorrespondence, that is (π • φ −1 1,0 • I 0 )(x 0 ) as shown by (4.1). This problem of determination of the inverse transformation for a given Lie transformation is not new, and several algorithms have been proposed (see e.g. [17,35]). We present here an approach based on [16].

(4.2)
In the context of our framework, Proposition 5 of [16] can be reformulated as the following Proposition 4.1 Let φ 1,t generated by H t defined in (4.2) and let φ −1 1,t generated by the τ -suspended vector field K t , then: Proof: We provide a proof for the sake of completeness. Let us denote the inverse of φ τ,t by η τ,t , i.e., for all (τ, t) ∈ R × R + , and x ∈ Ω, φ τ,t (η τ,t (x)) = x.
Taking the derivative with respect to τ yields We express the second term on the left-hand side in terms of the generator H t , and move it to the right-hand side to obtain which gives, after multiplying both sides with the inverse of (D The latter leads to Therefore, for each τ and each t, we deduce that η τ,t is generated by −(φ τ,t ) * H t , and thus substituting τ = 1 the proof is complete.
The following proposition, which is an obvious corollary of Theorem 2.9 and Proposition 4.1, where T Kt (1) is given by Definition 2.8, determines formally the inverse transformation for all t ≥ 0, and in particular the one in (4.1) for t = 0. Proposition 4.2 Let φ 1,t generated by a time-dependent vector field of the form H t = (1, H t ) T and let K t = −(φ 1,t ) * H t , then for every non-negative real t: This proposition allows us to compute the m th approximation of the solution of (1.1), which is in P ∞ dif f -correspondence to one of a n th averaged system by the following algorithm which defines the m th approximation of the n th type: (1). Compute the n th averaged system given by (3.12) and Proposition 3.1.
(2). Solve Lie's equations as explained in §3.2 until the order m, giving the m first W t,i 's terms of the τ -series of W .
(3). The initial condition x(0) of the exact system being given, truncate to the order m the K 0 -transformation (at time t = 0) taken at x(0), T K0 (1) · K 0 (x(0)), where K is loosely speaking the generator of the family of inverse transformations, and take this truncation plus x(0) as initial condition for the n th averaged system. (4). Compute numerically the solution x (n) of the n th averaged system, and add to x (n) (t) the truncation to the order m of the series T Wt (1) · W t (x (n) (t)), which gives finally the m th approximation of the n th type of the exact solution at time t.

Application to a problem in atmospheric chemistry
The original problem that motivated this work was to examine models of diurnal forcing in atmospheric dynamics and chemistry on long time-scales [8]. The day-to-night changes in the radiative heating and cooling of the planetary boundary layer -the lowest part of the atmosphere (1-2 km) -are very large. Still, one is often only interested in the slow, season-to-season or even year-to-year changes in the way this lower layer interacts with the underlying surface (land or ocean), on the one hand, and the free atmosphere above, on the other [8]. The diurnally averaged model we derived here provides insight into a similar problem, that of the basic chemistry of slow changes, from one day or week to the next, of a highly simplified system of photochemically active trace gases in the troposphere (i.e., the lower 10 km of the atmosphere). The system of two coupled ODEs we consider in this paper governs the concentration of the chemical species CO (carbon monoxide) and O3 (ozone) [20], This system belongs to the general class of ODEs systems given by (1.1).
In this system, the diurnal forcing is through the functions S i and Z i (see Figure 1 below for a typical example), which can have rather complicated shapes. The system (5.1) is a very simple model for air pollution in an urban environment. Changes in the chemistry at the day-night transitions and those in the emission of CO and O3 in a city during a day lead to S i and Z i which are only piece-wise smooth.
While it is not difficult to numerically integrate the system as it is by standard methods, the forcing carried by functions S i and Z i can have high frequencies and therefore require a time step too small compared to the total time of simulation, T max (say, several tens of years). This assumes that the goal is to obtain sufficiently smooth and therefore realistic numerical solutions.
We discuss in this section how well the numerical solutions of the averaged systems, at orders one and two, approximate those of the original one. An important part of the numerical work for this particular example is to carry out the transformations between the solutions of the original and averaged systems developed in §2. This includes the transformation of initial values of §4. Overall, the numerical simulations indicate that the correction performed by Lie transforms to the first order of the first averaged system M 2 or the second averaged system M 2 (2) obtained by applying results of §3, provides a good approximation to the original one. Furthermore, regularity of the averaged systems in time allows one to integrate them by simple methods, such as an Euler or a Runge-Kutta of order 4 (RK4 in the sequel). Last but not least, we can also use larger step sizes for the averaged systems than for the original one, which is the key to speed up simulations of long-term dynamical phenomena.

First and second averaged systems analysis
Let Ω := R 2 \{x 2 = 0}. Then the vector field Y associated to M 2 belongs to P ∞ (Ω) and Y ∈ C ∞ (Ω). Assume that there exists, for each t ≥ 0, a diffeomorphism φ 1,t ∈ P ∞ d ( Ω), generated by In order to demonstrate the numerical efficiency of the averaged systems, we present only an approximation of the family of diffeomorphisms (φ 1,t ) t∈R + up to order one.
For that, we compute the first corrector G 0,t = G 0,t (see Definition 2.8) given by integration of (2.22). Define G 0,1 (t, ξ) (resp. G 0,2 (t, ξ)) as the first (resp. second) component of G 0 (t, ξ). Then, from (5.1), where, for i ∈ {1, 2}, The constant parts C 0,1 and C 0,2 of (5.2) and (5.3) are given by assuming that for all ξ ∈ Ω, T 0 t 0 G 0 (ξ, s)dsdt = 0, which allows us to avoid "secular" terms at the first order, namely and from (5.3): Now, let x(0) be the initial condition of a solution x in system M 2. In order to compute the first approximation of the first type, following the procedure described at the end of §4, we have to truncate to order one the K 0 -transformation at time t = 0 of K 0 taken at x(0), i.e. T K0 (1) · K 0 (x(0)), where K is, loosely speaking, the generator of the family of inverse transformations. Thereby, we get a first approximation of the initial condition x(0) of the solution x which is P ∞ dif f -correspondent to x by the family (φ 1,t ) t∈R + .
So we have by Proposition 4.2, end of §4 and Theorem 7.1, that the first approximation of the first type x (1,1) (0) of the initial condition x(0) is with G 0 (0, x(0)) determined by the preceding constants evaluated at x(0). Note that the first index of the exponent (1, 1) in (5.7) indicates that we deal with the first type of approximation and note that the second superscript indicates that we make a truncation up to order one in the series T K0 (1) · K 0 (x(0)). An identical convention will be used in the following when we work with approximations of higher-order type. Next we numerically compute the solution z of M 2 through x (1,1) (0) and then we define the first pullback x (1) by truncating to "order one" the expression z(t) + T Gt (1) · G t (z(t)), which gives: We describe now the construction of the second pullback. By Proposition 3.1 with n = 2, the second averaged system M 2 (2) corresponds to the vector field Y Both direct computation by hand and symbolic manipulation software for Y T , yield for system (5.1) the following expression as function of ξ 1 and ξ 2 , namely, and: Denote by x (2) (0) the initial condition for M 2 (2) through which the solution x (2) is P ∞ dif fcorrespondent to x. We take as approximation of x (2) (0), following the preceding process and notations concerning the first averaged system and taking into account the invariance of the first corrector (cf (3.18)), the vector x (2,1) (0) is given by i.e. x (1,1) (0) defined in (5.7).
Then we numerically compute the solution v of M 2 (2) based on x (2,1) (0), and define the second pullback x (2) as the first approximation of the second type, i.e.: x (2) (t) := v(t) + G 0 (t, v(t)), for all t ∈ R + . (5.10) The choice of computing only the first approximation, based on the first and second averaged systems, allows us to compare the approximations of the original system M 2 by these averaged forms, the corrector being the same in each case, but applied to different solutions. This will be discussed in §5. 3.
The numerical experiments we present here correspond to choices of the forcing functions (S i , Z i ) which produce "broad" oscillations in the solutions of the original system. These are not very realistic but our goal here is to test how well the method performs for "large perturbations". The magnitude of perturbations is shown in Figure 1. In all the numerical simulations, Z 1 and S 1 were oscillatory but we used the constants Z 2 = 5.0 × 10 −2 and S 2 = 12.0 × 10 −2 .

Gain in numerical efficiency
The original system is integrated by standard methods (Euler, Runge-Kutta,...) with step size δt, which has to be small enough to resolve the oscillations induced by the forcing terms Z i and S i . The averaged system M 2 or M 2 (2) is integrated with step size t, which is typically ten times larger than δt.
A key point in numerical efficiency is the regularity of the solutions of the averaged system. Indeed, if the solutions of the averaged system are smooth enough in the sense that they do not show multiple oscillations over one forcing period, which is one day in our case, we can integrate the averaged system with a large step size t, without loss of regularity. This fact for system (5.1) is pointed out in Figure 2, for instance, where the solutions of the averaged systems M 2 and M 2 (2) , on the second component, are those which do not show multiple oscillations over one day. This fact is well known in the -dependent case where the drift described by the averaged system can be integrated with a step size chosen to be 1 times larger than for the non-averaged system [2]. In order to investigate how well the solutions of the averaged system corrected up to order one approximate that of the original one, we did the following: (i) For a given initial condition for the original system, x(0), we compute the 1 st -approximation of the 1 st -type x (1,1) (0) and the 1 st -approximation of the 2 nd -type x (2,1) (0) given by (5.7) and (5.9), of the initial conditions x(0) and x (2) (0), respectively.
(ii) Starting from these initial conditions, we compute x and x (2) by integrating M 2 and M 2 (2) with t as step size, by using a standard integrator. This costs less than integrating M 2 with a small enough step size that resolves the oscillating forcing terms.
Here it is sufficient to use the Euler method.
(iii) We do a simple linear interpolation on x and x (2) to obtain the solution on the temporal grid defined by δt.
(iv) According to (5.8) (resp. (5.10)) and the fact that the correctors are only composed of integrals, (5.2) and (5.3) are used for computing, on the grid defined by δt, the first and the second pullback based on solutions obtained in (ii). The numerical tests indicate that the solution of M 2, obtained by the procedure described through the items (i) to (iv), has accuracy which is close to that obtained by integrating M 2 itself by classical methods, such as a RK4 scheme. The CPU time linked with our procedure depends essentially on the one defined for solving the averaged systems M 2 or M 2 (2) with step size t, and on the one for the integral correction process. We observed that the use of the first pullback allows to obtain a gain in CPU time nearly equal to 75 percent when we use a large step size for t, while retaining good accuracy and regularity of the solutions, as shown in the following subsection where the tests and the related figures are for δt = 0.01 and t = 0.1.

Accuracy achieved by the first and second pullback
One of the main reasons for computing the second averaged system M 2 (2) were to obtain a better approximation. In this section we demonstrate numerically this fact for δt = 0.01 and t = 0.1. Note that the tests performed by RK4 on M 2 with t give no smooth solutions (not shown), whereas, as shown in Figures 3 and 4, this is not the case for solutions obtained by our method described in (i)-(iv) of §5.2.
Based on the method of approximation used here, in many cases, the second pullback is better than the first, and gives a solution which is close to the one obtained by RK4. In Figures 3 and 4, we observe that we have coincidence of the 1 st -, 2 nd -pullback and the full solution for the global trend in each component. A more accurate analysis shows that the 2 nd -pullback is a better approximation than the first as shown in Figure 5 and Figure  6 where the 2 nd -pullback is represented by dash, the 1 st -pullback by dots, and the solution computed with RK4 by solid line. This fact is made more evident by the analysis of the error in the supremum norm as shown in Figure 7, where the error produced by the 2 nd -pullback is under 6 × 10 −3 and in Figure 8 where the error with the 2 nd -pullback is represented by dash.         Nevertheless, we can easily imagine that for a system of dimension larger than two the CPU time required to calculate the 2 nd -pullback would increase with algebraic complexity of the 2 nd -averaged system (cf. §5.1), but as we shall see in §6 the 2 nd -pullback can be used for the problem of determination of periodic solutions in a general framework.
The tests presented in this study give satisfactory results, an accuracy of order 1 × 10 −3 for the second pullback. Moreover, we have to note that this accuracy is obtained with oscillations actually "far" from their averages (see Figure 1 again) producing oscillations, of magnitude about 1 × 10 −1 , on the full solution obtained by RK4, with a significant gain in CPU time compared to the standard method.
Therefore we can conclude that the Lie transforms averaging method developed here gives a rigorous setting for constructing the correctors and M 2 (2) to provide numerical approximations even for relatively large perturbations in time induced by the forcing terms. Furthermore, according to the tests performed in this work, the 2 nd -averaged system analysis seems to be more relevant than the first.       [4] with relaxed assumptions, for the first, second and third averaged systems, based on Brouwer degree theory. All these results are linked to a maximal size of perturbation 0 under which the existence of a hyperbolic fixed point p 0 of the first averaged system gives the existence of a unique periodic solution of the original system, revolving around p 0 . As was explained in §3.3, the analysis performed in this paper permits to get such results for higher-order averaging analysis based on explicit formulae (Proposition 3.1 for averaging and Theorem 2.9 for corrections) by revising the proof of classical results for first and second order averaged systems, given in the literature (e.g. [13]). For instance, if the k (k ≥ 1) first averaged systems vanish identically, then a proof of existence and unicity of a T -periodic solution in a -neighborhood of a hyperbolic fixed point of the (k + 1) th -averaged system can be given on the basis of principles given in the proof of theorem 4.1.1 of [13], essentially by showing that the Poincaré maps of nonautonomous and autonomous systems are -closed.
In the light of such considerations, a question arises naturally: do there exist such results for the -independent case?
A preliminary analysis can be done on time T -maps. Indeed, we can formulate the following lemma: Lemma 6.1. Let Y ∈ P r (Ω), and Z ∈ C r (Ω), (1 ≤ r ≤ ∞). Then Y and Z are P r dif fequivalent if and only if their time T -maps are C r -conjugate.
This lemma gives thereby a way to study sufficient conditions for the existence of solutions of the nonlinear functional equation (2.6), which will be investigated in a forthcoming paper. Note that the existence of time T -maps is realized under the assumption (λ).
In practice, if the conjugacy of time T -maps is achieved, we can localize a T -periodic orbit of the nonautonomous system Y . Indeed, suppose that the time T -map associated with an averaged system Z has a fixed point. Then by conjugacy, this fact still holds for the time T -map associated with the nonautonomous system. Making use of the proof of Lemma 6.1, we can observe that if (π • φ 1,t • I t ) t∈R + denotes the solution of (2.6) obtained by Lie transforms, then π • φ −1 1,o • I 0 is the diffeomorphism realizing the conjugacy between time T -maps associated with Y and Z. Therefore Proposition 4.2 and Theorem 2.9 allow us to compute an approximation of this diffeomorphism, which can be applied to a fixed point η associated with the system Z, leading to an approximation ξ 0 of an initial datum lying on the T -periodic orbit of the nonautonomous system which is in P ∞ dif f -correspondence with η. Thus, by a standard integrator based at ξ 0 , we can compute an approximation of this T -periodic solution. Such a method can be a useful tool for localizing T -periodic solutions of T -periodic nonautonomous dissipative systems, which are known to exist, in this case, as shown by classical results on the topic (cf. [27, p. 235]).
This process is relevant for our system (5.1), as proved by elementary analysis. Indeed, there exists a hyperbolic fixed point p 0 = (1.584, 0.431) T of M 2 and, as we can see on Figures 3,5 and Figures 4, 6, the first pullback solution is a good approximation of the periodic solution revolving around the fixed point.
Finally, note that the notion of equivalence considered in this paper (Definition 2.1) is the appropriate one from the numerical perspective described here in order to compute T -periodic solutions, in view of the fact that the main ingredient was to achieve a correspondence between fixed points (with all the periods) of an autonomous system and T -periodic orbits of a nonautonomous one.

Appendix. Solution of the time-dependent pullback problem via Lie transforms
Let p be an integer greater than or equal to one. Let A(τ, x) be a smooth vector field on R × R p expanded in powers of τ as and let H(τ, t, x) be a smooth vector field on R × R + × R p expanded as, H(τ, t, ξ) = H τ (t, ξ) = H τ,t (ξ) = n≥0 τ n n! H n (ξ, t) = n≥0 τ n n! H n,t (ξ), (7.2) where t is fixed. Denote by Ξ t the semiflow generated by H t , namely the semiflow of S Ht dξ dτ = H(τ, t, ξ) ; ξ ∈ R p , τ ∈ R.  with L H n−k,t expressing the Lie derivative with respect to H n−k,t (see e.g. [5,35], for details). A classical result is the following theorem which expresses the pullback of a vector field as a Lie transform which can be proved easily, by adapting, for instance, the work of [19] relying on Taylor series expansion and the key formula d dτ (Ξ τ,t ) * A τ = (Ξ τ,t ) * (L Hτ,t A τ + ∂ τ A τ ) (7.6) from the realm of differential geometry (e.g., [1,28]).
Theorem 7.1. Let τ ∈ R, as above. The pullback at time t of A τ by Ξ τ,t generated by H t is the Lie transform generated by H t of A, evaluated at τ , that is: (Ξ τ,t ) * A τ = L(H t )(τ ) · A, for all t ∈ R + .
Let B be another smooth vector field on R × R p with formal expansion One of the advantages of considering pullback as a Lie transform, is that the framework of the latter yields linear conditions for finding Ξ τ,t which satisfies (Ξ τ,t ) * A τ = B τ . Indeed, using the formal power series expansions (7.4) and (7.7), (Ξ τ,t ) * A τ = B τ leads to a sequence of recursive linear PDEs (e.g. [35]), namely, A (m) 0,t = B m , m ∈ Z + , t ∈ R + , (7.8) where the H n,t (n ≤ m) present in (7.2) are the unknowns contained in A (m) 0,t , determining by this way the generator H t of the Lie transform.
These equations are usually called Lie's equations.
As a consequence, we can state the following corollary of Theorem 7.1: Corollary 7.2. A necessary condition for the existence of a two-parameter family of diffeomorphisms (Ξ τ,t ) (τ,t)∈R×R + , generated by the one-parameter family of vector fields (H t ) t∈R + given by (7.2), such that (Ξ τ,t ) * A τ = B τ for all t ≥ 0 and all τ ∈ R, is that Lie's equations (7.8) are solvable for the unknowns H m,t for all (m, t) ∈ Z + × R + .