Analysis of a remarkable singularity in a nonlinear DDE

We investigate the dynamics of the nonlinear DDE (delay-differential equation) d2xdt2(t)+x(t-T)+x(t)3=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\hbox {d}^2x}{\hbox {d}t^2}(t)+x(t-T)+x(t)^3=0$$\end{document}, where T is the delay. For T=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=0$$\end{document}, this system is conservative and exhibits no limit cycles. For T>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0$$\end{document}, no matter how small T is, an infinite number of limit cycles exist, their amplitudes going to infinity in the limit as T approaches zero. We investigate this situation in three ways: (1) harmonic balance, (2) Melnikov’s integral, and (3) adding damping to regularize the singularity.

where T is the delay. Using Pontryagin's principle, Bhatt and Hsu [1] showed that the origin in this equation is linearly unstable for all values of T > 0. For T = 0 however, the origin is obviously Lyapunov stable. Thus, a stability change occurs as T changes from zero to any positive value, no matter how small. Associated with this change in stability is a remarkable bifurcation in which an infinite number of limit cycles exist for positive values of T in the neighborhood of T = 0, their amplitudes going to infinity in the limit as T approaches zero. We investigate this situation in three ways: 1. Harmonic balance, 2. Melnikov's integral, 3. Adding damping to regularize the singularity.

Harmonic balance
We seek an approximate solution to Eq. (1) in the form: Substituting Eq. (2) in Eq. (1), simplifying the trig, and equating to zero the coefficients of sin ωt and cos ωt, respectively, we obtain sin ωT = 0 and − ω 2 + cos ωT + 3 4 A 2 = 0 (3) Table 1 Limit cycle amplitudes for values of n in Eq. (4); harmonic balance and Melnikov's method for T = 0.3 Numerical integration of Eq. (1) using DDE23 in MATLAB shows limit cycles with amplitudes 12.31 and 33.56, which correspond to the approximate values 12.14 and 36.29 in Table 1. See Fig. 1. Examination of numerical results shows that limit cycles corresponding to n odd in Table 1 are found to be stable, whereas those corresponding to n even are unstable. In fact, initial condition (x, x ) = (26.681, 0) for t ≤ 0 leads to periodic motion with amplitude 12.31, while initial condition (x, x ) = (26.682, 0) for t ≤ 0 leads to periodic motion with amplitude 33.56, leading to the conclusion that there is an unstable periodic motion with amplitude approximately equal to 26.68, presumably corresponding to amplitude value 24.21 in Table 1.

Melnikov's integral
We begin by generalizing the discussion to a wider class of systems, returning to Eq. (1) later. We consider a conservative (Hamiltonian) system of the form: (1) using DDE23 in MAT-LAB. The initial condition for the upper figure is (x, x ) = (26.681, 0) for t ≤ 0 and leads to a periodic motion with amplitude 12.31, while the initial condition for the lower figure is (x, x ) = (26.682, 0) for t ≤ 0 and leads to periodic motion with amplitude 33.56, leading to the conclusion that there is an unstable periodic motion with amplitude approximately equal to 26.68, presumably corresponding to amplitude value 24.21 in Table 1 Note that Eq. (5) possesses the first integral H (x, y) = constant, since dH/dt = H xẋ + H yẏ = 0. Now we add a perturbation to the conservative system (5): where g 1 and g 2 are given functions of x and y.
For the system (6), the condition for one of the closed curves H (x, y) = constant to be preserved under the perturbation (6) turns out to be given by the vanishing of the following Melnikov integral: where represents the unperturbed closed curve H (x, y) = constant and whereẋ andẏ refer to time histories around in the unperturbed system. The derivation uses Green's theorem of the plane, and the result is approximate (see section 3.3 in [2]).
To apply the foregoing setup to Eq. (1), we write (1) in the following form: where x written without an argument stands for x(t).
That is, we consider Eq. (1) to be a perturbed Hamiltonian system (6) with Hamiltonian and with perturbations This choice for g 2 is based on the requirement that the perturbation be small compared to the other terms in Eq. (6). For small values of delay T , this choice of g 2 satisfies this requirement.
Thus, in our case the Melnikov integral condition where P is the period of the motion around in the unperturbed system, and where we have used the fact that: is the solution to Eqs. (5) with Hamiltonian (10) which turns out to be a Jacobian elliptic cn function, which may be written as where the parameters a 1 ,a 2 and k are related as follows (see section 2.2 in [2]): Thus, our Melnikov integral condition (12) simplifies to: where P = 4K (k)/a 2 , where K (k) is a complete elliptic integral of the first kind. In order to obtain an analytical approximation for this integral, we use the following expansions for the elliptic functions sn, cn and dn [3]: We take the first term in each of the expansions (16), (17), (18), whereupon the Melnikov integral condition (15) becomes: We are integrating over one full period, and thus, the second term will integrate to 0. The first term, sin 2 (πa 2 t/(2K )), is always positive and thus integrates to 0 only if the coefficient sin(πa 2 T /(2K )) is 0, i.e., Eq. (20) becomes: We are interested in the relationship between the amplitude a 1 and the delay T . The above gives an implicit relationship between a 1 and T since a 2 2 = 1 + a 2 1 and K is also determined by a 1 (through an elliptic integral). To make a much simpler explicit relationship, we will use the fact that we are in the regime of T 1, and in this parameter range we have empirically found that a 1 1. Then, from Eqs. (14) we can approximate a 2 ≈ a 1 , k 2 = a 2 1 /(2a 2 2 ) ≈ 1/2. This gives us where K = K (1/2) ≈ 1.854, giving the result: This result may be compared to the harmonic balance result of Eq. (4), which is These approximate analytical results may be checked by evaluating the Melnikov integral (15) numerically. For a fixed value of delay T , a value for the second integral in (15) may be computed in MATLAB once the amplitude a 1 is chosen. By varying a 1 , we obtained two plots, one with delay T = 0.05, and the other with T = 0.2, see Figs. 2 and 3. If we look at the zeros of both plots, it looks like they occur at integer multiplies of a certain amplitude. This agrees with the harmonic balance result of Eq. (4). Figure 4 compares the results obtained by numerically evaluating the Melnikov integral with those obtained by numerical integration of Eq. (1) using DDE23, and with those of harmonic balance for the first zero (corresponding to n = 1) for different values of delay.

