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

TL;DR

This paper introduces a word-series based method to derive high-order averaged equations for highly oscillatory delay differential equations, enabling better analysis and simulation of such systems.

Contribution

It extends word-series averaging techniques to delay differential equations, providing a new approach for high-order averaging in systems with delay and oscillations.

Findings

High-order averaged equations can be derived for delay differential equations.

The method simplifies the analysis of oscillatory systems with delay.

Potential applications in control and signal processing with delays.

Abstract

We show that, for appropriate combinations of the values of the delay and the forcing frequency, it is possible to obtain easily high-order averaged versions of periodically forced systems of delay differential equations with constant delay. Our approach is based on the use of word-series techniques to obtain high-order averaged equations for differential equations without delay.

Figures

Tables

Table 1. Table 1: Maximum errors in u 𝑢 u , with respect to the true solution, for the second-order and the third-order averaged solutions in the interval 0 ≤ t ≤ 2 0 𝑡 2 0\leq t\leq 2

	$\Omega =16\pi$	$\Omega =32\pi$	$\Omega =64\pi$	$\Omega =128\pi$	$\Omega =256\pi$	$\Omega =512\pi$
AS2	3.31(-4)	8.83(-5)	2.28(-5)	5.77(-6)	1.46(-6)	3.66(-7)
AS3	1.04(-4)	1.28(-5)	1.59(-6)	1.96(-7)	2.07(-8)	2.30(-9)

Equations

$$\frac{d}{dt}\xi =g(\xi ,\Omega t),t\ge 0,\xi (0)={\xi }_0\in {\mathbb{C}}^D.$$

$$g(\xi ,\theta )=\sum _{k\in \mathbb{Z}}\exp (ik\theta )g_k(\xi ),$$

$$g_{k_1\cdots k_n}(\xi )=g_{k_2\cdots k_n}^{\prime }(\xi )g_{k_1}(\xi ),$$

$$W_{\delta }(\xi )=\sum _{w\in \mathcal{W}}{\delta }_w{g}_w(\xi ),$$

$$\xi (t)=W_{\kappa (\Omega t)}(\Xi (t)),$$

$$\frac{d}{dt}\Xi =W_{\bar{\beta }}(\Xi ),\Xi (0)={\xi }_0\in {\mathbb{R}}^D.$$

$$\frac{d}{dt}\Xi =g_0(\Xi ),$$

$$\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{)}.$$

$$\frac{d}{dt}x(t)$$

$$x(t)$$

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

$$x^{(0)}(t)$$

$$x^{(\ell )}(t)$$

$$x^{(\ell )}(0)=x^{(\ell -1)}(\tau ),1\le \ell \le L.$$

$$\frac{d}{dt}{x}^{(\ell )}(t)=f(x^{(\ell )}(t),x^{(\ell -1)}(t),\Omega (t+(\ell -1)\tau )),0\le t\le \tau ,1\le \ell \le L.$$

$$\frac{d}{dt}{x}^{(\ell )}(t)=f(x^{(\ell )}(t),x^{(\ell -1)}(t),\Omega t),0\le t\le \tau ,1\le \ell \le L.$$

$$\xi =(\hat{t},x^{(0)},x^{(1)},\dots ,x^{(L)}),$$

$$\frac{d}{dt}\hat{t}=1,\frac{d}{dt}{x}^{(0)}=\dot{\phi }(\hat{t}-\tau ),0\le t\le \tau ,$$

$$\hat{t}(0)=0,x^{(0)}(0)=\phi (-\tau ).$$

$$X^{(\ell )}(0)=X^{(\ell -1)}(\tau ),1\le \ell \le L,$$

$$X(t)=X^{(\ell )}(t-(\ell -1)\tau ),-\tau \le t\le L\tau ,\ell =0,\dots ,L,$$

$$g_0=\left[\begin{array}{c}1\\ \dot{\phi }(\hat{t}-\tau )\\ f_0(x^{(1)},x^{(0)})\\ f_0(x^{(2)},x^{(1)})\\ ⋮\\ f_0(x^{(L)},x^{(L-1)})\end{array}\right],g_k=\left[\begin{array}{c}0\\ 0\\ f_k(x^{(1)},x^{(0)})\\ f_k(x^{(2)},x^{(1)})\\ ⋮\\ f_k(x^{(L)},x^{(L-1)})\end{array}\right],k\ne 0,$$

$$\frac{d}{dt}\hat{T}$$

$$\frac{d}{dt}{X}^{(0)}$$

$$\frac{d}{dt}{X}^{(\ell )}$$

$$\frac{d}{dt}X(t)$$

$$X(t)$$

$$\frac{d}{dt}X(t)$$

$$F_{2,1}$$

$$F_{2,2}$$

$$\frac{du}{dt}$$

$$\frac{dv}{dt}$$

