Stability criterion for time-delay 3-dimension damped Mathieu equation

  1. Yusry O. El - Dib 1

Department of Mathematics, Faculty of Education, Ain Shams University, Roxy, Egypt

  1. Corresponding author email

Associate Editor: Dr. Noor Danish Ahrar Mundari
Science and Engineering Applications 2016, 1, 76–88. doi:10.26705/SAEA.2016.1.21.76-88
Received 03 Nov 2016, Accepted 21 Nov 2016, Published 21 Nov 2016


In the present work, the investigation of the dynamics of three dimensions Mathieu equations containing periodic terms as well as the delayed parameters. Mathieu equations included the influence of damping terms and the coupled system involves both delayed and non-delay terms. The system is proposed as an extension of delayed two coupled Mathieu-type equations to higher dimensions, with emphasis on how resonance between the internal frequencies leads to a loss of stability. The method of multiple scales is used to examine the islands of stability near the resonance cases. The transition curves are analyzed using the method of harmonic balance, and we find we can use this method to easily predict the `resonance curves' from which bands of instability emanate. We note that the delayed higher dimension of the parametric excitation has a great interest and application to the design of nuclear accelerators.

Keywords: Three dimensions Mathieu equations; resonance curves; parametric excitation.


The stability of delay differential equations are used to variety of applications, including biology [1], population dynamics [2], control systems [3], manufacturing [4,5]. A fundamental and critical task is how finding the regions for which a system of delay differential equations is stable. Several methods have been proposed in the literature for studying the stability of delay differential equations with constant coefficients; however, the stability analysis of delay differential equations having time-depended coefficients, specially, in periodic configuration and time delayed is nontrivial. The closed-form stability chart for the delayed linear one-dimension Mathieu equation, defined as


is determined by T. Insperger and G. Stépán [6], and by Yunna Wu and Xu Xu [7], where the time delay is equal to the principal period 2 π, so it can also be viewed as a special resonant case of systems with optional principal period. where δ,α, β and τ are parameters: δ is the frequency-squared of the simple harmonic oscillator, α is the amplitude of the parametric resonance, β is the amplitude of delay, and τ is the delay period. Each term is caled by the small parameter ε .

A system of this type has been investigated by Morrison and Rand [8]. It was shown that the region of instability associated with 2:1 sub-harmonic resonance can be eliminated by choosing the delay T long enough.

Stability of damped Mathieu’s equation with time-periodic coefficients and time-delayed Mathieu’s equation has been considered by T. Insperger and G. Stépán [9] and by N. K. Garg et al [10]. Z. Ahsan et al [11] considered damped Mathieu equation with two different points delayed and defined by the following equation:


T. M. Morrison and R. H. Rand [8] investigate the dynamics of a delayed nonlinear one-dimension Mathieu equation:


Three different phenomena are combined in this system: 2:1 parametric resonance, cubic nonlinearity, and delay. The method of averaging at small ε is used to obtain a slow flow that is analyzed for stability and bifurcations.

The dynamics of a type of particle accelerator called a synchrotron, in which particles are made to move in nearly circular orbits of large radius. The stability of the transverse motion of such a rotating particle may be modeled as being governed by Mathieu’s equation. For a train of two such particles the equations of motion are coupled due to plasma interactions and resistive wall coupling effects [12]. A. Bernstein & R.H. Rand [13] has address investigation of coupled Parametrically Driven Modes in Synchrotron Dynamics. They studied a system consisting of a train of two particles which is modeled as two coupled nonlinear Mathieu equations with delay coupling. Recently, in (2016) the delay-coupled Mathieu equations in synchrotron dynamics has been addressed by A. Bernstein and R.H. Rand [14]. They investigate the dynamics of the model having two delay-coupled Mathieu equations having a single point delay:


They interested in the form of the above equations comes from an application in the design of nuclear accelerators. They used the two-time scales method [15] to study the dynamics for their coupled Mathieu equations. Only, one resonance case has been studied. In which, they detune off of the 2:1 sub-harmonic resonance by setting δ=1/4+ε δ1+o(ε2).The stability behavior has been determined from the equilibrium point at the origin.

In this work, we consider a generalized form of the delay-two coupled Mathieu equations considered by A. Bernstein and R.H. Rand [14]. In order to get a wide applications and for more generalization, the mathematical model has been extended in the dimension to become of 3-dimensional of the delayed-damping Mathieu equation.

The underlying mathematical problem of 3-dimension Mathieu equation, with weak viscous damping coefficients, and a single point delay is given below:


where ε is a non-zero small parameter, Ω is a frequency of the external excitation, t is an independent parameter, μjk and bjk scales the influence of the viscous damping state and the influence of delay state. Equation (5) represents a 3-coupled delayed Mathieu’s equation, which can be introduced in the vector extension of a standard Mathieu equation as shown below:


where A , B, Q and M are constant 3×3 matrices, with real entries. The vector

X(t)=(x1(t),x2(t),x3(t))T having three dependent variables on t, x1(t), x2(t) and x3(t).

2. Line of solution in the unperturbed pattern:

The 3-dimension system has three `natural' frequencies when the time dependent terms are switched off. The vector secondorder differential equation will reduce to


