Higher-Order Averaging, Formal Series and Numerical Integration II: The Quasi-Periodic Case

The paper considers non-autonomous oscillatory systems of ordinary differential equations with d≥1 non-resonant constant frequencies ω ₁,…,ω _d . Formal series like those used nowadays to analyze the properties of numerical integrators are employed to construct higher-order averaged systems and the required changes of variables. With the new approach, the averaged system and the change of variables consist of vector-valued functions that may be written down immediately and scalar coefficients that are universal in the sense that they do not depend on the specific system being averaged and may therefore be computed once and for all given ω ₁,…,ω _d . The new method may be applied to obtain a variety of averaged systems. In particular, we study the quasi-stroboscopic averaged system characterized by the property that the true oscillatory solution and the averaged solution coincide at the initial time. We show that quasi-stroboscopic averaging is a geometric procedure, because it is independent of the particular choice of co-ordinates used to write the given system. As a consequence, quasi-stroboscopic averaging of a canonical Hamiltonian (respectively, of a divergence-free) system results in a canonical (respectively, in a divergence-free) averaged system. We also study the averaging of a family of near-integrable systems where our approach may be used to construct explicitly d formal first integrals for both the given system and its quasi-stroboscopic averaged version. As an application we construct three first integrals of a system that arises as a nonlinear perturbation of five coupled harmonic oscillators with one slow frequency and four resonant fast frequencies.

1. Of course a problem with a resonant ω may be written in non-resonant form by lowering the number d of frequencies.

2. Alternatively it is also possible to let f depend smoothly on ε and expand in powers of ε as in [10].

3. In practice the change of variables is often expressed in alternative ways, for instance by using Lie series. With such alternatives, the recursion (7)–(10) has to be changed accordingly. However, the exact format of the change of variables and the recursion are not essential to the discussion that follows and we choose to focus on the form (4) for the sake of being definite. See in this connection Remark 2.16, which shows the relations between different averaged systems.

4. For the application of the idea of stroboscopic averaging to the construction of numerical integrators the reader is referred to [6, 7].

5. We work with real systems for the sake of being definite; all our results are valid for complex systems.

6. Repetitions among the u _ν are allowed. The order of the u _ν ’s is of no consequence: thus \([u_{1},u_{2}]_{\mathbf {k}}\) and \([u_{2},u_{1}]_{\mathbf {k}}\) are the same tree even if u ₁≠u ₂.

7. The effective computation of these coefficients will be considered in Sect. 4.

8. If A and B are sets, the notation B ^A means the set of all mappings ϕ:A→B.

9. This kind of initial-value problem was first considered in [5].

10. Recall that ω is assumed throughout to be non-resonant and that the resonant case may be rewritten in non-resonant form by lowering the value of d.

11. The method of characteristics shows that the solutions of (23) with are defined by

    $$\gamma(t, \theta) =\chi(\theta- t\omega) +\int _0^{t} \gamma\bigl(s, \theta+(s-t) \omega\bigr)*\beta\bigl(\theta+(s-t) \omega\bigr)\,\mathrm{d}s.$$

    The initial-value function \(\chi:\mathbb {T}^{d}\to \mathcal {G}\) has to vanish at \(\theta=\bf0\), but it is otherwise arbitrary.

12. The use of forests is essential in establishing the connection between formal series and Hopf algebras, see e.g. [4, 19]. That connection clarifies greatly the computation of the product ∗.

13. In fact (32) and (43) provide examples of the computation of an element in the Lie algebra—infinitesimal generator—in terms of the corresponding one-parameter semigroup.

14. This point is taken up in the next section.

15. Since the different references use different signs, we point out that here

    $$\{F,G\} = \sum_j \biggl( \frac{\partial F}{\partial p_i}\frac{\partial G}{\partial q_i} - \frac{\partial F}{\partial q_i}\frac{\partial G}{\partial p_i} \biggr).$$

16. The recursive formulae in Proposition 4.1 may be generalized to compute γ _u (t,θ) for arbitrary \(u\in \mathcal {T}\). Alternatively, one can compute the coefficients γ _u (t,θ) for all u by applying the explicit formulae given in [19] to find the coefficients α _u (\(u\in \mathcal {T}\)) of an arbitrary \(\alpha\in \hat { \mathcal {G}}\) in terms of the corresponding \(\alpha_{\mathbf {k}_{1}\cdots \mathbf {k}_{r}}\).

17. More precisely it is required that g _j and h are such that the transformed vector field f in (66) satisfies the requirements specified at the beginning of Sect. 2.