$$\frac{dU}{dt}$$

Full text

Mathematical Subject Classification (2010) 34C29

Keywords Delay differential equations, stroboscopic averaging, word series

1 Introduction

We show that, for appropriate combinations of the values of the delay and the forcing frequency, 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 [20, 14] 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, 15, 10, 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, 15, 10, 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.

2 Preliminaries

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

[TABLE]

Here the smooth function $g(\xi,\theta)$ is $2\pi$-periodic in its second argument $\theta\in\mathbb{R}$, with Fourier expansion

[TABLE]

and the parameter $\Omega\gg 1$ is the angular frequency of the fast periodic forcing. The corresponding period is $T=2\pi/\Omega$.

In the theory of word series the set of indices $\mathbb{Z}$ in (2) is seen as an (infinite) alphabet and for $n=1,2,\dots$, the elements of $\mathbb{Z}^{n}$ are regarded as words consisting of $n$ letters of this alphabet. There is also an empty word $\emptyset$ with $n=0$ letters. Given (1), to each word $w$ one associates the corresponding word basis function $g_{w}$, a map $\mathbb{C}^{D}\rightarrow\mathbb{C}^{D}$. By definition, for words with one letter $k\in\mathbb{Z}$, the word basis function $g_{w}(\xi)$ is the Fourier coefficient $g_{k}(\xi)$ in (2). For words with $n>1$ letters, we define recursively

[TABLE]

where $g^{\prime}_{k_{2}\cdots k_{n}}(\xi)$ is the Jacobian matrix of $g_{k_{2}\cdots k_{n}}(\xi)$. For the empty word $\emptyset$, the basis function is the identity map $\xi\mapsto\xi$.

The set of all words (including the empty word) is denoted by $\mathscr{W}$ and $\mathbb{C}^{\mathscr{W}}$ refers to the set of all mappings $\delta:\mathscr{W}\rightarrow\mathbb{C}$. For $\delta\in\mathbb{C}^{\mathscr{W}}$ and $w\in\mathscr{W}$, $\delta_{w}$ denotes the value of $\delta$ at $w$. To each $\delta\in\mathbb{C}^{\mathscr{W}}$ we associate the corresponding word series (relative to the mappings $g_{k}$ in (2)); this is the formal series

[TABLE]

whose terms are maps $\mathbb{C}^{D}\rightarrow\mathbb{C}^{D}$. The complex numbers $\delta_{w}$ are the coefficients of the series. As an example consider the case where $\delta_{\emptyset}=1$ and $\delta_{w}=0$ for each nonempty word $w$; in this case $W_{\delta}(\xi)=\xi$.

As proved in [7, 18] (see [19] for an alternative technique), there exist $\bar{\beta}\in\mathbb{C}^{\mathscr{W}}$ and a $2\pi$-periodic map $\theta\in\mathbb{R}\rightarrow\kappa(\theta)\in\mathbb{C}^{\mathscr{W}}$ such that the solution of (1) may be written as

[TABLE]

where $\Xi(t)$ is the solution of

[TABLE]

In this way, the time-dependent change of variables $\xi\mapsto\Xi$ 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 $\kappa_{w}(\theta)$ are such that, at $\theta=2k\pi$, $k=0,1,2,\dots$, we have $\kappa_{\emptyset}(\theta)=1$ and $\kappa_{w}(\theta)=0$ for each nonempty word $w$. This implies that, for those values of $\theta$, the transformation $\xi\mapsto\Xi$ is the identity map and it then follows that, in (4), $\xi(t)=\Xi(t)$ at the strobocopic times, $t=0,T,2T,\dots$ 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 $\xi\mapsto\Xi$ has 0 average over one period, rather than being the identity at $\theta=0$.

It is in order to point out that $\bar{\beta}$ and $\kappa(\theta)$ depend on $\Omega$, but are otherwise universal in the sense that they do not change with the dimension $D$ or with the choice of $g(\xi,\theta)$. 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 coeffcients. The values $\bar{\beta}_{w}$ and $\kappa_{w}(\theta)$ for $w\in\mathscr{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_{\bar{\beta}}(\Xi)$ 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_{\bar{\beta}}(\Xi)$ 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 $\Xi(t)$ coincides at stroboscopic times with the solution $\xi(t)$ except for an error of size $\mathcal{O}(1/\Omega^{n})$ as $\Omega\rightarrow\infty$.111To approximate $\xi$ 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) . The truncated averaged system with first order errors obtained by discarding contributions corresponding to words with two or more letters is found to be:

[TABLE]

as it may have been expected.

The truncated averaged system with second order errors $\mathcal{O}(1/\Omega^{2})$ is given by [7, 10, 18]

[TABLE]

