Multiple timing and spatial scaling for bifurcations

Scaling time and spatial domains often seem natural in ODEs, but unfortunately, its useful explicit form does not always agree with intuition. After reviewing the well-known methods, multiple timing and averaging, we will show that algebraic timelike variables may play a part in bifurcations. Part of this discussion is tied in with the problem of structural stability in the analysis of matrices, and another part is determined by bifurcation theory of nonlinear systems. Secondly, the theory of resonance manifolds for higher dimensional problems involving at least 2 angles will show the presence of unexpected small spatial domains that may emerge involving long timescales and containing interesting phenomena. A number of examples from mechanics are presented to demonstrate the theory.


Introduction
The theme of this paper is to discuss unexpected timelike and spatial variables in differential equations. Consider ordinary differential equations (ODEs) that contain a small positive parameter ε: F. Verhulst depending to some order smoothly on x and t for t 0 ≤ t < ∞ and on the parameter ε for 0 ≤ ε ≤ ε 0 and x ∈ D, D ⊂ R n ; the dot represents differentiation with respect to t. The smoothness implies that we can write the right-hand side as f (t, x, ε) = f (t, x, 0) + O(ε). A well-known example is the pendulum with oscillating support as displayed in Fig. 1.
A much more complicated problem is that of the rotating flywheel on a vibrating foundation displayed in Fig. 2. This problem is discussed in Sect. 8.
There are many more physical examples of such problems, see for instance [18], also [11] and [12].
In modelling the pendulum system of Fig. 1, the Equation of motion for the angle θ with the vertical yields after linearisation near θ = 0: x + (ω 2 + ε cos νt)x = 0. (2) Studying differential equations, in particular initial value problems, it seems natural to assume the presence of timelike variables t, εt, ε 2 t, . . . on which approximate solutions depend. In other mechanical problems, we have similar choices by scaling of spatial variables; as an example, the rotating flywheel problem is displayed in Fig. 2, it is discussed in Sect. 8. Contrasting with the approach of apriori guessing timelike and spatial variables are several other methods like averaging or renormalisation, normal form methods, where no apriori assumptions on the form of time-dependence or hidden spatial scales are made. Fig. 1 The pendulum has a vertical, harmonically oscillating support P, and it is described in Sect. 5 φ x excentric mass Fig. 2 A rotating flywheel has a small eccentric mass and is mounted on a spring. It can move through resonance but it can also move into a resonance domain and be locked. The analysis of the model can be found in Sect. 8 The multiple timing idea of assuming the presence of timelike variables was given already by Krylov and Bogoliubov in 1935 [13]; an application can be found in a paper by Kuzmak in 1959 [14]. After that, according to [15], the Kiev school of mathematics had no interest in multiple timing.
Much later multiple timing was studied in the 60 s and 70 s in [9], [4] and for instance [18]. In these papers, there is no reference to [13]. Spatial scaling in problems that did not present itself as singular perturbations with boundary layers came later. Small domains with different qualitative behaviour can be hidden in x-space D of equations like (1).
The scientific literature is rich on papers on the approximation of solutions of ODEs like (1), we can cite only a few of them. A critical comparison of aver-aging and multiple timing by a number of important examples can be found in [10]; this was the first paper to show weak points of multiple timing. In [30], the relation between averaging, multiple timing and the renormalisation method was discussed following [2], [3] and [16]. In [19], the asymptotic equivalence of the averaging method and multiple timing at first order in ε was established for standard variational equations likė with O(ε) error estimate on intervals of time of order 1/ε. See also the extensive discussions in [17] and [23].
The following concepts and description are based on [23] ch. 1. Asymptotic equivalence of methods would imply that, considering a solution of a differential equation x(t), expressionsx 1 (t) andx 2 (t) obtained by different methods would both represent an approximation of x(t) with error δ(ε) = o(1) as ε → 0 on the same interval of time (for instance of size 1/ε). Asymptotic expansions are not unique;x 1 (t) andx 2 (t) may be different but both acceptable approximations.
Following [23], we will indicate that an approximation with error δ(ε) is valid on an interval of size 1/ε. A more precise statement is that the error estimate is valid for t 0 ≤ εt ≤ t 0 + L with t 0 , L constants independent of ε.

