Word series high-order averaging of highly oscillatory differential equations with delay

We show that, when the delay is an integer multiple of the forcing period, it is possible to obtain easily high-order averaged systems of periodically forced systems of delay differential equations with constant delay.

The well-known theory of averaging [14, 15, 16, 17, 18, 19, 20] studies the reduction, by means of time-dependent changes of variables, of (nonautonomous) forced systems of differential equations to autonomous time-independent systems (averaged systems). Such a reduction is useful because autonomous systems are easier to analyze than their nonautonomous counterparts. In addition, the numerical integration of periodically or quasi-periodically forced systems may be expensive because typically integrators have to operate with step-sizes that are small with respect to the shortest period present in the forcing; in those cases, integrating an averaged version of the given system may be very advantageous.

In many occasions averaging is only carried out to first order. When the forcing is periodic this simply means replacing the right-hand side of the given oscillatory differential system by its average over one period of the forcing. In other cases, first order averaging is not enough and it is desirable to obtain averaged systems that provide more accurate approximations. For ordinary differential equations without delay the obtention of high-order averaging may be a daunting task, even with the help of a symbolic manipulator. The series of papers [6, 7, 8, 9, 10, 15, 18, 19] has developed a technique that, in the absence of delays, makes it possible to compute high-order averaged systems through simple recursions. The technique is based on the theory of word series [16], formal series [21] that have numerous applications in the fields of deterministic or stochastic dynamical systems [17] and numerical integration [1, 2].

The difficulties in obtaining high-order averaged systems are compounded if the system to be averaged has delays. In this paper we show that, for periodically forced differential systems with constant delay, it is possible to obtain high-order averaged systems by an application of the word series results in [6, 7, 8, 9, 10, 15, 18, 19]. The simple treatment presented here is only possible when the forcing period is a submultiple of the delay, a hypothesis whose scope is discussed in Section 3.

Section 2 recalls the word series averaging of systems of differential equations without delay. Section 3 contains the main idea and Section 4 describes an application that shows how third-order averaging succeeds where second-order averaging fails. The final Section 5 concludes.

In this section we summarize the word series averaging of differential equations without delay [18]. This summary will provide the foundation to address later the case with delay.

Let us consider highly oscillatory initial value problems of the form

ddtξ<!-- ? -->=g(ξ<!-- ? -->,Ω<!-- O -->t),t≥<!-- = -->0,ξ<!-- ? -->(0)=ξ<!-- ? -->0∈<!-- ? -->CD. $$\begin{array}{} \displaystyle \frac{d}{dt}\xi = g(\xi, \Omega t),\quad t\geq0,\qquad \xi(0) = \xi_0\in\mathbb{C}^D. \end{array}$$(1)

Here the smooth function g(ξ, θ) is 2π-periodic in its second argument θ ∈ ℝ, with Fourier expansion

g(ξ<!-- ? -->,θ<!-- ? -->)=∑<!-- ? -->k∈<!-- ? -->Zexp⁡<!-- ? -->(ikθ<!-- ? -->)gk(ξ<!-- ? -->), $$\begin{array}{} \displaystyle g(\xi,\theta) = \sum_{k\in\mathbb{Z}} \exp(ik\theta) g_k(\xi), \end{array}$$(2)

and the parameter Ω ≫ 1 is the angular frequency of the fast periodic forcing. The corresponding period is T = 2π/Ω.

In the theory of word series the set of indices ℤ in (2) is seen as an (infinite) alphabet and for n = 1, 2, …, the elements of ℤⁿ are regarded as words consisting of n letters of this alphabet. There is also an empty word ∅ with n = 0 letters. Given (1), to each word w one associates the corresponding word basis function g_w, a map ℂ^D → ℂ^D. By definition, for words with one letter k ∈ ℤ, the word basis function g_w(ξ) is the Fourier coefficient g_k(ξ) in (2). For words with n > 1 letters, we define recursively