Note that, $g_{0}$ is by definition the word-basis function associated with the one-letter words [math]. Furthermore, according to (3), the function $g_{0}^{\prime}g_{k}$ is the word basis function associated to the word $k0$ and similarly $g_{k}^{\prime}g_{0}$, $g_{k}^{\prime}g_{-k}$ 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 $\geq 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.

3 Highly oscillatory problems with delay

We now consider the $D$-dimensional system with constant delay $\tau>0$

[TABLE]

where $f(x,y,\theta)$ is $2\pi$-periodic in its third argument, with Fourier expansion

[TABLE]

and $\Omega$ is the angular frequency as above. Without loss of generality [22], it is assumed that the known function $\varphi$ that specifies the initial history is independent of $\Omega$. 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 $\tau$ is an integer multiple of the period $T$, or in other words that $t=\tau$ is a stroboscopic time. Let us discuss this hypothesis. In some applications, the value of $\tau$ is given and the interest is in studying the behaviour of the solution of (8)–(9) for large $\Omega$, but there is some freedom as to the choice of the specific value of $\Omega$. In those cases, the theory below applies by choosing $\Omega$ in such a way that $\tau\Omega/(2\pi)$ is an integer. Similarly, our theory applies to situations where the interest is in a given, large value of $\Omega$ and there is some freedom in the choice of $\tau$ to be used in the model. On the other hand, if the values of $\tau$ and $\Omega$ 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\leq t\leq L\tau$ 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\geq 0$.) We introduce the functions (see Fig. 1).

[TABLE]

and note that determining these functions is clearly equivalent to determining the solution of (8)–(9) in the interval $0\leq t\leq L\tau$. The $x^{(\ell)}(t)$ with $\ell>0$ satisfy the conditions

[TABLE]

and the differential equations

[TABLE]

Even though the relations (13) may be reminiscent of a two-point boundary value problem, we are facing here an initial value problem: $x^{(1)}(t)$ is determined by solving the differential equation with $\ell=1$ with initial condition $x^{(1)}(0)=\varphi(0)$, once $x^{(1)}(t)$ is known, $x^{(2)}(t)$ is determined by solving the differential equation with $\ell=2$ with initial condition $x^{(2)}(0)=x^{(1)}(\tau)$, etc.

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

[TABLE]

Note that this would be a system of the form (1) for the $D\times L$-dimensional vector $(x^{(1)},\dots,x^{(L)})$ if it were not for the fact that the right-hand side of the equation corresponding to $\ell=1$ contains $x^{(0)}$, 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 $\theta=\Omega t$, we proceed as follows. We introduce the $1+D(L+1)$-dimensional vector of unknown functions of the variable $t$

[TABLE]

and add to (14) the differential equations

[TABLE]

(the dot represents differentiation) subject to the initial conditions

[TABLE]

The solution of the initial value problem (15)–(16) is obviously $\widehat{t}=t$ and $x^{(0)}(t)=\varphi(t-\tau)$. 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^{(\ell)}$ is the averaged counterpart of $x^{(\ell)}$, we are interested in the solutions of the averaged system that satisfy (cf. (13))

[TABLE]

so that it is possible to define a continuous function $X$ in the interval $[-\tau,L\tau]$ by patching together the different $X^{(\ell)}$,

[TABLE]

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

3.2 The first-order averaged system

If $g$ denotes the right-hand side of the oscillatory system of differential equations (14)–(15) satisfied by $\xi$, the Fourier series for $g(\xi,\theta)$ (cf. (2)) has coefficients

[TABLE]

where the $f_{k}$ are the Fourier coefficients of $f$ as in (10). According to (6), the first-order averaged system of differential equation without delay in the interval $0\leq t\leq\tau$ is therefore given by (capital letters denote averaged dependent variables)

[TABLE]

This system is to be solved with the conditions $\widehat{T}(0)=0$, $X^{(0)}(0)=\varphi(-\tau)$ and (17). Clearly $\widehat{T}=t$, $X^{(0)}(t)=\varphi(t-\tau)$ and therefore the function $X$ in (18) satisfies the averaged delay problem (cf. (8)–(9))

[TABLE]

as one could have easily guessed. Note that $L$ does not feature in the averaged problem, as we pointed out above.

