Unbounded sequences of stable limit cycles in the delayed Duffing equation: an exact analysis

TL;DR

This paper provides an exact analysis demonstrating that the delayed Duffing equation has an infinite sequence of stable, rapidly oscillating periodic solutions when the delay parameter satisfies a specific inequality.

Contribution

It offers the first exact proof of unbounded stable limit cycles in the delayed Duffing equation, contrasting with previous approximate methods.

Findings

Infinite sequence of stable limit cycles identified

Stable solutions occur for delays satisfying T^2<1.5π^2

Analysis is exact, not approximate

Abstract

The delayed Duffing equation $\ddot{x}(t)+x(t-T)+x^3(t)=0$ is shown to possess an infinite and unbounded sequence of rapidly oscillating, asymptotically stable periodic solutions, for fixed delays such that $T^2<\tfrac{3}{2}\pi^2$. In contrast to several previous works which involved approximate solutions, the treatment here is exact.

Equations

$$\ddot{x}(t)+x(t-T)+x(t{)}^3=0,$$

$$p_n=2T/n$$

$$\ddot{x}(t)+(-1{)}^{n}x(t)+x(t{)}^3=0,n=1,2,3,\dots ,$$

$$H(t)=\frac{1}{2}{\dot{x}}^2+\frac{1}{2}(-1{)}^n{x}^2+\frac{1}{4}{x}^4.$$

$$H>0$$

$$0<x_n(0)=:A_n,{\dot{x}}_n(0)=0,$$

$$t={\int }_{x_n(0)}^{x_n(t)}\frac{dx}{\dot{x}(t)}=\pm {\int }_{x_n(t)}^{A_n}\frac{dx}{\sqrt{2E-(-1{)}^n{x}^2-x^4/2}}.$$

$$\frac{1}{4}{p}_n={\int }_0^{A_n}\frac{dx}{\sqrt{\left(2E-(-1{)}^n{x}^2-x^4/2\right)}}.$$

$$H_n=E=\frac{1}{2}(-1{)}^n{A}_n^2+\frac{1}{4}{A}_n^4$$

$$\frac{1}{4}{p}_n={\int }_0^{A_n}\frac{dx}{\sqrt{\left(A_n^2-x^2\right)\left((-1{)}^n+A_n^2/2+x^2/2\right)}}.$$

$$x_n(t)=A_n\mathrm{cn}({\omega }_{n}t,m_n).$$

$$m_n=\frac{A_n^2}{2(A_n^2+(-1{)}^n)}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{\omega }_n=\sqrt{A_n^2+(-1{)}^n}.$$

$$p_n=4K/{\omega }_n.$$

$$x_n(t-T)=(-1{)}^n{x}_n(t).$$

$$\ddot{x}(t)+x(t)+x(t{)}^3=0.$$

$$x_n(t)=A_n\mathrm{cn}({\omega }_{n}t,m_n),$$

$$m_n=\frac{A_n^2}{2(A_n^2+1)}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{\omega }_n=\sqrt{A_n^2+1}.$$

$$x_n(t-T)=x_n(t),\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{for all even }n.$$

$$\ddot{x}(t)-x(t)+x(t{)}^3=0.$$

$$H=\frac{1}{2}{\dot{x}}^2-\frac{1}{2}{x}^2+\frac{1}{4}{x}^4.$$

$$x(t)=-x(t-p/2),$$

$$\ddot{x}+x(t-p/2)+x^3=0.$$

$$x_n(t-T)=-x_n(t),\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{for all odd }n.$$

$$x_n(t)=A_n\mathrm{cn}({\omega }_{n}t,m_n),$$

$$m_n=\frac{A_n^2}{2(A_n^2-1)}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{\omega }_n=\sqrt{A_n^2-1}.$$

$$p\equiv p_n=2T/n,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}n=1,2,3,\dots ,$$

$$2T/n=p=4K(m(A_n))/\omega (A_n)$$