Adding damping to regularize the singularity
We have seen that in the case of Eq. (1), even infinitesimal delay gives rise to effective negative damping and growing oscillations. Accordingly, we expect that if damping is added to Eq. (1), as in the case of the following DDE: then if α is held fixed and delay T is increased from 0, there will be a point at which the equilibrium at the Fig. 5 A plot of T crit from (28) as a function of α. Note that at small α the function is approximately given by T = α. Limit cycles exist above this line, but not below it origin will make a transition from stable to unstable. Supposing that such a transition is a Hopf bifurcation, we linearize Eq. (25) by dropping the x 3 term, and then set x = exp iωt, giving the real and imaginary parts: Squaring and adding (26) and (27) and using (26) again yields the critical delay for a Hopf bifurcation: A plot of T crit as a function of α can be seen in Fig. 5.
In addition to this Hopf bifurcation, it turns out that additional limit cycles can occur in this system by being born in a fold (also known as a saddle node of cycles). In order to see this, we again use the method of harmonic balance. Assuming an approximate solution of the form x(t) = A cos ωt, we substitute into Eq. (25), simplify the trig, and equate to zero the coefficients of sin ωt and cos ωt, respectively, giving: Suppose the value of T is fixed and α is started from a high value. The first equation of (29) can be viewed in terms of two functions of the variable ω; one is the straight line αω and the second is the sinusoid sin ωT . See Fig. 6.
If α > T , then the two curves intersect only at the trivial point ω = 0 and there is no limit cycle. This corresponds to curve a in Fig. 6. Now consider the  Fig. 6), and upon being slightly lower still, the curves develop a non-trivial intersection (curve c in Fig. 6). This means that a sinusoidal response with frequency ω is a possible state of the system. Once the frequency is specified, the amplitude of the motions gets determined by the second equation in (29). We thus have a limit cycle with amplitude A and frequency ω.
As we lower α still further, the non-trivial intersection point between the straight line and the sinusoid will shift rightwards. Ultimately, the two graphs will touch at a second point (curve d in Fig. 6) and a new pair of limit cycles will get born there, since further lowering α will turn the tangency into a pair of intersections, one Fig. 7 Plot of the bifurcation curves using the tangency condition (32), (33). The number of limit cycles in the various regions is shown in the first few cases (0,1,3,5) Fig. 8 Frequency of the limit cycles which appear as α is varied keeping the delay T fixed. Note the exponential decrease of the alpha-axis scale as one moves toward the right. Note also that the Hopf bifurcation occurs close to α = T = 0.1, cf. Fig. 5. Thereafter, as the damping decreases, the successive limit cycles keep emerging in folds. At zero damping, an infinitude of LCs will appear, corresponding to the singularity discussed in this article corresponding to a stable limit cycle and the other to an unstable one. Thus, we can say that a saddle-node bifurcation of cycles is occurring there. The points of tangency are given by the relation which imply that at the nth intersection point where β n 's are the solutions of tan x = x. A bifurcation diagram using these relations is shown in Fig. 7.
In this figure, the line on the left is the Hopf bifurcation. The uppermost line on the right is the first saddle-node bifurcation of cycles, while the next lines going down correspond to the subsequent saddle-node bifurcations. These predictions are in agreement with numerical simulation results.
The frequency of limit cycles is shown in Fig. 8 as the damping parameter alpha is decreased towards zero for fixed delay T. The smallest frequency omega at a fixed alpha corresponds to a limit cycle that is born in a Hopf bifurcation, while the higher frequency limit cycles are born in fold bifurcations.
A comparison between theory and simulation performed on Matlab using the routine DDE23 is presented in Table 2. This table shows the birth of a limit cycle (LC) as α is lowered to a value smaller than T and its branching out into multiple cycles as α is lowered further. It also shows the eigenvalues characterizing the stability of the origin (NRP negative real part, DNE does not exist.) In some regions of the space, 3 and 5 LCs are seen in the harmonic balance. In this case, the first one and then the subsequent alternate ones are observed numerically with the intermediate amplitudes acting as separatrix. It is also seen that the LC is born when the real part of the eigenvalue at the origin crosses from negative to positive, the archetypal signature of a Hopf bifurcation

Conclusions
We have investigated the occurrence of limit cycles in the delay-differential equation: Besides numerical integration, we used harmonic balance and Melnikov's integral to obtain analytic approximations for this system. All of these approaches support the conclusion that this system exhibits an infinite number of limit cycles for positive values of T in the neighborhood of T = 0, their amplitudes going to infinity in the limit as T approaches zero.
This work may be extended by applying center manifold reduction [2,4,5] and Hopf bifurcation calculation, which would yield not only approximation of the limit cycles but also their stability.
We thank Professor Gabor Stepan for noting that the phenomenon discussed in this paper may be related to delay equations with two time delays. See [6] where an example is given in which stability depends on whether the ratio of delays is rational or irrational.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.