Let the solution of (7) having the exponential form so that


where ω is called the eigenvalue represents the frequency of the wave-train solution of (7), which satisfies the following characteristic equation:


where Ak is the sum of all kth-order principal minors of the matrix A and I denoted the identity matrix[16] . This is a cubic equation in . ω2 Therefore the necessary and sufficient condition for stability is that ωj2 must be real and positive. It is easily verified from (9) that this restriction for the stability implies the following conditions:


From elementary algebra, the three distinct roots (ωj2; j=1,2,3) will be real as the discriminant for the cubic polynomial (9) has a positive value [17]. This requires that


According to Cardano's formula [17], the three distinct eigenvalues of (9) are given by


where ηj represents the three cubic roots of unity as


In the light of the exponential solution as assumed in (4) we may write


where ; λj;j =1,2,3 are arbitrary constants of integration and λj* is the complex conjugate of λj. Here x̃ , ỹ and z̃ are the unknown elements of an eigenvector, which to be found from equation (7). The latter is a system of linear equations with zero on the right side. It will have nonzero solutions if and only if the determinant of the system is equal to zero. These solutions can be found and put in the vector in the following vector form:


where the constant vector Rj is given by


3. The analysis through the method of multiple time -scales

Although the solutions and properties of a one dimension Mathieu equation is well known [15], there is no general analytical solutions available for the 3-dimension Mathieu equation (6). Therefore we shall discuss the stability of (6) using asymptotic expansion treatment. For turn on the effect of the variable coefficient in equation (6), we need to use a perturbation technique. We shall apply the well-known of the multiple scale method [15]. This method enables us to discuss the stability of the problem (6).

On applying the method of multiple scales we may use the scale T0 ,T1 such that Tn = εnt, n= 0,1, 2,... T0 is the variable appropriate to fast variable, and T1 is the slow variable. In addition, the delay time can scaled as τnnτ . The differential operator can now expressed as the partial derivative expansions:


Assuming that ω1222 and ω32 are both real and positive, then the dependent vector variable X(t) can be expanded in the form:


where the vector X0(T0,T1) has been found to be


while the perturbed vector X1(T0,T1) is given by


where c.c represents the complex conjugate of preceding terms. Since λj(T1-ετ) is unknown function which can be expanded about the un-delayed variable. The most general way to express it in terms of T1 is with a Taylor series:


Consequently, (19) becomes


The above equation contains secular vectors that are proportional to the factor (e±iωjT0). Before eliminating these secular vectors we are in need to distinguish between several possible combinations of Ω, ω12 and ω3 . These cases of Ω near ωj or near the combinations (ω1±ω2)/2 ,(ω1±ω3)/2 and (ω2±ω3)/2 are known as resonance cases. The non-resonant case arises when Ω away from ω123 and (ωp±ωq)/2;p,q=1,1,3;p≠q .

3.1 The non-resonance case

The elimination of the source of secular vectors from the vector equation (21), in the non-resonant case, leads to


This is a first-order vector differential equation represent the amplitude equation. In order to transform the vector equation to a scalar case, one can multiply both sides of (22) from the left by the complex vector


and use the following normalized condition:


Then the vector equation (22) will comes in the following scalar form:


where (SjM Rj) , (SjQ Rj) and (SjB Rj) are complex scalar quantities. The complex vector Sj can be split into the real part and the imaginary part. This can be done by help of the properties of inverse matrices [16]. It can easily to write Sj=SjRe+iSjIm the upper super script “Re” and “Im” denote the real and imaginary parts,


The amplitude equation (25) represents a scalar ordinary first-order differential equation with complex coefficient. Solution of this equation may be found in the form