$$\displaystyle\begin{split}p=\frac{\gamma}{\sqrt{\pi}}\Big{(}&A^{-1}-\left(\tfrac{1}{2}+4\pi^{2}/\gamma^{2}\right)A^{-3}+\left(\tfrac{1}{2}+6\pi^{2}/\gamma^{2}\right)A^{-5}-\left(\tfrac{5}{8}+9\pi^{2}/\gamma^{2}\right)A^{-7}+\\ &+\left(\tfrac{85}{96}+14\pi^{2}/\gamma^{2}\right)A^{-9}-\left(\tfrac{87}{64}+\tfrac{903}{40}\pi^{2}/\gamma^{2}\right)A^{-11}\Big{)}+\mathcal{O}\left(A^{-13}\right)\,.\end{split}$$

$$\displaystyle\begin{split}A=\frac{\gamma}{\sqrt{\pi}}\Big{(}&p^{-1}-\pi\left(\tfrac{1}{2}\gamma^{2}+4\pi^{2}\right)\gamma^{-4}p-2\pi^{4}\left(\gamma^{2}+16\pi^{2}\right)\gamma^{-8}p^{3}-\\ &-8\pi^{7}\left(3\gamma^{2}+56\pi^{2}\right)\gamma^{-12}p^{5}+\tfrac{1}{96}\pi^{4}\left(\gamma^{8}-36\,864\,\gamma^{2}\pi^{6}-737\,280\,\pi^{8}\right)\gamma^{-16}p^{7}+\\ &+\tfrac{1}{960}\pi^{5}\left(5\gamma^{10}+328\gamma^{8}\pi^{2}-6\,758\,400\,\gamma^{2}\pi^{8}-140\,574\,720\,\pi^{10}\right)\gamma^{-20}p^{9}\Big{)}+\\ &+\mathcal{O}\left(p^{11}\right)\,.\end{split}$$

$$\displaystyle\begin{split}p=\frac{\gamma}{\sqrt{\pi}}\Big{(}&A^{-1}+\left(\tfrac{1}{2}+4\pi^{2}/\gamma^{2}\right)A^{-3}+\left(\tfrac{1}{2}+6\pi^{2}/\gamma^{2}\right)A^{-5}+\left(\tfrac{5}{8}+9\pi^{2}/\gamma^{2}\right)A^{-7}+\\ &+\left(\tfrac{85}{96}+14\pi^{2}/\gamma^{2}\right)A^{-9}+\left(\tfrac{87}{64}+\tfrac{903}{40}\pi^{2}/\gamma^{2}\right)A^{-11}\Big{)}+\mathcal{O}\left(A^{-13}\right)\,.\end{split}$$

Full text

1 Introduction

This work concerns a differential-delay equation (DDE) known as the delayed Duffing equation

[TABLE]

where $T>0$ is the time delay. The existence of an infinite number of stable limit cycles, i.e. of asymptotically stable periodic solutions, in this DDE was first suggested in a paper by Wahi and Chatterjee [WaCha04]. Formally and to leading order, they performed the method of averaging and obtained a slow flow that predicted infinitely many stable limit cycles. In their DDE, the time delay was fixed at $T=1$. Mitra&al [MiChaBa17] studied the same DDE with an added linear stiffness. By assuming an approximate solution in harmonic form $x(t)=A\sin(\omega t)$, they claimed that the system exhibits an infinite number of stable limit cycles for any value of the time delay $T$. In a paper by Davidow&al [DaShaRa17], the same claim was supported by a) harmonic balance, b) Melnikov’s integral with Jacobi elliptic functions, and c) the introduction of damping.

Strictly speaking, however all these works on the delayed Duffing equation (1.1) were restricted to small amplitudes of the limit cycles. In our work, we present an exact treatment of (1.1), in the limit of unboundedly large amplitudes. In particular, the previously studied infinite sequences of “stable limit cycles” lose stability, eventually, for delays $T$ such that $T^{2}>\tfrac{3}{2}\pi^{2}$.

Section 2 gives a brief account of the numerical integration method used for our simulations. In section 3 we study exact periodic solutions $x_{n}(t)$ of a slightly generalized Duffing ordinary differential equation (ODE), with vanishing time delay $T=0$; see (3.1). We show how the non-delay ODE solutions $x_{n}(t)$ of minimal (or fundamental) periods $p_{n}$ lift to exact solutions of the original delayed Duffing DDE (1.1) with positive delay $T>0$, provided their minimal periods

[TABLE]

