Analysis of a remarkable singularity in a nonlinear DDE
AAnalysis of a remarkable singularity in a nonlinear DDE
Matthew Davidow ∗ , B. Shayak ∗∗ and Richard H. Rand ∗∗∗ ∗ Center for Applied Mathematics, Cornell University, Ithaca, NY ∗∗ Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY ∗∗∗
Department of Mathematics and Sibley School of Mechanical and Aerospace Engineering,Cornell University, Ithaca, NY
In this work we investigate the dynamics of the nonlinear DDE (delay-differential equation) d xdt + x ( t − T ) + x = 0 (1)where T is the delay. Using Pontryagin’s Principle, Bhatt and Hsu [1] showed that the originin this equation is linearly unstable for all values of T >
0. For T = 0 however, the origin isobviously Liapunov stable. Thus a stability change occurs as T changes from zero to any positivevalue, no matter how small. Associated with this change in stability is a remarkable bifurcationin 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. We seek an approximate solution to eq.(1) in the form: x ( t ) = A cos ωt (2)Substituting eq.(2) in eq.(1), simplifying the trig, and equating to zero the coefficients of sin ωt and cos ωt respectively, we obtainsin ωT = 0 and − ω + cos ωT + 34 A = 0 (3)1 a r X i v : . [ m a t h . D S ] J a n able 1: Limit cycle amplitudes A for values of n in eq.(4), for T = 0 . n A ωT = nπ for n =1,2,3, · · · , whereupon the second gives A = 2 √ (cid:115) n π T ± , n =1,2,3, · · · (4)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 amplitudes12.31 and 33.56, which correspond to the approximate values 12.14 and 36.29 in Table 1. SeeFig.1. Presumably the reason we do not see limit cycles with the other amplitudes listed inTable 1 is that they are unstable. In fact, initial condition (x,x’)=(26.681,0) for t ≤ t ≤ We begin by generalizing the discussion to a wider class of systems, returning to eq.(1) later. Weconsider a conservative (Hamiltonian) system of the form: dxdt = ∂H∂y , dydt = − ∂H∂x (5)Note that eq.(5) possesses the first integral H ( x, y ) = constant, since dH/dt = H x ˙ x + H y ˙ y = 0.Now we add a perturbation to the conservative system (5): dxdt = ∂H∂y + g , dydt = − ∂H∂x + g (6)where g and g are given functions of x and y .igure 1: Numerical integration of eq.(1) using DDE23 in MATLAB shows that initial condition(x,x’)=(26.682,0) for t ≤ t ≤ H ( x, y ) = constant to be preservedunder the perturbation (6) turns out to be given by the vanishing of the following Melnikovintegral: (cid:73) Γ ( g ˙ y − g ˙ x ) dt = 0 (7)where Γ represents the unperturbed closed curve H ( x, y ) = constant and where ˙ x and ˙ y refer totime histories around Γ in the unperturbed system. The derivation uses Green’s Theorem of thePlane, 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:˙ x = y (8)˙ y = − x − x + ( x − x ( t − T )) (9)where x written without an argument stands for x ( t ). That is we consider eq.(1) to be a perturbedHamiltonian system (6) with Hamiltonian H ( x, y ) = 12 y + 12 x + 14 x (10)and with perturbations g = 0 and g = x − x ( t − T ) (11)Thus in our case the Melnikov integral condition (7) becomes (cid:90) P − ( x ( t ) − x ( t − T )) ˙ x ( t ) dt = (cid:90) P x ( t − T )) ˙ x ( t ) dt = 0 (12)here P is the period of the motion around Γ in the unperturbed system, and where we haveused the fact that: (cid:90) P − x ˙ x ( t ) dt = x ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P = x ( P ) − x (0) x ( P ) = x (0) because x ( t ) is periodic with period P . Here x ( t ) is the solution toeqs.(5) with Hamiltonian (10) which turns out to be a Jacobian elliptic cn function, which maybe written as x = a cn( a t, k ) , (13)where the parameters a , a and k are related as follows (see section 2.2 in [2]): a = a + 1 , k = a a ) . (14)Thus our Melnikov integral condition (12) simplifies to: (cid:90) P cn( a ( t − T ) , k ) ddt (cn( a t, k ) dt = (cid:90) P cn( a ( t − T ) , k ) sn( a t, k ) dn( a t, k ) dt = 0 (15)where P = 4 K ( k ) /a , 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 expansionsfor the elliptic functions sn, cn and dn [3]:sn( z, k ) = 2 πkK ∞ (cid:88) n =0 q n +1 / sin((2 n + 1) G )1 − q n +1 (16)cn( z, k ) = 2 πkK ∞ (cid:88) n =0 q n +1 / sin((2 n + 1) G )1 + q n +1 (17)dn( z, k ) = π/ (2 K ) + 2 πK ∞ (cid:88) n =0 q n cos(2 nG )1 + q n (18)where G = πz/ (2 K ( k )), q = e − πK (cid:48) ( k ) /K ( k ) and K (cid:48) ( k ) = K ( √ − k ). We take the first term ineach of the expansions (16),(17),(18), whereupon the Melnikov integral condition (15) becomes: (cid:90) P cos( πa ( t − T ) / (2 K )) sin( πa t ) / (2 K )) dt = 0 (19)Expanding the cosine term gives (cid:90) P [sin ( πa t/ (2 K ))sin( πa T / (2 K ))+sin( πa t/ (2 K ))cos( πa t/ (2 K ))cos( πa T / (2 K ))] dt = 0 (20)We are integrating over one full period, and thus the second term will integrate to 0. Thefirst term, sin ( πa t/ (2 K )), is always positive and thus integrates to 0 only if the coefficientsin( πa T / (2 K )) is 0, i.e. eq.(20) becomes:sin( πa T / (2 K )) = 0 (21)igure 2: Melnikov Integrals at T=0.05We are interested in the relationship between the amplitude a and the delay T . The abovegives an implicit relationship between a and T since a = 1 + a and K is also determined by a (through an elliptic integral). To make a much simpler explicit relationship we will use thefact that we are in the regime of T <<
1, and in this parameter range we have empirically foundthat a >>
1. Then from eqs.(14) we can approximate a ≈ a , k = a / (2 a ) ≈ / πa T / (2 K (1 / ⇒ a = 2 Kn/T (22)where K = K (1 / ≈ . a ≈ . n/T. (23)This result may be compared to the Harmonic Balance result of eq.(4), which is a ≈ (2 π/ √ n/T ≈ . n/T. (24)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 computedin MATLAB once the amplitude a is chosen. By varying a we obtained two plots, one withdelay T = 0 .
05, and the other with T = 0 .
2, see Figs.2 and 3. If we look at the zeros of bothplots, it looks like they occur at integer multiplies of a certain amplitude. This agrees with theHarmonic Balance result of eq.(4). Fig.4 compares the numerical results with those of HarmonicBalance in a plot of the first zero (corresponding to n = 1) for different values of delay.igure 3: Melnikov Integrals at T=0.2Figure 4: Comparison of limit cycle amplitudes obtained numerically versus analytically. Nu-merical values correspond to the first zero of the Melnikov integral (12), while analytical valuesare those obtained by Harmonic Balance, eq.(24), for n =1. Adding damping to regularize the singularity