On systems of differential equations with extrinsic oscillation

We present a numerical scheme for an efficient discretization of nonlinear systems 
 of differential equations subject to highly oscillatory perturbations. This 
 method is superior to standard ODE numerical solvers in the presence 
 of high frequency forcing terms, and is based on asymptotic expansions 
 of the solution in inverse powers of the oscillatory parameter 
 $\omega$, featuring modulated Fourier series in the expansion 
coefficients. Analysis of numerical stability and numerical examples are included.

1. Introduction. In this paper we are concerned with systems of ordinary differential equations (ODEs) of the form y ′ (t) = h(y(t)) + g ω (t)f (y(t)), where y(t) : R → R d , f (y), h(y) : R d → R d are analytic functions, and the scalar term g ω (t) can be expressed as a modulated Fourier expansion (MFE), that is Observe that we allow the coefficients a m (t) to depend on t. Therefore it is possible to think of the function g ω (t) as periodic in the variable ωt, but not necessarily in t, which allows us to consider a wide range of different perturbations. Modulated Fourier series have been used in the field of Geometric integration and numerical analysis of Hamiltonian and oscillatory ODEs, see [7,6,15], and also in connection with the Heterogeneous Multiscale Method, [24]. This is a very general framework for ODEs with extrinsic oscillation, and we note two particular cases which are of importance: firstly, if the function f (y) is An alternative is given by perturbation theory, albeit not in a completely standard form. The setting (1) belongs to the general framework of the so-called parametric perturbation (or parametric modulation) of differential equations. In the mathematical literature this has been a recurring topic in the field of perturbation theory, see for instance the classical reference [3], [25] or [17]. However, most standard methods, such as averaging, are developed for systems where the perturbation is multiplied by a small parameter ε. In this context, the general idea is that the solution of the unperturbed system plus corrections in powers of ε may yield a good approximation to the solution of the perturbed problem.
In our case the perturbation is not necessarily small, but a somewhat related idea can be applied. Intuitively, for large values of ω, the exceedingly fast oscillations of g ω (t) will produce massive cancellation between positive and negative parts of the forcing term g ω (t)f (y(t)), and will therefore give a small contribution, roughly speaking. More rigorously, if we consider the unperturbed system z ′ (t) = h(z(t)), z(0) = y 0 , then nonlinear variation of constants [16] allows us to connect the two solutions where Φ solves the so-called variational equation The matrix Φ may not be analytically available in general, but the important fact is that if the integrand is smooth enough and g ω (t) is a trigonometric function (see the examples cited before), then integration by parts gives t 0 Φ(t − s)f (y(s))g ω (s)ds = O(1/ω), ω → ∞.
This motivates the fundamental ansatz that we propose later on, that is, that the solution y(t) admits an expansion in inverse powers of the oscillatory parameter ω. The terms in this expansion can be computed explicitly in a quite general setting, and they adopt the form of modulated Fourier expansions themselves. The expansion can be seen therefore as a correction (in inverse powers of ω this time) of the solution of the unperturbed system.
We note that another possibility is presented in [8,9], namely the use of Filon quadrature (plus waveform relaxation) to approximate efficiently the oscillatory integral present in (5). This method is numerically effective, but it has the disadvantage of time-stepping, which can be expensive. As we shall see, if one expands the solution in inverse powers of ω with the right ansatz , the terms are computed beforehand and then the desired value of t can be substituted therein, without subdividing the interval of integration.
We remark that one cannot expect this approach to be satisfactory in general when the system exhibits chaotic behaviour, for example. It is known that a parametric perturbation can be used to take such a system into or out of chaotic behaviour, or to alter the chaotic states of the original problem, see for instance [26], but the very high sensitivity of the problem to changes in the data will typically render useless any method based on perturbation theory. The paper is organised as follows: in Sections 2.1 and 2.2 we present the basic features of our method and how to construct the expansion explicitly. Next, in Section 3 we analyse the bandwidth of the different terms in the expansion, depending on the original bandwidth of the forcing term g ω (t). In Section 4 we pay special attention to the stability properties of the algorithm, which essentially depends, as expected, on the behaviour of the linearisation of the system around the nonoscillatory base function. We conclude the paper with several examples to illustrate the performance of the method.
2.1. Construction. Our ansatz when constructing asymptotic-numerical solvers is that the solution y(t) admits an asymptotic expansion in inverse powers of the oscillatory parameter ω: where ψ s (t) depend on ω, but ψ s (t) = O(1), ω ≫ 1, for s ∈ Z + . Given the structure of the original ODE, a reasonable assumption is that each ψ s (t) in (6), except when s = 0, has itself the form of a modulated Fourier expansion: Furthermore, ψ 0 (t) = p 0,0 (t), i.e. p 0,m (t) ≡ 0 when m = 0, since there are no positive powers of ω present in the equation. In order to satisfy the initial condition of the differential equation, we impose ψ 0 (0) = y(0) = y 0 , which means that ψ s (0) = 0 for s ≥ 1, or equivalently Thus our ansatz reads Before we continue, we point out the difference between this ansatz and the one proposed for second order ODEs with forcing terms in [10]. That problem clearly is a particular case of (1), and the reason why the structure of the solution is apparently different is that in that case we are only interested in the first component of the solution, and then in the first level p 1,m (t) should be 0 when m = 0 in order to avoid positive powers of ω when taking the second derivative. Once we work with all the components of y(t), the ansatz should be the one presented before. In [10], this is consistent with the fact that y(t) shows small oscillations with amplitude O(ω −2 ) whereas in the case of y ′ (t) this is O(ω −1 ). See also the structure of the term ψ 1 (t) in the first example later on.
Differentiation term by term in (8) gives, formally, Since h, f : R d → R d are analytic functions, we can expand them into Taylor series around ψ 0 , assuming that this term will give the main contribution: Here h n is an n-tensor related to the n-th derivative of h at ψ 0 , etc. In general we can write Note that each h n (ψ 0 , θ, . . . , θ) is linear in each of the θ k s. It is clear that if the functions h and f are not analytic, but say C r , then we can still consider the first few terms of the expansion and adapt everything accordingly.
In the sequel, we will use the following notation: where and with the standard notation for multi-indices |k| = k 1 + k 2 + . . . + k n . A similar formula applies to the function f . Furthermore, we can express this expansion in the following way: k∈In,s l∈Kn,m h n (p 0,0 , p k1,l1 , · · · , p kn,ln ), where K n,m = {(l 1 , . . . , l n ) ∈ Z n : |l| = m}.
Observe that for simplicity of notation, we have omitted the dependence on t of the different terms p s,m in the expansion.
Now we identify coefficients in two levels: first we consider orders of magnitude (inverse powers of ω), and then frequencies (values of m) within each order of magnitude. The first level (corresponding to s = 0) is clear: Separation of frequencies yields a differential equation for the coefficient p 0,0 (t): together with the initial condition p 0,0 (0) = y(0) = y 0 , and additionally Observe that the ODE for p 0,0 is nonoscillatory. Hence, even if it is not solvable explicitly (due to the nonlinear terms h and f ), it can be efficiently solved using standard numerical methods. Note also that the components p 1,m (t) depend on t in general, even if the coefficients a m (t) are independent of t (that is, even if the forcing term g ω (t) is periodic and has a classical Fourier expansion).
For s ≥ 1 and m ∈ Z, we have where and similarly with b s,m [f ]. Note that once again we have omitted the dependence on t for brevity. Observe that for s ≥ 1 we have on the one hand the ODE with initial condition since we have imposed that ψ s (0) = 0, and on the other hand the recursion for m = 0. This is the pattern that we will find for each value of s ≥ 1. From a computational perspective, an alternative to this scheme would be to solve a system of differential-algebraic equations (DAEs) involving all the terms p s,m (t) up to the desired value of s. In this context, the analysis of bandwidth presented in Section 3 is relevant to determine the number of terms that we need within each level.