are integer fractions of the double delay $2T$. In particular we show how the more and more rapidly oscillating periodic solutions $x_{n}(t)$ develop unbounded amplitudes $A_{n}\nearrow\infty$, for $n\rightarrow\infty$. In section 4 we indicate how to determine the amplitudes $A_{n}$ of the lifted solutions $x_{n}$ , numerically and by series expansions for $n\rightarrow\infty$. Section 5 recalls our stability results from [Fie&al19]. These mathematical results basically assert local asymptotic stability of the solutions $x_{n}(t)$, for any fixed positive delay $T$ such that $T^{2}<\tfrac{3}{2}\pi^{2}$ and for sufficiently large odd $n=1,3,5,\dots$. They also show instability, for sufficiently large even $n=2,4,6,\dots$. For full mathematical details, including added linear stiffness, we refer to [Fie&al19]. We conclude with numerical illustrations of our results, in section 6, and a short summary 7.

Acknowledgment. Just as the more mathematically inclined account in [Fie&al19], the present work has originated at the International Conference on Structural Nonlinear Dynamics and Diagnosis 2018, in memoriam Ali Nayfeh, at Tangier, Morocco. We are deeply indebted to Mohamed Belhaq, Abderrahim Azouani, to all organizers, and to all helpers of this outstanding conference series. They indeed keep providing a unique platform of inspiration and highest level scientific exchange, over so many years, to the benefit of all participants. This work was partially supported by DFG/Germany through SFB 910 project A4. Authors RHR, BS and SMS gratefully acknowledge support by the National Science Foundation under grant number CMMI-1634664.

2 Numerical integration

For zero delays, $T=0$, the delayed Duffing DDE (1.1) reduces to a non-delayed ordinary differential equation (ODE) known as the classical Duffing equation. The equation is conservative and hence exhibits a continuum of periodic orbits, rather than any asymptotically stable limit cycles.

Even for arbitrarily small fixed positive delays, $T>0$, in contrast, approximate analysis and numerical simulations suggest that an infinite number of stable limit cycles may coexist, their amplitudes going to infinity [DaShaRa17].

Figure 2.1 shows the time history (a) and phase plane (b) of the first three stable limit cycles obtained by numerical integration of the delayed Duffing DDE (1.1), for $T=0.3$ and with different initial conditions. The numerical integrations in the present work were performed using the Python library pydelay for DDEs [Flu11]. The integrator is based on the Bogacki-Shampine method [BoSha89]. The maximal step size used to produce the plots in the present work was fixed at $\Delta t=10^{-4}$. See section 6 for further numerical examples.

3 Lifting periodic solutions from the non-delayed to the delayed Duffing equation

In this section we show the existence of infinitely many rapidly oscillating periodic solutions of specific periods $p$ in the delayed Duffing DDE  (1.1). Our approach is based on a lift of certain periodic solutions of the ordinary non-delayed Duffing ODE (3.1) below, with minimal (or, fundamental) period $p$, to periodic solutions of the delayed Duffing DDE (1.1) with time delay $T$. We will show this remarkable fact for minimal periods $p$ which are integer fractions of the doubled delay $2T=np$ ; see claim (1.2). We first recall some elementary facts on the non-delayed Duffing ODE, in subsection 3.1. We separately address the cases of even and odd fractions $n$ in subsections 3.2 and 3.3, respectively.

3.1 General Duffing equation

We consider the following two general forms of the classical Duffing ODE [KoBr11]:

[TABLE]

The time-independent Hamiltonian energy of (3.1) takes the form

[TABLE]

See Figure 3.1. For even $n$ the Hamiltonian is always positive; see Figure 3.1a,c. For odd $n$, however, see Figure 3.1b,d: the Hamiltonian is either strictly negative (green), identically zero (blue) or strictly positive (red), depending on the ODE initial conditions. Note how single trajectories in the $(x,\dot{x})$-plane are point symmetric to the origin, if and only if the positive energy condition

[TABLE]

is satisfied. We assume this restriction to hold throughout our further analysis.

For $H>0$, we may time-shift solutions $(x_{n}(t),\dot{x}_{n}(t))$ of (3.1) such that the initial conditions

[TABLE]

are satisfied. In particular, $A_{n}=\max|x_{n}(t)|$ is the amplitude of the solution $x_{n}$. For odd $n$, note how our positivity condition (3.3) requires an amplitude $A_{n}>\sqrt{2}$ in (3.2); see the red curve in Figure 3.1b,d, outside the blue figure-8 shaped separatrix loops. The periodic closed curves fill the part of the phase space $(x,\dot{x})$ where $H>0$. Each periodic orbit corresponds to specific initial conditions and possesses a specific minimal period.

