## Abstract

We present a method to compute efficiently solutions of systems of ordinary differential equations (ODEs) that possess highly oscillatory forcing terms. This approach is based on asymptotic expansions in inverse powers of the oscillatory parameter, and features two fundamental advantages with respect to standard numerical ODE solvers: first, the construction of the numerical solution is more efficient when the system is highly oscillatory, and, second, the cost of the computation is essentially independent of the oscillatory parameter. Numerical examples are provided, featuring the Van der Pol and Duffing oscillators and motivated by problems in electronic engineering.

## 1. General setting

In this paper, we are concerned with second-order ordinary differential equations (ODEs) with highly oscillatory forcing terms, more precisely equations of the form
1.1
for *t*≥0, where the forcing term *g*_{ω}(*t*) can be expressed as a modulated Fourier expansion (MFE), that is
1.2

We will further assume that *R*(*y*) and *S*(*y*) are analytic, which ensures the existence and uniqueness of the solution *y*(*t*). These conditions can be weakened to cater for less smooth data if needed, and the scheme can be adapted in a straightforward way.

This setting includes some differential equations with important applications, in particular the *Van der Pol oscillator*
1.3
where *g*_{ω}(*t*) is of the form (1.2) and *μ*>0 is given. The standard *forced Van der Pol oscillator* is given by , which is clearly a special case of equation (1.3), with *α*_{−1}=−*α*_{1}=*i**A*/2.

Another important example belonging to this type of differential equation is the *Duffing oscillator*
1.4
where *k*>0 is the damping constant, and *b*>0 corresponds to the so-called hard spring case and *b*<0 to the soft spring case.

There are two examples of forcing terms of importance in electronic engineering
1.5
where *ω*_{1}≫*ω*_{2}≫1 and *c*_{1},*c*_{2}≠0 are constants. The first example represents a simple sinusoidal signal (possibly highly oscillatory), whereas the second one corresponds to an AM signal (if *c*_{1}=0, we have a double-sideband suppressed carrier AM modulation). The presence of two different frequencies is motivated by the fact that RF communication circuits are marked by the presence of signals with widely varying time scales. In particular, for modulated signals, low-frequency information (in this case given by *ω*_{2}) is superimposed on a high-frequency carrier (given by *ω*_{1}), so that aerials of practical dimensions can be employed. In addition, different signals can be modulated onto carriers of different frequencies, thereby enabling a large number of radio transmitters to transmit at the same time.

In both the Van der Pol and Duffing equations, if the forcing term *g*_{ω}(*t*) has period *T*>0, the existence of a non-constant *T*-periodic solution to the forced equation is known; see, for instance, Farkas (1994). These two equations have been extensively studied, notably in the context of singular perturbation theory, where the damping parameter is supposed to be small; see, for instance, Jordan & Smith (2007) or the classical reference Bogoliubov & Mitropolsky (1961). Both equations have been widely used as well in the modelling of electronic circuits; for instance, in Hilborn (2000), Pulch (2005) and Volos *et al.* (2007).

In this paper, we investigate the properties and computation of solutions of this type of equation when the forcing term is highly oscillatory, that is, when *ω*≫1. Similar to what happens in the general case (Condon *et al.*2009*a*–*c*), the oscillatory nature of the solution imposes a very small stepsize on standard numerical methods for ODEs, thereby rendering them exceedingly expensive.

Our approach is a combination of asymptotic and numerical techniques: asymptotic expansion in inverse powers of the oscillatory parameter *ω* provides a convenient and fast-converging representation of the solution (especially for large values of *ω*), while numerical discretization of *non-oscillatory* differential equations generates expansion coefficients in an efficient way.

In §§2 and 3, we present the construction of our asymptotic-numerical solvers, and the explicit derivation of the first few terms. As we shall see, the first few terms in our asymptotic expansion of the solution of the forced ODE preserve the bandwidth of the original input, which is important for efficiency issues. However, as we progress to higher order terms, we find that the bandwidth of the solution of the ODE increases, a phenomenon that we call *blossoming*. This is an unavoidable consequence of the nonlinearity of the differential equation, but we are able to quantify this increase in the number of frequencies in §4.

Before we commence with the theory, we display in figure 1 the limit cycle and the trajectories of the unforced (solid line) and forced (dashed line) oscillators. Looking at the limit cycle, it is clear that forcing induces oscillations. However, a closer look at the trajectories indicates that these oscillations follow a pattern. While they are hardly visible in *y*(*t*), the variable *y*′(*t*) exhibits significant oscillations. This observation is critical to our analysis.