where ∧j is a complex constant of integration. Since the stability occurs for negative real values of the exponential given in (28), then, one can conclude that the system is stable at the non-resonance case whence


This condition reads that the influence of the viscous damping matrix M plays a stabilizing role as well as the amplitude of the delayed matrix B .

3.2 The resonant case of Ω near ωj

Introducing the detuning parameter σ in equation (21) to convert the small divisor term into secular term as follows:


and write


Using (31), the small-divisor term arising from exp[±i(ω-2Ω)T0] in equation (21) can be transformed into a secular term. Then, remove the source of secular terms. At this stage, the solvability condition given by (25) will be modified to become


with its complex conjugate. This is the amplitude equation governed the stability behavior at the resonance case. It’s a first-order deferential equation with complex coefficients and having parametric term associated with the complex conjugate of the variable λj(T1) .The use the complex conjugate of the solvability condition (32) consist a coupled system of the first-order in the variable λj(T1 and its complex conjugate. Solution of this system can be sought in the form


with complex unknown constant γj which can be calculated using the condition of the nontrivial solution to be


where Sj* is the complex conjugate of the vector Sj .The eigenvalues for the characteristic exponent is found from the following characteristic equation:


where Hj is a square matrix given by


Clearly, both the elements on the principal diagonal of the matrix Hj are a complex conjugate for each other as well as the nonprincipal diagonal. Stability criteria of matrices are related to the stability criteria of polynomials given by Routh-Hurwitz [18] . In which the polynomial must have negative real parts. A necessary and sufficient condition for the stability of a square matrices with real (complex conjugate) entries is performed by Michael Y. Li and Liancheng Wang [19]. Thus, if all the eigenvalues, of the above characteristic equation, have negative real parts then the stability arises whence


The first condition reduces to


This is the same condition that arises at the non-resonance case. Therefore, the critical condition for stability is the second one in (37). This condition can be arranged in powers of the detuning parameter σ as


This condition has two zeros, namely, σ1 and σ2 with (σ12) . Thus the instability is found at the resonance case whence the detuning parameter σ lies inside the open interval (σ21). In terms of the frequency Ω , the transition curves separating stable region from unstable one corresponding to


where the constants k1 and k0 are given below:


The region lies between the two curves represent the unstable (resonance) case, which impeded through the stable (the non-resonant) case

3.3 The resonant case of Ω near (ω12)/2

Here, we shall consider the positive sign of ωq;p≠q and p,q = 1,2,3 , while the negative one can be obtained for replacing the sign in final results. We express the nearness of Ω to (ωpq)/2 by introducing the detuning parameter δ such that


Accordingly, we have


At this end, the secular terms appear in equation (22) can be rearranged to introducing the following two solvability conditions:


This is a coupled system with complex entries. Its solution may be sought in the following form:


where γp,q is a complex constants and the characteristic exponent Ξ at this case is given by


The above characteristic equation is a quadratic in Ξ . It has two different complex roots Ξ1 and . Ξ2 In an ordinary differential equation with complex coefficients, the trivial solution is asymptotically stable if and only if all roots of the corresponding characteristic equation have negative real parts. Since the characteristic function is a polynomial, the well-known Routh-Hurwitz criterion [18] can be used in order to determine the negativity of the real parts of the roots Ξ1 and Ξ2 for characteristic equation (46). This criterion imposes the following conditions for the stabilization at the present case:


First condition of the above stability criteria (47) depends on the influence of the amplitude of the parametric force Q and the amplitude of the delay terms B as well as the delayed parameter τ. In addition the matrix M , of the damping coefficients, is included. But the influence of the detuning parameter δ has been included, only, in the second condition. The second condition has been arranged in powers in detuning parameter δ . This condition has two zeros δ1 and δ2 (say, δ12 ). Then one can concluded that the stability behavior is found at the present resonance case whence


where δ1 and δ2 are given by the well-known low of the quadratic polynomial equation. The transition curves which separating resonance region from non-resonance one corresponding to:


Clearly the resonance region lies between the two curves given by (48). The regions outside these curves represent the stable case. Similar results can be obtained for the case of Ω near (ωpq)/2 by changing the sign of ωq in the above analysis.

The real constants fl and gl ; l= 0,1,2 that appear in characteristic equation (45) are given below:


4. Conclusion

From the preceding analysis, it is quite clear that we have studied the linear instability analysis of dynamics of three dimensions Mathieu equations containing periodic terms as well as the delayed parameters has been don. These three equations are formulated through a matrix representation. In order to clarify the problem, we have introduced a brief introduction, in Section 1, that shows a little history on the Mathieu equations and their properties and applications. In Section 2, Line of solution in the unperturbed system in the absence of the small parameters ε. The solution of the vector second-order differential equation leads to an eigen-values governed by a characteristic equation of cubic degree in ω2. The influence of the small parameter ε is introduced in Section 3. The method of two-time scales is used. The stability analyses studied in subsections 3.1, 3.2 and 3.3. A first-order vector differential equation represent the amplitude equation in the non-resonance case as well as in the resonance cases have been imposed. These vector differential equations have been transformed to a scalar differential equations. The presence of the delay parameter τ leads to produce a complex matrix of type 3×3 .The formula of the inverse of complex matrix is used. Stability criteria of matrices are applied and conditions of stability imposed. According to the Floquet theorem, the transition curves that separate the stable from the unstable regions are theoretically derived in both the resonance case of Ω near ωj and the case of Ω near (ωpq)/2,ωq ; p≠q and p,q= 1,2,3.


  1. Bocharov, G.A., Rihan, F.A.: Numerical modelling in biosciencesusing delay differential equations. J.Comput.Appl. Math. (2000), 125(1–2), 183–199)
  2. Kuang, Y.: Delay Differential Equations:With Applications in Population Dynamics. Academic Press, SanDiego (1993)
  3. Richard, J.P.: Time-delay systems: an overview of some recent advances and open problems. Automatica, (2003) 39(10),1667–1694
  4. Balachandran, B.: Nonlinear dynamics of milling processes. Philos. Trans. R. Soc. A, (2001).359(1781), 793–819
  5. Insperger, T., Stépán, G.: Stability of the milling process. Period. Polytech. Mech. Eng., (2000) 44(1), 47–57.
  6. T. Insperger and G. Stépán, Proc. R. Soc. Lond. A (2002) 458, 1989-1998.
  7. Yunna Wu and Xu Xu, THEORETICAL & APPLIED MECHANICS LETTERS (2013) 3, 063007-1- 063007-3.
  8. T. M. Morrison and R. H. Rand, Nonlinear Dynamics (2007) 50:341–352
  9. T. Insperger and G. Stépán, Journal of Dynamic Systems, Measurement, and Control (2003), Vol. 125 166-171.
  10. N. K. Garg, B. P. Mann, N. H. Kim and M. H. Kurdi, Journal of Dynamic Systems, Measurement, and Control(2007), Vol. 129,125-135.
  11. Z. Ahsan, A. Sadath ,T. K. Uchida and C. P. Vyasarayani, Nonlinear Dynamics (2015) 82:1893–1904.
  12. A. Bernstein, & R.H. Rand. Coupled Parametrically Driven Modes in Synchrotron Dynamics. Chapter 8, pp.107-112 in Nonlinear Dynamics, Volume 1: Proceedings of the 33rd IMAC, A Conference and Exposition on Structural Dynamics, G. Kerschen, editor, Springer (2016).
  13. A. Bernstein and R.H. Rand. Coupled Parametrically Driven Modes in Synchrotron Dynamics. Presented at Society for Experimental Mechanics (SEM) IMAC XXXIII Conference on Structural Dynamics February 2-5, (2015), Orlando FL
  14. A. Bernstein and R.H.R, Journal of Applied Nonlinear Dynamics (2016),(to appear)
  15. A. H. Nayfeh and D. T. Mook, Nonlinear Oscillations, Wiely-Interscience, New York, (1979).
  16. K. B. Petersen and M. S. Pedersen, The Matrix Cookbook[ ],15, (2012).
  17. Internet references Wikipedia:Cubic function-Cardano’s method.
  18. S. S. Hasan, Automatic Control Systems, Katson Publisher, Delhi, 2008
  19. Michael Y. Li and Liancheng Wang , J. Mathematical analysis and applications, (1998),225, 249-264.

© 2016 Yusry O. El - Dib; licensee Payam Publishing Pvt. Lt..
This is an Open Access article under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The license is subject to the Science and Engineering Applications terms and conditions: (

Keep In Touch

RSS Feed

Subscribe to our Latest Research Articles RSS Feed.


Follow the Science and Engineering Applications


Twitter: @SAEA

Facebook: Facebook

Back to Article List