Set-up of the paper
To anticipate simple timelike variables as t and εt is not a bad idea in the cases that the set of solutions, or if you wish the dynamics of the problem is studied in a structural stable setting. We mean by this that there are no qualitative changes in the behaviour of the solutions; later we will be more precise about structural stability. However, if there is a qualitative change in the dynamics, we may find the presence of unusual or unexpected timelike variables and spatial domains.
A serious point is then that in research, we will be especially interested in structural changes like the emergence of periodic solutions, stability changes, the presence of tipping points in dynamics, etc. Anticipating the timelike variables of a problem like εt, ε 2 t, etc. conflicts with having an open mind about the possible outcome of research.
In Sects. 2 and 3, we will review multiple timing and averaging, and we explain the need and presence of algebraic timelike variables like ε q t with q a rational number in Sect. 4. Such variables arise in stability problems after linearisation at an equilibrium and in standard bifurcations that are found in many applications. In Sects. 4 and 5, we show that structural stability problems of matrices produce already unexpected timelike variables in linear ODEs.
Analysing nonlinear perturbation problems, we have to use bifurcation theory, when discussing for instance the Van der Pol-equation the use of Hopf bifurcation is natural. Bifurcation theory is a large field, we discuss only a few prominent cases.
A different topic starts with Sect. 7. In equations like (1) describing oscillatory behaviour, there may arise local resonance manifolds in x-space that are characterised by small spatial size and unexpected timelike variables. The presence of resonance manifolds is not obvious and requires analysis. These problems are tied in with passage through resonance and capture into resonance. Interestingly, such problems are found in both conservative and dissipative problems. In Sect. 8, we consider dissipative ODEs, in Sect. 9 Hamiltonian systems. A few examples show that we have to develop fairly high order approximations to characterise the dynamics. A short discussion of other methods and some conclusions are given in Sect. 11.

The multiple timescale method
Many small ε parameter problems are studied using timescales like t, εt, ε 2 t and in general ε n t with n ∈ N. A simple but typical form of multiple timing runs as follows. Consider the equatioṅ with f (t, x) T -periodic in t, the initial value x(0) is given. We will look for solutions of the form with τ = εt, the dots represent the higher order expansion terms. As the unknown functions x 0 , x 1 , . . . are supposed to depend on two variables, we have to transform the differential operator; we have to first order in ε: Using this differential operator and the expansion, we find Suppose we can Taylor-expand the function f to a certain order, collecting then the terms of order 1 and ε, we find the simple partial differential equations The first equation produces with A(τ ) still an unknown function; A will be determined in the next step. For x 1 we find by integration The function B(τ ) is unknown and has to satisfy B(0) = 0. If we are looking for bounded solutions of Eq. (3), or even for periodic solutions, the integral has to be bounded. This is called the secularity condition. We can achieve this by determining A(τ ) such that Assuming that f (t, x) has a Fourier expansion, this is a natural condition as it means that the 'constant' term of the expansion vanishes. The determination of A(τ ) implies that satisfying the secularity condition corresponds with averaging the function f (t, x). The idea of secularity conditions can be traced to the end of the 18th century, for instance in the writings of Lagrange and Laplace (see [23]).

Example 1 Consider the Van der Pol-equation
x We will look for solutions of the form (4) in t and τ = εt; this leads to the well-known first-order result: x 0 (t, τ ) = re If we have initially r = 2, the first-order approximation is periodic. It has been shown that this first-order term represents an O(ε) asymptotic approximation, valid on the timescale 1/ε; see for instance [23] or [27]. Multiple timing and averaging (see next section) produce the same first-order approximations in this problem.