The exact periodic solutions of the Duffing ODE $(x_{n}(t),\dot{x}_{n}(t))$ of (3.1) are easily determined. Indeed the energy $H\equiv E$ is identically constant. Solving (3.2) for $\dot{x}$ and classical separation of variables therefore lead to the elliptic integrals

[TABLE]

Here we have substituted the initial condition (3.4) for $x_{n}(0)$. The minimal (fundamental) period $p_{n}$ can be determined as the special case $t=p_{n}/4$, where symmetry implies $x_{n}(t)=0$:

[TABLE]

Evaluating the invariant Hamiltonian $H_{n}\equiv E$ at the initial condition (3.4) provides the energy

[TABLE]

and the explicit elliptic integral

[TABLE]

The elliptic integral (3.5) allows us to express the exact periodic solution of the general Duffing ODE (3.1) in terms of Jacobi elliptic function as

[TABLE]

Here $\mathrm{cn}$ denotes the Jacobi elliptic cosine function. The arguments $A_{n}$, $\omega_{n}$ and $0<m_{n}<1$ are the amplitude, the angular frequency, and the elliptic modulus, respectively. The frequency $\omega_{n}$ and the modulus $m_{n}$ in the solution (3.9) are related to the amplitude $A_{n}$ such that

[TABLE]

The minimal period (3.8) can be expressed in terms of the complete elliptic integral of the first kind $K\equiv K(m_{n})$ as

[TABLE]

See [Rand94]. Figure 3.2 indicates the relation between amplitude and frequency for the general Duffing ODE (3.1). The two black curves are obtained from the second equation of (3.10), and they correspond to the relation between amplitude and frequency of the periodic solutions (3.9) in the non-delayed Duffing ODE (3.1), for $n$ odd (upper curve) and $n$ even (lower curve). Each point represents a periodic orbit of the general Duffing ODE (3.1). In the phase plane, each of the black curves therefore indicates a foliation by periodic solutions. For the delayed Duffing equation (1.1), the same periodic solutions $x_{n}$ of minimal period $p_{n}=2T/n$ on the upper curve ($n$ odd) will turn out locally asymptotically stable, for $T^{2}<\tfrac{3}{2}\pi^{2}$ and large $n$, while large $n$ of even parity (lower curve) always turn out linearly unstable; see Theorems 5.1, 5.2 below.

Our lift construction from solutions of the non-delayed Duffing ODE (3.1) to the delayed Duffing ODE (1.1) is based on two interpretations of the mathematical expression $x(t-T)$. On the one hand, $x(t-T)$ represents a delay, as in (1.1). The same expression, on the other hand, represents a periodic solution when equated to $\pm x(t)$ by

[TABLE]

Here $2T$ represents any (not necessarily minimal) period of the periodic solution $x(t)$. Indeed, any positive energy solution of the Duffing ODE (3.1) is periodic and will automatically satisfy the periodicity condition (3.12), for some $T>0$. Upon substitution of the periodicity condition (3.12), however, the non-delayed Duffing ODE (3.1) produces the delayed Duffing DDE (1.1), where now the (half) period $T$ represents the delay. Thus any periodic solution of the Duffing ODE (3.1) with periodicity condition (3.12) lifts to a periodic solution of the DDE (1.1), for that choice of the delay $T$. The marked points (red) on the two black curves in Figure 3.2, for example, correspond to periodic solutions of the delayed Duffing DDE (1.1), with delay $T=3$.

Actually, the non-delayed ODE Duffing equation (3.1) possesses an uncountable continuum of periodic orbits, foliating the phase plane. The number of periodic orbits $x_{n}$ which satisfy the periodicity condition (3.12), however, is (at most) countable. In particular, our lift construction from the non-delayed Duffing ODE (3.1) to the delayed Duffing DDE (1.1) restricts the allowable points on the curves in Figure 3.2 to a countable set and therefore produces only a countable set of periodic solutions for the delayed Duffing DDE. We do not claim that our lift construction covers all possible periodic solutions of the DDE (1.1); in section 5 we will see indications of additional periodic solutions which cannot be obtained by our lift.