## 2. An asymptotic-numerical solver

We seek to represent *y*(*t*) as an asymptotic series in inverse powers of the oscillatory parameter *ω*
2.1

This representation should be understood in an asymptotic sense, that is,
2.2
and in general no assumption is made about its convergence for fixed values of *ω*. Furthermore, because of the structure of the forcing term *g*_{ω}(*t*), we make the ansatz that each term *ψ*_{r}(*t*) has the form of a modulated Fourier series
2.3

As already explained in Condon *et al.* (2009*b*), MFEs provide a natural framework for solving differential equations with highly oscillatory forcing terms. This type of expansion has already been used in the context of Hamiltonian systems (see Hairer *et al.* 2006, ch. XIII; Cohen *et al.* 2003), and also as the basic ingredient in the heterogeneous multiscale method (see Sanz-Serna 2009).

It is important to observe that we need to impose *p*_{0,m}(*t*)≡0 and *p*_{1,m}(*t*)≡0 for *m*≠0, since otherwise differentiation with respect to *t* would produce positive powers of *ω*, which are not present in the original equation (1.1). For this reason, we modify the previous ansatz to the following:
2.4

Note also that in this setting the oscillations in *y*(*t*) have amplitude of order 1/*ω*^{2}, whereas in *y*′(*t*) they are of order 1/*ω*, consistent with the behaviour that can be observed in the previous example.

This approach is clearly related to perturbation theory, in the sense that we express the solution *y*(*t*) as a base function *ψ*_{0}(*t*) plus corrections in terms of inverse powers of *ω*, which will be small if *ω* is large. The difference between the perturbed and unperturbed solutions *y*(*t*)−*p*_{0,0}(*t*) can be written using variation of constants (see Hairer *et al.* 1994). For fixed values of *ω*, a standard stability analysis can be carried out using a linearization of the solution about *p*_{0,0}(*t*). We refer the reader to Condon *et al.* (2009*b*) for more details.

Following the theory presented in Condon *et al.* (2009*c*), we proceed to expand formally the different terms in the ODE and identify those multiplying equal powers of *ω*. For the time being we disregard convergence issues, since we intend to verify the validity of the ansatz. We first observe that term-by-term differentiation in equation (2.4) gives
2.5and
2.6

Since we have assumed that *R*(*y*) and *S*(*y*) are analytic, we can expand in Taylor series about the function *p*_{0,0}(*t*)
2.7
and analogously for *S*(*y*), where
2.8
and
with the standard notation for multi-indices |** l**|=

*l*

_{1}+

*l*

_{2}+⋯+

*l*

_{n}. A similar formula applies to

*S*(

*y*).

For simplicity of notation, in the sequel we suppress the dependence on *t* of the different terms in the expansion.

### (a) Separation of orders of magnitude

We now attempt to separate the 𝒪(*ω*^{−r}) term for *r*≥0 in the differential equation. The corresponding term from *y*′′(*t*) follows from equation (2.6), while the term *S*(*y*) yields a contribution given by equation (2.7). Finally, the 𝒪(*ω*^{−r}) term from the product −*R*(*y*)*y*′ can be extracted from equations (2.5) and (2.7), and after some manipulations the outcome is as follows:
2.9

Note that the last term is nil for *r*=2. Putting the three contributions together, we obtain the term.

### (b) Separation of frequencies

Next we separate the different frequencies within each term. First, we observe that where

Note that, unlike the multi-indices in , in this case the components can also be non-positive integers. Likewise, and

Combining everything, we obtain where

This can be somewhat simplified, because
the last line being justified by the fact that *l*_{n+1}*p*_{1,ln+1}=0. Moreover, since for *j*∈{0,1} we have *p*_{j,ln+1}′=0 (unless *l*_{n+1}=0), we obtain
Therefore, separating frequencies, we have for every
2.10

Further simplification is possible, using the following results.

### Proposition 2.1.

*For every r*≥1 *and* *, it is true that*
2.11

### Proposition 2.2.

*For every r*≥1 *and**, it is true that*
2.12

The proofs of both propositions are relegated to appendix A.

If we substitute equations (2.11) and (2.12) into equation (2.10), the outcome is 2.13

An important observation is that in each step we obtain from equation (2.13) both an ODE for *p*_{r,0}(*t*), which is non-oscillatory since there is no dependence on *ω*, and a recursion for *p*_{r+2,m}(*t*), *m*≠0. More precisely, since *p*_{r,0}(*t*) terms on the right feature only for *n*=1, we have
2.14
where
2.15
is the linearization of the original ODE *y*′′−*R*(*y*)*y*′+*S*(*y*)=0 about *y*=*p*_{0,0}.

Likewise, for *m*≠0, we have a recursion for the *p*_{r,m}(*t*) terms
2.16

The initial conditions for the ODE (2.14) are determined by imposing
2.17
consequently the rest of the *p*_{r,m}(*t*) coefficients, together with their derivatives, should be equal to 0 when *t*=0, that is,
2.18
and
2.19

In the next section, we present the first few terms of the expansion, computed using the differential equation and the recursion presented above.

## 3. Construction of the asymptotic expansion

### (a) The zeroth term

When *r*=0, we have, directly from the differential equation,
together with the initial conditions (2.17). It is also possible to show that
3.1
directly in terms of the modulated Fourier coefficients of the forcing term. We thus deduce that
3.2
in accordance with equation (2.18). These initial conditions will be used to solve the ODE for *p*_{1,0}(*t*), which is given by the analysis of the terms.

### (b) The first term

We now look at the 𝒪(*ω*^{−1}) terms and separate scales. First we have a differential equation for *p*_{1,0} from equation (2.14), namely , where is given in equation (2.15), and a recursion for the next level
3.3
Moreover, it follows from equation (2.19) that
3.4

### (c) The second term

The differential equation for *p*_{2,0}(*t*) reads

However, note that, because *p*_{1,m}(*t*)≡0 when *m*≠0, we have

The recursion for *p*_{4,m}(*t*), *m*≠0, after some similar simplification, reads
3.5

### (d) The third term

There is an important reason to consider the case *r*=3 in detail, and it is related to the blossoming phenomenon that we consider in the next section. With the same simplifications that we applied before, the ODE for *p*_{3,0}(*t*) reads

The recursion for *p*_{5,m}(*t*), *m*≠0, can be deduced from equation (2.16), setting *r*=3. Again using the fact that *p*_{1,m}(*t*)≡0 except when *m*=0, the expressions can be considerably simplified to yield

It is clear that the process can be continued, at the price of increasingly more complicated algebra. This, however, can be overcome with the use of a symbolic algebra package, and it is also worth noticing that in most important examples the functions *R*(*y*) and *S*(*y*) are quite simple (e.g. polynomials in *y*), and hence some terms in the previous expressions are identically 0.

## 4. Band-limited input and blossoming

The stage *r*=3 in the previous computations has a special significance. We note that the forcing terms in equation (1.5) share a common feature, namely that they are clearly *band limited*, since the number of frequencies is finite. In the case of the AM signal, this follows from elementary trigonometric identities, and the spectrum of the forcing term contains the carrier frequency *ω*_{1} (unless we implement a suppressed carrier modulation) and the two sidebands *ω*_{1}±*ω*_{2}, together with the negative frequencies −*ω*_{1} and −*ω*_{1}±*ω*_{2}.

If the original oscillator is *band limited*, i.e. there exists such that *α*_{m}=0 for |*m*|≥*ϱ*+1, then it is clear from our narrative that *p*_{2,m},*p*_{3,m},*p*_{4,m}=0 for |*m*|≥*ϱ*+1. In other words, the relevant MFEs stay band limited with the same original bandwidth *ϱ*. However, things change with regard to *p*_{5,m}. The fact that *p*_{2,m} is band limited implies that
This leads to the non-empty range of summation when
hence |*m*|≤2*ϱ*. In other words, *p*_{5,m} is band limited of bandwidth 2*ϱ*: we call this phenomenon blossoming and note that it has obvious implications in the programming of the method.

The bandwidth of linear systems is the same as that of the highly oscillatory input (i.e. there is no blossoming). It is known, however, that nonlinearity might interfere with bandwidth, and our analysis quantifies this phenomenon for equations of type (1.1). Of course, the bandwidth is likely to blossom further as *r* increases, indicating that resonance shifts energy between frequencies. What is remarkable is that all this occurs only once we hit *p*_{5,m}. Consequently, as long as we do not go beyond *r*=4, disregarding error of 𝒪(*ω*^{−5}), we retain the original bandwidth of the forcing term. If *ω* is large enough, this is likely to be sufficient for virtually all cases of interest.

Let us denote by *θ*_{r} the bandwidth of the 𝒪(*ω*^{−r}) term; therefore,

Before we state a general theorem, let us acquire basic intuition in manipulating the relevant expressions. In order to compute *θ*_{r+1}, we consider terms of the form
for *n*∈{1,…,*r*−1} and *n*∈{1,…,*r*}, respectively. Clearly, it is enough to consider only the terms in , since we wish to maximize the bandwidth.

Consider *r*=5 as a first example, so *n*∈{1,…,5}. Because of symmetry, we might assume without loss of generality that the entries of are ordered monotonically.

*n*=1. The only possible term is*p*_{5,m}; hence, in this case*θ*_{5}=2*ϱ*.*n*=2. We have two monotone choices:=(1,4), which results in*k**θ*_{1}+*θ*_{4}=*ϱ*, and=(2,3), which gives*k**θ*_{2}+*θ*_{3}=2*ϱ*.*n*=3. Now there are two monotone possibilities, (1,1,3) and (1,2,2), with maximal bandwidths of 2*θ*_{1}+*θ*_{3}=*ϱ*and*θ*_{1}+2*θ*_{2}=2*ϱ*, respectively.*n*=4. Just a single monotone choice, (1,1,1,2), leading to 3*θ*_{1}+*θ*_{2}=*ϱ*.*n*=5. Only one 5-tuple, (1,1,1,1,1), and the maximal bandwidth is 5*θ*_{1}=0.

Therefore, *θ*_{6}=2*ϱ*, the maximum over all possible choices.

A similar analysis of monotone sequences in the case *r*=6 yields *θ*_{7}=3*ϱ*. This bandwidth is obtained from *n*=3 and ** k**=(2,2,2). This observation is crucial in the proof of the following general theorem.

## Theorem 4.1.

*It is true that*
4.1

## Proof.

The theorem is certainly true for *r*≤5. We continue by induction on *r*. Thus, we assume that it is true up to *r*≥2 and wish to prove it for *r*+1.

Given *r*≥2, the formula for *p*_{r+1,m} is a linear combination of terms of the form
and their derivatives—as the derivative does not change the bandwidth, we can disregard differentiation in this context. The bandwidth is provided by the largest such that
where , |** k**|=

*r*, and |

**|=**

*l**m*. Each

*k*

_{j}can contribute at most the bandwidth

*θ*

_{j}; hence, altogether the bandwidth of

*ϖ*

_{k,l}is at most the sum of all

*θ*

_{kj}, for 1≤

*j*≤

*n*. As a consequence, 4.2

First we observe that *θ*_{r+1}≥*θ*_{r} when *r*≥2. This can be deduced by considering *n*=1, whence *ϖ*_{r,m}=*p*_{r,m}. Next, we prove that

If , then let us choose and . Since *θ*_{2}=*ϱ*, we deduce that the bandwidth of *ϖ*_{k,l} is maximized by ; therefore, . Likewise, if , then we choose and ** k**=(1,2,2,…,2). Since

*θ*

_{1}=0 and

*θ*

_{2}=

*ϱ*, it follows again that .

In order to prove the reverse inequality
let , *n*∈{1,2,…,*r*}. We may assume without loss of generality that none of the *k*_{j}s equals 1. The reason is as follows. Suppose that *s* of the *k*_{j}s are 1 and denote by the vector that we obtain after excising these terms. Since *θ*_{1}=0, the maximal bandwidths of *ϖ*_{k,l} are the same as the maximal bandwidth of for some . But, since , the term has already featured while forming *p*_{r−s+1,m} for some . Now, we already know that *θ*_{r+1}≥*θ*_{r}≥⋯≥*θ*_{r−s}; hence, this term can be disregarded and we can assume without loss of generality that *k*_{j}≠1, *j*=1,2,…,*n*.

We write
where *β*_{1}=⋯=*β*_{n1}=2, *γ*_{1},…,*γ*_{n2} are even, and *δ*_{1},…,*δ*_{n3} are odd, , with , . We thus have *n*=*n*_{1}+*n*_{2}+*n*_{3} and
(Thus, *n*_{3} is necessarily of the same parity as *r*.) Moreover, by induction,
Now,
and we deduce that

Now we need to maximize the quantity on the right-hand side for because of equation (4.2). Recalling that *r* and *n*_{3} must have the same parity, we distinguish between even and odd *r*. If , then ; therefore,
Likewise, for , we have and again

Taking all this together completes the proof of the theorem. ▪

The implications of this result when programming the method are clear, since it is possible to determine in advance how many terms *p*_{r,m}(*t*) we need for any given *r*≥2.

## 5. Examples

We consider first the Van der Pol oscillator. In this case, we have *R*(*y*)=*μ*(1−*y*^{2}) and *S*(*y*)=*y*, and the base equation is

Using the same notation as before, see equation (2.15), we have, for *r*≥1,
and

We take the initial values *y*(0)=0 and *y*′(0)=1, and the forcing term . Thus, we have , and *α*_{m}≡0 for *m*≠±1. The initial values for the system of ODEs follow from equations (3.2), (3.4) and (2.19),

Moreover, from equation (3.1), we have consequently,

Also, from equation (3.3), so

We will use all this information to assemble the numerical solver up to order 3. In table 1, we display the errors in the approximation of the solution *y*(*t*) and its derivative *y*′(*t*) for different values of *ω*, setting *μ*=1 and using the first few terms in the asymptotic expansion. The non-oscillatory systems of ODEs for the coefficients *p*_{1,0}, *p*_{2,0} and *p*_{3,0} are solved numerically using the ode45 routine in Matlab, and we compare the results with the solution of the original ODE, using relative and absolute tolerance equal to 10^{−12}. The notation used for the errors is

The next example is an application to the Duffing equation with damping

We take *k*=1/2, *a*=1, *b*=−1/3 and a forcing term that is an AM signal
with *c*_{1}=40, *c*_{2}=20 and frequencies *ω*_{1}=1000 and *ω*_{2}=100. In order to construct the asymptotic expansion in a modulated Fourier series, we define *ω*:=(*ω*_{1},*ω*_{2}), which is the greatest common divisor of the two frequencies. Additionally, let *m*_{1}=*ω*_{1}/*ω* and *m*_{2}=*ω*_{2}/*ω*.

In this case, the base equation is
and, for *r*≥1,

Moreover,

We can easily work out the initial values
and
where *m*_{+}=*m*_{1}+*m*_{2} and *m*_{−}=*m*_{1}−*m*_{2}. The other non-zero terms are
and
so

Similarly, and therefore,

Figures 2 and 3 display the errors in the approximation of the solution *y*(*t*) and its derivative *y*′(*t*), using a different number of terms in the asymptotic expansion.

## 6. Conclusions and further research

We have presented a combined asymptotic-numerical method to solve efficiently second-order differential equations with highly oscillatory forcing terms. The approach is based on using asymptotic expansions in inverse powers of the oscillatory parameter *ω* together with MFEs. With the aid of a computer algebra package such as Maple, it is possible to compute all the terms in this type of expansions to high accuracy.

A key feature of this approach is that, unlike classical numerical algorithms for ODEs, the performance of this method improves in the presence of high oscillation, that is, when *ω* is large. This is a consequence of the asymptotic methodology that we have used, instead of applying the Taylor expansion of the solution.

We have presented numerical examples based on two equations that are very relevant in applications, the forced Van der Pol and Duffing oscillators. This is just one possible application of this type of asymptotic-numerical solvers. Other scenarios that are currently being analysed are ODEs where the coefficients depend on *ω* (a situation that includes important examples such as the inverted pendulum) and differential-algebraic equations, which are relevant in the modelling of electronic circuits.

## Acknowledgements

A.D. acknowledges financial support from the Spanish Ministry of Education under the programme of postdoctoral grants (Programa de becas postdoctorales) and project MTM2006-09050. The material is based upon works supported by Science Foundation Ireland under Principal Investigator grant no. 05/IN.1/I18.

## Appendix A. Two propositions in §2.2

In this appendix, we present the proofs of the two propositions that we used before.

## Proposition A.1.

*For every r*≥1 and*, it is true that*

## Proof.

By direct differentiation, the right-hand side is equal to
However, because of symmetry, for any 1≤*q*≤*n*,
therefore
In the last summation, we let *s*=*k*_{1}+⋯+*k*_{n}. Since *s*+*k*_{n+1}=*r* and *k*_{j}≥1, we deduce that *s*∈{*n*,*n*+1,…,*r*−1}. Moreover, *k*_{n+1}=*r*−*s* and
The last step follows because, letting *n*∈{1,2,…,*r*} and *s*=*r* and noting that *p*_{0,ln+1}≠0 only for *l*_{n+1}=0, we recover the first sum. The proposition follows. ▪

## Proposition A.2.

*For every r*≥1 and*, it is true that*

## Proof.

Similar to the proof of the previous proposition. We let and |** k**|=

*s*; hence,

*k*

_{n+1}=

*r*−

*s*+1, while

*s*∈{

*n*,

*n*+1,…,

*r*}. In other words, Therefore, shifting the index

*n*, Finally, for

*n*=1, we have , ; therefore,

Using the underlying symmetry, it is true that The lemma follows by straightforward substitution. ▪

## Footnotes

- Received September 14, 2009.
- Accepted December 4, 2009.

- © 2010 The Royal Society