gk1⋯<!-- ? -->kn(ξ<!-- ? -->)=gk2⋯<!-- ? -->kn′<!-- ' -->(ξ<!-- ? -->)gk1(ξ<!-- ? -->), $$\begin{array}{} \displaystyle g_{k_1\cdots k_n}(\xi) = g^\prime_{k_2\cdots k_n}(\xi)g_{k_1}(\xi), \end{array}$$(3)

where gk2⋯<!-- ? -->kn′<!-- ' -->(ξ<!-- ? -->) $\begin{array}{} \displaystyle g^\prime_{k_2\cdots k_n}(\xi) \end{array}$ is the Jacobian matrix of g_k2⋯kn(ξ). For the empty word ∅, the basis function is the identity map ξ ↦ ξ.

The set of all words (including the empty word) is denoted by 𝒲 and ℂ^𝒲 refers to the set of all mappings δ : 𝒲 → ℂ. For δ ∈ ℂ^𝒲 and w ∈ 𝒲, δ_w denotes the value of δ at w. To each δ ∈ ℂ^𝒲 we associate the corresponding word series (relative to the mappings g_k in (2)); this is the formal series

Wδ<!-- d -->(ξ<!-- ? -->)=∑<!-- ? -->w∈<!-- ? -->Wδ<!-- d -->wgw(ξ<!-- ? -->), $$\begin{array}{} \displaystyle W_\delta(\xi) = \sum_{w\in\mathscr{W}}\delta_w g_w(\xi), \end{array}$$

whose terms are maps ℂ^D → ℂ^D. The complex numbers δ_w are the coefficients of the series. As an example, consider the case where δ_∅ = 1 and δ_w = 0 for each nonempty word w; in this case W_δ(ξ) = ξ.

As proved in [7, 18] (see [19] for an alternative technique), there exist β̄ ∈ ℂ^𝒲 and a 2π-periodic map θ ∈ ℝ → κ(θ) ∈ ℂ^𝒲 such that the solution of (1) may be written as

ξ<!-- ? -->(t)=Wκ<!-- ? -->(Ω<!-- O -->t)(Ξ<!-- ? -->(t)), $$\begin{array}{} \displaystyle \xi(t) = W_{\kappa(\Omega t)}(\Xi(t)), \end{array}$$(4)

where Ξ(t) is the solution of

ddtΞ<!-- ? -->=Wβ<!-- ? -->¯<!-- ? -->(Ξ<!-- ? -->),Ξ<!-- ? -->(0)=ξ<!-- ? -->0∈<!-- ? -->RD. $$\begin{array}{} \displaystyle \frac{d}{dt}\Xi = W_{\bar \beta}(\Xi),\qquad \Xi(0) = \xi_0\in\mathbb{R}^D. \end{array}$$(5)

In this way, the time-dependent change of variables ξ ↦ Ξ in (4) transforms the given highly oscillatory initial value problem into the initial value problem (5) where there is no periodic forcing. Therefore (5) provides an averaged version of (1). In addition, the coefficients κ_w(θ) are such that, at θ = 2kπ, k = 0, 1, 2, …, we have κ_∅(θ) = 1 and κ_w(θ) = 0 for each nonempty word w. This implies that, for those values of θ, the transformation ξ ↦ Ξ is the identity map and it then follows that, in (4), ξ(t) = Ξ(t) at the stroboscopic times, t = 0, T, 2T, … We then say that the change of variables is stroboscopic and that (5) is obtained from (1) through stroboscopic averaging. There are of course alternative forms of averaging; for instance one may impose the condition that the periodic change of variables ξ ↦ Ξ has 0 average over one period, rather than being the identity at θ = 0.

It is in order to point out that β̄ and κ(θ) depend on Ω, but are otherwise universal in the sense that they do not change with the dimension D or with the choice of g(ξ, θ). Such a universality implies that they may be computed once and for all; averaging a new differential system requires the computation of new basis functions but not of new coefficients. The values β̄_w and κ_w(θ) for w ∈ 𝒲 may be computed easily by recursion with respect to the number of letters in the word w; the reader is referred to [7, 10, 18] for details.