The origin of timelike variables t and εt
Consider Eq. (1) in the expanded forṁ As this is supposed to be a perturbation problem, we should be able to solve the 'unperturbed' probleṁ to obtain the 'unperturbed' solution where C is an n-dimensional constant of integration. Apply variation of constants (Lagrange) by putting For y we obtain the equation: a so-called variational or slowly varying system. Note that variation of constants is in general easier to apply if the unperturbed problem (6) is linear. Averaging or using the multiple timescale method produces a transformed (normal form) equation: a slowly varying equation forȳ. We have transformed x → y →ȳ without giving the details of the process and until this point, no approximation has been applied.
Omitting the O(ε 2 ) terms to start the approximation process, and dividing the equation forȳ by ε, we note that the 'natural' first-order timelike variable forȳ is τ = εt. Because of the transformation x → y, t is the zeroorder time variable for the original perturbation problem in x and so we have the variables t, εt. The only assumptions until now are smoothness of the vector functions on a suitable domain and the possibility of inversion of the variation of constants relations. We conclude that at first order, t and εt are natural variables.

Averaging techniques
To explain the emergence of unexpected timelike variables, it is necessary to discuss briefly second-order averaging and so we have to begin with first order. Suppose that the variational system (7) has a right-hand side that is T -periodic in t. We will consider the averaged vector field where we keep y fixed during integration. Omitting the O(ε 2 ) terms in system (8), we have the approximating systeṁ We can derive the error estimate that if y(0) =ȳ(0) we have |y(t) −ȳ(t)| = O(ε) on the timescale 1/ε. For a proof see [23] or [26].
To obtain a second-order approximation is much more work, we will present the idea without all the details. Consider eq. (7) in the form: Put h 0 (y) for the average of h(t, y) over t, in general we use the superscript 0 for an averaged vector field. We introduce the n × n Jacobian matrix Dg(y, t) (differentiation with respect to y only) and the vector field We compute the vector field with average F 10 . Consider the equation: then the expression w(t)+εu 1 (t, w) approximates y(t) with error O(ε 2 ) on the timescale 1/ε. For a proof and discussions see [23]. We gave explicitly these expressions to show that we made no assumptions on relevant timelike variables. An O(ε) approximation produces timelike variable εt; a subsequent timelike variable will be introduced by solving the variational Eq. (11). We will see examples later.
behaviour of the dynamical system studied. If the system at hand experiences qualitative changes when the parameters of the system pass certain critical values, we call them bifurcations. This can involve many causes; it can be stability changes, emergence or destruction of solutions, transitions to tori, chaos in its different forms and other phenomena.
As in many research problems the unperturbed problem (6) is linear, we need basic results from matrix theory. In addition, when we will study a special solution like an equilibrium of Eqs. (10) or (11), we will usually linearise near this solution, for instance to determine stability. This also asks for matrix theory.
A second cause of qualitative changes is the bifurcations of a nonlinear part of the vector field when a parameter varies. Standard cases are co-dimension 1 bifurcations like for instance pitchfork and saddlenode; see Sect. 6.
A different cause of qualitative changes will be if in phase-space we find local behaviour as encountered in boundary layers that is very different from the global behaviour. The occurrence of such local changes can be quite unexpected, see Sect. 7.
All this knowledge will help to avoid making incorrect apriori assumptions on timelike and spatial variables.

