A stroboscopic averaging algorithm for highly oscillatory delay problems

TL;DR

This paper introduces a multiscale numerical method called the stroboscopic averaging method (SAM) for efficiently solving highly oscillatory delay differential equations with fast periodic forcing, achieving high accuracy with reduced computational effort.

Contribution

The paper presents a novel heterogenous multiscale method (SAM) that provides accurate approximations for delay differential equations under fast periodic forcing, with error bounds independent of the forcing frequency.

Findings

Achieves (H^2+1/\u03a9^2) error bounds.

Computational effort scales like H^{-1}, independent of forcing frequency.

Uniform accuracy across a range of high-frequency forcings.

Abstract

We propose and analyze a heterogenous multiscale method for the efficient integration of constant-delay differential equations subject to fast periodic forcing. The stroboscopic averaging method (SAM) suggested here may provide approximations with $\(\mathcal{O}(H^2+1/\Omega^2)\)$ errors with a computational effort that grows like $\(H^{-1}\)$ (the inverse of the stepsize), uniformly in the forcing frequency Omega.

Figures

Tables

Table 1. Table 1: SAM Algorithm

Compute $X_1$

Micro-integration

$u_{0,0}=\phi \left(0\right)$ % initial value

For $\nu =0:{\nu }_{max}$, $v_{0,\nu }=\phi \left(-\tau +\nu h\right)$, end % history

For $\nu =0:{\nu }_{max}-1$, $u_{0,\nu +1}=u_{0,\nu }+hf(u_{0,\nu },v_{0,\nu },t_0+\nu h,\Omega \nu h;\Omega )$, end % Euler

$F_0=\left(u_{0,{\nu }_{max}}-u_{0,0}\right)/T$ % slope for macro-step

Macro-integration

$t_0=0$, $X_0=\phi \left(0\right)$, $t_1=t_0+H$, $X_1=X_0+H{F}_0$ % Euler macro at n=0

Compute $X_2,\dots ,X_N$

For $n=1:N-1$

Micro-integration

$u_{n,0}=X_n$ % initial value

For $\nu =-{\nu }_{max}:{\nu }_{max}$, $v_{n,\nu }=\phi \left(-\tau +nH+\nu h\right)$, end % history

For $\nu =0:{\nu }_{max}-1$, $u_{n,\nu +1}=u_{n,\nu }+hf(u_{n,\nu },v_{n,\nu },t_n+\nu h,\Omega \nu h;\Omega )$

$u_{n,-\nu -1}=u_{n,-\nu }-hf(u_{n,-\nu },v_{n,-\nu },t_n-\nu h,-\Omega \nu h;\Omega )$, end % Euler

$F_n=\left(u_{n,{\nu }_{max}}-u_{n,-{\nu }_{max}}\right)/\left(2T\right)$ % slope for macro-step

Macro-integration

$t_{n+1}=t_n+H$, $X_{n+1}=X_n+\left(3H/2\right)F_n-\left(H/2\right)F_{n-1}$ % Adams-Bashforth2 macro

end

Compute $X_{N+1}$

Micro-integration

$u_{N,0}=X_N$ % initial value

For $\nu =1:{\nu }_{max}$, $u_{0,-\nu }=\phi \left(-\nu h\right)$, end % values taken from history

For $\nu =-{\nu }_{max}:{\nu }_{max}$, $v_{N,\nu }=u_{0,\nu }$, end % save for later history

For $\nu =0:{\nu }_{max}-1$, $u_{N,\nu +1}=u_{N,\nu }+hf(u_{N,\nu },v_{N,\nu },t_N+\nu h,\Omega \nu h;\Omega )$

$u_{N,-\nu -1}=u_{N,-\nu }-hf(u_{N,-\nu },v_{N,-\nu },t_N-\nu h,-\Omega \nu h;\Omega )$, end % Euler

$F_N=\left(u_{N,{\nu }_{max}}-u_{N,0}\right)/T$ % slope for macro-step

Macro-integration

$t_{N+1}=t_N+H$, $X_{N+1}=X_N+H{F}_N$ % Euler macro at n=N

Compute $X_{N+2},\dots$

For $n=N+1:\lfloor t_{max}/H\rfloor$

Micro-integration

$u_{n,0}=X_n$ % initial value

For $\nu =-{\nu }_{max}:{\nu }_{max}$, $v_{n,\nu }=u_{n-N,\nu }$, end % save for later history

For $\nu =0:{\nu }_{max}-1$, $u_{n,\nu +1}=u_{n,\nu }+hf(u_{n,\nu },v_{n,\nu },t_n+\nu h,\Omega \nu h;\Omega )$

$u_{n,-\nu -1}=u_{n,-\nu }-hf(u_{n,-\nu },v_{n,-\nu },t_n-\nu h,-\Omega \nu h;\Omega )$, end % Euler

$F_n=\left(u_{n,{\nu }_{max}}-u_{n,-{\nu }_{max}}\right)/\left(2T\right)$ % slope for macro-step

Macro-integration

$t_{n+1}=t_n+H$, $X_{n+1}=X_n+\left(3H/2\right)F_n-\left(H/2\right)F_{n-1}$ % Adams-Bashforth2 macro

end

Table 2. Table 2: Errors in x 1 subscript 𝑥 1 x_{1} for SAM with respect to the averaged solution for problem ( 1 ). Step points are not stroboscopic times

N	$\Omega =25$	$\Omega =50$	$\Omega =100$	$\Omega =200$	$\Omega =400$	$\Omega =800$	$\Omega =1600$	$\Omega =3200$
1	6.28(-2)	3.42(-2)	1.71(-2)	7.87(-3)	3.14(-3)	1.66(-3)	2.04(-3)	2.30(-3)
2	***	7.66(-3)	3.74(-3)	1.66(-3)	8.27(-4)	7.56(-4)	7.20(-4)	7.02(-4)
4	***	***	1.11(-3)	4.45(-4)	2.60(-4)	2.20(-4)	1.99(-4)	1.88(-4)
8	***	***	***	1.80(-4)	6.35(-5)	5.57(-5)	5.06(-5)	4.77(-5)
16	***	***	***	***	3.20(-5)	1.27(-5)	1.22(-5)	1.18(-5)
32	***	***	***	***	***	6.31(-6)	2.81(-6)	2.85(-6)
64	***	***	***	***	***	***	1.36(-6)	6.46(-7)
128	***	***	***	***	***	***	***	3.22(-7)

Table 3. Table 3: Errors in x 1 subscript 𝑥 1 x_{1} for SAM with respect to the averaged solution for problem ( 1 ). Step points are stroboscopic times

N	$\Omega =8\pi$	$\Omega =16\pi$	$\Omega =32\pi$	$\Omega =64\pi$	$\Omega =128\pi$	$\Omega =256\pi$	$\Omega =512\pi$	$\Omega =1024\pi$
1	6.25(-2)	3.40(-2)	1.70(-2)	7.82(-3)	3.11(-3)	1.66(-3)	2.04(-3)	2.30(-3)
2	***	7.62(-3)	3.72(-3)	1.65(-3)	8.26(-4)	7.56(-4)	7.20(-4)	7.02(-4)
4	***	***	1.11(-3)	4.42(-4)	2.59(-4)	2.20(-4)	1.99(-4)	1.88(-4)
8	***	***	***	1.78(-4)	6.34(-5)	5.57(-5)	5.06(-5)	4.77(-5)
16	***	***	***	***	3.16(-5)	1.26(-5)	1.22(-5)	1.18(-5)
32	***	***	***	***	***	6.24(-6)	2.80(-6)	2.85(-6)
64	***	***	***	***	***	***	1.35(-6)	6.47(-7)
128	***	***	***	***	***	***	***	3.18(-7)