2.2.
Explicit form of the asymptotic expansion. In this section we derive explicitly the first few terms of our asymptotic expansion. We first note that there is an important simplification in the above formulas when the number of a r (t) terms that are different from 0 is finite (in other words when the input function g ω (t) is band limited). As an example, we illustrate the case where the perturbation is of the form g ω (t) = µ sin ωt. In that case it is clear that a −1 (t) = iµ/2, a 1 (t) = −iµ/2 and a m (t) ≡ 0 if |m| = 1. Thus, the original bandwidth is ̺ = 1. Similar results hold for the case g ω (t) = µ cos ωt.
2.2.1. The zeroth term. As explained before, the term p 0,0 (t) obeys the following differential equation: , together with the initial condition p 0,0 (0) = y(0) = y 0 . Equating terms with the same frequency, we obtain additionally With g ω (t) = µ sin ωt, the previous ODE can be reduced to and also together with p 1,m (t) ≡ 0 when |m| > 1.

2.2.2.
The first term. When s = 1 we obtain Here and similarly for f . Additionally we have the initial condition which follows from ψ 1 (0) = 0. Furthermore, we get When g ω (t) = µ sin ωt, then the differential equation for p 1,0 (t) (23) reads because of the symmetry of the coefficients p 1,−1 and p 1,1 . Since we impose ψ 1 (0) = 0, the initial condition for p 1,0 (t) is Putting all the information together, and taking into account (22), we obtain Furthermore, we can compute p 2,m for m = 0 from (24). According to Section 3, the predicted bandwidth is equal to 2̺ = 2, and therefore we only need to compute the coefficients p 2,±1 and p 2,±2 . Thus Because of symmetry of the p 1,m (t) coefficients, we deduce that and clearly p 2,−2 (t) = p 2,2 (t). Therefore, we have The equation and initial conditions for p 2,0 are obtained when analysing the O(ω −2 ) terms.