Structural stability of matrices
Part of the (classical) material in this section can be found in [27], in particular appendix 15.3. Before formulating the theory, we consider as an introduction the Mathieu equation (2) in its fundamental 1 : 2resonance with a slight detuning of the frequency ω. The equation models the pendulum motion with oscillating support shown in Fig. 1. Near the vertical axis linearisation and replacing θ by x leads to the equation: with ω 2 = 1 + εa + ε 2 b; a and b are free parameters independent of ε. See also Fig. 1.
We summarise the second-order approximation analysis in [27]. We apply Lagrange variation of constants x,ẋ → y 1 , y 2 to Eq. (12): x = y 1 cos t + y 2 sin t,ẋ = −y 1 sin t + y 2 cos t. The slowly varying equations for y = (y 1 , y 2 ) are after averaging of the formẏ = A(ε)y; The trivial solution is stable if |a| > 1/2, unstable if |a| < 1/2. For a = ±1/2, we have a first-order approximation of the curves separating stability and instability domains, see Fig. 3. The matrix A(ε) is singular if a = ±1/2. The Floquet tongues are bounded by the bifurcation curves where the transition from unstable to stable solutions takes place in (ω 2 , ε)-parameter space. What happens at the tongue boundary, for instance if ω 2 = 1 + εa with a = 1/2 ? In this case we have to first order: , a typical degenerate matrix from bifurcation theory. Following [23] or [27] we perform second-order averaging following Sect. 3.1 to find as perturbation of A 1 : The disc rotates with constant frequency , its foundation produces at point P small vibrations of the form ε cos ω 0 t We find for the eigenvalues of A(ε) to this second order of approximation produces a more precise location of the Floquet tongue.
Near the boundary of the Floquet tongue we have that λ 2 = O(ε 3 ); ii is remarkable that the timescale ε 3 2 t plays a part in this problem. The timescales characterising the flow near the Floquet tongue are derived from second order of the eigenvalues: The presence of the timelike variable ε 3 2 t was noted for the Mathieu equation earlier in [2], using the renormalisation method.

The rotating disc
A heavy disc is rotating on a vertical shaft. The shaft is fixed at its suspension point P, but the centre of the disc is free to make small vibrations in the horizontal directions, see Fig. 4. The point of suspension is elastically attached to the foundation z = 0. As a first approximation, we assume that the suspension point oscillates harmonically in the vertical direction. Following [20], the equations of motion are as follows: with α dependent on the inertial moments and inverse proportional to the rotational speed , for the vibration of the foundation we have a harmonic function. Damping is added with positive coefficient κ.
In [20], the case κ = 0, no damping, is analysed first. The frequencies in the case ε = 0 are as follows: We have a so-called sum resonance if ω 1 + ω 2 = 2 √ 1 + α 2 = 2η. Introducing detuning as in Eq. (12), we obtain the same timelike variables characterising the dynamics. For stability boundaries for the position on the vertical axis, we find to first order in ε: For the standard calculations on Mathieu equations see also [26]. Mechanical rotation effects in combination with the parametric resonance produces with κ = 0 an instability tongue depicted in Fig. 5. A remarkable result is found if we add small damping; if κ > 0, we find by eigenvalue analysis of the matrices the stability boundaries to order ε: As this result is valid for arbitrary positive κ, the resulting boundaries differ essentially from the stability boundaries given by (15). The domain of instability becomes actually larger when damping is introduced. Mathematically, this phenomenon is caused by the structural instability of matrices as explained in [8]. Physically, it can be understood from the fact that damping introduces an extra coupling between the 2 degrees of freedom of the rotating disc. See also Fig. 5. Instability caused by damping is an important phenomenon in 2 or more degrees-of-freedom systems with rotating components. For introductory surveys see [11] and [12].

