#### 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
(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.
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:// ecommons.library.cornell.edu/handle/1813/28989 ( 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 )