In the following we will further detail the lift construction (3.12) which is based on the known exact periodic solutions (3.9) of the general Duffing ODE (3.1). We consider the two cases, $n$ even and $n$ odd, separately.

3.2 Even $n$

For even $n$, the general Duffing ODE (3.1) reduces to the single-well case

[TABLE]

By (3.9), the exact periodic solutions are expressed as

[TABLE]

where now (3.10) becomes

[TABLE]

According to (3.11), minimal periods $p$ decrease monotonically from $p=2\pi$, at amplitude $A=0$, to $p=0$, for unbounded amplitudes $A\nearrow\infty$; see Figure 3.2.

To perform the lift from the Duffing ODE (3.13) to the Duffing DDE (1.1), we fix a time delay T (black dot $T=2$ in Figure 3.3), a priori, such that $T<2\pi$. Then we can always find a solution (3.9) to (3.13) with minimal period $p_{2}=T$; see solution $x_{2}(t)$ in Figure 3.3. If we shift the curve of $x_{2}(t)$ to the right by $T$ we obtain a new curve $x_{2}(t-T)$ that coincides with $x_{2}(t)=x_{2}(t-T)$; see Figure 3.3b. We can also find another solution $x_{4}(t)$, of larger amplitude $A_{4}>A_{2}$, whose minimal period is $p_{4}=T/2$. Shifting by $T$ we obtain a new curve $x_{4}(t-T)$ that coincides with $x_{4}(t)=x_{4}(t-T)$, see Figure 3.3b again. In the same manner, we can find infinitely many periodic solutions $x_{n}(t)$ with minimal periods $p_{n}=2T/n$, for $n=2,4,6,\dots$. After time shift by their shared (non-minimal) period $T$ we obtain

[TABLE]

Substituting (3.15) into (3.13) lifts all those ODE Duffing solutions $x_{n}(t)$ to the delayed Duffing DDE (1.1), for fixed delay $T<2\pi$. Note how $p_{n}\searrow 0$ implies unbounded amplitudes $A_{n}\nearrow\infty$, for $n\rightarrow\infty$. As the amplitudes $A_{n}$ of the periodic solutions of the Duffing ODE (3.13) increase to infinity, the minimal periods $p_{n}$ decrease to zero. Thus we obtain an unbounded sequence of more and more rapidly oscillating periodic solutions, with minimal periods $T,T/2,T/3,\dots$, which are also periodic with (non-minimal) period $T$. This proves our claim (1.2), for even $n$.

3.3 Odd $n$

For odd $n$, the general Duffing ODE (3.1) reduces to the double-well case

[TABLE]

Any solution conserves the Hamiltonian energy

[TABLE]

We recall how the phase portrait of the double-well Duffing ODE (3.16) is characterized by a figure-8 shaped separatrix $H=0$; see the blue curve in Figure 3.1b,d. For positive energy $H>0$, the (red) solutions of (3.16) oscillate around the exterior of the separatrix. Again, minimal periods $p$ decrease monotonically: this time from $p=\infty$, at the separatrix amplitude $A=\sqrt{2}$, to $p=0$, for $A\nearrow\infty$.

Since each level of positive energy $H>0$ consists of a single periodic orbit $(x,\dot{x})$, with odd force law, the time taken to travel from any point $(x,\dot{x})$ on a level set to its antipode $(-x,-\dot{x})$ is half its minimal period, $p/2$. Indeed this fact holds for any odd force law, by time reversibility of the oscillator. Therefore, every solution of the double-well Duffing ODE (3.16) with positive energy $H$ and minimal period $p$ satisfies the oddness symmetry

[TABLE]

for all $t$.

To perform the lift from the double-well Duffing ODE (3.16) to the delayed Duffing DDE (1.1), we now fix any time delay $T>0$ (black dot $T=2$ in Figure 3.4), this time without any further constraint. For $p_{1}:=2T$, the delay $T$ coincides with half the minimal period of the solution $x_{1}(t)$ of the non-delayed double-well Duffing ODE (3.16). The oddness symmetry (3.18) at half period $p_{1}/2=T$ therefore implies that $x_{1}(t)$ also solves our original delayed Duffing DDE (1.1),

[TABLE]