Remark on weak coupling and damping
It is interesting to consider the case of a weaker coupling of the rotating system by putting α → εα. One can also put κ → εκ. With these assumptions, the stability matrix requires a second-order averaging calculation, and the resulting 4 × 4 matrix contains elements with terms mixed of ε and ε 2 . To first order, we have the situation of Fig. 5 (left figure) without damping. It is not known if adding the terms O(ε 2 ) may produce qualitative and quantitative changes. We start with Eq. (1) of the formẋ = f (x, t, ε) to find an equilibrium or special solution ψ(t) and study the behaviour of this solution as the parameters are changing; we need to calculate eigenvalues, Lyapunov exponents or characteristic multipliers. As we saw in the Mathieu equation and will also see in the examples later on, bifurcation phenomena in ODEs lead by averaging and local linearisation to studying systems of the form: We assume that we can expand The n×n-matrices A n are independent of ε. Usually, we have A 0 derived from the unperturbed problem, A 1 is produced by perturbation methods, and sometimes we will have some knowledge about higher order terms. An important question is then if the eigenvalues of A 0 and A 0 + ε A 1 are in a sense typical for the eigenvalues of A(ε). This question is determined by the structural stability of the matrices and whether eigenvalues are single or multiple. Failure of structural stability and the presence of multiple eigenvalues is characteristic for bifurcation phenomena. We give a definition: A n × n matrix is called structurally stable if it is nonsingular and all eigenvalues have nonzero real part. If we have a zero eigenvalue or purely imaginary eigenvalues, we can expect bifurcations. As we shall see later on, multiple eigenvalues may affect the form of the expansions and timelike variables.
Multiple eigenvalues of the matrix A 0 or A 0 + ε A 1 may produce eigenvalues of the order ε q with q a rational number. Consequently timelike variables like ε q t play a part. We can predict the occurrence of such algebraic timescales from the actual eigenvalues. We add an example.
Example 2 Consider a system that can be obtained from linearisation near an equilibrium of a chemical reaction equation or an interacting system in population dynamics: with constants a, b > 0. The eigenvalues of the coefficient matrix are as follows: −ε, ±ε 3/2 √ ab with timelike variables εt, ε 3/2 t.

Classical results
Results for timelike variables from matrix expansions are essential for a sound analysis of our problems.
We summarise a few 19th century results referring to [27] appendix 15.3 for references. Consider the matrix expansion with A 0 structurally stable: 2. According to Newton-Puisieux: If λ 0 is multiple, fractional powers of ε are possible in the expansion of the eigenvalues.

Example 3 Newton-Puisieux expansion
Consider again the equationẋ = A(ε)x with for the matrix A(ε): A 0 is the matrix with zero elements so we consider A 1 . The characteristic equation of ε A 1 produces 4 equal eigenvalues ε. The matrix ε A 1 + ε 2 A 2 has the characteristic equation: The eigenvalues are as follows: Again we find timelike variables with a fractional exponent of ε.

Co-dimension 1 bifurcations
The theory of bifurcations is a large and well-researched topic. Bifurcations, qualitative changes in the dynamics, are often found when analysing variational equations. This leads often to interesting phenomena in a nonlinear setting where certain parameter values are identified that may cause qualitative changes We consider here the simplest but often occurring case where we have one or two parameters involved. The examples we discuss are low dimensional, so-called codimension 1 bifurcations. Such bifurcations may arise in subsystems of ODEs after first-or second-order averaging.
A simple example is the saddle-node bifurcation described by: If ab < 0, there is no critical point; suppose ab > 0, then we have the critical points x 0 = ± √ a/b. Their stability and local behaviour with time is described by the coefficient of linearisation −2bx 0 = ∓2 √ ab near the critical points. If for instance a = ε 2 , b = ε, the leading timelike variable is ε 3/2 t.
Consider the system inspired by the pitchfork bifurcation: Three critical points (equilibria) are (x 0 , y 0 ) = (0, 0), (0, ± √ ε). Linearisation near the critical points (o, y 0 )) produces: x = ε 2 y − 3εy 2 0 y + . . . ,ẏ = εx, where the dots represent the neglected nonlinear terms. We have the characteristic eigenvalue equations and timelike variables near the critical points: In general, we expect in regions where bifurcations occur and for higher dimensional problems the presence of timelike variables of the form ε q t with q a positive rational number. As we will see, an example of the pitchfork bifurcation is found for the amplitude in the Van der Polequation.

Example 4 Consider the Van der Pol-equation in the following form:
x If parameter a < 0, the oscillations will be damped, and there is no periodic solution. With a > 0, we have after first-order averaginġ If parameter a starts at a negative value and we let it increase, it will pass through zero and for a > 0 a periodic solution emerges with amplitude 2 √ a by a pitchfork bifurcation. It starts with a being small, say a = ε. This is the situation where we have the timelike variable ε 3/2 t. We improve the reasoning by writing Eq. (22) as: x + x = −εẋ x 2 + ε 2ẋ .
In amplitude-phase variables r, φ the variational system to O(ε) becomes: First-order averaging produces: Computing the quantities D f and u 1 in the notation of Sect. 3.1, we find an O(ε 2 ) contribution for the phase and the amplitude; for the amplitude we have: We find that the the amplitude of the periodic solution grows as 2 √ ε.