2.2.3.
The second term. When s = 2 we obtain Here and similarly for f , together with and When g ω (t) = µ sin ωt then we have the following ODE for p 2,0 (t), in accordance with (27): where we have used a can be simplified using the symmetry of the coefficients: Furthermore, we have the initial condition Equation (28) gives the coefficients p 3,m (t). Notice that now the bandwidth is 3̺ = 3, so we only need these terms when |m| ≤ 3.
It is clear that the process can be iterated, at the price of increasingly cumbersome expressions. Nevertheless, we observe that in many relevant cases the functions h and f are quite simple, for example multivariate polynomials of low degree, and hence many terms involving high order derivatives vanish identically. We omit any further steps for brevity.

Bandwidth and blossoming.
There is an important observation to make regarding the bandwidth of the different functions ψ s (t). Due to the nonlinearity of the original equation, the number of nonzero frequencies (i.e. values of m) increases as we move to higher values of s. We call this phenomenon blossoming, and it is relevant when discussing efficiency issues, since it quantifies the number of terms that we need to compute in each step of the algorithm.
Let us suppose that there exists ̺ ∈ N such that a m (t) ≡ 0 when |m| ≥ ̺ + 1. In other words, we suppose that the forcing term g ω (t) is band limited. Let θ s be the maximum bandwidth of the term ψ s (t), that is θ s := max{m ∈ Z : p s,|m| = 0}.
It can be checked that the first few values are For the general case, we can prove the following result: Theorem 3.1. For s ≥ 0, the maximum bandwidth θ s of the term ψ s (t) is Proof. We shall use formula (21) We note that differentiation does not alter the bandwidth, and also that the is different from 0 in general, since for example it contains the element f s (p 0,0 , p 1,̺ , . . . , p 1,̺ ).
If we try m > (s + 1)̺, then m − r > s̺, and according to the result given before we get b s,m−r [f ](t) ≡ 0. A similar argument can be used with negative frequencies, setting m = −(s + 1)̺ and r = −̺. Therefore, the maximum total bandwidth of ψ s+1 (t, ω) is indeed (s + 1)̺.
It is important to observe that this result corresponds to the maximal possible bandwidth, and that the actual one could well be smaller. For instance, if the function f is constant (which corresponds to the case of a system of ODEs with an oscillatory forcing term), then b s,m [f ](t) is identically 0 for s ≥ 1. Therefore, we would have The corresponding result is the following: Proof. We shall use formula (30) and proceed along the lines of the proof of Theorem 4.1 in [10]. We recall that we have a combination of terms of the form By induction, with equality when n 2 = 0, which corresponds to the case of the multi-index { (1, 1, . . . , 1)}.
Further reduction of the actual bandwidth can occur in some situations, particularly if higher order derivatives of h and f vanish identically. This is the case if these functions are (multivariate) polynomials, which is typical in many applications. Compare with the analysis in [10], where blossoming appears at a slower rate (in the function y(t)), though consistently with the results presented above. 4. Stability. As a direct consequence of the construction of the method in the previous subsections, one obtains that y(t) − ψ 0 (t) = O(1/ω), where y(t) is the solution of the perturbed system and ψ 0 (t) is the solution of the unperturbed one.
More precisely, suppose that y(t) is a solution of the perturbed system and consider the unperturbed system Then, if we write y(t) = ψ 0 (t) + w(t), assuming that w(t) = O(1/ω) for large ω, we obtain to leading order where A(t) and B(t) are respectively the Jacobian matrices of h and f evaluated at ψ 0 (t). This is nothing but a local linearisation of the solution around ψ 0 . Now we split the matrices A(t) + B(t)g ω (t) = U + V (t), where U is a constant matrix, and compare (31) with the system with trivial solution z(t) ≡ 0. Using standard variation of constants it follows that where Φ(t) = e tU is the fundamental matrix of the system and F (s) = V (s)w(s) + g ω (s)f (ψ 0 (s)).
It is clear from (32) that w(t) represents a deviation from the zero solution z(t), and therefore the behaviour of w(t) as t > 0 is related to the stability of this zero solution. This in turn is governed by the eigenvalues of the fundamental matrix U , see for instance [25,Ch. 6]. More explicitly, we can state the following result. • all the eigenvalues of the matrix U , say λ k , 1 ≤ k ≤ d, satisfy that ℜλ k ≤ 0, and those eigenvalues with zero real part are simple, and • it is true that for any t > 0 there exist constants c 1 , c 2 > 0 such that then the zero solution solution z(t) is stable in the sense of Lyapunov, and w(t) is bounded.
Proof. The proof follows the lines of the one given in [25, Th. 6.2]. Since Because of the first condition, we have Φ < c 3 , for a certain constant c 3 , and therefore Applying Gronwall's inequality, we get It is possible to obtain a stronger result if we impose strict negativity of the real part of all the eigenvalues λ k of the matrix U . In that case, one has the bound for suitable positive constants c 4 and ν, see [25]. Thus the zero solution z(t) is asymptotically stable and w(t) will tend exponentially fast to 0 with t. This fact, together with large values of ω will make our method very effective.
In the particular case where the system has a linear part, that is for some constant matrix M ∈ C n×n , then clearly M = U and stability is determined by the eigenvalues of the linear part of the system. If the conditions of the theorem are not met, we expect the solution w(t) to grow unboundedly in t, and then the difference between y(t) and ψ 0 (t) can be very large. A typical example of this situation occurs when the system exhibits chaotic behaviour (like in the case of the Lorenz system), and as such is very sensitive to small perturbations in the data. However, we remark that from a computational perspective, if the size of the eigenvalues is moderate, it may happen that the difference between the perturbed solution y(t) and the unperturbed one ψ 0 (t) is small enough to be acceptable when ω is large. Some examples further on will illustrate this last point.