The solutions $x$ and $X$ of (8)–(9) and (19)–(20) differ at stroboscopic times by $\mathcal{O}(1/\Omega)$ as we prove next. That $x(t)-X(t)=\mathcal{O}(1/\Omega)$ at stroboscopic times $t\leq\tau$ follows by considering that, for those values of $t$, $x(t)=x^{(1)}(t)$, $X(t)=X^{(1)}(t)$ and $x^{(1)}(t)-X^{(1)}(t)=\mathcal{O}(1/\Omega)$ because of the properties of stroboscopic averaging of ordinary differential equations without delay. For stroboscopic times $t$ in the interval $\tau<t\leq 2\tau$ the situation is slightly more complicated because in addition to the error $\mathcal{O}(1/\Omega)$ arising from truncating (5) there is an additional source of error: in this interval $x(t)=x^{(2)}(t+\tau)$, $X(t)=X^{(2)}(t+\tau)$ but $x^{(2)}(t)$ and $X^{(2)}(t)$ do not share the same initial value at $t=0$. However the difference of initial values is $x^{(2)}(0)-X^{(2)}(0)=x^{(1)}(\tau)-X^{(1)}(\tau)$, which we know is of size $\mathcal{O}(1/\Omega)$. It follows that $x(t)-X(t)=\mathcal{O}(1/\Omega)$ also at stroboscopic times $\tau<t\leq 2\tau$. The iteration of this argument leads to the conclusion we seek, i.e. $x(t)-X(t)=\mathcal{O}(1/\Omega)$ whenever $t$ is a stroboscopic time. Of course the constant implied in the $\mathcal{O}$ notation of course depends on $t$ (and in general will grow as $t$ increases).

3.3 Second-order averaged system

If rather than using the truncated averaged system (6), we use (7), the procedure described above leads, after some simple algebra, to an averaged problem given by (20) and a delay differential equation of the form

[TABLE]

with

[TABLE]

and

[TABLE]

valid for $t\geq\tau$. In these expressions $\partial_{x}f_{k}(x,y)$ and $\partial_{y}f_{k}(x,y)$ denote the Jacobian matrices of $f_{k}(x,y)$ with respect to $x$ and $y$ respectively. By arguing as in the first order case, one proves that at stroboscopic times the difference between the solution $x(t)$ of (8)–(9) the solution $X$ of (20)–(22) is $\mathcal{O}(1/\Omega^{2})$.

3.4 Third- and higher order averaged systems

It is still feasible to obtain in closed form the third-order averaged system for delay problems starting from the corresponding expression for the case without delay given in [7, 10, 18]. The expression of the averaged delay system one obtains in this way is lengthy and will not be reproduced here. In fact for order 3 or higher, the best approach is to find recursively the word basis functions for the system (14)–(15) for the specific $f$ of interest.

For an averaged system of order $n$, the right-hand side of the system has $n$ expressions corresponding to the intervals $0\leq t<\tau$, …, $(n-2)\tau\leq t<(n-1)\tau$, and $t\geq(n-1)\tau$. (Recall that the second-order averaged system displayed above has two expressions.)

4 An example

As an illustration, we consider the system

[TABLE]

where $\alpha$ and $\beta$ are parameters. There are two periodic forcing terms; the forcing $A\sin(\omega t)$ is slow (i.e. the frequency $\omega$ is of moderate size) and $B\sin(\Omega t)$, $\Omega\gg 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(\omega t)$ as $A\sin(\omega\hat{t})$, where $\hat{t}$ is a new dependent variable defined by the differential equation $(d/dt)\hat{t}=1$ and the initial condition $\hat{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

[TABLE]

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

[TABLE]

for $t\geq\tau$. (As pointed out above, in general, the third-order averaged system has an analytic expression in $\tau\leq t<2\tau$ and a different analytic expression in $t\geq 2\tau$; for the problem at hand those two expressions happen to coincide.) The second-order averaged system may be retrieved from the last displays by omitting the terms that have an $\Omega^{2}$ 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 $\alpha=2.5$, $\beta=2$, $A=0.1$, $\omega=0.1$, $B=2.0$, and $\tau=0.5$, with constant history $u(t)=0.5$, $v(t)=2.0$, $-\tau\leq t\leq 0$ (an equilibrium of the unforced system), Table 1 presents the maximum error in $u$ at stroboscopic times in the short interval $0\leq t\leq 2$; in agreement with the theory, the error behavior is $\mathcal{O}(1/\Omega^{2})$ for second-order averaging and $\mathcal{O}(1/\Omega^{3})$ for third-order averaging.

In Fig. 2, the parameters are $\alpha=2.5$, $\beta=2$, $A=0.2$, $\omega=0.2$, $B=2.0$, $\Omega=4\pi$ and $\tau=0.5$, with constant history $u(t)=2.0$, $v(t)=0.5$, $-\tau\leq t\leq 0$ (again an equilibrium of the unforced system). The integration is perfomed in the interval $0\leq t\leq 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.

5 Conclusion

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 $\Xi$ 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 [22, 5].

Acknowledgements. J.M.S. was supported by project MTM2016-77660-P(AEI/ FEDER, UE) funded by MINECO (Spain). B. Z. is supported by the National Natural Science Foundation of China (Grant No. 11771438) and the Postdoctoral Fund Project of China (Grant No. 2018M641506) and in addition by the National Center for Mathematics and Interdisciplinary Sciences, CAS.