Resonance manifolds
To obtain variational equations can be more difficult if the unperturbed system ε = 0 of Eq. (1) is nonlinear.
In such cases, we may obtain a formulation like: In general, the Fourier expansion of X (x, φ) will contain combinations of the angles (φ 1 , . . . , φ m ). As we will see, in systems of the form (25) we have to account for the presence of resonance manifolds. Problems of this type arise both in dissipative and in Hamiltonian systems; see for instance [23] or [27] and references there. Also, in these problems, higher order algebraic timescales and asymptotically small domains are natural. We start with a simple example.
Example 5 Consider the system: The vector field X (x, φ) is scalar in this case and contains 2 angles, φ 1 and the combination angle 2φ 1 − φ 2 . Putting χ = 2φ 1 − φ 2 for the combination angle we have: As φ 1 (t) = 2t +φ 1 (0), this angle is timelike. The angle χ is timelike except in a neighbourhood of x = ±2. Averaging over the angles φ 1 and χ outside a neighbourhood of x = ±2, we find: A neighbourhood of x = ±2 will be called a resonance domain, x = 2 a resonance manifold. Note that the approximate solution starts in x = 1 and will increase To determine the size of the resonance domain and the local dynamics, we rescale: with δ(ε) → 0 as ε → 0 (δ is a small parameter to be determined). From system (26) we find: The 2 equations are balanced if δ(ε) = √ ε; in the theory of singular perturbations, this is called a significant degeneration of the system, see for the theory [27]. With this assumption, the small size of the resonance domain near x = 2 is √ ε. After omitting the higher order terms and averaging over φ 1 , we havė so that by differentiation we find the forced pendulum equation for χ in the resonance domain: To first order the timescale of the dynamics in the resonance domain and manifold is O( √ ε), the timelike variable √ εt, the error of the first approximation will also be O( √ ε), see Fig. 6. Note that this resonance domain is in a sense hidden in system (26).
We will see that the size of resonance manifolds and the timescale found in this simple problem are typical for much more complicated problems. Even for dissipative systems of the form (25), the first-order approximation of the solutions in a resonance domain will always be conservative; dissipation may be shifted to second order (Fig. 7).

Resonance manifolds in dissipative systems
A typical example from [27], example 12.11, describes a slightly eccentric flywheel, see Fig. 2. in the Introduction. The vertical displacement x of the flywheel and its rotation angle φ are given bÿ It turns out that there exists a resonance domain; the domain is of size O(ε 1 2 ), the timelike variable of the dynamics is √ εt in the resonance domain. A difference with conservative systems is the possibility of locking into resonance for various initial conditions in this mechanical problem where we have both flywheel and spring oscillating. For details, we refer to [27] and the references there.

Example 6
A simpler example is the 3-dimensional system: The function f (x) is smooth. The equation for x contains 2 angles, φ 1 and χ = 2φ 1 − φ 2 . The angle φ 1 is clearly timelike for any value of x; we can average over φ 1 . We have for χ : so χ is also timelike. Averaging over both angles we findẋ = 0 and the trivial dynamics x(t) = x(0)+O(ε). One can check this result numerically for instance for the choices f (x) = x and f (x) = sin x. For ε = 0.1, the error stays below 0.1.
A less trivial dynamics is shown in the following example.
The function f (x) is smooth, f (1) = 1. Again the angle φ 1 is clearly timelike and so is φ 2 , we can average over φ 1 , φ 2 to find the approximate equatioṅ producing O(ε) approximations outside resonance domains.
We have for χ : to find to first order the system: To balance the equations, we choose δ(ε) = √ ε leading to the first-order equation for χ : The timelike variable in the resonance domain is again √ εt, and the approximations at this order have error O( √ ε). Critical points (equilibria) of the pendulum equation Eq. (32) are found for solutions of sin χ = 0. If χ = 0, the solutions are to first order neutrally (Lyapunov) stable; for χ = π , we have instability. As the second equilibrium is a saddle point, the instability persists to all orders of ε. It is interesting the repeat the computation for ε = 0.01, see Fig. 8.