Analogously, we can perform the lift from the non-delayed double-well Duffing ODE (3.16) to the delayed Duffing DDE (1.1), for any odd $n=1,3,5,\dots$, as follows. Let $x_{n}$ denotes the ODE solution of (3.16) with minimal period $p_{n}:=2T/n$. Then oddness symmetry (3.18) implies

[TABLE]

Substitution into (3.16) implies that $x_{n}(t)$ also solves (1.1). See Figure 3.4a,b for illustrations of the cases $n=1,3$. Note how $p_{n}\searrow 0$ implies unbounded amplitudes $\sqrt{2}<A_{n}\nearrow\infty$, for $n\rightarrow\infty$. Thus we obtain an unbounded sequence of more and more rapidly oscillating periodic solutions to the delayed Duffing DDE (1.1), with minimal periods $2T,2T/3,2T/5,\dots$, which are also periodic with (non-minimal) period $2T$. This proves our claim (1.2), for odd $n$.

By (3.9), the exact periodic solutions $x_{n}(t)$ are expressed as Jacobi elliptic functions

[TABLE]

where (3.10) becomes

[TABLE]

As we have mentioned in subsection 3.1, the positivity and symmetry condition $H>0$ becomes equivalent to $A>\sqrt{2}$.

Figure 3.5 schematically illustrate the lift from the non-delayed Duffing ODE (3.1) to the delayed Duffing DDE (1.1), for both even and odd $n$. This lift will be used in the next section to numerically determine the amplitudes $A_{n}$ of the lifted, rapidly oscillating periodic solutions of the DDE (1.1).

4 Amplitudes

We sketch two practical approaches to determine the amplitudes $A_{n}$ of the rapidly oscillating periodic solutions $x_{n}(t)$ in the delayed Duffing equation (1.1). One approach is essentially numerical; the other approach is analytic, based on an exact series expansion at $n=\infty$ and at infinite amplitude.

The amplitudes $A_{n}$ of the lifted solutions $x_{n}(t)$ arise from the closed curves $H>0$ in the non-delayed Duffing ODE (3.1) with specific values

[TABLE]

of their minimal period. See (1.2) and section 3 for details.

Substitution of (4.1) into the explicit elliptic integral (3.11) provides the implicit equation

[TABLE]

for $A_{n}$, given $T$ and $n$. Here the functions $m(A_{n})$ and $\omega(A_{n})$ are specified in (3.10); we have suppressed explicit dependence on the parity of $n$ in this abbreviated notation.

For high precision numerical solutions $A_{n}$ of (4.2) we rely on the Python-based Newton solver fsolve. The Newton-method requires initial approximations for the desired solution $A_{n}$; for initial guesses we use the formal expansions in [DaShaRa17], Eq. (4). The complete elliptic integral $K(m)$ in (4.2) is evaluated using the Python-based quadrature quad. The integration is performed using a Clenshaw-Curtis method which uses Chebyshev moments. For $T=3$, for example, the reference amplitudes $A_{n}$ corresponding to the red marked points in Figure 3.2 are found to be $A_{1}=1.74566491$…, $A_{2}=2.16089536$…, $A_{3}=3.90053028$…, and $A_{4}=4.79499435$….

Note that time delays $T$ and $\bar{T}$ share the same reference amplitudes, if the relation $T/n=\bar{T}/\bar{n}$ holds. Here $n$ and $\bar{n}$ are required to be both odd, or both even. For example the amplitude $A_{n}=A_{n}(T)$, for $n=1$ and $T=0.1$, coincides with the amplitude $A_{\bar{n}}(\bar{T})$, for $\bar{n}=3$ and $\bar{T}=0.3$.

Our second approach is analytic in nature. We start from (4.2) with an exact Taylor expansion of $p(A):=4K(m(A))/\omega(A)$, at $A=\infty$, with respect to $1/A$. For even $n$ the functions $m(A)$ and $\omega(A)$ have been specified in (3.14). Up to errors of order 13 in $1/A$ we obtain

[TABLE]

Here $\gamma:=\Gamma(1/4)^{2}$ denotes the square of the Euler Gamma-function, evaluated at 1/4. Note $p=0$ at $A=\infty$. Inverting the above series provides an expansion of the inverse function $A(p)$. Specifically, the Taylor expansion of $A$ as a function of $1/p$ at $p=0$, up to errors of order 11 in $1/p$ reads