5.
Examples. In this section we present several examples that illustrate the construction and properties of the expansion that we have presented in the previous sections. The first and third examples (where the term h(y) is constant) were already considered in [11] in an equivalent form, but here notation and computations are considerably simplified. The second example is a second order ODE, but features a nonconstant term f (y). We finally include an example with a 3 × 3 Lorenz-type system.
In all cases we will compare the approximation given by the first few terms of the asymptotic-numerical solver with the exact solution (which is either analytically available or computed numerically with standard Matlab routines up to prescribed accuracy). We will normally use the standard ODE solver ode45 in Matlab, with absolute and relative tolerance equal to 10 −12 .
We use the notation for the errors, taken componentwise. We stress that the values of ω that we use are much smaller than the ones normally present in applications. This restriction is essentially imposed by the fact that the comparison with the exact solution should be reliable and affordable. Increasing ω will benefit the asymptotic-numerical solver, since the approximation with a fixed number of terms will be more accurate, and the computational cost will be roughly similar. 5.1. A linear system. As a first example, we consider a forced harmonic oscillator with damping. This system is modelled by a simple second order ODE: x ′′ (t) + bx ′ (t) + kx(t) = µ cos ωt, where b is the damping coefficient, k the spring constant and we have set the mass m to 1 for simplicity. We introduce a forcing term with amplitude µ and frequency ω, and we suppose that ω ≫ ω 0 , where ω 0 is the natural frequency of the unperturbed oscillator in the underdamped case. In matrix form: thereby using our notation In this case, since the system is already linear, the matrix U is directly with eigenvalues Since both b, k > 0, the real part of both eigenvalues is always negative and we have asymptotic stability according to Section 4. The construction of the asymptotic expansion is particularly simple in this case, since we have after brief computation: Actually, because of the function f is constant, we have from (19) and (20)  In particular, this means that p ′ 1,0 = U p 1,0 , p s,0 (0) = 0, which leads to p 1,0 ≡ 0. Hence we conclude that in this case the first term is Note that the O(1/ω) term is 0 for the first component of the solution. In other words, x(t) is superimposed with tiny oscillations of amplitude O(1/ω 2 ), whereas in the case of the derivative x ′ (t) these oscillations are of order O(1/ω). This is intuitively consistent with what can be observed in Figure 1, and also with the ansatz taken in [10]. Analogously, and since the bandwidth in this example is θ 2 = 1, see Section 3, then we know that p 2,m ≡ 0 if |m| > 1. Furthermore hence ψ 2 (t) = p 2,0 + µ −1 b cos ωt. In Figures 1 and 2 we plot the solution of the perturbed system with parameters k = 4.2, b = 0.6, µ = 0.8, ω = 100 and initial values x(0) = x ′ (0) = 1/2, and the errors when we compare with the first terms of the approximation, respectively.