Hamiltonian resonance
The following results are based on [21], [23] and [27]. Consider the two degrees-of-freedom (dof) Hamiltonian in local coordinates with Taylor-expansion: with H k homogeneous polynomials of degree k (= 2, 3, . . .) in positions and momenta ( p, q). H 2 takes the standard form with the integers m, n positive and in most cases relative prime. The phase-flow in a neighbourhood of the origin takes place on compact manifolds parametrised by the Hamiltonian (energy) integral. Resonance domains can be found but because of the recurrence properties of the phase-flow capture into resonance is not possible. Near the origin of phase-space Hamiltonian (33) was obtained by rescaling q = εq, p = εp, dividing by ε 2 in the Hamiltonian and leaving out the bars.
Most of the attention in the literature went to the primary resonance 1 : 2 and to the secondary resonances 1 : 1 and 1 : 3, see [23]. In these resonance cases, the dominant part of the phase-flow is characterised by the timescales t, εt, ε 2 t and the time intervals of validity of approximation 1/ε and 1/ε 2 . This picture changes drastically for higher order resonance.
The higher order normal form The frequency cases where m +n ≥ 5 are called higher order resonances by definition. To study these resonances, we have to compute higher order normal forms; this involves intervals of time longer than 1/ε, even longer than 1/ε 2 . In the Hamiltonian normal form, the first resonant term, involving not only actions but also angles, arrives from H m+n at O(ε m+n−2 ).
The first basic approach to higher order resonance was given in [21] with applications in [22]. In [25], an improvement of the estimates has been given, together with a number of applications, for instance the elastic pendulum. Introducing action-angle variables the normal (averaged) form is obtained by nearidentical transformation and will look like with resonance combination angle with E 0 , E 1 constants. for the combination angle we havė To compute the termH 4 , we can use second-order averaging. The Hamiltonian and equations of motion truncated after ε 2 -terms produce an O(ε) approximation of the solutions on the timescale 1/ε 2 .
To determine the term to order O(ε m+n−2 ) is in general a lot of work.
The phase-flow of higher order resonance for 2 dof We summarise. Consider for 2 dof time-independent Hamiltonians near stable equilibrium the higher order resonances as defined by m + n ≥ 5. When averaging we usually transform to amplitude-phase variables before averaging over t. We can use instead of the phases ψ 1 , ψ 2 the timelike variables φ 1 = mt + ψ 1 , φ 2 = nt +ψ 2 . A resonance domain with resonance manifold M exists if we have solutions of: It turns out there are two domains in phase-space where the dynamics is very different and is characterised by different timescales. We have the results:  Fig. 9 The Poincaré map for the 1 : 6-resonance of the elastic pendulum (ε = 0.75, large for illustration purposes). In the resonance domain, the saddles are connected by heteroclinic cycles and inside the cycles are 6 centre fixed points, see [25]. It is shown in [25] that for Hamiltonians derived from a potential, we have α = 0, and so the combination angle χ = nφ 1 − mφ 2 . For the elastic pendulum, after the first-order 4 : 2-resonance, the higher order 4 : 1resonance is the most prominent one with resonance manifold of size O(ε 1 2 ) and time interval of interaction O(ε − 5 2 ); for a Poincaré map of the 1 : 6-resonance of the elastic pendulum see Fig. 9. It was constructed by numerical integration for fixed energy, the 2 positions q 1 , q 2 are on the axes with the map arising from the transversal points when passing recurrently the plane p 2 = 0.

