Analysis of a remarkable singularity in a nonlinear DDE

Nonlinear Dynamics, Jul 2017

We investigate the dynamics of the nonlinear DDE (delay-differential equation) \(\frac{\hbox {d}^2x}{\hbox {d}t^2}(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 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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Analysis of a remarkable singularity in a nonlinear DDE

Analysis of a remarkable singularity in a nonlinear DDE Matthew Davidow 0 1 2 B. Shayak 0 1 2 Richard H. Rand 0 1 2 0 R. H. Rand ( 1 B. Shayak Sibley School of Mechanical and Aerospace Engineering, Cornell University , Ithaca, NY , USA 2 M. Davidow Center for Applied Mathematics, Cornell University , Ithaca, NY , USA We investigate the dynamics of the nonlinear DDE (delay-differential equation) 2 ddt2x (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 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. Delay-differential equation; Harmonic balance; Melnikov's integral 1 Introduction In this work, we investigate the dynamics of the nonlinear DDE (delay-differential equation) dd2t x2 + x (t − T ) + x 3 = 0 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. 2 Harmonic balance x (t ) = A cos ωt 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 + 43 A2 = 0 ( 1 ) ( 2 ) ( 3 ) Melnikov The first of these gives ωT = nπ for n 1, 2, 3, . . ., whereupon the second gives 2 A = √3 n2π 2 T 2 ± 1, n = 1, 2, 3, . . . 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. 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. 3 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: dx ∂ H dt = ∂ y , Note that Eq. ( 5 ) possesses the first integral H (x , y) = constant, since dH/dt = Hx x˙ + Hy y˙ = 0. Now we add a perturbation to the conservative system ( 5 ): dx ∂ H dt = ∂ y + g1, dy ∂ H dt = − ∂ x + g2 where g1 and g2 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: (g1 y˙ − g2x˙) dt = 0 where represents the unperturbed closed curve H (x , y) = constant and where x˙ and y˙ 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: ( 5 ) x˙ = y y˙ = −x − x 3 + (x − x (t − T )) ( 6 ) (7) (8) (9) H (x , y) = 21 y2 + 21 x 2 + 41 x 4 and with perturbations g1 = 0 and g2 = x − x (t − T ) This choice for g2 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 g2 satisfies this requirement. Thus, in our case the Melnikov integral condition (7) becomes 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 0 P where P is the period of the motion around in the unperturbed system, and where we have used the fact that: −x x˙(t )dt = x (t )2 P 2 0 = x ( P)2 − x (0)2 2 = 0 Note that x ( P) = x (0) because x (t ) is periodic with period P. Here x (t ) is the solution to Eqs. ( 5 ) with Hamiltonian (10) which turns out to be a Jacobian elliptic cn function, which may be written as x = a1cn(a2t, k), where the parameters a1,a2 and k are related as follows (see section 2.2 in [2]): a22 = a12 + 1, k2 a12 = 2(1 + a1 ) Thus, our Melnikov integral condition (12) simplifies to: 0 P We are integrating over one full period, and thus, the second term will integrate to 0. The first term, sin2(π a2t /(2K )), is always positive and thus integrates to 0 only if the coefficient sin(π a2T /(2K )) is 0, i.e., Eq. (20) becomes: sin(π a2T /(2K )) = 0 We are interested in the relationship between the amplitude a1 and the delay T . The above gives an implicit relationship between a1 and T since a22 = 1 + a12 and K is also determined by a1 (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 a1 1. Then, from Eqs. (14) we can approximate a2 ≈ a1, k2 = a12/(2a22) ≈ 1/2. (10) (11) (12) (13) (14) 0 0 P P (16) (17) (18) (19) (21) where K = K (1/2) ≈ 1.854, giving the result: a1 ≈ 3.71 n/ T . (22) (23) 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 a1 is chosen. By varying a1, 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. 4 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: dd2t x2 + α ddxt + x (t − T ) + x 3 = 0 (25) a1 ≈ (2π/ √3)n/ T ≈ 3.63 n/ T . 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: −ω2 + cos ωT = 0 αω − sin ωT = 0 Squaring and adding (26) and (27) and using (26) again yields the critical delay for a Hopf bifurcation: Tcrit = √2 arccos −α + 2 √α4 + 4 2 −α2 + A plot of Tcrit 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: sin ωT = αω and − ω2 + cos ωT + 43 A2 = 0 (29) 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 which imply that at the nth intersection point 1 ω = T βn αn = T cos βn (32) (33) 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. 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 sin ωT = αω T cos ωT = α (30) (31) α 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.1 Eigenvalues 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 5 Conclusions We have investigated the occurrence of limit cycles in the delay-differential equation: (34) d2 x dt 2 + x (t − T ) + x 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 (, 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. 1. Bhatt , S.J. , Hsu , C.S.: Stability criteria for second-order dynamical systems with time lag . J. Appl. Mech . 33 ( 1 ), 113 - 118 ( 1966 ) 2. Rand , R.H. : Lecture Notes in Nonlinear Vibrations. Published on-line by The Internet -First University Press. http:// ( 2012 ) 3. Byrd , P.F. , Friedman , M.D.: Handbook of Elliptic Integrals for Engineers and Scientists, 2nd edn . Springer, Berlin ( 1971 ) 4. Stepan , G.: Great delay in a predator-prey model . Nonlinear Anal. Theory Methods Appl . 10 , 913 - 929 ( 1986 ) 5. Kalmar-Nagy , T. , Stepan , G. , Moon , F.C. : Subcritical Hopf bifurcation in the delay equation model for machine tool vibrations . Nonlinear Dyn . 26 , 121 - 142 ( 2001 ) 6. Zhang , L. , Stepan , G.: Exact stability chart of an elastic beam subjected to delayed feedback . J. Sound Vib . 367 , 219 - 232 ( 2016 )

This is a preview of a remote PDF:

Matthew Davidow, B. Shayak, Richard H. Rand. Analysis of a remarkable singularity in a nonlinear DDE, Nonlinear Dynamics, 2017, 1-7, DOI: 10.1007/s11071-017-3663-2