In general, the formal series W_β̄(Ξ) in the averaged system (5) and the formal series in the change of variables (4) do not converge and have to be truncated and, of course, the truncation introduces an error, see [8, 9]. For our purposes here we just mention that, if W_β̄(Ξ) is truncated by eliminating all the terms in the series that correspond to words with more than n letters, one obtains an initial value problem whose solution Ξ(t) coincides at stroboscopic times with the solution ξ(t) except for an error of size 𝓞(1/Ωⁿ) as Ω → ∞.

To approximate ξ at times that are not stroboscopic one needs to apply to the solution of the truncated averaged system a change of variables obtained by truncating the series in (4).

ddtΞ<!-- ? -->=g0(Ξ<!-- ? -->), $$\begin{array}{} \displaystyle \frac{d}{dt} \Xi = g_0(\Xi), \end{array}$$(6)

as it may have been expected.

The truncated averaged system with second-order errors 𝓞(1/Ω²) is given by [7, 10, 18]

ddtΞ<!-- ? -->=g0(Ξ<!-- ? -->)+∑<!-- ? -->k≠<!-- ? -->0ikΩ<!-- O -->(g0′<!-- ' -->(Ξ<!-- ? -->)gk(Ξ<!-- ? -->)−<!-- - -->gk′<!-- ' -->(Ξ<!-- ? -->)g0(Ξ<!-- ? -->)+gk′<!-- ' -->(Ξ<!-- ? -->)g−<!-- - -->k(Ξ<!-- ? -->)). $$\begin{array}{} \displaystyle \frac{d}{dt} \Xi = g_0(\Xi)+\sum_{k\neq 0} \frac{i}{k\Omega} \Big( g_0^\prime(\Xi)g_k(\Xi)- g_k^\prime(\Xi)g_0(\Xi) +g_k^\prime(\Xi)g_{-k}(\Xi) \Big). \end{array}$$(7)

Note that, g₀ is by definition the word basis function associated with the one-letter word 0. Furthermore, according to (3), the function g0′<!-- ' -->gk $\begin{array}{} \displaystyle g_0^\prime g_k \end{array}$ is the word basis function associated to the word k0 and similarly gk′<!-- ' -->g0,gk′<!-- ' -->g−<!-- - -->k $\begin{array}{} \displaystyle g_k^\prime g_0, g_k^\prime g_{-k} \end{array}$ are associated to 0k and –kk. Thus the right-hand side of (7) is indeed a word series (whose coefficients vanish for all words with 3 or more letters).

The closed form expression for the third-order averaged system for the problem (1) with general g may be seen in [7, 10, 18]. For order ≥ 4 it is not practical to find the general expression of the averaged system and then to apply it to the specific instance of (1) at hand. One should rather find recursively the word basis functions corresponding to the g of interest, perhaps with the help of a symbolic manipulator, see [15, 19] for additional details.

We now consider the D-dimensional system with constant delay τ > 0

ddtx(t)=f(x(t),x(t−<!-- - -->τ<!-- t -->),Ω<!-- O -->t),t≥<!-- = -->0 $$\begin{array}{} \displaystyle \frac{d}{dt}x(t)=f(x(t),x(t-\tau),\Omega t),\qquad t\ge 0 \end{array}$$(8)

x(t)=φ<!-- f -->(t),−<!-- - -->τ<!-- t -->≤<!-- = -->t≤<!-- = -->0, $$\begin{array}{} \displaystyle x(t)=\varphi(t),\qquad -\tau\le t\le 0, \end{array}$$(9)

where f(x, y, θ) is 2π-periodic in its third argument, with Fourier expansion

f(x,y,θ<!-- ? -->)=∑<!-- ? -->k∈<!-- ? -->Zexp⁡<!-- ? -->(ikθ<!-- ? -->)fk(x,y), $$\begin{array}{} \displaystyle f(x,y,\theta) = \sum_{k\in\mathbb{Z}} \exp(ik\theta)f_k(x,y), \end{array}$$(10)

and Ω is the angular frequency as above. {It is assumed that the known function φ that specifies the initial history is independent of Ω; this involves no loss of generality because the general case may be reduced to the Ω-independent case by using the new dependent variable x(t)-Φ(t), where Φ(t) coincides with ϕ(t) for –τ ≤ t ≤ 0. Reference [22] contains a list of references to problems of this form that appear in different applications.

The developments that follow are based on the standing hypothesis that the delay τ is an integer multiple of the period T, or in other words that t = τ is a stroboscopic time. Let us discuss this hypothesis. In some applications, the value of τ is given and the interest is in studying the behaviour of the solution of (8)–(9) for large Ω, but there is some freedom as to the choice of the specific value of Ω. In those cases, the theory below applies by choosing Ω in such a way that τΩ/2π is an integer. Similarly, our theory applies to situations where the interest is in a given, large value of Ω and there is some freedom in the choice of τ to be used in the model. On the other hand, if the values of τ and Ω are fixed, the material in this paper does not apply. In those cases finding high-order averaged systems may be extremely complicated (see e.g. the computations in [22]).

3.1

The averaging procedure

Assume that (8)–(9) is to be studied in a bounded interval of the form 0 ≤ t ≤ Lτ for a suitable integer L > 0. (This hypothesis is made at this stage for mathematical convenience to avoid systems of infinitely many differential equations; the averaged systems that we will find will be valid for all t ≥ 0.) We introduce the functions (see Fig. 1).

Fig. 1

x(0)(t)=φ<!-- f -->(t−<!-- - -->τ<!-- t -->),0≤<!-- = -->t≤<!-- = -->τ<!-- t -->, $$\begin{array}{} \displaystyle x^{(0)}(t)=\varphi(t-\tau), \qquad 0\leq t \leq \tau, \end{array}$$(11)

x(ℓ<!-- l -->)(t)=x(t+(ℓ<!-- l -->−<!-- - -->1)τ<!-- t -->),0≤<!-- = -->t≤<!-- = -->τ<!-- t -->,ℓ<!-- l -->=1,…<!-- ? -->,L, $$\begin{array}{} \displaystyle x^{(\ell)}(t)=x(t+(\ell-1)\tau), \qquad 0\leq t \leq \tau,\qquad \ell=1, \ldots, L, \end{array}$$(12)

and note that determining these functions is clearly equivalent to determining the solution of (8)–(9) in the interval 0 ≤ t ≤ Lτ. The x^(ℓ)(t) with ℓ > 0 satisfy the conditions

x(ℓ<!-- l -->)(0)=x(ℓ<!-- l -->−<!-- - -->1)(τ<!-- t -->),1≤<!-- = -->ℓ<!-- l -->≤<!-- = -->L. $$\begin{array}{} \displaystyle x^{(\ell)}(0)=x^{(\ell-1)}(\tau),\qquad 1 \leq \ell \leq L. \end{array}$$(13)

and the differential equations

ddtx(ℓ<!-- l -->)(t)=f(x(ℓ<!-- l -->)(t),x(ℓ<!-- l -->−<!-- - -->1)(t),Ω<!-- O -->(t+(ℓ<!-- l -->−<!-- - -->1)τ<!-- t -->)),0≤<!-- = -->t≤<!-- = -->τ<!-- t -->,1≤<!-- = -->ℓ<!-- l -->≤<!-- = -->L. $$\begin{array}{} \displaystyle \frac{d}{dt}x^{(\ell)}(t)=f(x^{(\ell)}(t),x^{(\ell-1)}(t),\Omega (t+(\ell-1)\tau)),\quad 0\leq t \leq \tau, \: 1 \leq \ell \leq L. \end{array}$$

Even though the relations (13) may be reminiscent of a two-point boundary value problem, we are facing here an initial value problem: x⁽¹⁾(t) is determined by solving the differential equation with ℓ = 1 with initial condition x⁽¹⁾(0) = φ(0), once x⁽¹⁾(t) is known, x⁽²⁾(t) is determined by solving the differential equation with ℓ = 2 with initial condition x⁽²⁾(0) = x⁽¹⁾(τ), etc.

By taking into account the periodicity of f(x, y, θ) as a function of θ and our standing hypothesis that τ is a multiple of the period, the last display may be written as

ddtx(ℓ<!-- l -->)(t)=f(x(ℓ<!-- l -->)(t),x(ℓ<!-- l -->−<!-- - -->1)(t),Ω<!-- O -->t),0≤<!-- = -->t≤<!-- = -->τ<!-- t -->,1≤<!-- = -->ℓ<!-- l -->≤<!-- = -->L. $$\begin{array}{} \displaystyle \frac{d}{dt}x^{(\ell)}(t)= f(x^{(\ell)}(t),x^{(\ell-1)}(t),\Omega t),\quad 0\leq t \leq \tau, \: 1 \leq \ell \leq L. \end{array}$$(14)

Note that this would be a system of the form (1) for the D × L-dimensional vector (x⁽¹⁾, …, x^(L)) if it were not for the fact that the right-hand side of the equation corresponding to ℓ = 1 contains x⁽⁰⁾, a known function of t given in (11). In order to have a system of the form (1), where the right-hand side only depends on t through the combination θ = Ωt, we proceed as follows. We introduce the 1 + D(L + 1)-dimensional vector of unknown functions of the variable t

ξ<!-- ? -->=(t^<!-- ^ -->,x(0),x(1),…<!-- ? -->,x(L)), $$\begin{array}{} \displaystyle \xi=(\widehat{t}, x^{(0)},x^{(1)}, \dots, x^{(L)}), \end{array}$$

and add to (14) the differential equations

ddtt^<!-- ^ -->=1,ddtx(0)=φ<!-- f -->˙<!-- ? -->(t^<!-- ^ -->−<!-- - -->τ<!-- t -->),0≤<!-- = -->t≤<!-- = -->τ<!-- t -->, $$\begin{array}{} \displaystyle \frac{d}{dt} \widehat{t} = 1,\qquad \frac{d}{dt} x^{(0)} = \dot \varphi(\widehat{t}-\tau),\quad 0\leq t\leq \tau, \end{array}$$(15)

(the dot represents differentiation) subject to the initial conditions

t^<!-- ^ -->(0)=0,x(0)(0)=φ<!-- f -->(−<!-- - -->τ<!-- t -->). $$\begin{array}{} \displaystyle \widehat{t}(0) = 0,\qquad x^{(0)}(0) = \varphi(-\tau). \end{array}$$(16)

The solution of the initial value problem (15)–(16) is obviously t̂ = t and x⁽⁰⁾(t) = φ(t – τ). Now (15) in tandem with (14) is a 1 + D(L + 1)-dimensional system of the form (1) that may be stroboscopically averaged by following the procedure described in the previous section. If X^(ℓ) is the averaged counterpart of x^(ℓ), we are interested in the solutions of the averaged system that satisfy (cf. (13))

X(ℓ<!-- l -->)(0)=X(ℓ<!-- l -->−<!-- - -->1)(τ<!-- t -->),1≤<!-- = -->ℓ<!-- l -->≤<!-- = -->L, $$\begin{array}{} \displaystyle X^{(\ell)}(0)=X^{(\ell-1)}(\tau),\qquad 1 \leq \ell \leq L, \end{array}$$(17)

so that it is possible to define a continuous function X in the interval [–τ, Lτ] by patching together the different X^(ℓ),

X(t)=X(ℓ<!-- l -->)(t−<!-- - -->(ℓ<!-- l -->−<!-- - -->1)τ<!-- t -->),−<!-- - -->τ<!-- t -->≤<!-- = -->t≤<!-- = -->Lτ<!-- t -->,ℓ<!-- l -->=0,…<!-- ? -->,L, $$\begin{array}{} \displaystyle X(t) = X^{(\ell)}(t-(\ell-1)\tau),\quad -\tau\leq t \leq L\tau, \quad \ell = 0, \dots, L, \end{array}$$(18)

thus undoing the process that we used to move from the solution x of (8)–(9) to the functions x^(ℓ) in (11)–(12), see Fig. 1. In this way the averaged delay equation is obtained by patching the equations for the X^(ℓ). The whole procedure will be clear after we present the simple case of the first-order averaged system.

As an illustration, we consider the system

dudt=α<!-- a -->1+vβ<!-- ? -->−<!-- - -->u(t−<!-- - -->τ<!-- t -->)+Asin⁡<!-- ? -->(ω<!-- ? -->t)+Bsin⁡<!-- ? -->(Ω<!-- O -->t),dvdt=α<!-- a -->1+uβ<!-- ? -->−<!-- - -->v(t−<!-- - -->τ<!-- t -->), $$\begin{array}{} \begin{split} \displaystyle \frac{du}{dt}&=&\frac{\alpha}{1+v^{\beta}}-u(t-\tau)+A\sin(\omega t)+B\sin(\Omega t),\\ \frac{dv}{dt}&=&\frac{\alpha}{1+u^{\beta}}-v(t-\tau), \end{split} \end{array}$$

where α and β are parameters. There are two periodic forcing terms; the forcing A sin(ωt) is slow (i.e. the frequency ω is of moderate size) and B sin(Ωt), Ω ≫ 1, is fast. When there is no forcing (A = B = 0), the system represents a delayed genetic toggle switch, a synthetic gene regulatory network [12]. The forced system was studied in [11] as an instance of the emergence of vibrational resonance [13, 15], i.e. the enhancement, due to the presence of the fast forcing, of the response of the system to the slow forcing.

In order to have a system of the form (8) it is necessary to rewrite A sin(ωt) as A sin(ωt̂), where t̂ is a new dependent variable defined by the differential equation (d/dt) t̂ = 1 and the initial condition t̂ = 0 at t = 0 (this is of course the technique used above in (15)–(16) to deal with the slow dependence on t of the right-hand side of (14)). Thus only the fast forcing is averaged. After computing recursively the required word basis functions and coefficients, the third-order averaged system is found to be given by

dUdt=α<!-- a -->1+Vβ<!-- ? -->−<!-- - -->U(t−<!-- - -->τ<!-- t -->)+Asin⁡<!-- ? -->(ω<!-- ? -->t),dVdt=α<!-- a -->1+Uβ<!-- ? -->−<!-- - -->V(t−<!-- - -->τ<!-- t -->)−<!-- - -->BΩ<!-- O -->α<!-- a -->β<!-- ? -->Uβ<!-- ? -->−<!-- - -->1(1+Uβ<!-- ? -->)2+B2Ω<!-- O -->23α<!-- a -->β<!-- ? -->Uβ<!-- ? -->−<!-- - -->2(Uβ<!-- ? -->−<!-- - -->β<!-- ? -->+β<!-- ? -->Uβ<!-- ? -->+1)4(1+Uβ<!-- ? -->)3, $$\begin{array}{} \begin{split} \displaystyle \frac{dU}{dt}&=&\frac{\alpha}{1+V^{\beta}}-U(t-\tau)+A\sin(\omega t),\\ \nonumber \frac{dV}{dt}&=&\frac{\alpha}{1+U^{\beta}}-V(t-\tau)-\frac{B}{\Omega}\frac{\alpha\beta U^{\beta-1}}{(1+U^{\beta})^2} +\frac{B^2}{\Omega^2}\frac{3\alpha\beta U^{\beta-2}(U^{\beta}-\beta+\beta U^{\beta}+1)}{4(1+U^{\beta})^3}, \end{split} \end{array}$$

for 0 ≤ t < τ, and

dUdt=α<!-- a -->1+Vβ<!-- ? -->−<!-- - -->U(t−<!-- - -->τ<!-- t -->)+Asin⁡<!-- ? -->(ω<!-- ? -->t)−<!-- - -->BΩ<!-- O -->,dVdt=α<!-- a -->1+Uβ<!-- ? -->−<!-- - -->V(t−<!-- - -->τ<!-- t -->)−<!-- - -->BΩ<!-- O -->α<!-- a -->β<!-- ? -->Uβ<!-- ? -->−<!-- - -->1(1+Uβ<!-- ? -->)2+B2Ω<!-- O -->23α<!-- a -->β<!-- ? -->Uβ<!-- ? -->−<!-- - -->2(Uβ<!-- ? -->−<!-- - -->β<!-- ? -->+β<!-- ? -->Uβ<!-- ? -->+1)4(1+Uβ<!-- ? -->)3, $$\begin{array}{} \begin{split} \displaystyle \frac{dU}{dt}&=&\frac{\alpha}{1+V^{\beta}}-U(t-\tau)+A\sin(\omega t)-\frac{B}{\Omega},\\ \nonumber \frac{dV}{dt}&=&\frac{\alpha}{1+U^{\beta}}-V(t-\tau)-\frac{B}{\Omega}\frac{\alpha\beta U^{\beta-1}}{(1+U^{\beta})^2} +\frac{B^2}{\Omega^2}\frac{3\alpha\beta U^{\beta-2}(U^{\beta}-\beta+\beta U^{\beta}+1)}{4(1+U^{\beta})^3}, \end{split} \end{array}$$

for t ≥ τ. (As pointed out above, in general, the third-order averaged system has an analytic expression in τ ≤ t < 2τ and a different analytic expression in t ≥ 2τ; for the problem at hand those two expressions happen to coincide. Also note that, for the problem at hand, the past history φ happens not to feature in the averaged system.) The second-order averaged system may be retrieved from the last displays by omitting the terms that have an Ω² factor in the denominator.

We have integrated numerically the true dynamics and the second- and third-order averaged systems. The integrations were carried out with the Matlab function dde23 with relative and absolute tolerances 10^–8 and 10^–10 respectively. For the choice α = 2.5, β = 2, A = 0.1, ω = 0.1, B = 2.0, and τ = 0.5, with constant history u(t) = 0.5, v(t) = 2.0, –τ ≤ t ≤ 0 (an equilibrium of the unforced system), Table 1 presents the maximum error in u at stroboscopic times in the short interval 0 ≤ t ≤ 2; in agreement with the theory, the error behavior is 𝓞(1/Ω²) for second-order averaging and 𝓞(1/Ω³) for third-order averaging.

Table 1

Maximum errors in u, with respect to the true solution, for the second-order and the third-order averaged solutions in the interval 0 ≤ t ≤ 2. The notation 3.31(–4) means 3.31 × 10^–4, etc.

In Fig. 2, the parameters are α = 2.5, β = 2, A = 0.2, ω = 0.2, B = 2.0, Ω = 4π and τ = 0.5, with constant history u(t) = 2.0, v(t) = 0.5, –τ ≤ t ≤ 0 (again an equilibrium of the unforced system). The integration is perfomed in the interval 0 ≤ t ≤ 100. The numerical integration of the oscillatory problem needed 261.0 seconds, a quantity that has to be compared with the 2.1 and 2.6 seconds required to integrate the second- and third-order averaged systems respectively. In studies like those performed in [11], where the oscillatory systems has to be integrated in long time intervals for many different choices of the values of the parameters, the advantage of averaging is then clear.

Fig. 2

We have showed that, when the delay is a multiple of the forcing period, it is possible to extend to periodically forced, constant delay problems, the word series approach to the systematic derivation of high-order averaged systems. We have presented an example where the new technique has been applied to a system of interest in connection with the phenomenon of vibrational resonance.

To conclude let us point out that, for ordinary differential equations, it is possible to compute numerically a stroboscopically averaged solution Ξ without the explicit knowledge of the corresponding averaged system; the information on the averaged system required by the integrator is derived by numerically simulating the oscillatory system (1) [3, 4]. Such techniques have been extended to delay differential equations [5, 22].