A m : n toy problem
A well-known model for orbits in axi-symmetric galaxies is the family of Hénon-Heiles potential problems, see [7]. The applicability of this potential is much more general as symmetry of the potential, think of pendulum equations, occurs in many mechanical problems. A simplification was studied numerically in [5]. Using the notation of [5], the Hamiltonian is as follows: The equations of motion are as follows: x + m 2 x = εy 2 ,ÿ + n 2 y = 2εx y.
So, the x-normal mode y =ẏ = 0 is an exact harmonic solution. Using amplitude-phase variables x = r 1 cos(mt + φ 1 ),ẋ = −mr 1 sin(mt + φ 1 ), y = r 2 cos(nt + φ 2 ),ẏ = −nr 2 sin(nt + φ 2 ) we find after second-order averaging: For the combination angle χ = nφ 1 − mφ 2 , we find to O(ε 2 ): Putting the right-hand side of (42) zero and using the energy integral we find for the location of a resonance manifold M the values: So, for the existence of a resonance manifold, we have the requirement: From the Hamiltonian at higher order resonance (35), we have apart from the r 1 , r 2 values (43) the values χ = 0, π. In physical space, this leads by elimination of the goniometric expressions to polynomial relations between x(t) and y(t).
Examples are shown in Fig. 10 for the m : n = 4 : More details can be found in [22].

Discussion and conclusions
Considering the scientific literature, one observes that the use of asymptotic series to approximate solutions of differential equations takes all kind of different forms: averaging, multiple timing, harmonic balance, renormalisation, WKBJ, etc., see [30]. Averaging is the only method with explicit error estimates and intervals of validity for first-and higher order approximations. Multiple timing is, with the right conditions, correct at first order but has counterexamples for higher order. Harmonic balance is a method without any foundation or justification, see [24], ch. 9. The choice of a particular method seems to be often a matter of taste. In this respect, it is very important for well-founded research to have comparative and unifying studies as [19], [16], [17], [23] and [2], [3], to name a few.

Conclusions
• Initial value problems with a small parameter may involve timelike variables as t, εt, in general ε n t, n = 1, 2, . . . . In the interesting case of qualitative changes, tipping points and bifurcations, algebraic timelike variables may arise. This occurs already in linear problems with eigenvalues causing structural instability. • Oscillating systems where we can identify 2 or more angles can contain resonance manifolds. We have seen such problems in dissipative and in Hamiltonian systems. The quantitative description of these resonance manifolds require again higher order algebraic timescales and asymptotically small domains • It is essential to use approximation methods that do not anticipate the timelike variables that are relevant for the approximations. These variables present themselves naturally in the course of the analysis for normal form methods and averaging.
In engineering, there is a trend to consider exclusively numerical methods to solve problems. This is useful for a number of isolated practical problems where one looks for numbers only and not for theoretical insight. Think of the design of a building structure of given dimensions. See also the comments in [6].
For general theoretical insight, it is important to have the use of both analytical and numerical, computational methods. Note that validation of results by numerics is possible only for isolated cases, general validation of results and methods require mathematical analy- sis. A recent, rather simple example of this approach is [1], where Neimark-Sacker bifurcation is analysed for interaction of a self-excited and a parametrically excited oscillator. It turns out that its simple formulation is deceptive, numerical explorations show various phenomena. These results inspire the asymptotic analysis that yield values of the parameters that produce families of quasi-periodic solutions that are organised on tori surrounding periodic solutions. In its turn, the analysis suggests more numerical illustrations.
Advanced numerical methods and high speed computers are available for research problems; the following strategy can be useful: A hybrid strategy 1. Start with a few numerical explorations.

Identify parameters and use asymptotic expansions
to obtain more general insight.

Use numerical bifurcation programs like Auto and
Matcont to extend general insight. 4. Summarise the information in pictures typical for the dynamics.
request. The numerics and output of the figures were obtained by Matcont routine ode78 under Matlab.

Conflict of interest The author has no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.