Analysis of a remarkable singularity in a nonlinear DDE

In this work we investigate the dynamics of the nonlinear DDE (delay-differential equation) x''(t)+x(t-T)+x(t)^3=0 where T is the delay. For T=0 this system is conservative and exhibits no limit cycles. For T>0, no matter how small, 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.


Introduction
In this work we investigate the dynamics of the nonlinear DDE (delay-differential equation) 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 Liapunov 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.

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 The first of these gives ωT = nπ for n=1,2,3,· · ·, whereupon the second gives where the upper sign refers to n odd, and the lower sign refers to n even. For example, when T = 0.3, Table 1 gives values for amplitudes of limit cycles for given values of n, from eq.(4). 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. Presumably the reason we do not see limit cycles with the other amplitudes listed in Table 1 is that they are unstable. In fact, initial condition (x,x')=(26.681,0) for t ≤ 0 leads to periodic motion with amplitude 12.14, while initial condition (x,x')=(26.682,0) for t ≤ 0 leads to periodic motion with amplitude 36.29, 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: 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. (1) using DDE23 in MATLAB shows that initial condition (x,x')=(26.682,0) for t ≤ 0 leads to periodic motion with amplitude 36.29 while initial condition (x,x')=(26.681,0) for t ≤ 0 leads to periodic motion with amplitude 12.14, 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.
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 g 1 = 0 and Thus in our case the Melnikov integral condition (7) becomes 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: 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]: where G = πz/(2K(k)), q = e −πK (k)/K(k) and K (k) = K( √ 1 − k 2 ). We take the first term in each of the expansions (16),(17),(18), whereupon the Melnikov integral condition (15) becomes: Expanding the cosine term gives P 0 [sin 2 (πa 2 t/(2K))sin(πa 2 T /(2K)) +sin(πa 2 t/(2K))cos(πa 2 t/(2K))cos(πa 2 T /(2K))]dt = 0 (20) 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 sin(πa 1 T /(2K(1/2)) = 0 ⇒ a 1 = 2Kn/T 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). Fig.4 compares the numerical results with those of Harmonic Balance in a plot of the first zero (corresponding to n = 1) for different values of delay. 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 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: A plot of T crit as a function of α can be seen in Fig.5. Figure 5: A plot of T crit from (28) as a function of α. Note that at small α the function is like the identity T = α. Limit cycles exist above this line, but not below it.
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: sin ωT = αω and − ω 2 + cos ωT + 3 4 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 situation as α is lowered. When it becomes equal to T there is a tangency at the origin (curve b in 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 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 n th intersection point where β n 's are the solutions of tan x = x.
A bifurcation diagram using these relations is shown in Fig.7. Figure 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).
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. A comparison between theory and simulation performed on Matlab using the routine DDE23, is presented in Table 1. 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 is also seen that the LC is born when the real part of the eigenvalue crosses from negative to positive.

Conclusions
We have investigated the occurrence of limit cycles in the delay-differential equation: Besides numerical integration, we used three different approximate analytic approaches to study 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.