5.2.
A model for an injection-locked frequency divider. Next we consider a nonlinear example, given by the following system, which is used in [21] and [2]: and A, V dd , C, L and R are parameters of the system. A periodic perturbation can be introduced as follows: After normalisation and scaling, this system can be written as where α = L R 2 C , β = LA RC , µ = LB RC , Φ(t) = β + µ sin ωt, ω = ΩL R Hence, using our notation, we have together with g ω (t) = µ sin ωt, and thus a −1 (t) = iµ/2 and a 1 (t) = −iµ/2. The first term of the expansion, following (15), solves the system p ′ 0,0 = h(p 0,0 ), since a 0 ≡ 0, and we also get from (16) together with p 1,m ≡ 0 when |m| > 1. From (23), the coefficient p 1,0 satisfies the ODE , because of the parity of the coefficients p 1,±1 . Here Furthermore, the initial condition is Putting everything together, we have ψ 1 (t) = p 1,0 − µf (p 0,0 ) cos ωt.
The term ψ 2 (t) can be obtained from (26) substituting all the data corresponding to the functions h and f in this case.
We will vary the parameters C and Ω in the examples, but C is always positive and small and Ω large enough so that the scaled frequency ω = ΩL/R is large. We observe that the matrix U in this case is with eigenvalues With the parameters given (in particular because L ≪ 1), the term (β + 1) 2 − 4α is negative, and we have complex eigenvalues. This is consistent with the observed oscillatory behaviour of the solutions, see Figure 3. The real part of the eigenvalues is given by 1 2 (β − 1), and for stability we need β − 1 < 0, see Section 4. Thus which holds roughly when C ≥ 10 −7 .
In Figure 3 we plot the solution of the unperturbed system and the absolute errors when using the first terms of the expansions with C = 10 −6 , Ω = 2π × 10 6 (which gives ω = 451.73 after scaling) and initial values u(0) = v(0) = 1/2. Similar results, with smaller errors, are obtained if we consider larger values of Ω. The  Top row, errors in u(t) using the zeroth term (left), using up to the first term of the approximation (centre) and using up to the second term of the approximation (right). Bottom row, same for v(t). systems of ODEs for the coefficients p 0,0 , p 1,0 and p 2,0 can be solved by standard numerical methods, since they are not oscillatory. In the example, the Matlab routine ode45 was used.
If we use C = 10 −8 then the eigenvalues of the matrix U are λ ± = 8.49 ± 47.53i, therefore we do not expect numerical stability in this case. Indeed, Figure 5 shows that no significant improvement is obtained when adding more terms in the expansion. However, it is worth noticing that the errors are quite small on the whole interval of integration, so the numerical solution might be acceptable, depending on the required accuracy.

