CUBO A Mathematical Journal Vol.14, No¯ 03, (103–113). October 2012 An Elementary Study of a Class of Dynamic Systems with Two Time Delays Akio Matsumoto 1 Department of Economics, Chuo University, 742-1, Higashi-Nakano, Hachioji, Tokyo, 192-0393, Japan email: akiom@tamacc.chuo-u.ac.jp Ferenc Szidarovszky Department of Systems and Industrial Engineering, University of Arizona, Tucson, 85721-0020, USA. email: szidar@sie.arizona.edu ABSTRACT An elementary analysis is developed to determine the stability region of a certain class of ordinary differential equations with two delays. Our analysis is based on determining stability switches first where an eigenvalue is pure complex, and then checking the conditions for stability loss or stability gain. In the case of both stability losses and stability gains Hopf bifurcation occurs giving the possibility of the birth of limit cycles. RESUMEN Se realiza un análisis básico para determinar la estabilidad de la región de una cierta clase de ecuaciones diferenciales ordinaras con dos retrasos. Nuestro análisis se basa en la determinación de switches de estabilidad, en primer lugar cuando un autovalor es complejo puro, y luego revisando las condiciones para la pérdida o ganancia de estabilidad. En el caso de ambas pérdidas de estabilidad y ganancias de estabilidad, se obtiene la bifurcación de Hopf dando la posibilidad del nacimiento de ciclos ĺımites. Keywords and Phrases: dynamic systems, time delays, stabiliy analysis. 2010 AMS Mathematics Subject Classification: 34K20, 37C75 1 The authors highly appreciate financial supports from the Japan Society for the Promotion of Science (Grant- in-Aid for Scientific Research (C) 24530202) and Chuo University (Grant for Special Research and Joint Research Grant 0981). 104 Akio Matsumoto and Ferenc Szidarovszky CUBO 14, 3 (2012) 1 Introduction Dynamic models with time delays have many applications in many fields of quantitative sciences (see for example, Cushing (1977) and Invernizzi and Medio (1991)). The case of a single delay is well established in the literature (Hayes (1950) and Burger (1956)), however the presence of multiple delays makes analysis much more complicated. Sufficient and necessary conditions were derived for several classes of models giving a complete description of the stability region (Hale (1979), Hale and Huang (1993) and Piotrowska (2007)). In this paper a special class of dynamic systems is considered which are governed by delay differential equations with two delays. It is well known (Hayes (1950) and Cooke and Grossman (1982)) that stability can be lost or gained on a curve of stability switches, where an eigenvalue is pure complex. We will therefore determine these curves and then by bifurcation analysis char- acterize those segments where stability is gained or lost. In this way the stability region can be completely described. This paper is the continuation of our previous work (Matsumoto and Szidarovszky (2011)) where an elementary analysis was presented with a single delay. The paper is organized in the following way. Section 2 determines the curves where stability switches are possible and characterizes those segments where stability is lost or gained in the nonsymmetric cases. Section 3 discusses the symmetric case and Section 4 concludes the paper. 2 Stability Switches and Stability Region We will examine the asymptotical stability of the delay differential equation ẋ(t) + Ax(t − τ1) + Bx(t − τ2) = 0 (2.1) where A and B are positive constants. The characteristic equation can be obtained by looking for the solution in the exponential form αeλt. By substitution, αλeλt + Aαeλ(t−τ1) + Bαeλ(t−τ2) = 0 or λ + Ae−λτ1 + Be−λτ2 = 0. (2.2) Introduce the new variables ω = A A + B , 1 − ω = B A + B , λ̄ = λ A + B γ1 = τ1(A + B) and γ2 = τ2(A + B) to reduce equation (2.2) to the following: λ̄ + ωe−λ̄γ1 + (1 − ω)e−λ̄γ2 = 0. (2.3) CUBO 14, 3 (2012) An Elementary Study of a Class of Dynamic Systems ... 105 Because of symmetry we can assume that ω ≥ 1/2.In order to find the stability region in the (γ1, γ2) plane we will first characterize the cases when an eigenvalue is pure complex, that is, when λ̄ = iυ. We can assume that υ > 0, since if λ̄ is an eigenvalue, its complex conjugate is also an eigenvalue. Substituting λ̄ = iυ into equation (2.3) we have ιυ + ωe−iυγ1 + (1 − ω)e−iυγ2 = 0. In the special case of γ1 = 0, the equation becomes ιυ + ω + (1 − ω)e−iυγ2 = 0. The real and imaginary parts imply that ω + (1 − ω) cos(υγ2) = 0 υ − (1 − ω) sin(υγ2) = 0. We can assume first ω > 1/2, so from the first equation cos(υγ2) = − ω 1 − ω < −1 so no stability switch is possible. If ω = 1/2, then cos(υγ2) = −1 implying that sin(υγ2) = 0 and so υ = 0 showing that there is no pure complex root. Hence for γ1 = 0 the system is asymptotically stable with all γ2 ≥ 0. Assume now that γ1 > 0, γ2 ≥ 0. The real and imaginary parts give two equations: ω cos(υγ1) + (1 − ω) cos(υγ2) = 0 (2.4) and υ − ω sin(υγ1) − (1 − ω) sin(υγ2) = 0. (2.5) We consider the case of ω > 1/2 first and the symmetric case of ω = 1/2 will be discussed later. Introduce the variables x = sin(υγ1) and y = sin(υγ2), then (2.4) implies that ω2(1 − x2) = (1 − ω)2(1 − y2) or − ω2x2 + (1 − ω)2y2 = 1 − 2ω. (2.6) From (2.5), υ − ωx − (1 − ω)y = 0 106 Akio Matsumoto and Ferenc Szidarovszky CUBO 14, 3 (2012) implying that y = υ − ωx 1 − ω (2.7) Combining (2.6) and (2.7) yields −ω2x2 + (1 − ω)2 ( υ − ωx 1 − ω )2 = 1 − 2ω from which we can conclude that x = υ2 + 2ω − 1 2υω (2.8) and then from (2.7), y = υ2 − 2ω + 1 2υ(1 − ω) . (2.9) Equations (2.8) and (2.9) provide a parameterized curve in the (γ1, γ2) plane: sin(υγ1) = υ2 + 2ω − 1 2υω and sin(υγ2) = υ2 − 2ω + 1 2υ(1 − ω) . (2.10) In order to guarantee feasibility we have to satisfy − 1 ≤ υ2 + 2ω − 1 2υω ≤ 1 (2.11) and − 1 ≤ υ2 − 2ω + 1 2υ(1 − ω) ≤ 1. (2.12) Simple calculation shows that with ω > 1/2 these relations hold if and only if 2ω − 1 ≤ υ ≤ 1. From (2.10) we have four cases for γ1 and γ2, since γ1 = 1 υ { sin−1 ( υ2 + 2ω − 1 2υω ) + 2kπ } or γ1 = 1 υ { π − sin−1 ( υ2 + 2ω − 1 2υω ) + 2kπ } (k = 0, 1, 2, ...) and similarly γ2 = 1 υ { sin−1 ( υ2 − 2ω + 1 2υ(1 − ω) ) + 2nπ } or γ2 = 1 υ { π − sin−1 ( υ2 − 2ω + 1 2υ(1 − ω) ) + 2nπ } (n = 0, 1, 2, ...). CUBO 14, 3 (2012) An Elementary Study of a Class of Dynamic Systems ... 107 However from (2.4) we can see that cos(υγ1) and cos(υγ2) must have different signs, so we have only two possibilities: L1(k, n) :              γ1 = 1 υ { sin−1 ( υ2 + 2ω − 1 2υω ) + 2kπ } γ2 = 1 υ { π − sin−1 ( υ2 − 2ω + 1 2υ(1 − ω) ) + 2nπ } (2.13) and L2(k, n) :              γ1 = 1 υ { π − sin−1 ( υ2 + 2ω − 1 2υω ) + 2kπ } γ2 = 1 υ { sin−1 ( υ2 − 2ω + 1 2υ(1 − ω) ) + 2nπ } (2.14) For each υ ∈ [2ω − 1, 1] these equations determine the values of γ1 and γ2. At the initial point υ = 2ω − 1, we have υ2 + 2ω − 1 2υω = 1 and υ2 − 2ω + 1 2υ(1 − ω) = −1 and if υ = 1, then υ2 + 2ω − 1 2υω = 1 and υ2 − 2ω + 1 2υ(1 − ω) = 1. Therefore the starting point and end point of L1(k, n) are given as γs1 = 1 2ω − 1 ( π 2 + 2kπ ) , γs2 = 1 2ω − 1 ( 3π 2 + 2nπ ) and γe1 = π 2 + 2kπ and γe2 = π 2 + 2nπ. Similarly, the starting and end points of L2(k, n) are as follows: γS1 = 1 2ω − 1 ( π 2 + 2kπ ) , γS2 = 1 2ω − 1 ( − π 2 + 2nπ ) and γE1 = π 2 + 2kπ and γE2 = π 2 + 2nπ. Figure 1 illustrates the loci L1(k, n) and L2(k, n) of the corresponding points (γ1, γ2), when υ increases from 2ω − 1 to unity. The parameter value ω = 0.8 is selected. The red curves show L1(0, n) and the blue curves show L2(0, n) with n = 0, 1, 2, .... Notice that γ S 2 is infeasible at n = 0 and only the segment of L2(0, 0) between υ = √ 2ω − 1 and υ = 1 is feasible. With fixed value of k, L1(k, n) and L2(k, n) have the same end point, however the starting point of L1(k, n) is the same as that of L2(k, n + 1). Therefore the segments L1(k, n) and L2(k, n) with fixed k form a continuous curve with n = 0, 1, 2, .... 108 Akio Matsumoto and Ferenc Szidarovszky CUBO 14, 3 (2012) L1H0,0L L1H0,1L L1H0,2L L2H0,0L L2H0,1L L2H0,2L E0 E1 E2 S0 S1 Γ1 m Γ1 M Γ1 3Π 2 H2Ω-1L 7Π 2 H2Ω-1L Γ2 Figure 1. Partition curve in the (γ1, γ2) plane with fixing k = 0. Consider first the segment L1(k, n). Since ( υ2 − 2ω + 1 ) /(2υ(1 − ω)) is strictly increasing in υ, γ2 is strictly decreasing in υ. By differentiation and substitution of equation (2.4), we have ∂γ1 ∂υ ∣ ∣ ∣ ∣ L1 = − 1 υ2 ( sin−1 ( υ 2 +2ω−1 2υω ) + 2kπ ) + 1 υ √ 1− ( υ 2 +2ω−1 2υω ) 2 2υ(2υω)−(υ 2 +2ω−1)2ω 22υ2ω2 = − 1 υ2 υγ1 + 1 υ cos(υγ1) υ2 − 2ω + 1 2υ2ω = − 1 υ2 (υγ1 + tan(υγ2)) . (2.15) Consider next segment L2(k, n), similarly to (2.15) we can shown that ∂γ1 ∂υ ∣ ∣ ∣ ∣ L2 = − 1 υ2 (υγ1 + tan(υγ2)) which is the same as in L1(k, n), since from (2.14), cos(υγ1) < 0. Similarly ∂γ2 ∂υ ∣ ∣ ∣ ∣ L2 = − 1 υ2 (υγ2 + tan(υγ1)) (2.16) where we used again equation (2.4). In order to visualize the curves L1(k, n) and L2(k, n), we change the coordinates (γ1, γ2) to (υγ1, υγ2) to get the transformed segments: ℓ1(k, n) :              υγ1 = sin −1 ( υ2 + 2ω − 1 2υω ) + 2kπ υγ2 = π − sin −1 ( υ2 − 2ω + 1 2υ(1 − ω) ) + 2nπ (2.17) CUBO 14, 3 (2012) An Elementary Study of a Class of Dynamic Systems ... 109 and ℓ2(k, n) :              υγ1 = π − sin −1 ( υ2 + 2ω − 1 2υω ) + 2kπ υγ2 = sin −1 ( υ2 − 2ω + 1 2υ(1 − ω) ) + 2nπ (2.18) They also form a continuous curve with each fixed value of k, and they are periodic in both directions υγ1 and υγ2. Figure 2 shows them with k = 0 where the curves ℓ1(0, n) are shown in red color while the curves ℓ2(0, n) with blue color. l1H0,0L l1H0,1L l2H0,0L l2H0,1L Π�2xm x M ΝΓ1 - Π 2 0 Π 2 Π 3 Π 2 2 Π 5 Π 2 3 Π 7 Π 2 ΝΓ2 Figure 2. Partition curve in the (υγ1, υγ2) plane with fixing k = 0 We will next examine the directions of the stability switches on the different segments of the curves L1(k, n) and L2(k, n). We fix the value of γ2 and select γ1 as the bifurcation parameter, so the eigenvalues are functions of γ1 : λ̄ = λ(γ1). By differentiating the characteristic equation (2.3) implicitly with respect to γ1 we have dλ̄ dγ1 + ωe−λ̄γ1(− dλ̄ dγ1 γ1 − λ̄) + (1 − ω)e −λ̄γ2 ( − dλ̄ dγ1 γ2 ) = 0 implying that dλ̄ dγ1 = λ̄ωe−λ̄γ1 1 − ωγ1e −λ̄γ1 − (1 − ω)γ2e −λ̄γ2 (2.19) From (2.3) we have (1 − ω)e−λ̄γ2 = −λ̄ − ωe−λ̄γ1, so dλ̄ dγ1 = λ̄ωe−λ̄γ1 1 + λ̄γ2 + ω(γ2 − γ1)e −λ̄γ1 110 Akio Matsumoto and Ferenc Szidarovszky CUBO 14, 3 (2012) If λ̄ = ιυ, then dλ̄ dγ1 = iυω(cos(υγ1) − i sin(υγ1)) 1 + iυγ2 + ω(γ2 − γ1)(cos(υγ1) − i sin(υγ1)) and the real part of this expression has the same sign as υω sin(υγ1)[1 + ω(γ2 − γ1) cos(υγ1)] + υω cos(υγ1)[υγ2 − ω(γ2 − γ1) sin(υγ1)] = υω [sin(υγ1) + υγ2 cos(υγ1)] Hence Re ( dλ̄ dγ1 ) R 0 if and only if sin(υγ1) + υγ2 cos(υγ1) R 0 Consider first the case of crossing any segment L1(k, n) from the left. Here υγ1 ∈ (0, π/2], so both sin(υγ1) and cos(υγ2) are positive. Hence stability is lost everywhere on any segment of L1(k, n). Consider the case when crossing the segments of L2(k, n) from the left. Here υγ1 ∈ [π/2, π], so cos(υγ1) < 0. Combining (2.16) and the conditions for the sign of Re[dλ̄/dγ1], we have that Re ( dλ̄ dγ1 ) R 0 if and only if ∂γ2 ∂υ R 0. That is, stability is lost when γ2 increases in υ and stability is gained when γ2 decreases in υ. We can also show that at any intercept with L1(k, n) or L2(k, n) the complex root is single. Otherwise λ = iυ would satisfy both equations λ + ωe−λγ1 + (1 − ω)e−λγ2 = 0 and 1 − ωγ1e −λγ1 − (1 − ω)γ2e −λγ2 = 0, from which we have e−λγ1 = 1 + λγ2 (γ1 − γ2)ω and e−λγ2 = −1 − λγ1 (γ1 − γ2)(1 − ω) . By substituting λ = iυ and comparing the real and imaginary parts yield sin(υγ1) + υγ2 cos(υγ1) = sin(υγ2) + υγ1 cos(υγ2) = 0. Therefore this intercept is at an extremum in υ of a segment L1(k, n) and also at an extremum of a segment L2(k̄, n̄) which is impossible. For each γ2 > 0, define m(γ2) = min γ1 {(γ1, γ2) ∈ L1(k, n) ∪ L2(k, n), k, n ≥ 0} (2.20) At γ1 = 0 the system is asymptotically stable with all γ2 > 0. With fixed value of γ2 by increasing the value of γ1 the first intercept with m(γ2) should be a stability loss, since there is no stability switch before. Then by increasing the value of γ1 further, the next intercept is either a stability CUBO 14, 3 (2012) An Elementary Study of a Class of Dynamic Systems ... 111 gain or a stability loss. In the first case the equilibrium becomes asymptotically stable. In the second case the equilibrium remains unstable, which will not change even if the next intercept is an stability gain, since the real part of only one eigenvalue becomes negative. Consider next a point (γ∗ 1 , γ∗ 2 ) with γ∗ 1 , γ∗ 2 > 0 which is not located on any curve L1(k, n) or L2(k, n), and consider the horizontal line γ2 = γ ∗ 2 and its segment with γ1 ∈ (0, γ∗1). If it has no stability switch, then the equilibrium is asymptotically stable. This is the case even if the number of stability losses equals the number of stability gains, otherwise the equilibrium is unstable. The stability region is shown as the shaded region in Figure 1. Notice that this is the same result which was obtained earlier by Hale and Huang (1993) by using different approach. 3 The Symmetric Case Assume next that ω = 1/2. Then equations (2.4) and (2.5) imply that cos(υγ1) + cos(υγ2) = 0 υ − 1 2 (sin(υγ1) + sin(υγ2)) = 0 and the curves L1(k, n) and L2(k, n) are simplified as follows: L1(k, n) :          γ1 = 1 υ ( sin−1(υ) + 2kπ ) γ2 = 1 υ ( π − sin−1(υ) + 2nπ ) (3.1) and L2(k, n) :          γ1 = 1 υ ( π − sin−1(υ) + 2kπ ) γ2 = 1 υ ( sin−1(υ) + 2nπ ) (3.2) which are shown in Figure 3. The same argument as shown above for the nonsymmetric case can be applied here as well to show that stability region is left of L1(0, 0) and below L2(0, 0), where the shape of the stability region differs from that of the nonsymmetric case. It is illustrated in Figure 3 by the shaded domain. 112 Akio Matsumoto and Ferenc Szidarovszky CUBO 14, 3 (2012) l1H0,0L l1H0,1L l2H0,0L l2H0,1L Π�2xm x M ΝΓ1 - Π 2 0 Π 2 Π 3 Π 2 2 Π 5 Π 2 3 Π 7 Π 2 ΝΓ2 Figure 3. Partition curve in the (γ1, γ2) plane with ω = 1 2 Notice that at each segment of ℓ2(k, n) there are at most two intercepts with the υγ2 = − tan(υγ1) curve, so the same holds for L2(k, n). At every other point Re[dλ̄/dγ1] 6= 0, so at these points Hopf bifurcation occurs giving the possibility of the birth of limit cycles. 4 Conclusions Ordinary differential equation were examined with two delays. After finding the possible stability switches, their curves were determined. Hopf bifurcation was used to find segments with stability losses and stability gains. The boundary of the stability region are the γ2 = 0, γ1 = 0 and a continuous curve consisting of certain portions of the segments L1(0, n) and L2(0, n). All other points on the curves L1(k, n) and L2(k, n) for k ≥ 1 do not lead to actual stability switches, since the system is already unstable. Received: July 2011. Revised: June 2012. References [1] Burger, E. (1956), On the Stability of Certain Economic Systems. Econometrica, 24(4), 488-493. [2] Cooke, K. L. and Z. Grossman (1982), Discrete Delay, Distributed Delay and Stability Switches. J. of Math. Analysis and Appl., 86, 592-627. [3] Cushing, J. M. (1977), Integro-Differential Equations and Delay Models in Population Dynam- ics. Springer-Verlag, Berlin/Heidelberg/New York. [4] Hayes, N. D. (1950), Roots of the Transcendental Equation Associated with a Certain Difference-Differential Equation. J. of the London Math. Society, 25, 226-232. CUBO 14, 3 (2012) An Elementary Study of a Class of Dynamic Systems ... 113 [5] Hale, J. (1979), Nonlinear Oscillations in Equations with Delays. In Nonlinear Oscillations in Biology (K. C. Hoppenstadt, ed.). Lectures in Applied Mathematics, 17, 157-185. [6] Hale, J. and W. Huang (1993), Global Geometry of the Stable Regions for Two Delay Differ- ential Equations. J. of Math. Analysis and Appl., 178, 344-362. [7] Invernizzi, S. and A. Medio (1991), On Lags and Chaos in Economic Dynamic Models. Journal of Math. Econ., 20, 521-550. [8] Matsumoto, A. and F. Szidarovszky (2011), An Elementary Study of a Class of Dynamic Sys- tems with Single Delay. Mimeo: DP 161 (http://www2.chuo-u.ac.jp/keizaiken/discuss.htm), Institute of Economic Research, Chuo University, Tokyo, Japan. [9] Piotrowska, M. (2007), A Remark on the ODE with Two Discrete Delays. Journal of Math. Analysis and Appl., 329, 664-676.