18. The exponent d _j in the scaling factor for each conserved quantity is determined by the scaling of the initial condition in tandem by the degree, as polynomials in the solution components, of the Poisson brackets discarded when truncating.

P. Chartier has received financial support from INRIA through the associated team MIMOL. A. Murua and J.M. Sanz-Serna have been supported by projects MTM2010-18246-C03-03 and MTM2010-18246-C03-01, respectively, from Ministerio de Ciencia e Innovación.

Authors and Affiliations

1. INRIA Rennes and ENS Cachan Bretagne, Campus Ker-Lann, av. Robert Schumann, 35170, Bruz, France

   P. Chartier

2. Konputazio Zientziak eta A. A. Saila, Informatika Fakultatea, UPV/EHU, 20018, Donostia–San Sebastián, Spain

   A. Murua

3. Departamento de Matemática Aplicada e IMUVA, Facultad de Ciencias, Universidad de Valladolid, Valladolid, Spain

   J. M. Sanz-Serna

Corresponding author

Correspondence to J. M. Sanz-Serna.

Communicated by Ernst Hairer.

An erratum to this article is available at http://dx.doi.org/10.1007/s10208-016-9311-2.

Appendix: Multiplication of Coefficients

In this appendix we give details of the product ∗ in \(\mathcal {G}\). We hope that this will increase the readability of the paper; for a more complete treatment the reader is referred to [15], Chap. III.

Given \(\delta\in \mathbb {C}^{ \mathcal {T}\cup \{\emptyset \}}\) with δ _∅=1 (that is \(\delta\in \mathcal {G}\), so that B(δ,y) is a near-identity B-series), and an arbitrary B-series B(η,y), the composition B(η,B(δ,y)) is a B-series of the form B(ζ,y), where ζ _∅=η _∅ δ _∅ and for \(u\in \mathcal {T}\) the coefficient ζ _u =(δ∗η) _u , \(u\in \mathcal {T}\) equals η _∅ δ _u +η _u δ _∅=η _∅ δ _u +η _u plus a sum of products of the form \(\eta_{u_{0}}\delta_{u_{1}} \cdots\delta_{u_{m}}\) with |u ₀|+|u ₁|+⋯+|u _m |=|u|, \(u_{0},\ldots,u_{m} \in{ \mathcal {T}}\). For instance, for the tree in the last row of Table 1 we have (we use brackets rather than subscripts for typographical convenience):

(80)

Multiplication of the formula (80) by \(\exp(i(\mathbf {k}+\mathbf {l}+\mathbf {m})\cdot\theta)\) establishes the validity of (52) at

In the general case, the coefficient (δ∗η) _u , \(u\in \mathcal {T}\) can be written as

(81)

where, for each fixed \(\delta\in \mathcal {G}\), \(L_{\delta}: \mathcal{V}\rightarrow\mathcal{V}\) is a linear map (defined below) in the vector space \(\mathcal{V}=\mathrm{span}( \mathcal {T}\cup \{\emptyset \})\) of linear combinations (with complex coefficients) of rooted trees, and η (originally defined as a map \(\mathcal {T}\cup \{\emptyset \}\rightarrow \mathbb {C}\)) has been extended by linearity to \(\mathcal{V}\) (i.e. \(\eta(\sum_{j}\lambda_{j} u_{j}) = \sum_{j} \lambda_{j} \eta_{u_{j}}\) for each linear combination ∑ _j λ _j u _j of trees).

In order to define the linear map L _δ , we first extend the notation \([u_{1} \cdots u_{m}]_{\mathbf {k}}\) (used to recursively represent trees in \(\mathcal {T}\)), to define for each m≥1 and each multi-index \(\mathbf {k}\in \mathbb {Z}^{d}\), an m-linear, commutative m-ary operation on \(\mathcal{V}\) by setting

Now, L _δ (u) is defined recursively by

(82)

One can check that for any \(u \in \mathcal {T}\), L _δ (u)=δ(u) ∅+u+v where v a sum of rooted trees with fewer vertices than u.

Observe that, pictorially, in the right-hand side of (80), the tree associated with the right factor η ranges in the set of trees obtained by pruning some edges in the tree in the left-hand side (‘pruning’ includes complete uprooting, as in the first term, or no pruning at all, as in the last). In each term, the trees associated with the left factor δ are precisely the parts that have been chopped off. It can be shown that (81)–(82) imply that the same pruning procedure is valid to write the formula for (δ∗η) _u for arbitrary \(u \in \mathcal {T}\).

Reprints and permissions