Table 4. Table 4: Errors in x 1 subscript 𝑥 1 x_{1} for SAM with respect to the true oscillatory solution for problem ( 1 ). Step points are stroboscopic times

N	$\Omega =8\pi$	$\Omega =16\pi$	$\Omega =32\pi$	$\Omega =64\pi$	$\Omega =128\pi$	$\Omega =256\pi$	$\Omega =512\pi$	$\Omega =1024\pi$
1	6.10(-2)	3.30(-2)	1.66(-2)	7.74(-3)	3.09(-3)	1.66(-3)	2.04(-3)	2.30(-3)
2	***	6.65(-3)	3.41(-3)	1.56(-3)	8.31(-4)	7.57(-4)	7.20(-4)	7.02(-4)
4	***	***	7.95(-4)	3.56(-4)	2.63(-4)	2.21(-4)	1.99(-4)	1.88(-4)
8	***	***	***	9.25(-5)	6.62(-5)	5.64(-5)	5.07(-5)	4.77(-5)
16	***	***	***	***	1.50(-5)	1.34(-5)	1.23(-5)	1.18(-5)
32	***	***	***	***	***	3.03(-6)	2.95(-6)	2.88(-6)
64	***	***	***	***	***	***	6.44(-7)	6.76(-7)
128	***	***	***	***	***	***	***	1.43(-7)

Table 5. Table 5: Errors between SAM solution and the averaged solution for problem ( 32 ). (H2) holds

N	$\Omega =8\pi$	$\Omega =16\pi$	$\Omega =32\pi$	$\Omega =64\pi$	$\Omega =128\pi$	$\Omega =256\pi$	$\Omega =512\pi$
1	3.08(-2)	3.00(-2)	2.97(-2)	2.95(-2)	2.95(-2)	2.95(-2)	2.94(-2)
2	***	8.41(-3)	8.36(-3)	8.35(-3)	8.34(-3)	8.33(-3)	8.33(-3)
4	***	***	2.21(-3)	2.26(-3)	2.28(-3)	2.29(-3)	2.29(-3)
8	***	***	***	5.64(-4)	5.84(-4)	5.91(-4)	5.94(-4)
16	***	***	***	***	1.42(-4)	1.48(-4)	1.50(-4)
32	***	***	***	***	***	3.57(-5)	3.73(-5)
64	***	***	***	***	***	***	8.94(-6)

Table 6. Table 6: Errors in x 𝑥 x between SAM solution and the averaged solution for problem ( 32 ). Neither (H1) nor (H2) hold

N	$\Omega =8\pi +\frac{\pi }{64}$	$\Omega =16\pi +\frac{\pi }{32}$	$\Omega =32\pi +\frac{\pi }{16}$	$\Omega =64\pi +\frac{\pi }{8}$	$\Omega =128\pi +\frac{\pi }{4}$	$\Omega =256\pi +\frac{\pi }{2}$	$\Omega =512\pi +\pi$
1	3.08(-2)	3.00(-2)	2.97(-2)	2.95(-2)	2.95(-2)	2.95(-2)	2.95(-2)
2	***	8.43(-3)	8.38(-3)	8.36(-3)	8.36(-3)	8.36(-3)	8.36(-3)
4	***	***	2.23(-3)	2.28(-3)	2.30(-3)	2.31(-3)	2.32(-3)
8	***	***	***	5.78(-4)	6.00(-4)	6.12(-4)	6.23(-4)
16	***	***	***	***	1.58(-4)	1.69(-4)	1.79(-4)
32	***	***	***	***	***	5.63(-5)	6.60(-5)
64	***	***	***	***	***	***	3.76(-5)

Table 7. Table 7: Errors in x 1 subscript 𝑥 1 x_{1} for SAM with respect to the averaged solution for problem ( 33 ).

N	$\Omega =8\pi$	$\Omega =16\pi$	$\Omega =32\pi$	$\Omega =64\pi$	$\Omega =128\pi$	$\Omega =256\pi$	$\Omega =512\pi$	$\Omega =1024\pi$
1	4.10(-2)	4.09(-2)	4.08(-2)	4.08(-2)	4.08(-2)	4.07(-2)	4.07(-2)	4.07(-2)
2	***	8.08(-3)	7.89(-3)	7.81(-3)	7.77(-3)	7.76(-3)	7.75(-3)	7.75(-3)
4	***	***	1.78(-3)	1.70(-3)	1.67(-3)	1.65(-3)	1.64(-3)	1.64(-3)
8	***	***	***	4.29(-4)	4.08(-4)	3.99(-4)	3.99(-4)	4.06(-4)
16	***	***	***	***	1.06(-4)	1.01(-4)	1.04(-4)	1.07(-4)
32	***	***	***	***	***	2.62(-5)	2.50(-5)	2.68(-5)
64	***	***	***	***	***	***	6.53(-6)	6.47(-6)
128	***	***	***	***	***	***	***	1.63(-6)

Equations

$${\dot{x}}_1(t)$$

$${\dot{x}}_2(t)$$

$$\dot{x}(t)$$

$$y(t)$$

$$x(t)=\phi (t),-\tau \le t\le 0.$$

$$\dot{X}(t)$$

$$X(t)$$

$${\dot{X}}_1(t)$$

$${\dot{X}}_2(t)$$

$$X_{n+1}=X_n+\frac{3}{2}H{F}_n-\frac{1}{2}H{F}_{n-1},$$

$$X_1=X_0+H{F}_0,X_{N+1}=X_N+H{F}_N.$$

$$u_{n,\nu +1}=u_{n,\nu }+hf(u_{n,\nu },v_{n,\nu },t_n+\nu h,\Omega \nu h;\Omega ),\nu =0,1,\dots ,{\nu }_{max}-1,$$

$$u_{n,-\nu -1}=u_{n,-\nu }-hf(u_{n,-\nu },v_{n,-\nu },t_n-\nu h,-\Omega \nu h;\Omega ),\nu =0,1,\dots ,{\nu }_{max}-1.$$

$$F_n=\frac{u_{n,{\nu }_{max}}-u_{n,-{\nu }_{max}}}{2T}$$

$$F_0=\frac{u_{0,{\nu }_{max}}-u_{0,0}}{T},F_N=\frac{u_{N,{\nu }_{max}}-u_{N,0}}{T},$$

$$\dot{x}(t)=f(x(t),y(t),\Omega t).$$

$$f(x,y,\theta )=\sum _{k\in \mathbb{Z}}\exp (ik\theta )f_k(x,y).$$