[TABLE]

Inserting $p=2T/n$ readily provides Taylor expansions of $A$ with respect to $n$, in the limit of large $n\rightarrow\infty$ and for any fixed delay $T>0$. Alternatively, of course, we may consider $n$ fixed and read (4.4) as an expansion with respect to small delays $T>0$, or with respect to any small combination of $T/n$.

For odd $n$, the analogous expansions have to be based on the functions $m(A)$ and $\omega(A)$ specified in (3.21). With the same notation as above we obtain

[TABLE]

[TABLE]

Comparing the even and odd cases, we observe how their sign patterns are related by the complex linear transformation $p\mapsto\mathrm{i}p,\ A\mapsto\mathrm{i}A$. This is in agreement with a scaling of the Duffing ODE.

We emphasize that all Taylor expansions (4.3)–(4.6) are convergent and hence can be performed up to any order. Worries like secular terms and other nuisances ubiquitous in formal asymptotics, disappear. In summary, analytic expansions work best for small $T/n$, e.g. for large $n$, where numerical methods face increasing difficulties. The numerical approach, on the other hand, is the method of choice for larger $T/n$, e.g. for small $n$.

5 Stability

In this section we summarize results from [Fie&al19] on local asymptotic stability and instability of the rapidly oscillating periodic solutions $x_{n}(t),\ n=1,2,3,\dots$, of the delayed Duffing DDE (1.1), as constructed in section 3. We recall how the ODE solutions $x_{n}$ of (3.1) with positive energy $H$ are uniquely determined by their minimal periods $p_{n}=2T/n$, where $T>0$ denotes the delay in (1.1); see (1.2) and (3.9)–(3.11).

To be precise we recall that a periodic reference orbit $x_{*}$ is called stable limit cycle, or also locally asymptotically stable, if any other solution $x(t)$, which starts sufficiently nearby, remains near the set $x_{*}$ and converges to that set, for $t\rightarrow\infty$. A sufficient (but not necessary) condition for local asymptotic stability is linear asymptotic stability. In other words, all Floquet (alias Lyapunov) exponents $\eta$ of the periodic orbit $x_{*}$ possess strictly negative real part (except for the algebraically simple trivial exponent $\eta=0$). We speak of linear instability, in contrast, if $x_{*}$ possesses any Floquet (alias Lyapunov) exponent with strictly positive real part. Deeper results on unstable manifolds then imply nonlinear instability. In fact, there exists a solution $x(t)$ which is defined for all $t\leq 0$ and converges to $x_{*}$ in backwards time $t\rightarrow-\infty$.

The stability results of [Fie&al19] specialize to our present context as follows.

Theorem 5.1**.**

Let $n$ be odd and assume

[TABLE]

Moreover assume that $n\geq n_{0}(T)$ is chosen large enough.

Then the periodic orbit $x_{n}$ of the delayed Duffing equation (3.1) is asymptotically stable, both linearly and locally.

Theorem 5.2**.**

Let $n$ be even, $T>0$, and assume $n\geq n_{0}(T)$ is chosen large enough.

Then the periodic orbit $x_{n}$ of the delayed Duffing equation (3.1) is linearly and nonlinearly unstable.

For the leading Floquet exponent $\eta$, i.e. the nontrivial exponent with real part closest to zero, the precise asymptotics

[TABLE]

has been derived, for even and odd $n\rightarrow\infty$.

Towards the stability boundary $T^{2}=\tfrac{3}{2}\pi^{2}$ of Theorem 5.1, the periodic orbits $x_{n}$ with odd $n$ lose stability, and undergo a torus bifurcation of Neimark-Sacker-Sell type. In particular, rational rotation numbers on the bifurcating torus will indicate periodic orbits of the delayed Duffing DDE (3.1) which are not lifts of the ODE Duffing orbits $x_{n}$ studied in the present paper.

We caution the reader that Floquet theory for delay differential equations is not an entirely trivial matter. Therefore we only illustrate our stability results in the next section, numerically. For detailed mathematical proofs we have to refer to [Fie&al19].

6 Discussion