5.
3. An expcos oscillator. A more complicated example features a linear part plus a perturbation of the form g ω (t) = e µ cos ωt . For example, in the modelling of diode circuits with inductive loads, we would find an equation of the form where L, R, C, I s and V T are parameters. We will take the values L = 10 −4 , R = 100, C = 10 −6 , I s = 10 −12 and V T = 0.0259. The forcing term is g ω (t) = µ cos ωt, with large ω. The constant term −I s L/C can be added in the O(1) level in a straightforward way. Thus the resulting system is = −L/RC −L/C 1 0 where β = I s L/C. The properties of similar types of oscillator have been analysed in [8] and [9]. The relevant fact is that the function g ω (t) can be expanded in Fourier series using modified Bessel functions, see [1, Eq. 9.6.34]: and the asymptotic behaviour of the modified Bessel functions for large orders guarantees convergence for fixed values of µ and t. It is clear that the coefficients are a m = I m (µ) for m ∈ Z, using the fact that for integer orders I m (µ) = I −m (µ), see [1, Eq. 9.6.6]. The base equation follows from (15): and also The differential equation for p 0,0 cannot be analytically solved because of the nonlinearity coming from the function f , but it is nonoscillatory and therefore amenable to numerical solution using standard methods.
The differential equation for p 1,0 apparently involves an infinite number of terms: the last sum being 0 because of the symmetry of the modified Bessel functions with respect to the order. Hence: which implies that p 1,0 ≡ 0. Therefore, Observe that this last sum converges as well due to the decay of the modified Bessel functions, so its numerical implementation is not problematic.
The second term ψ 2 (t) can be computed from the general setting, although the different sums involving modified Bessel functions that appear can be quite expensive to evaluate. In particular, observe that µ = 1/V T , so if V T = 0.0259 then µ is moderately large and the convergence of the series (36) can be slow. As a compensation, we note that our expansion has two important advantages, namely that its cost is essentially independent of ω (whereas any standard numerical method will need to reduce the stepsize considerably when ω grows), and that large values of ω will yield a more accurate expansion with the same number of terms.

5.4.
A Lorenz-type system. A final example is provided by the 3×3 Lorenz-type system, see [4]: where σ, ρ and β are given parameters. Usual values of the parameters are σ = 10, ρ = 28 and β = 8/3, for which the system exhibits chaotic behaviour and develops a strange attractor.
We note that in this particular example p 2,2 ≡ 0, because of the structure of f and the matrix J[f ].
In Figures 6 and 7 we show the errors when taking the zeroth, the zeroth plus first and the zeroth plus first plus second term of the asymptotic solver to approximate the solution of the perturbed Lorenz system. The systems of ODEs for p 0,0 and p 2,0 have been solved numerically using the ode45 routine. We have taken σ = 10, β = 8/3, ω = 500 and k = 0.1. Figure 6 illustrates the case ρ = 8 and Figure 7 the case ρ = 28. This last value is typical of chaotic behaviour, and consequently, the approximation based on this perturbation method is not that satisfactory. We note, however, that even if the error builds up when t grows, see Figure 7, these errors may be acceptable for small values of t.  Figure 7. Absolute errors in the approximation of the solution of the Lorenz system with σ = 10, ρ = 28, β = 8/3, ω = 500 and k = 0.1. The first row represents (from left to right) the error up to zeroth, the first and the second approximant for x(t). The second and third rows show the same result for y(t) and z(t).