$$\{\begin{array}{rl}x& =X+\frac{1}{\Omega }\sum _{k\ne 0}\frac{\exp (ik\Omega t)-1}{ik}{f}_k(X,Y),t\ge 0,\\ x& =X,-\tau \le t<0,\end{array},$$

$$\dot{X}$$

$$f_k=f_k(X,Y),f_{k\tau }=f_k(Y,Z),k\in \mathbb{Z}.$$

$$\dot{X}=f_0+\mathcal{O}(1/\Omega ),$$

$$f_0(X,Y)=\frac{1}{2\pi }{\int }_0^{2\pi }f(X,Y,\theta )d\theta .$$

$$X$$

$$[f_i,f_j]=\frac{\partial f_j}{\partial X}{f}_i-\frac{\partial f_i}{\partial X}{f}_j,i,j\in \mathbb{Z}.$$

$$\dot{X}(t)$$

$$Y(t)$$

$$Z(t)$$

$$X(t)$$

$$F^{(1)}(X,Y,t;\Omega )=f_0+\sum _{k>0}\frac{i}{k\Omega }\left([f_k-f_{-k},f_0]+[f_{-k},f_k]\right)-\sum _{k\ne 0}\frac{i}{k\Omega }\frac{\partial f_k}{\partial Y}\dot{\phi }(t-\tau ),$$

$$F^{(2)}(X,Y,Z;\Omega )$$

$$\underset{0\le t^{(k)}\le t_{max}}{\max }\Vert x(t^{(k)})-X(t^{(k)})\Vert =\mathcal{O}(\frac{1}{{\Omega }^2}),\Omega \to \infty .$$

$$F_n=\frac{u_n(T)-u_n(-T)}{2T},$$

Full text

1 Introduction

We propose and analyze a heterogenous multiscale method (HMM) (E, 2003; E & Engquist, 2003) for the efficient integration of constant-delay differential equations subject to fast periodic forcing. The stroboscopic averaging method (SAM) suggested here may provide approximations with $\mathcal{O}(H^{2}+1/\Omega^{2})$ errors with a computational effort that grows like $H^{-1}$ (the inverse of the stepsize), uniformly in the forcing frequency $\Omega$.

Delay differential equations with fast periodic forcing appear in a number of recent contributions to the nonlinear Physics literature. As shown by Daza et al. (2017), under fast periodic forcing, the Delayed Action Oscillator (Boutle et al. (2007)) that describes El Niño phenomenon may generate basins of attraction with the Wada property, i.e. each point on the boundary of one of the basins is actually on the boundary of all basins. The system

[TABLE]

($\tau$, $\alpha,\beta,\omega,A,B$ are constants) describes, for $A=0$, $B=0$ a time-delayed genetic toggle switch, a synthetic gene-regulatory network (Gardner et al., 2000). Studied by Daza et al. (2013a) is the phenomenon of vibrational resonance (Landa & McClintock, 2000), i.e. the way in which the presence of the high-frequency forcing $B\sin(\Omega t)$ enhances the response of the system to the low-frequency forcing $A\sin(\omega t)$. Jeevarathinam et al. (2011) and Daza et al. (2013b) investigate in a similar way the forced, delayed Duffing oscillator. A new kind of nonlinear resonance of periodically forced delay systems has recently been described by Coccolo et al. (2018).

The numerical integration of highly oscillatory differential equations with or without delay may be a very demanding task, as standard methods typically have to employ timesteps smaller than the periods present in the solution. For systems without delay, the literature contains many suggestions of numerical schemes specially designed to operate in the highly-oscillatory scenario; many of them are reviewed in Hairer et al. (2002). Perhaps counterintuitively, some of those methodologies take advantage of the large frequency and their efficiency actually increases with $\Omega$ (Iserles & Nørsett, 2005). On the other hand, schemes for highly oscillatory problems may suffer from unexpected instabilities and inaccuracies (Calvo & Sanz-Serna, 2009).

The algorithm suggested here is based on ideas presented, for systems without delay, by Calvo et al. (2011a) and Calvo et al. (2011b). In these references, given an oscillatory problem, a stroboscopically averaged problem is introduced such that, at the stroboscopic times $t^{(k)}=kT$, $T=2\pi/\Omega$, $k=0,1,\dots$, its solution $X(t)$ (approximately) coincides with the true oscillatory solution $x$. The stroboscopically averaged problem does not include rapidly varying forcing terms and therefore, if available in closed form, may be integrated numerically without much ado. The algorithms in Calvo et al. (2011a) and Calvo et al. (2011b) compute numerically values of $X$, without demanding that the user finds analytically the expression of the averaged system. More precisely, the algorithms only require evaluations of the right-hand side of the originally given oscillatory problem. The solution $X$ is advanced with a standard integrator (the macro-integrator) with a stepsize $H$ that, for a target accuracy, may be chosen to be independent of $\Omega$. When the macro-integrator requires a value $F$ of the slope $\dot{X}$, $F$ is found by numerical differentiation of a micro-solution $u$, i.e. a solution of the originally given oscillatory problem. While the micro-integrations to find $u$ are performed with stepsizes $h$ that are submultiples of the (small) period $T$, the corresponding computational cost does not increase as $\Omega\rightarrow\infty$, because $u$ is only required in windows of width $mT$, $m$ a small integer.

The extension of the material in Calvo et al. (2011a) and Calvo et al. (2011b) to systems with delay is far from trivial. A first difficulty stems from the well-known fact that, in the delay scenario, regardless of the smoothness of the equation, solutions may be non-smooth at points $t$ that are integer multiples of the (constant) delay. Therefore, the algorithm presented here has to make special provision for that lack of smoothness. In addition, the analysis of the algorithm (but, as emphasized above, not the algorithm itself) is built on the knowledge of the stroboscopically averaged systems. While the construction of a stroboscopically averaged system with errors $x(t^{(k)})-X(t^{(k)})=\mathcal{O}(1/\Omega)$ is not difficult, here we aim at errors $x(t^{(k)})-X(t^{(k)})=\mathcal{O}(1/\Omega^{2})$ and this requires much additional analysis. The classical reference Lehman & Weibel (1999) only considers zero-mean, $\mathcal{O}(1/\Omega)$ averaging.

In Section 2 we present the main ideas of the paper and a detailed description of the algorithm. Due to the difficulties imposed by the lack of smoothness in the solution, the algorithm uses low-order methods: the second-order Adams-Bashforth formula as a macro-integrator and Euler’s rule as a micro-integrator. Section 3 contains the construction of the stroboscopically averaged system with $\mathcal{O}(1/\Omega^{2})$ accuracy. Sections 4 and 5 are devoted to the analysis of the SAM algorithm. In the first of these, we assume that the micro-integrations carried out in the algorithm are performed exactly. Under suitable hypotheses, the errors in SAM are $\mathcal{O}(H^{2}+1/\Omega^{2})$. The effect of the errors in the micro-integration is studied in Section 5: it is shown that, with a computational cost that grows like $1/H$, SAM may yield errors of size $\mathcal{O}(H^{2}+1/\Omega^{2})$. The $H^{2}$ (second order) behaviour of the error may come as a surprise, because micro-integrations are performed by Euler’s rule; of key importance here is a superconvergence result (see the bound in (27)) for the Euler solution of oscillatory problems when the integration is carried out over a whole number of periods. The last two sections report numerical experiments that, on the one hand, confirm the theoretical expectations and, on the other, show the advantage of SAM when compared with a direct numerical integration of the oscillatory problem. Speed-ups larger than three order of magnitude are reported at the end of Section 7.

2 The stroboscopic averaging method (SAM)

This section motivates and describes the SAM algorithm.

2.1 Motivation

We consider highly oscillatory delay differential systems of the form

[TABLE]

where the solution $x$ is defined for $t\geq-\tau$ and takes values in $\mathbb{R}^{D}$, the function $f(x,y,t,\theta;\Omega):\mathbb{R}^{D}\times\mathbb{R}^{D}\times\mathbb{R}\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{D}$ depends $2\pi$–periodically on its fourth argument $\theta$, $\tau>0$ is the (constant) delay, and the frequency $\Omega$ is seen as a parameter, $\Omega\gg 1$. Note that $f$ depends explicitly on $t$ through its third and fourth arguments; the fourth is the rapidly rotating phase $\theta=\Omega t$ and the third corresponds to a slow (i.e. $\Omega$-independent) dependence on $t$ (see the toggle-switch equations above). The values of $x$ on the interval $[-\tau,0]$ are prescribed through an $\Omega$*–independent111Assuming that $\varphi$ does not depend on $\Omega$ implies no loss of generality, as the general case may be reduced to the $\Omega$-independent case by introducing a new dependent variable $x(t)-\Phi(t)$, where $\Phi(t)$ coincides with $\varphi(t)$ for $-\tau\leq t\leq 0$. * function $\varphi$:

[TABLE]

It is well known that, regardless of the smoothness of $f$ and $\varphi$, the function $x(t)$ will typically not be differentiable at $t=0$ and that in (2) $\dot{x}(0)$ has to be understood as a right derivative. Furthermore the discontinuity of $\dot{x}(t)$ at $t=0$ will lead to the discontinuity of $\ddot{x}(t)$ at $t=\tau$, etc.

We assume that, at the stroboscopic times $t^{(k)}=kT$, where $T=2\pi/\Omega$ is the period and $k=0,1,\dots$, the solution $x(t)$ of the oscillatory delay problem (2)–(3) may be approximated (in a sense to be made precise later) by the solution $X(t)$ of an averaged problem

[TABLE]

where the value of the function $F$ may depend on $X(t)$, on the history $X(s)$, $-\tau\leq s<t$, on the slow time $t$ and on $\Omega$ but is independent of the fastly varying phase $\theta=\Omega t$.

We illustrate the situation in the particular case of the toggle-switch problem (1). Fig. 1 displays, in a short time interval, a solution of the given oscillatory system and the corresponding averaged solution, found by solving the stroboscopically averaged system

[TABLE]

where $Y_{i}(t)=X_{i}(t-\tau)$ and $1_{\{t\geq\tau\}}(t)$ is the (indicator) function that takes the value 1 for $t\geq\tau$ and vanishes for $t<\tau$ (see the next section for the derivation of this averaged system). Note that the slow time-dependent forcing $A\sin(\omega t)$ has not been averaged out and that, due to the presence of $1_{\{t\geq\tau\}}(t)$, the right-hand side of the averaged system is discontinuous at $t=\tau$ (this discontinuity manifests itself in Fig. 1 through a discontinuity in the slope of $X$, not to be confused with the discontinuity at $t=0$ typically present, as mentioned above, in solutions of delay problems). The notion of stroboscopic averaging studied in detail in Chartier et al. (2012) is far from new. However standard treatments of the theory of averaging favour alternative techniques, specially the zero-mean approach where the functions necessary to express the required change of variables are chosen so as to have zero mean over one period of the oscillation. In stroboscopic averaging the freedom available in the choice of change of variables is used to impose that the old and new variables coincide at stroboscopic times. This is advantageous for the numerical methods studied here. The zero-mean approach may be better for analytic purposes at it usually leads to simpler high-order averaged systems. If zero-mean averaging had been used in Fig. 1, the averaged solution would had been located halfway between the maxima and the minima of the oscillatory solution.

The study of the vibrational resonance of (1) requires to simulate over long time intervals (the interval $0\leq t\leq 300$ is used in Daza et al. (2013a)) for many choices of the values of the constants $\alpha,\beta,\omega,A,B,\tau$ and the parameter $\Omega$; the presence of the fast-frequency oscillations makes such a task costly. It is then of interest to simulate, if possible, averaged systems like (4) rather than highly-oscillatory models like (2). However, obtaining $F$ analytically may be difficult or even impossible and we wish to have a numerical method that approximates $X$ by using only $f$. SAM is such a technique.

The idea behind SAM is as follows. Let $x$ and $X$ be respectively the oscillatory and averaged solutions corresponding to the same $\varphi$. At fixed $t$, $\dot{x}$ and $\dot{X}$ may differ substantially. However difference quotients such as $(x(t+T)-x(t))/T$ or $(x(t+T)-x(t-T))/(2T)$ may provide a good approximation to the slope $\dot{X}(t)$ (see Fig. 1). As other heterogeneous multiscale methods (see e.g. E (2003); E & Engquist (2003); E et al. (2007); Engquist & Tsai (2005); Li et al. (2007); Ariel et al. (2009); Sanz-Serna (2009); Calvo & Sanz-Serna (2010)), the algorithm includes macro-integrations and micro-integrations. Macro-integations are used to advance $X$ over macro-steps of length $H$ larger than the period $T$. The necessary slopes $\dot{X}$ are obtained by forming difference quotients of auxiliary oscillatory solutions found by micro-integrations with small steps $h$.

2.2 The algorithm

Let us now describe the algorithm.

2.2.1 Macro-integration

Choose a positive integer $N$ and define the macro-stepsize $H=\tau/N$. If the solution is sought in an interval $0\leq t\leq t_{max}$, SAM generates approximations $X_{n}$ to $X(t_{n})$, $t_{n}=nH$, $n=0,1,\dots,\lfloor t_{max}/H\rfloor$ by using the second-order Adams-Bashforth formula (macro-integrator)

[TABLE]

starting from $X_{0}=x(0)=X(0)=\varphi(0)$; here $F_{n}$ is an approximation to $\dot{X}(t_{n})$ obtained by numerical differentiation of the micro-solution. The formula is not used if $n=0$ and $n=N$, where it would be inconsistent in view of the jump discontinuities of $\dot{X}$ at $t=0$ and $t=\tau$ noted above. For $n=0$ and $n=N$ we use Euler’s rule, i.e. we set

[TABLE]

2.2.2 Micro-integration

If $\nu_{max}$ is a positive integer, the micro-stepsize $h$ is chosen to be $T/\nu_{max}$ (recall that $T=2\pi/\Omega$ denotes the period). We use Euler’s rule, starting from $u_{n,0}=X_{n}$, first to integrate forward the oscillatory problem (2) over one period

[TABLE]

and then to integrate backward over one period

[TABLE]

Here $v$ denotes the past values for $u$ given by $v_{n,\nu}=u_{n-N,\nu}$ if $n>N$, and $v_{n,\nu}=\varphi(-\tau+nH+\nu h)$ for $n<N$ (if $n=N$ $v_{N,\nu}=u_{0,\nu}$ for $\nu\geq 0$ and $v_{N,\nu}=\varphi(\nu h)$ for $\nu<0$).

It is crucial to observe the values $\Omega\nu h$ and $-\Omega\nu h$ used for the fast argument $\theta$ in (8) and (9) respectively; each micro-integration starts from $\theta=0$ rather than from the value $\theta=\Omega t_{n}$ which may perhaps have been expected. The reason for this is that in stroboscopic averaging, the resulting averaged system changes with the initial value of the phase $\theta$; to work with one and the same averaged system micro-integrations have to start from the value $\theta=0$ that the phase takes at the initial point of the interval $[0,t_{max}]$ (see Calvo et al. (2011a), Calvo et al. (2011b) for a detailed discussion).

The slopes $F_{n}$ to be used in (6) or (7) are given by the central difference formula,

[TABLE]

if $n\neq 0$ and $n\neq N$, while for $n=0$ and $n=N$, we use the forward difference formula

[TABLE]

due to the discontinuity of $\dot{X}(t)$ at $t=0$ and $t=\tau$. A detailed description of the algorithm is provided in Table 1.

Fig. 2 may help to better understand the procedure. The upper time axis corresponds to the macro-integration; all the information needed to obtain $X_{n+1}$ is $X_{n}$, $F_{n}$ and $F_{n-1}$ ($F_{n-1}$ is actually not required for $n=0,N$). The value of $F_{n}$ is derived by numerical differentiation of the micro-solution and passed to the macro-integrator to compute $X_{n+1}$. For fixed $n$, the computation of the Euler micro-solution $u_{n,\nu}$ uses the past values $v_{n,\nu}$ and the initial datum $u_{n,0}=X_{n}$ (note that $X_{n}$ is the most recent vector found in the macro-integration).

SAM only operates with macro-stepsizes $H$ that are submultiples of the delay $\tau$; this restriction is imposed to enforce that $t=\tau$ be a step-point to better deal with the discontinuity in slope at $t=\tau$ (see Fig. 1). In general the quotient $\tau/T$ will not be a whole number and then the step points $t_{n}$ will not be stroboscopic times; this is the case in the simulation in Fig. 1.

We emphasize that, given $N$ and $\nu_{max}$, the complexity of the algorithm is independent of $\Omega$. When $\Omega$ increases, the micro-stepsize $h$ decreases to cater for the more rapid variation of the oscillations but the window of width $2T$ (or $T$) for each micro-integration becomes correspondingly narrower.

Finally we point out that we have tested several alternative algorithms. For instance we alternatively performed the micro-integrations with the Adams-Bashforth second order method, or we used second-order forward approximations for $F_{0}$, $F_{N}$. While those modifications improve the accuracy of the results for small stepsizes, experiments reveal that they are not always beneficial for large stepsizes; therefore we shall not be concerned with them here.

3 The averaged system

In this section and the three that follow, we assume that the system (2) is of the particular form

[TABLE]

When comparing with the general format in (2) we note that now $f(x,y,\theta)$ has three arguments rather than five. The case where $f$ includes a slow explicit dependence on $t$, i.e. $f=f(x,y,t,\theta)$, may be trivially reduced to (12) by adding a component $x_{D+1}$ to the state vector $x\in\mathbb{R}^{D}$ and setting $\dot{x}_{D+1}=1$. The case where $f$ depends on $\Omega$, i.e. $f=f(x,y,\theta;\Omega)$, is taken up in the final section. We assume that $f(x,y,\theta)$ and the initial function $\varphi$ in (3) are of class $C^{3}$, and that the solution $x$ exists in the interval $[0,t_{max}]$ where $t_{max}$ is a constant (i.e. does not change with the parameter $\Omega$).

Using an approach similar to that in Chartier et al. (2010), Chartier et al. (2012), Chartier et al. (2015), we use a Fourier decomposition

[TABLE]

The coefficients $f_{k}(x,y)$ satisfy $f_{k}\equiv f_{-k}^{*}$ because the problem is real.

It is easily seen that, under the preceding hypotheses, $x$ undergoes oscillations of frequency $\Omega$ and amplitude $\mathcal{O}(1/\Omega)$ as $\Omega\rightarrow\infty$. To reduce the amplitude of the oscillations to $\mathcal{O}(1/\Omega^{2})$, we consider the near identity change of variables

[TABLE]

where we note that $x$ and $X$ coincide at stroboscopic times i.e. the change is stroboscopic. In what follows, we use the notation $z(t)=y(t-\tau),~{}t\geq\tau$, $Y(t)=X(t-\tau),t\geq 0$, $Z(t)=Y(t-\tau),~{}t\geq\tau$.

The proof of the following result is a straightforward but very lengthy exercise on Taylor expansion.

Lemma 3.1**.**

The change of variables (13) transforms the system (12) into

[TABLE]

for $0\leq t<\tau$, and into

[TABLE]

for $t\geq\tau$. Here

[TABLE]

Note that system in (14)–(15) is of the (expected) form

[TABLE]

where $f_{0}$ is the average

[TABLE]

By suppressing the remainder in (16), we obtain the averaged system $\dot{X}=f_{0}(X,Y)$ with $\mathcal{O}(1/\Omega)$ error; this is not sufficiently accurate for the values of $\Omega$ of interest and we take the averaging procedure to higher order. To do so, we perform an additional stroboscopic change of variables chosen so as to annihilate the oscillatory parts of the $\mathcal{O}(1/\Omega)$ terms in (14)–(15). Specifically we take:

[TABLE]

for $0\leq t<\tau$ and

[TABLE]

for $t\geq\tau$, where $f_{k}=f_{k}(\tilde{X},\tilde{Y})$ and $f_{k\tau}=f_{k}(\tilde{X},\tilde{Y})$ ($\tilde{Y}(t)=\tilde{X}(t-\tau)$, $\tilde{Z}(t)=\tilde{Y}(t-\tau)$).

Taking the last displays to (14)–(15) and discarding the $\mathcal{O}(1/\Omega^{2})$ remainder results in the averaged system (18)–(3.2) below, where, for notational simplicity, we have suppressed the tildes over $X$, $Y$ and $Z$ and $[\cdot,\cdot]$ is the Lie-Jacobi bracket or commutator defined by

[TABLE]

The averaged solution $X$ will obviously approximate $x$ at stroboscopic times with $\mathcal{O}(1/\Omega^{2})$ errors (arising from discarding the $\mathcal{O}(1/\Omega^{2})$ remainder); this implies that, at least for $\Omega$ sufficiently large, $X$ exists in the interval $[0,t_{max}]$. We then have proved the following result.

Theorem 3.2**.**

Consider the averaged problem

[TABLE]

with $F(X,Y,Z,t;\Omega)=F^{(1)}(X,Y,t;\Omega)$,

[TABLE]

for $0\leq t<\tau$ and $F(X,Y,Z,t;\Omega)=F^{(2)}(X,Y,Z;\Omega)$,

[TABLE]

for $t\geq\tau$. For $\Omega$ sufficiently large, the averaged solution $X$ exists in the interval $[0,t_{max}]$. Furthermore, for the approximation at stroboscopic times, we may write

[TABLE]

We point out that, for $t\geq\tau$, the dependence of $\dot{X}(t)$ on the past values $X(s)$, $-\tau\leq s<t$ is through both $Y(t)=X(t-\tau)$ and $Z(t)=X(t-2\tau)$; this is to be compared with the situation for the original oscillatory system (2), which does not include the delay $2\tau$. The functions $F^{(1)}$ and $F^{(2)}$ are of class $C^{2}$, but of course $F$ is discontinuous at $t=\tau$. By implication, $X(t)$ will be of class $C^{3}$ except at $t=0,\tau$ (where $\dot{X}(t)$ has a jump discontinuity), at $t=2\tau$ (where the second derivative jumps) and $t=3\tau$ (where the third derivative jumps).

In the particular case where $f$ does not depend on the delayed argument $y$, so that we are dealing with an ordinary differential system, (19) and (3.2) are in agreement with formula (64) of Chartier et al. (2012).

4 Error analysis: exact micro-integration

In this section we investigate the global error of the algorithm under the assumption that the micro-integration is exact, so that the macro-integration and the numerical differentiations performed to find the slopes $F_{n}$ are the only sources of error. This scenario is of course relevant when the micro-step $h$ is chosen to be very small. The errors due to the Euler micro-integration will be studied in the next section.

In order to avoid misunderstandings, we state that ‘exact micro-integration’ has to be understood as follows. Consider e.g. the computation of $X_{2}$, …, $X_{N}$ in Table 1; in this section we assume that

[TABLE]

where $u_{n}$ solves the problem $\dot{u}_{n}=f(u_{n},v_{n},\Omega t)$, $u_{n}(0)=X_{n}$, $v_{n}(t)=\varphi(-\tau+nH+t)$. Of course similar modifications of the algorithm in Table 1 have to be carried out for the computation of $F_{n}$ for other values of $n$.

We begin by proving a stability bound for the macro-integrator. We consider a sequence $\{\hat{X}_{n}\}$ of vectors in $\mathbb{R}^{D}$ such that $\hat{X}_{-n}=\varphi(-nH)$, $n=1,\dots,N$, and furthermore satisfy the macro-integration equations with residuals $\{\hat{\rho}_{n}\}$, i.e.

[TABLE]

Furthermore we consider a second sequence $\{\tilde{X}_{n}\}$ with $\tilde{X}_{-n}=\varphi(-nH)$, $n=1,\dots,N$, satisfying the macro-integration equations above with residuals $\{\tilde{\rho}_{n}\}$, rather than $\{\hat{\rho}_{n}\}$.

Proposition 4.1**.**

To each bounded set $B\subset\mathbb{R}^{D}$, there corresponds a constant $C>0$, independent of $H$ and $\Omega$, such that for any sequences $\{\hat{X}_{n}\}$ and $\{\tilde{X}_{n}\}$ as above contained in $B$,

[TABLE]

Proof 4.2**.**

From the hypotheses in the preceding section, $F$ is a Lipschitz continuous function of $X$, $Y$ and $Z$ with a Lipschitz constant that is uniform as $t$ varies in the interval $0\leq t\leq t_{max}$ and $X$, $Y$ and $Z$ vary in $B$. The stability bound is then proved in a standard way by recurrence.

In our next result we investigate the consistency of the formulas (10)–(11) used to recover the slopes $F_{n}$. There are four cases corresponding to the four successive blocks in Table 1.

Proposition 4.3**.**

The following results hold:

If the function $u_{n}$ solves the problem

[TABLE]

then, with $Y_{n}=\varphi(-\tau+t_{n})$,

[TABLE] 2. 2.

If $u_{n}$ is as in the preceding item, then

[TABLE]

where $f_{k}$, $k\in\mathbb{Z}$ are evaluated at $(X_{n},Y_{n})$, $Y_{n}=\varphi(-\tau+t_{n})$. 3. 3.

If $u_{n}$ and $v_{n}$ satisfy

[TABLE]

with $w_{n}$ a continuously differentiable function, then

[TABLE] 4. 4.

If $u_{n}$ and $v_{n}$ are as in 3., then

[TABLE]

where $f_{k}$, $k\in\mathbb{Z}$ are evaluated at $(X_{n},Y_{n})$ and $f_{k\tau}$ stands for $f_{k}(Y_{n},Z_{n})$.

Proof 4.4**.**

Here we prove the fourth case; the other proofs are similar (and slightly simpler). For simplicity the subindex $n$ is dropped. We rewrite the equation for $u$ in integral form

[TABLE]

and use Picard’s iteration. Clearly for $-T\leq t\leq T$,

[TABLE]

We take this equality to the integral equation above and integrate with respect to $s$ to find

[TABLE]

where

[TABLE]

Then,

[TABLE]

By Taylor expanding $f$ at $X,Y$, we obtain, for $-T\leq t\leq T$:

[TABLE]

The proof concludes by evaluating this expression at $t=\pm T$ and taking those values to the left-hand side of (24).

According to this result, (21) and (23) provide approximations with $\mathcal{O}(1/\Omega)$ errors to (19) and (3.2) (evaluated at $X=X_{n}$, $Y=Y_{n}$, $t=t_{n}$) respectively, as expected from the use of forward differencing. Similarly (22) approximates (19) (at $X=X_{n}$, $Y=Y_{n}$, $t=t_{n}$) with $\mathcal{O}(1/\Omega^{2})$ error, as expected of central differences. However, when comparing (24) with (3.2) (at $X=X_{n}$, $Y=Y_{n}$, $Z=Z_{n}$, $t=t_{n}$), we observe that the last sum in (3.2) has a factor $\exp(ik\Omega\tau)$, which is not present in the last sum in (24) and therefore the error is, in general, only $\mathcal{O}(1/\Omega)$. To achieve $\mathcal{O}(1/\Omega^{2})$ errors we may assume that the functions $f_{k}(X,Y)$, $k\neq 0$, are independent of the second argument $Y$, i.e. the delay argument $y$ only appears in $f$ through $f_{0}$ (this is the case in (1)). Alternatively, we may assume that, for all $k\neq 0$, $\exp(ik\Omega\tau)=1$, i.e. that $\tau$ is an integer multiple of the period $T=2\pi/\Omega$.

We are now ready to give the main result of this section.

Theorem 4.5**.**

Assume that the problem (12), with the smoothness assumptions stated in the preceding section, is integrated with SAM with exact micro-integrations. In addition assume that one of the following hypotheses holds:

• •

(H1) The oscillatory Fourier coefficients $f_{k}$, $k\neq 0$ of $f$ do not depend on the delayed argument $y$.

• •

(H2) The delay $\tau$ is an integer multiple of the period $T=2\pi/\Omega$.

Then there exist constants $\Omega_{0}$, $H_{0}$ and $K$ such that, for $\Omega\geq\Omega_{0}$ , $H\leq H_{0}$, $0\leq t_{n}=nH\leq t_{max}$, the difference between the numerical solution and the solution of the averaged problem has the bound

[TABLE]

Proof 4.6**.**

We apply Proposition 4.1 with $B$ taken as a large ball containing the trajectory $X(t)$, $0\leq t\leq t_{max}$ in its interior; the vectors $X(t_{n})$ play the role of $\tilde{X}_{n}$ and the vectors $X_{n}$ play the role of $\hat{X}_{n}$. It is clear that each $\tilde{X}_{n}$ is in $B$. For $H$ sufficiently small and $\Omega$ sufficiently large the same will be true for each $\hat{X}_{n}$; this is established by means of a standard argument by contradiction using (25).

Each residual $\tilde{\rho}_{n}$ is $\mathcal{O}(H^{2})$ with the exceptions of $\tilde{\rho}_{0}$, $\tilde{\rho}_{N}$ and $\tilde{\rho}_{2N}$; these are only $\mathcal{O}(H)$ because the first two correspond to Euler steps and in the third there is a jump discontinuity in the second time derivative of the averaged solution. According to Proposition 4.3, the residuals $\hat{\rho}_{n}$ for the numerical solution are $\mathcal{O}(1/\Omega^{2})$, with the exceptions of $\hat{\rho}_{0}$ and $\hat{\rho}_{N}$, which are of size $\mathcal{O}(1/\Omega)$. Taking these results to Proposition 4.1, we get a global error bound of the desired form.

Remark 1. It is clear that the bound in (25) may be replaced by one of the form $K^{\prime}(H^{2}+1/\Omega^{2})$. We prefer the form (25), as it relates to three sources of error: the macro-integration $H^{2}$ error, the error $H/\Omega$ arising from differencing at [math], $\tau$, $2\tau$, and the error from second-order differencing at all other step points $t_{n}$.

Remark 2. The discrepancy between (24) and (3.2) that leads to the introduction of (H1) and (H2) stems from the fact that the values of the phase $\theta$ at $t_{n}$ and $t_{n-N}=t_{n}-\tau$ are in general different in the oscillatory problem but the same in SAM. The hypothesis (H1) holds in most applications; in fact in all the papers mentioned in the introduction, the delay system is obtained by adding to the right-hand side of an oscillatory ordinary differential system a linear feedback term $My$, with $M$ a constant matrix. (H2) is relevant in those studies where there is freedom in choosing the exact value of the large frequency $\Omega$.

Remark 3. If (H1) and (H2) do not hold, the same proof yields for our algorithm a bound of the form $K(H^{2}+1/\Omega)$ under the assumption that $f$ and $\varphi$ are $C^{2}$ functions. Numerical experiments reported in Section 6 reveal that in that case the bound cannot be improved to $K(H^{2}+1/\Omega^{2})$. However if (H1) and (H2) do not hold, errors of size $K(H^{2}+1/\Omega)$ may be obtained by means of a simpler algorithm based on applying forward differences at each step point $t_{n}$; obviously that alternative algorithm does not require the backward integration legs (9).

5 Error analysis: micro-integration errors

We now take into account the errors introduced by the Euler micro-integration. We begin with an auxiliary result. Note the improved error bound at the end of the integration interval.

Proposition 5.1**.**

Consider the application of Euler’s rule with constant step size $h=T/\nu_{max}$ to integrate in the interval $0\leq t\leq T$ the initial value problem $\dot{u}=f(u,v,\Omega t)$, $u(0)=X$, where $v$ is a known $C^{1}$ function. Denote by $u_{\nu}$ the Euler solution at $t=\nu h$. There are constants $C$, $\Omega_{0}$ and $h_{0}$ such that for $h\leq h_{0}$, $\Omega\geq\Omega_{0}$, the following bounds hold:

[TABLE]

Proof 5.2**.**

A standard error bound for Euler’s rule is

[TABLE]

where $L$ is the Lispchitz constant of $f$ with respect to $u$ in a neighbourhood of the solution and $M$ is an upper bound for $\|(1/2)\ddot{u}(t)\|$, $0\leq t\leq T$. In the present circumstances we have to take into account that, as $\Omega\rightarrow\infty$, the length $T=2\pi/\Omega$ of the integration interval decreases and $M$ grows like $\Omega$, because

[TABLE]

From the elementary inequality $(\exp(LT)-1)/L\leq T\exp(LT)$ and the standard bound, we have

[TABLE]

and therefore (26) holds.

By adding all the Euler equations we find

[TABLE]

and from (26),

[TABLE]

a relation that has to be compared with

[TABLE]

The bound (27) will be established if we show that the quadrature sum in (29) approximates the integral in (30) with errors of size $\mathcal{O}(hT)$. To this end we decompose the function being integrated as

[TABLE]

It is easily seen (for instance by expanding in a Fourier series) that the total time derivative $(d/dt)f_{2}$ remains bounded as $\Omega\rightarrow\infty$; elementary results then show that the quadrature of $f_{2}$ has errors of the desired size $\mathcal{O}(hT)$. On the other hand, the time derivative of $f_{1}$ grows like $\Omega$, and quadrature errors of size $\mathcal{O}(h)$ may be feared. Fortunately the quadrature for $f_{1}$ is actually exact, because one checks by an explicit computation that it is exact for each Fourier mode $f_{k}(X,v(0))\exp(ik\Omega s)$.

The third bound (28) is a trivial consequence of (27).

It goes without saying that the corresponding result holds for backward integrations with $-T\leq t\leq 0$. The following theorem provides the main result of this paper.

Theorem 5.3**.**

Assume that the problem (12), with the smoothness assumptions stated in the preceding sections, is integrated with SAM. In addition assume that one of the following hypotheses holds:

• •

(H1) The oscillatory Fourier coefficients $f_{k}$, $k\neq 0$ of $f$ do not depend on the delayed argument $y$.

• •

(H2) The delay $\tau$ is an integer multiple of the period $T=2\pi/\Omega$.

Then there exist constants $\Omega_{0}$, $H_{0}$, $h_{0}$ and $K$ such that, for $\Omega\geq\Omega_{0}$ , $H\leq H_{0}$, $h\leq h_{0}$, $0\leq t_{n}=nH\leq t_{max}$, the difference between the numerical solution and the solution of the averaged problem may be bounded as follows

[TABLE]

In particular, if the grids are refined in such a way that $h$ is taken proportional to $H/\Omega$ (i.e. $\nu_{max}$ is taken proportional to $N$), then the bound becomes

[TABLE]

Proof 5.4**.**

We argue as in Theorem 4.5. Now in the residual for the numerical solution $\{X_{n}\}$ we have to taken into account the micro-integration error. For $n\leq N$, the bound (28) (and the corresponding bound for the backward integration) show that the micro-integration adds a term of size $\mathcal{O}(h)$ to the residual. For $n>N$ the situation is slightly more complicated, because the algorithm uses past values $v_{n,\nu}$ that are themselves affected by micro-integration errors. However the stability of the micro-integrator guarantees that even when those errors are taken into account an estimate like (28) holds.

We recall that taking $h$ proportional to $H/\Omega$ makes the complexity of the algorithm independent of $\Omega$.

6 Numerical experiments

We now present experiments that illustrate the preceding theory.

We first integrated with SAM the toggle switch problem (1) with $\alpha=2.5$, $\beta=2$, $A=0.1$, $\omega=0.1$, $B=4.0$, and $\tau=0.5$. The function $\varphi$ is constant: $x_{1}(t)=0.5$, $x_{2}(t)=2.0$, $-\tau\leq t\leq 0$, which corresponds to the system staying at an equilibrium up to time $t=0$ and then switching on the slow and fast oscillatory forcing terms. For the macro-stepsize we set $H=\tau/N$, $N=1,2,4,\dots$ and the micro-stepsize was chosen as $h=T/(2N)$ (note that for the coarsest macro-step size $N=1$, there are only two Euler steps per period). Simulations took place in the interval $0\leq t\leq 2$ that includes the locations $t=0$, $t=\tau$ and $t=2\tau$ where the growth in the error of the SAM solution is more pronounced due to the lack of smoothness. SAM has no difficulty in coping with the longer intervals required in practical simulations, but we have chosen a short interval because on longer intervals it may be very expensive to obtain a sufficiently accurate reference solution to measure errors; see in this connection the computational times quoted at the end of Section 7.

Table 2 reports, for varying $\Omega$ and $N$, the maximum error, over $0\leq t\leq 2$, in the $X_{1}$ component of the SAM solution with respect to the averaged solution obtained by integrating (5) (this integration was carried out with the Matlab function dde23 with relative and absolute tolerances $10^{-8}$ and $10^{-10}$ respectively). The combinations of $N$ and $\Omega$ leading to values of $H$ not significantly larger than $T$ were not attempted, as the HMM idea does not make sense for them. Note that here $\tau/T$ is irrational and therefore the step points $t_{n}$ are not stroboscopic times. Fig. 3 displays the errors in Table 2 as functions of $N$; for clarity not all values of $\Omega$ are included. By looking at the columns of the table (or at each of the four solid lines in the figure), we see that the error behaves as $N^{-2}$, i.e as $H^{2}$, except at the bottom of each column, where the behaviour is as $N^{-1}$. This is of course the behaviour of the bound in (31). Errors along rows saturate if $\Omega$ is very large; for those values one just observes the error in the macro-integration. This behaviour is also seen in the figure by comparing points corresponding to the same value of $N$ and varying $\Omega$. Along the main diagonal of the table, errors approximately divide by four, which is also in agreement with the bound (31). In the figure this corresponds to observing the behaviour of the right-most point of each of the solid lines.

Table 3 differs from Table 2 in that now $\Omega$ is taken from the sequence $8\pi$, $16\pi$, … that consists of values not very different from those in Table 2. In fact the errors in Table 3 are very similar to those in Table 2. However for the sequence $8\pi$, $16\pi$, …  the step points are stroboscopic times and it makes sense to compare the SAM solution with the true oscillatory solution. The results are reported in Table 4. From Theorems 3.2 and 5.3 the errors with respect to the true solution possess a bound of the form (31) and this is consistent with the data in the table.

In order to check numerically that the hypotheses (H1)–(H2) are necessary to ensure a bound of the form (31), we have considered the simple scalar equation

[TABLE]

this is an academic example where (H1) does not hold (recall that in all the systems from the literature cited (H1) holds). The averaged version is

[TABLE]

for $0\leq t<\tau$ and

[TABLE]

for $t\geq\tau$.

We used $\tau=0.5$, $H=\tau/N$, $h=T/(5N)$, a constant $\varphi=0.1$ and, as before, measured errors in the maximum norm for $0\leq t\leq 2$. Table 5222For typographic reasons only, this table and the next have one column less than those presented before. There is nothing unexpected in the results corresponding to the omitted frequency $\Omega=1024\pi$ (or $\Omega=1024\pi+2\pi$). gives errors when $\Omega$ is taken from the sequence $8\pi$, $16\pi$, …, so that the periods $T$ are $1/4$, $1/8$, …Hypothesis (H1) does not hold, but (H2) does, so that Theorem 5.3 may be applied. In fact an $N^{-2}$ (or equivalently $\Omega^{-2}$) behaviour is seen along the diagonals of the table.

We next slightly changed the frequencies and used the sequence $8\pi+\pi/64$, $16\pi+\pi/32$, …(this represents in increase of less than $0.2\%$ in each frequency). Neither (H1) nor (H2) are fulfilled and Theorem 5.3 may not be applied. The results in Table 6 show that, for $\Omega$ large the second order behaviour along the main diagonal is lost, indicating that the error does not behave as in the bound (31). Note the substantial difference between Tables 5 and 6 at $N=64$, in spite of the very small relative change in the value of $\Omega$.

7 Extensions

We finally consider the application of SAM to problems that are not of the form (12). The number of variants that may arise is very high and we restrict the attention to reporting numerical results for a case study. The corresponding analysis may be carried out by adapting the proofs given in the preceding sections.

We study again the toggle switch problem, but now in an alternative asymptotic regime. The system is given by

[TABLE]

where $\hat{B}$ is a constant and the other symbols are as before. As $\Omega\rightarrow\infty$, the variable $x_{1}$ undergoes oscillations of frequency $\Omega$ and $\mathcal{O}(1)$ amplitude, which, for $\Omega$ large, makes the direct numerical integration of the system more expensive than that of (1) (the amplitude there is $\mathcal{O}(1/\Omega)$). For an analytic treatment, we begin by performing, for $t\geq 0$, the preliminary stroboscopic change of variables

[TABLE]

which differs from the identity in $\mathcal{O}(1)$ quantities. This leads to:

[TABLE]

The highly oscillatory forcing has been reduced from $\mathcal{O}(\Omega)$ to $\mathcal{O}(1)$ and, in principle, it is possible to average (34) by the techniques used to deal with (12). Unfortunately, finding the required Fourier coefficients in closed form does not appear to be possible in general. In the particular case where $\beta=2$, the [math]-Fourier coefficient may be found by evaluating the relevant integral (17) with the help of the residue theorem. This leads to the averaged system of the form $\dot{X}=f_{0}$ explicitly given by

[TABLE]

with $M=X_{1}^{2}+2\hat{B}X_{1}-1$, $N=M^{2}+4(X_{1}+\hat{B})^{2}$, whose solutions approximate $x$ with errors of size at least $\mathcal{O}(1/\Omega)$ at stroboscopic times. However, since (34) is even in $\Omega$, in this particular case the errors are actually $\mathcal{O}(1/\Omega^{2})$.

Table 7 presents the errors in the SAM solution measured with respect to the solution of (35). The experiments have $\hat{B}=0.1$; all other details are as in the toggle switch simulations in the preceding section. The table shows that the performance of SAM is very similar to that encountered in problem (1). We also measured errors with respect to the oscillatory solution and found that they are very close to those reported here, i.e. the situation is similar to that seen when comparing Tables 3 and 4.

Finally we mention that, for the choice of constants considered here, the integration (in the interval $0\leq t\leq 2$) of the oscillatory problem for $\Omega=1024\pi$ with dde23 in a laptop computer took more than 9,000 seconds. The corresponding SAM solution with the smallest value of $H$ took approximately one second. Since as pointed out before, the study of vibrational resonance requires integrations in time intervals two orders of magnitude larger than $0\leq t\leq 2$, for many choices of the values of the constants that appear in the model, it is clear that a direct numerical integration of the oscillatory problem is not feasible for large values of $\Omega$.

Acknowledgements. J.M. Sanz-Serna has been supported by projects MTM2013-46553-C3-1-P from Ministerio de Economía y Comercio, and MTM2016-77660-P(AEI/FEDER, UE) from Ministerio de Economía, Industria y Competitividad, Spain. Beibei Zhu has been supported by the National Natural Science Foundation of China (Grant No. 11371357 and No. 11771438). She is grateful to Universidad Carlos III de Madrid for hosting the stay in Spain that made this work possible and to the Chinese Scholarship Council for providing the necessary funds. The authors are thankful to M. A. F. Sanjuán and A. Daza for bringing to their attention the vibrational resonance phenomenon, the toggle switch problem and other highly-oscillatory systems with delay.