Figure 6.1 plots two solutions of the delayed Duffing equation (1.1) with delay $T=0.5$: a numerical solution $x(t)$ (blue), and the lifted exact solution $x_{1}(t)$ (red) specified in (3.9). The minimal period $p_{1}$ of $x_{1}(t)$ coincides with $2T$; see (1.2). Figure 6.1 contains the time history (a) and the phase plane (b). The green curve denotes the initial history function

[TABLE]

for $-T<t<0$ and with initial amplitude $A=4.3$. The values of $m$ and $\omega$ are obtained from (3.10) with $n=1$. Note how $x(t)$ is a solution of the non-delayed Duffing ODE (3.1) with minimal period $p=p(A)=1.7972608\dots$. However, the initial history function $x(t)$ is not a solution of the delayed Duffing DDE (1.1), because $T=0.5$ is not an integer multiple of the larger ODE period $p=1.7972608\dots$ . Therefore the simulated solution $x(t)$ of the delayed Duffing DDE (3.1) (blue), is not periodic.

Instead, the simulated solution (blue), with initial amplitude $A=4.3$, approaches the exact periodic solution $x_{1}(t)$ (red) of minimal period $2T$ and with amplitude $A_{1}=7.5139958\dots$ . Indeed, the black curve indicates the history function, for $100-T<t<100$, of the final state of the blue solution $x(t)$ at $t=100$. The stability result of Theorem 5.1 only asserts local convergence to $x_{n}$ for large odd $n$, but not for $n=1$. The convergence to $x_{1}$ indicates how that stability result might actually extend, all the way, down to the smallest possible choice $n=1$. Moreover, “local” attraction to $x_{1}$ holds sway over quite a distance, down to an initial amplitude $A=4.3$ significantly smaller than the asymptotic amplitude $A_{1}=7.5139958\dots$ of $x_{1}$.

Figure 6.2 compares two lifted exact periodic solutions, $x_{1}(t)$ (red) and $x_{2}(t)$ (teal). Two numerical solutions of the delayed Duffing DDE (1.1) for $T=0.5$ are included, which arise from the two initial history functions (6.1) with initial amplitudes $A=1.42$ (blue) and $A=14.77$ (violet), respectively. The reference amplitudes corresponding to the exact $n=1$ (red) and $n=2$ (teal) periodic solutions (3.9) are $A_{1}=7.5139958\dots$ and $A_{2}=14.7834172\dots$, respectively. Figure 6.2 indicates how both simulated solutions (blue and violet) approach the exact stable limit cycle $x_{1}$ (red); see Theorem 5.1. Also note how the simulated solution with initial condition $A=14.77$ (violet) starts very close to the exact, but linearly unstable, periodic solution (teal) of $A_{2}=14.7834172\dots$, but eventually diverges as time $t$ increases. See Theorem 5.2. This indicates the presence of a heteroclinic orbit $x(t)$, from $x_{2}$ to $x_{1}$, which is defined for all positive and negative times $t$ and converges to $x_{2}$, for decreasing $t\searrow-\infty$, and to $x_{1}$, for increasing $t\nearrow+\infty$.

Our periodicity Ansatz requires half minimal periods $p/2=T/n$ to be integer fractions $n=1,2,3,\dots$ of the delay $T$. Of course we have to caution the reader that there may be many periodic solutions of the DDE (1.1) which are not captured by this Ansatz.

7 Conclusion

In this work we showed how the Duffing equation (1.1) with time delay $T$ possesses an unbounded sequence of infinitely many rapidly oscillating periodic solutions $x_{n}(t),\ n=1,2,3,\dots$ .

Each solution $x_{n}$ arises from a periodic solution $x_{n}(t)$ of the non-delayed classical Duffing equation (3.1) with minimal period $p_{n}=2T/n$. In particular, the classical non-delayed Duffing oscillator provides an unbounded sequence of exact periodic solutions of the delayed Duffing equation. Based on the Hamiltonian energy of the classical Duffing equation, and standard Jacobi elliptic integrals, we have also derived high-precision reference amplitudes of these periodic solutions $x_{n}$.

For delays $T$ such that $0<T^{2}<\tfrac{3}{2}\pi^{2}$, and for odd $n$ large enough, the solutions $x_{n}$ are locally asymptotically stable limit cycles. For large even $n$, in contrast, the solutions $x_{n}$ are linearly and nonlinearly unstable.

We have illustrated our results with numerical simulations, for low $n=1,2$.