Shilnikov-type dynamics in three-dimensional piecewise smooth maps

Chaos, Solitons and Fractals 133 (2020) 109655

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Shilnikov-type dynamics in three-dimensional piecewise smooth maps Indrava Roy a,∗, Mahashweta Patra b, Soumitro Banerjee c a

Institute of Mathematical Sciences, HBNI, Chennai, India Indian Institute of Technology Madras, India c Indian Institute of Science Education and Research, Kolkata, India b

a r t i c l e

i n f o

Article history: Received 4 September 2019 Revised 20 December 2019 Accepted 27 January 2020

Keywords: Three dimensional piecewise smooth map Shilnikov chaos Stable manifold and unstable manifold

a b s t r a c t We show the existence of Shilnikov-type dynamics and bifurcation behaviour in general discrete threedimensional piecewise smooth maps and give analytical results for the occurence of such dynamical behaviour. Our main example in fact shows a ‘two-sided’ Shilnikov dynamics, i.e. simultaneous looping and homoclinic intersection of the one-dimensional eigenmanifolds of fixed points on both sides of the border. We also present two complementary methods to analyse the return time of an orbit to the border: one based on recursion and another based on complex interpolation.

1. Introduction The phenomenon of chaos is observed in many nonlinear deterministic systems in both experimental and computer-simulation contexts. In 1965, L. Shilnikov (also written Šil’nikov) [14] showed that in a continuous-time dynamical system, if the real eigenvalue of an equilibrium point has a larger magnitude than the real part of the complex conjugate eigenvalues, then there are horseshoes present in the return maps defined near the homoclinic orbit. Such behavior normally appears when a parameter is varied towards the homoclinic condition associated with a saddle focus. The criteria are called Shilnikov criteria and the resulting dynamics is called Shilnikov chaos. Shilnikov attractors can be found in many different systems, including Rössler system, Arenodo-Coullet-Tresser systems, and Rosenzweig-Mac-Arthur system. Shilnikov chaos is observed in many physical systems, including a single mode lasers with feedback [1] and with a saturable absorber [5], the BelousovZhabotinskii reaction, a glow discharge plasma [3], an optically bistable device [10], a multimode laser [18] and some other systems. Theoretical and experimental study of a chaotic attractor in a discrete time model, appearing via the mechanism described by Shilnikov is shown in a CO2 laser in [11]. Similar situations may also occur in maps. Attractors of a spiral type that appear in accordance with this scenario in discrete dynamical systems are called discrete Shilnikov attractors. First exam-

∗

Corresponding author. E-mail address: [email protected] (I. Roy).

https://doi.org/10.1016/j.chaos.2020.109655 0960-0779/© 2020 Elsevier Ltd. All rights reserved.

© 2020 Elsevier Ltd. All rights reserved.

ples of such attractors were found in 3-D generalized Hénon map [20]. Pisarchik et al. [10] have reported the first experimental observation of the discrete behavior of Shilnikov chaos. Zhou et al. [19] have shown that there exist only three kinds of chaos—homoclinic chaos, heteroclinic chaos, and a combination of homoclinic and heteroclinic chaos. They construct a new chaotic system of quadratic polynomial ordinary differential equations (ODE) in three dimensions, which has a single equilibrium point. They rigorously prove that this system satisfies all conditions stated in the Shilnikov theorem [17], which clearly reveals its chaos formation mechanism and implies the existence of Smale horseshoes. Silva [15] has given a brief introduction of Shilnikov’s method to detect analytically the presence of chaos in continuous autonomous systems. As an application to a piecewise linear system they took Chua’s circuit and have shown the existence of a homoclinic orbit from a saddle focus as well as heteroclinic orbit between two saddle foci. The mechanism of formation of Shilnikov chaos has not yet been investigated in piecewise smooth discrete dynamical systems. The questions we address in this paper are: What are the criteria for Shilnikov chaos to occur in a piecewise smooth discretetime dynamical system? What sort of theoretical conditions can be given that can guarantee existence or non-existence of transverse homoclinic intersections, which will help create more efficient computer programs that can search for such intersections? Since fixed points of saddle-focus type do not appear in a 2-D piecewise linear map, to investigate the Shilnikov phenomenon in a piecewise smooth system we take a 3-D piecewise linear normal form map [4,9,12] and answer the questions mentioned above. We

2

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

derive an necessary and sufficient condition analytically for the occurrence of a homoclinic intersection, thereby, the occurrence of a chaotic orbit. This condition is obtained keeping in view that checking the existence of a homoclinic intersection requires a computer simulation; the implementation of the algorithm is rendered more convenient by our analytic condition. In particular, it provides a finite range of iterations which should be checked, outside of which no homoclinic intersection on first return is possible. We also show a motivating numerical example which exhibits Shilnikov-type behaviour, and show the utility of our methods with the example. The two-sided Shilnikov dynamics, (i.e., Shilnikovtype behaviour exhibited simultaneously by fixed points on both sides of the border) that this example shows is an interesting feature that occurs in three-dimensional piecewise smooth systems and may well be unique to them. The paper is organized as follows. The description of the normal form for a three-dimensional piecewise smooth system is given in Section 2. A motivating example is shown for a two-sided Shilnikov-type dynamics in Section 3. Section 4 is the main theoretical part of the paper which gives two complementary methods to analyse the existence of transverse homoclinic intersections. The last section contains discussions of the results obtained and concluding remarks. 2. 3-Dimensional piecewise smooth maps: system description The piecewise linear approximation of a general piecewise smooth 3D system evaluated in a close neighborhood of the border, called the ‘normal form’ map [4,6,9,12] is given by



Xn+1 = Fμ (Xn ) =

Al Xn + μC, if xn ≤ 0 Ar Xn + μC, if xn ≥ 0

(2.1)

where Xn = (xn , yn , zn ) ∈ R3 , C = (1, 0, 0 ) ∈ R3 and μ is a realvalued parameter. The phase space of this map is divided by the ‘border’ Xb : {x = 0} into two regions L := {(x, y, z ) ∈ R3 : x < 0} and R := {(x, y, z ) ∈ R3 : x > 0}. We shall frequently refer to L and R as the ‘left’ side and the ‘right’ side of the border, respectively. In each region, the dynamics is governed by an affine map and the equations are continuous across Xb . Al and Ar are real valued 3 × 3 matrices



τl Al = −σl δl

1 0 0



0 1 0



τr and Ar = −σr δr

1 0 0

0 1 0



If λ1 , λ2 and λ3 are the eigenvalues of the Jacobian matrix of the original PWS map evaluated at a fixed point placed on the left side close to the border, then the parameters of the matrix Al are simply the trace τl = λ1 + λ2 + λ3 , the second trace σl = λ1 λ2 + λ2 λ3 + λ3 λ1 and the determinant δl = λ1 λ2 λ3 . The parameters of the matrix Ar depends, in a similar manner, on the eigenvalues of the Jacobian matrix computed at a fixed point located on the right side. The fixed points of the system in both sides of the boundary are given by



L∗ =

 ∗

R =

μ(−σl + δl ) μδl , 1 − τl + σl − δl 1 − τl + σl − δl 1 − τl + σl − δl μ



,

μ(−σr + δr ) μδr μ , , 1 − τr + σr − δr 1 − τr + σr − δr 1 − τr + σr − δr



If L∗ lies on the left side of the border, it is called admissible, otherwise it is called virtual. Similarly, R∗ is admissible if it lies on the right side of the border and virtual otherwise. We will assume the generic condition μ = 0 for our discussions.

3. Motivating example for Shilnikov-type dynamics and bifurcation scenario Let us consider the system given by Eq. (2.1) for the following parameter values:

τl = 1.0, σl = −0.25, δl = 0.3, τr = 0.58, σr = 0.38, δr = −1.27, μ = 1.0 For these values, the left fixed point is admissible, located at L∗ = (−1.8182, −1, −0.5455 ). The right fixed point is also admissible and located at (0.4831, − 0.7971, − 0.6135 ). The matrix Al has a positive unstable eigenvalue λ1,l = 1.3499,

and a pair of stable complex eigenvalues λ2,l , λ2,l = −0.1749 ± 0.4378i with absolute value |λ2,l | = 0.4714. Therefore the left fixed point is of saddle-focus type. The matrix Ar governing the dynamics on the right side of the border x = 0, has a negative stable eigenvalue λ1,r = −0.8251 and pair of unstable complex eigenvalues λ2,r , λ2,r = 0.7025 ± 1.0226i with absolute value 1.2407. Therefore, R∗ is a flip saddle. At this setting of parameter values, there exists a transverse homoclinic intersection between the 1-dimensional stable manifold Sr and the 2-dimensional unstable plane Ur of R∗ (see Fig. 1(d)), resulting in a chaotic attractor which is stable under small perturbations for τ r < 0.62, see Fig. 1(a). The Lyapunov exponents and bifurcation diagrams are shown in Fig. 2. We also observe that the unstable manifold Ul for the left fixed point L∗ grows until it hits the border (say at B∗l ) and passes to the right side, and then loops back towards the left side due to the unstable nature of R∗ . When Ul crosses the border again, the first iterate that crosses the border lies ‘above’ the stable manifold Sl , i.e., in the same side of B∗l with respect to Sl . Thus the dynamics mimics the classical Shilnikov looping behaviour of the unstable manifold Ul , and a Shilnikov-type bifurcation occurs for these parameter values when τ r increases from 0.58 to 0.62. At the critical value τr = 0.62, the unstable manifold Ul has a transverse homoclinic intersection with the stable plane Sl , which results in the sudden vanishing of the chaotic attractor and appearance of a chaotic saddle. See Fig. 1(c) which shows the transverse homoclinic intersection of Ul and Sl . Note that due to the transverse intersection of Sr and Ur , chaotic dynamics is still present in the system, but there is no attractor. The basin of attraction is plotted in the planes z = constant in Fig. 3, for the values z = −1.1, 0 and 1.0 at τr = 0.58 and τr = 0.63. On the other hand, increasing τ l from 1 to 1.05 causes the homoclinic intersection of Ur and Sr to be broken, and thus the chaotic dynamics also ceases to occur. Thus near the critical values τr = 0.62 and τl = 1.05, Shilnikov-type dynamics occur on both sides of the border. This kind of dynamics is most likely unique to piecewise smooth systems due to its inherent asymmetry. 4. Conditions for existence of homoclinic intersections A homoclinic intersection between stable and unstable manifold implies an infinite number of intersections, therefore a horseshoe structure is born. We can give a general condition for occurrence of chaos through the occurrence of homoclinic intersection. Let L∗ be an admissible fixed point which has an unstable eigenvalue and a pair of stable eigenvalues. It can be both real or complex. In Fig. 4 we have shown a schematic diagram to illustrate the algorithm for finding the conditions of occurrence of homoclinic intersection. Grey plane shows the border. L∗ is shown in black asterisk, stable manifold (SL ) is shown in red, unstable eigenvector is shown in blue, iterations of the fold points are shown in blue points. If the iterations of the unstable eigenvector or the fold points of L∗ crosses the stable manifold, there is homoclinic intersection between stable manifold and unstable manifold of L∗ .

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

3

Fig. 1. (a) Chaotic attractor; (b), (c) Shilnikov-type dynamics and transverse homoclinic intersection between Sl (black) and Ul (blue) at (b) τr = 0.58 and (c) τr = 0.62. (d) The homoclinic intersection of Ur (blue) and Sr (red). (e) is the zoomed version of figure (d). The rest of the parameters are the same as in the beginning of Section 3; the green and red stars show the location of L∗ and R∗ , respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

If the intersection point is on the same side of the fixed point then there is a homoclinic intersection, and therefore, chaotic dynamics must occur.

Fig. 2. (a). Lyapunov exponent plot, and (b) Bifurcation diagram for the same parameter values as in the beginning of Section 3.

The algorithm for finding whether there is homoclinic intersection or not is as follows: 1. Calculate the unstable eigenvector. Calculate the point P0 where the unstable eigenvector crosses the border x = 0. Consider that point as the initial point. 2. Calculate the nth iteration of this initial point i.e., Pn , where n is the minimum number of iterations needed for the unstable manifold to cross the border again (assumed that the iterates of P0 , i.e., P1 must have crossed the border). 3. Calculate Pn+1 . 4. Calculate the equation of the stable eigenplane. Check whether Pn+1 and Pn are on the same side of the stable eigenplane equation or not. 5. Calculate the intersection point of the stable eigenplane and the line going through the Pn and Pn+1 . Check whether the intersection point is in the same side of the fixed point or not.

Here, to calculate the Pn , i.e., the nth iteration, we have taken two complementary approaches. The first one is based on a recursion method following [2] (see also Saha and Banerjee [13]. The recursion stops as soon as the orbit crosses the border, at which point it is easily checked whether there exists a transverse homoclinic intersection with the 2-dimensional stable manifold of L∗ . The second method is then used to provide upper bounds for the number of iterations that need to be checked for the border crossing and subsequent homoclinic intersection. This method is based on a complex interpolation scheme whereby we define fractional iterations of our system in Eq. (2.1). This interpolation is then used to get a transcendental equation in the positive reals; the least positive solution of this equation then provides the border return time. In this algorithm we assume that the orbit will eventually return to the border and the algorithm to check for homoclinic intersection is valid only under this assumption. 4.1. Recursive method Here in this subsection, we are introducing a recursive method. Using this method, we can calculate Pn directly from P0 without calculating P1 , P2 , . . . given the condition that all the iterations are on the same side of the border. In this method we do not need to calculate all the iterations P1 , P2 , . . . , to calculate Pn . Thus the calculation becomes more efficient. Let us consider a situation where a period-1 fixed point L∗ = (xl , yl , zl ) is admissible. The unstable eigenvector crosses the x = 0

4

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

Fig. 3. (a), (c), (e). Basins of attraction for the planes z = −1.1, 0, 1.1 and τr = 0.58, (b), (d), (f) Basins of attraction for the planes z = −1.1, 0, 1.1 and τr = 0.63. The xcoordinates and the y-coordinates both lie in the interval [−1, 1]. color bar shows the vector norm of the orbit point after 10 0 0 iterations, normalized to lie between 0 and 1. Other parameter values are same as in the beginning of Section 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

plane at P0 = (0, C1 , C2 ) and the image of this point is the 1st fold point of the unstable manifold of L∗ and its coordinate is v P1 = (x1 , y1 , z1 ) = (C1 + μ, C2 , 0 ), where, C1 = yl − v21 xl , and C2 = 11 v31 zl − v xl . Here V1 = (v11 , v21 , v31 ) is the unstable eigenvector as11 sociated with Al . We shall assume that both matrices Al and Ar are non-singular and do not have an eigenvalue 1. Further, we suppose that C1 + μ > 0, so that P1 lies on the right side of the border. Now we need to calculate the minimum number of iterations needed to cross the border again. Let’s suppose n iterations are

needed to come again to the left side. We shall refer to n as the border return time.

⎛ ⎞ 1

Pn = F n (P0 ) = Anr (P0 ) + μ(Ar − I )−1 (Anr − I )⎝0⎠

(4.1)

0 Suppose that δ r = 0. We need to calculate the nth power of a matrix A. We use a recursion method similar to the one used by

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

5

Note that the recursion relations can be solved explicitly using standard methods (for the 2-dimensional case see [2], appendix B) and we omit it here. One can also use a recursion scheme to compute the matrix (Ar − I )−1 (Anr − I ) efficiently, analogous to the recursion method given in [13], Appendix E. With this recursion scheme, we can thus calculate P1 , P2 , . . . , Pn efficiently until the xcoordinate of Pn is negative. So, we get the minimum number of iterations n at which the orbit crosses the border again. Equation of the stable eigenplane: L∗ is a saddle fixed point which has one unstable eigenvector and two stable eigenvectors. Equation of the stable eigenplane containing the two stable eigenvectors in vector notation, where V2 and V3 are the two stable eigenvectors is

r · (V2 × V3 ) = 0 In Cartesian co-ordinates, we can express this equation as

a·x+b·y+c·z+d = 0 Fig. 4. Schematic diagram to illustrate the algorithm. Grey plane shows the border. L∗ is shown in black asterisk, stable manifold (SL ) is shown in red, unstable eigenvector is shown in blue, iterations of the fold points are shown in blue points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Avrutin et al. [2] and calculate An for n ≥ 1 as



An =

 =

an dn gn

bn en hn

cn fn in



an δr .an−2 − σr .an−1 δr .an−1

an−1 δr .an−3 − σr .an−2 δr .an−2

an−2 δr .an−4 − σr .an−3 δr .an−3



(4.2) which follows the recursive equation given by

an = τr an−1 − σr an−2 + δr an−3

(4.3)

where the initial conditions are

a0 = 1, a−1 = 0, a−2 = 0, a−3 =

1

δr

Assuming all the iterations are on the right side, after nth iteration the coordinates are given by the following expression:

 

1 Pn = F n (P0 ) = Anr (P0 ) + μ(Ar − I )−1 (Anr − I ) 0 0



=

an dn gn



+μ

 =

 Pn =

an dn gn

bn en hn a11 a21 a31 bn en hn

cn fn in

  x0 y0 z0

a12 a22 a32

a13 a23 a33

cn fn in

x0 y0 z0



an − 1 dn gn



 

+μ

a11 a21 a31

bn en − 1 hn a12 a22 a32

 

cn fn in − 1 a13 a23 a33



1 0 0



μ(a11 (−1 + an ) + a12 dn + a13 gn ) + an x0 + bn y0 + cn z0 μ(a21 (−1 + an ) + a22 dn + a23 gn ) + dn x0 + en y0 + fn z0 μ(a31 (−1 + an ) + a32 dn + a33 gn ) + gn x0 + hn y0 + in z0

where (Ar − I )−1 =

a11 a21 a31

a12 a22 a32

a13 a23 a33

In this subsection we present an interpolation method which considers the positive real exponents of a non-singular real matrix A, i.e., matrices of the form At for t ∈ R+ . Note that such exponents are defined using the complex logarithm of A and the result takes values in the algebra of matrices with complex entries. The matrix At is defined by the following formula:

At := exp(t ln A ) In general, this is a multi-valued function due to the complex logarithm. However, in this paper we shall always restrict ourselves to the principal value of the complex logarithm to get a unique matrix At . Therefore, unless otherwise stated ln (z) will denote the principal value of the complex logarithm.  z1 0 0 For a 3 × 3-diagonal matrix D = 0 z2 0 with non-zero 0 0 z3 complex entries zi , i = 1, 2, 3, the matrix ln (D) is easily defined as

ln(D ) :=



(4.4)



4.2. Complex interpolation and estimation of border return time

ln(z1 ) 0 0

0 ln(z2 ) 0



0 0 ln(z3 )

Therefore, the matrix Dt for t ∈ R+ is given by





As the eigenplane passes through the fixed point (xl , yl , zl ), the constant term d can be calculated as (−a · xl − b · yl − c · zl ). We then check whether the points Pn+1 and Pn are on opposite sides of the stable eigenplane or not. If they are on opposite sides, we infer that the unstable manifold must have gone through the stable manifold, and therefore there is homoclinic intersection resulting in the occurrence of chaotic dynamics. Similar calculation can be carried out for the fixed point at R∗ . In the next subsection we give a method based on complex interpolation of the discrete system given by Eq. (2.1) to estimate the border return time.



an − 1 dn gn

where (a, b, c ) := V2 × V3

(4.5)

t

D :=

exp(t ln(z1 )) 0 0

0 exp(t ln(z2 )) 0



0 0 exp(t ln(z3 ))

Remark 4.1. In case of a matrix A having repeated eigenvalues, given its block-diagonal Jordan normal form, one can also define the matrix ln (A) using the so-called Jordan-Chevalley decomposition of A into a sum of two matrices S and N, where S is a semisimple matrix and N a nilpotent matrix. For more details, we refer the reader to [7], chapter 4. In this paper however, we restrict ourselves to the case of distinct eigenvalues only.

6

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

Turning to the system given by Eq. (2.1), we wish to diagonalize the matrices Al and Ar , and then use the above formulae to define the continuous complex interpolation of the system. Let Pl (resp. Pr ) be the complex base change matrix for Al (resp. Ar ). Let us assume that neither Al nor Ar has repeated roots, so that the matrices Dl = Pl Al Pl−1 and Dr = Pr Ar Pr−1 are diagonal (with complex entries). This excludes the case of block-diagonal matrices, although much of the analysis goes through once the matrices ln (Dl ) (or ln (Dr )) have been defined using the Jordan-Chevalley decomposition, as in Remark (4.1). We also assume as before that the matrices Al and Ar do not have any eigenvalue equal to 0 or 1. As a consequence, the matrices Al , Ar , Al − I and Ar − I are invertible. For t ∈ R+ , we set Atl := Pl−1 Dtl Pl and similarly Atr := Pr−1 Dtr Pr . We use the notation (z) to denote the real part of any m-tuple of complex numbers z ∈ Cm . Definition 4.2. The complex interpolation of the discrete system in Eq. (2.1) is given by the following formula: for any t ∈ [0, 1] and a starting point X0 = (x0 , y0 , z0 ) ∈ R3 such that y0 + μ = 0,



X (t ) ≡ Fμ (X0 ) := t

4.2.1. Estimating the border return time In this section we give some results regarding upper bounds for the border return time as well as algebraic conditions for not returning to the border. This is important in view of implementing the algorithm based on the recursive Eq. (4.2) in Section 4.1. Indeed, in order to check the transverse intersection condition by implementing the algorithm, one must determine in advance an

Atl X0 + μ(Al − I )−1 (Atl − I )C, if x0 < 0 or (x0 = 0 and y0 < −μ) Atr X0 + μ(Ar − I )−1 (Atr − I )C, if x0 > 0 or (x0 = 0 and y0 > −μ)

The system given by Eq. (4.6) yields a continuous, piecewise smooth curve (X(t)) (the real part of the complex function X(t)) in R3 for t ∈ [0, 1]. The definition for the complex interpolation for any t ∈ R+ is then straightforward to deduce: {t }

meant to imply that it must intersect with the stable manifold infinitely many times. In particular, Fμt does not follow the semi-group composition law Fμt+s = Fμt ◦ Fμs in general. 2. We have nevertheless found via numerical simulations that in many cases, the companion orbits (or companion stable or unstable manifolds) do capture a lot of global dynamics of the system in Eq. (2.1), see Fig. (6) below in which a companion orbit as well as the companion unstable manifold of Ul for the motivating example in Section 3 is given.

t

Fμt (X0 ) := Fμ (Fμ (X0 )) where t is the greatest integer less than or equal to t and {t} is the fractional part of t. It follows immediately from the definition that the complex interpolation Fμt agrees with the system in Eq. (2.1) whenever t is an integer. Since the first iterate of points on the line y0 = −μ lie again on the Y − Z plane, it is not evident whether to use the formula for the left side or the right side. Hence, we have excluded the line y0 = −μ in the Y − Z plane so that there is no ambiguity in the definition of the interpolation. However, the interpolation formulae can be extended to this case according to the situation where the second iterate of X0 under Fμ has x-coordinate positive or negative, equivalently, whether z0 > −μ or z0 < −μ. That leaves us with the case X0 = (0, −μ, −μ ) whose second iterate also lies on the Y − Z plane. For this point one should consider on which side of the border the third iterate lies. This “direction” for the third iterate for X0 can be decided simply by the condition μ > 0 or μ < 0 (see Fig. 5 for a schematic diagram). In this way one can extend the interpolation scheme to every point in the Y − Z plane.

(4.6)

upper bound for the number of iterations n since the border return time equation is transcendental in nature (see Eq. (4.8)) and therefore has no general solution in closed form. Thus, one has to rely on a computer program to find the least positive solution for such equations. For such a program, hard-coding a fixed number of iterates to check the border crossing condition for different parameter values is clearly not a reliable method, and this is where upper bounds for various parameter regions help us determine the number of iterates to check whether a transverse homoclinic intersection exists after the orbit crosses the border. Let us consider the system in Eq. (2.1) under the change of coordinates corresponding to the eigenbasis of Al . Let us denote the eigenvalues of Al by (λ1 , λ2 , λ3 ). Let Pl = ( pi j )i, j=1,2,3 be the change of basis matrix with respect to the given eigenvalues, and let Pl−1 be given by a matrix (ρi j )i, j=1,2,3 . Thus, the new co-ordinates, which shall be denoted in vector form by χ = (χ1 , χ2 , χ3 ) and the old (standard Euclidean) co-ordinates denoted by X = (x, y, z )

Definition 4.3. We shall call the interpolation X(t) the companion orbit of the dynamical system (2.1) starting at X0 . If X0 lies on the stable (resp. unstable) manifold of L∗ or R∗ , we shall call X(t) the companion stable (resp. companion unstable) manifold. Since the companion orbit must agree with the actual orbit at integer points, the following lemma is immediate: Lemma 4.4. The orbit of a point under Fμ is bounded if and only if the corresponding companion orbit Fμt is bounded for all t ∈ R+ . The complex interpolation function Fμt will be used to estimate the border return time in the next subsection. A few remarks are in order. Remark 4.5. 1. The terminology ‘companion stable’ does not imply that the curve X(t) itself is stable under iterations of Fμ , rather it is

Fig. 5. A schematic diagram showing the direction of orbits. The red and blue points lying on the Y − Z plane outside the line y0 = −μ take only one iteration to land in either the left (L) or right (R) region. The green point lying on y0 = −μ takes in general two iterations to move towards L or R, while the purple point (0, −μ, −μ ) takes three iterations to do so. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

7

mentioned, the least positive real solution t0 of Eq. (4.9) gives the crossing iteration n = t0 , i.e., the nth iteration lies on the left side of the border while the (n + 1 )th iteration passes to the right side of the border. Note that the above equation interpolates the actual orbit of X0 under Fμ correctly only for t ≤ t0 , after which the corresponding equation for the matrix Ar must be used. If the matrix Al has all positive real eigenvalues, then of course the expression within the big brackets in Eq. (4.9) is already real. However, if any of the eigenvalues are negative or complex, then there is a non-zero complex part of the interpolation, which we nevertheless ignore. Remark 4.6. A special case arises when the right fixed point R∗ has a flip unstable direction Ur and the point X0 is the intersection point of Ur with the border x = 0. In that case, the first iterate will never land on the left side because it flips to the ‘opposite side’ of the flip direction with respect to the orientation of Ur on the right side of the border, and the next iterate then lands on Ur again but now on the left side of the border (since the direction is unstable). Therefore, in this case the starting point X0 has non-zero x-coordinate (in the standard Euclidean basis) and the Eq. (4.9) has to be modified accordingly. This modification is nonetheless straightforward and we omit the details here. The border return time is then estimated as t0 + 2, where t0 is the least positive solution to the modified equation. The corresponding function f(t) has the same general form as in the positive eigenvalue case (see Eqs. (4.8) and (4.10) below).

Fig. 6. (a) An example of a companion orbit, and (b) The companion unstable manifold for Ul , for the motivating example of Section 3. In (b), the black dots along the companion orbit represent the actual orbit points, while the green dot at the center is the right fixed point R∗ . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

f (t ) := α1 λt1 + 2α2 cos(t θ0 + δ )r0t + α3 = 0,

are related by the change of base matrix Pl :

χ=

Pl−1 X

  χ1 x −1 y ⇔ χ2 = Pl χ3 z

Under the χ -coordinates, for a starting point ξ0 = (ξ01 , ξ02 , ξ03 ) ∈ R3 , the complex interpolation equation takes

a particularly simple form. For i = 1, 2, 3, we have:

χi (t ) = λti ξ0i +

λ −1 μ λi − 1 i t i

(4.7)

where μi = μρ1i for i = 1, 2, 3. To compute the border return time for a point X0 = (0, y0 , z0 ) in the old-coordinate system for which y0 + μ < 0 (which means the first iterate lands on the left side of the border), we derive the following equation. Let ξ0 = Pl−1 X0 . Then, the least positive real solution of the equation below gives the border return time.

( p11 χ1 (t ) + p12 χ2 (t ) + p13 χ3 (t )) = 0 When expressed in the old Eq. (4.8) takes the following form:

co-ordinate

  λt − 1 μ1 λt1 (ρ12 y0 + ρ13 z0 ) + 1 λ1 − 1   λt2 − 1 t + p12 λ2 (ρ22 y0 + ρ23 z0 ) + μ λ2 − 1 2   λt3 − 1 t + p13 λ3 (ρ32 y0 + ρ33 z0 ) + μ =0 λ3 − 1 3

4.2.2. Case A: one positive unstable eigenvalue and a pair of stable complex conjugate eigenvalues: Let us analyse Eq. (4.9) in more detail for the special case when λ1 > 1, and λ2 = λ3 = r0 eiθ0 , r0 > 0, θ0 ∈ (−π , π ] are a pair of complex eigenvalues with absolute value less than 1. Eq. (4.9) then takes the following general form:

(4.8) system,

the

 p11

t≥0

(4.10)

where α1 , α2 , α3 , δ ∈ R. Our goal is to give upper bounds for the smallest positive root of f(t). Note that the cosine term makes f(t) an oscillating function. It has two enveloping curves given by

α1 λt1 + 2α2 r0t + α3 , f (t ) = α1 λt1 − 2α2 r0t + α3 f + (t ) =

and

−

(4.11)

The condition f − (t ) ≤ f (t ) ≤ f + (t ) is then satisfied for all t ≥ 0. The functions of the kind f ± (t) are sometimes called exponential polynomials or generalized Dirichlet polynomials [8,16]. We thus get the following general form of the enveloping curves given in Eq. (4.14):

g(t ) := a1 κ1t + a2 κ2t + a3 ,

t ≥ 0;

κ1 > 1 > κ2 > 0,

and ai ∈ R for i = 1, 2, 3.

(4.12)

We are interested in the location of zeros of g(t) in the interval (0, ∞). Note that for such functions, Descartes’ rule of signs is applicable, and we have the following nice theorem, sometimes called the lost cousin of the Fundamental Theorem of algebra: Theorem 4.7. [8] Let σ (g) be the number of sign changes in the sequence (a1 , a2 , a3 ). Then the number of zeros Z(g) of g(t) (counted with multiplicity) in the interval (0, ∞), is bounded above by σ (g). Moreover, the difference σ (g) − Z (g) is always an even number, and Z (g ) ≥ Z (g) − 1.

(4.9)

It is straighforward to check that t = 0 is always a solution to Eq. (4.9), which corresponds to the starting position X0 . As already

The following corollaries are immediate. Corollary 4.8. With the same assumptions as in 4.7, the number of zeros of g(t) is at most 2.

8

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

Corollary 4.9. With the same assumptions as in 4.7, if σ (g) = 0 then g(t) has no solution. Note that our original goal was to find roots of the function f(t) of Eq. (4.10). Since f ± (t) bound f(t), the following lemma is obvious: Lemma 4.10. If f(t) has at least one solution for t > 0, then either f + (t ) or f − (t ) must have at least one solution for t > 0. The least positive solution of f(t) is therefore bounded above by the largest positive solution of f ± (t). Therefore to find an upper bound for the least positive solution of f(t), it suffices to find an upper bound for the largest solution of f ± (t), and in turn it suffices to bound the solutions of the general form g(t). Assuming that at least one solution of g(t) exists in the interval (0, ∞), we give elementary estimates for the upper bound in the following lemma. Lemma 4.11. If Z(g) > 0, there is no solution of g(t) in the region t ≥ t0 , where

t0 =

ln((|a2 | + |a3 | )/|a1 | ) . ln(κ1 )

Proof. We treat separately the following two cases: 1. σ (g) = 1: we have either (a) a1 , a2 > 0, a3 < 0, or (b) a1 > 0, a2 , a3 < 0, since other cases can be reduced to these two. 2. σ (g) = 2: a1 , a3 > 0, a2 < 0. Case 1 (a): Since we have assumed that at least one solution exists, we must have |a3 | > a1 . In this case, we also have |a3 | > |a3 | − a2 κ2t . Therefore there is no solution in the region t > t0a , ln(|a3 |/a1 ) , as in ln(κ1 ) a1 , t0a is positive.

where t0a =

this region |a3 | ≤ a1 κ1t . Since κ 1 > 1

and |a3 | > Case 1 (b): One can use similar arguments as in the previous ln((|a2 |+|a3 | )/a1 ) case to show that there is no solution for t > t0b = . ln(κ ) 1

Case (2): Since a2 < 0 and κ 1 < 1, we have g (t ) = a1 ln(κ1 )κ1t + a2 ln(κ2 )κ2t > 0 for all t > 0. Therefore g can have at most 1 positive solution. Since σ (g) = 2, we must have Z (g) = 0 by 4.7, which contradicts our assumption that Z(g) > 0. Therefore the case σ (g) = 2 cannot arise for g(t). Therefore we get that there is no solution for g(t) for t > ln((|a2 |+|a3 | )/|a1 | ) max{t0a , t0b } = t0b =  ln(κ ) 1

Using elementary arguments, one can also give upper bounds for solutions of the border return time Eq. (4.9) in many more cases, e.g., when all three eigenvalues are real. However, we do not claim that any of these upper bounds are tight. We also note that in case g(t) has exactly one solution, the least positive solution can also be bounded below by the (unique) extremum point t∗ of g(t) given by

t∗ =

ln(−a2 ln(κ2 )/a1 ln(κ1 )) ln(κ1 /κ2 )

Note that this may or may not be a positive number. Now we can summarize the discussions above and state a necessary condition for the occurence of chaos through transverse homoclinic intersection on the first return of an unstable manifold Ur for R∗ , following a Shilnikov-type dynamical behaviour. The condition is also sufficient for the existence of homoclinic intersections, if one drops the assumption that the intersection takes place on the first return. Theorem 4.12. Suppose that the left-side fixed point L∗ is admissible and has eigenvalues (λ1 , λ2 , λ3 ) where λ1 is real with λ1 > 1, and λ2 = λ3 is a pair of complex-conjugate eigenvalues with absolute value r0 < 1. Assume also that the right-side fixed-point R∗ is admissible with one unstable flip eigenvalue and two stable eigenvalues;

denote the associated unstable manifold by Ur and the stable manifold as Sr . Take X0 = (0, y0 , z0 ) to be the intersection of Ur with the border x = 0 with y0 + μ < 0, and consider the corresponding function f(t) of Eq. (4.9) (see also Remark (4.6)), with upper and lower bounding curves f ± (t), with either f + (t ) or f − (t ) having at least one solution (say f + (t )). Then, a necessary condition for a transverse homoclinic intersection on first return to occur between Ur and Sr , is given by:

((Pk − R∗ ) · nSr ).((Pk+1 − R∗ ) · nSr ) < 0 for some k ≤ t0 + 2 ln((2|α2 | + |α3 | )/|α1 | ) wheret0 = . ln(λ1 ) Here Pk is the kth iterate of X0 , nSr is the normal vector to the 2dimensional stable eigenplane of R∗ and v · w denotes the dot product of two vectors v, w ∈ R3 (v is the transpose of v). 4.2.3. Case B: one stable negative eigenvalue and a pair of unstable complex conjugate eigenvalues: We now turn to another case where one of the eigenvalues λ1 is negative with absolute value strictly less than 1, and a pair of complex conjugate eigenvalues λ2 , λ3 with absolute value strictly greater than 1, i.e. we have λ1 < 1, and λ2 = λ3 = r0 eiθ0 , r0 > 0, θ0 ∈ (−π , π ] are a pair of complex eigenvalues with r0 > 1. Eq. (4.9) then takes the following general form:

f (t ) := α1 cos(t π )λt1 + 2α2 cos(t θ0 + δ )r0t + α3 = 0,

t≥0 (4.13)

where α1 , α2 , α3 , δ ∈ R. Again our objective is to give upper bounds for the smallest positive root of f(t). We shall only treat the following special case;

α1 < 0, α2 > 0, α3 > 0, θ0 ∈ (0, π /2 ),

f (1 ) > 0

( ∗ ).

Note that due to the oscillatory coefficient for r0 and the fact that r0 > 1, the function f(t) has infinitely many zeroes in the positive real line. To estimate the least positive zero, we consider the two bounding curves of f(t) as before:

α1 cos(t π )λt1 + 2α2 r0t + α3 , f (t ) = α1 cos(t π )λt1 − 2α2 r0t + α3 f + (t ) =

and

−

(4.14)

To deduce an upper bound for the least positive solution of f(t), we shall use the geometry of the curve f(t) and the bounding curves f ± (t). The essential idea is that since f − (t ) is strictly negative for t large enough due to the term involving λ1 dying out, if we consider a point of intersection of f(t) and f − (t ) for t = t0 large enough, the points f (t0 + s ), 0 ≤ s ≤ 1 will also lie below the xaxis. This is due to the assumption that θ 0 lies between 0 and π /2, which means that the term cos((t0 + s )θ0 + δ ) does not change sign and is strictly increasing for s ∈ [0, 1], as at the intersection point we must have t0 θ0 + δ = (2k + 1 )π for some k ∈ N. This also rules out the situation where there is a positive local maximum of f(t) between t0 and t0 + 1, because the curve f + (t ) must intersect f(t) before the local maxima, but at such a point cos((t0 + s )θ0 + δ ) must take the value +1 (see Fig. 7 for an example). This argument is formalized in the following lemma. Lemma 4.13. With the assumptions in (∗ ), we have: t

1. if r00 > α3 /2α2 for some t0 ∈ (0, ∞), then f − (t0 ) < 0. 2. Suppose that t0 ∈ (1, ∞) is such that t0 θ0 + δ = (2k + 1 )π for some k ∈ N, and

r0t0 >

|α1 |λ1 + α3 . 2α2 cos(θ0 )

Then we have,

f (t0 + s ) < 0

for all s ∈ [0, 1].

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

9

Proof. The first statement is an immediate consequence of ∗ ∗ Lemma 4.13, since we have   chosen t such that f(t ) < 0. On the other hand, since n0 = t ∗

lies between t∗ and t ∗ + 1, from part

(2) of Lemma 4.13, we see that f(n0 ) < 0. By assumption (∗ ), f(1) is positive, hence the sequence {f(n)} has changed sign at least once.  Therefore Lemma 4.14 not only gives an estimate for upper bound of the least positive solution of f(t), but also that of sign change in the sequence {f(n)}, which is what we are really after since we have a discrete dynamical system. Observe that even if f(t) crosses the x-axis at some point t = T , due to the oscillatory nature of f(t) with increasing amplitudes as t increases,   the sequence {f(n)} may not change sign between T and

T . However,

making the careful choice T = t ∗ forces the sequence to change  sign before T .

Fig. 7. The plot of the functions f(t) (blue curve), f + (t ) (red curve), f − (t ) (green curve), and the threshold value t∗ (shown by a red dot on the x-axis). The green diamonds are x-coordinates of the actual orbit points for the dynamical system. Note that since the interpolation is only valid until the orbit crosses the border, the orbit points lie on f(t) only until the first crossing of the x-axis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Proof. To prove (1), we note that since cos(t0 π ) ≥ −1, we have:

f − (t0 ) = α1 cos(t0 π )λt10 − 2α2 r0t0 + α3 ≤ α1 λt10 − 2α2 r0t0 + α3 t

t

Since we have r00 > α3 /2α2 , or equivalently −2α2 r00 + α3 < 0, and since α 1 < 0 from the assumption (∗ ), we get the result. To prove (2), we use the fact that cos (sθ 0 ) ≥ cos (θ 0 ) and r0s ≥ 1∀s ∈ [0, 1]. We get,

α1 cos((t0 + s )π )λt0 +s + 2α2 cos(t0 θ0 + δ + sθ0 )r0t0 +s + α3 < α1 cos((t0 + s )π )λt0 +s − 2α2 cos(sθ0 )r0t0 +s + α3 < α1 cos((t0 + s )π )λt0 +s r s (|α1 |λ1 + α3 ) − 2α2 cos(sθ0 ) 0 + α3 2α2 cos(θ0 ) ≤ α1 cos((t0 + s )π )λt0 +s − (|α1 |λ1 + α3 ) + α3 = α1 (cos((t0 + s )π )λt10 +s + λ1 )

f (t0 + s ) =

< 0, since t0 > 1, λ1 < 1 and α 1 < 0.



Let us consider the region where both the conditions in Lemma 4.13 are satisfied, i.e. t > t1 where,



t1 = max





1 1 α3 ln , ln ln(r0 ) 2α2 ln(r0 )



  |α1 |λ1 + α3 ,1 , 2α2 cos(θ0 )

and suppose that k∗ is the positive odd integer given in terms of t1 and the known parameters as follows:

k∗ =

   θ0t1 + δ  odd    π 

Now we can state a theorem similar to that of Theorem 4.12 concerning the Shilnikov-type dynamics in this case. Theorem 4.15. Suppose that the left-side fixed point L∗ is admissible and has eigenvalues (λl1 , λl2 , λl3 ) where λ1 is real with λl1 > 1,

and λl2 = λl 3 is a pair of complex-conjugate eigenvalues with absolute value r0l < 1. Assume also that the right-side fixed-point R∗ is admissible with one stable flip eigenvalue −1 < λr1 < 0 and a pair of a

pair of unstable complex-conjugate eigenvalues λr2 = λr 3 with absolute value strictly greater than 1. Denote the unstable manifold for L∗ by Ul and the stable manifold of L∗ as Sl . Take X0 = (0, y0 , z0 ) to be the intersection of Ul with the border x = 0 with y0 + μ > 0, and consider the corresponding function f(t) of Eq. (4.13). Then, a necessary condition for a transverse homoclinic intersection on first return to occur between Ul and Sl , is given by:

 

((Pk − L∗ ) · nSl ).((Pk+1 − L∗ ) · nSl ) < 0

for some k ≤ t ∗ ,

where t∗ is determined by Lemma 4.14. As before, Pk is the kth iterate of X0 , nSl is the normal vector to the 2-dimensional stable eigenplane of L∗ . 4.2.4. Discussion of the motivating example Recalling the motivating example from Section 3, we see that it satisfies the assumptions for the left and right fixed-points as well as the assumptions for the eigenvalues in Theorem 4.15. Moreover, we have the following values for the various parameters in Eq. 4.13:

α1 = −0.0633, α2 = 0.285, α3 = 0.483, θ0 = 0.968, δ = −2.396. The

border

intersection

point

of

Ul

is

given

by

(0, −0.364, −0.141 ), thus we also have y0 + μ > 0, which im-

plies f(1) > 0. Hence all the assumptions in (∗ ) are satisfied. The threshold value t∗ is computed to be t ∗ = 5.716 which turns out to be a good estimate for the least positive solution in this case. See Fig. 7 which shows the functions f(t), f ± (t), as well the actual orbit points.

 odd

where x

denotes the least positive odd integer greater than or

equal to x. The choice of

t θ0 + δ = k ∗ π ,

k∗

provides the solution of the equation

t > t1

Lemma 4.14. The least positive solution of f(t) must be strictly less ∗ −δ than t ∗ = k π . Moreover, the sequence { f (n )}n∈N must change sign θ 0

 

at least once, on or before n0 = t ∗ .

5. Conclusion In the continuous case, Shilnikov bifurcation occurs because the 1-dimensional unstable manifold loops back and intersects the 2dimensional stable manifold. We have shown that a similar phenomenon occurs in a discrete three-dimensional piecewise smooth system, albeit the looping of the 1-dimensional unstable manifold occurs due to the nature of the fixed point on the other side of

10

I. Roy, M. Patra and S. Banerjee / Chaos, Solitons and Fractals 133 (2020) 109655

the border. Therefore we call this Shilnikov-type dynamics and the resulting chaos as Shilnikov-type chaos. In this paper we have also derived analytical conditions for the occurrence of a transverse homoclinic intersection, and therefore the occurrence of a chaotic dynamics. It is well-known that such transverse homoclinic intersections cause chaotic behaviour via the presence of Smale horseshoes. More precisely, we employ two analytical methods to give necessary conditions for the occurrence of a homoclinic orbit: one that uses recursion and another that uses complex interpolation. The two methods are complementary to each other and are meant to be used together in practice. As far as we know, this method seems to be new in the case of discrete piecewise smooth systems. In particular, since the interpolation technique provides an explicitly-defined continuous curve that can be used to probe the dynamics of the system, we expect that it will be useful to study other dynamical properties of such systems as well. We also present an example, which shows a ‘two-sided Shilnikov dynamics’, where looping behaviour of 1-dimensional eigenmanifolds occurs for both fixed points on either side of the border, which then intersect transversally their respective 2dimensional eigenmanifolds. The resulting chaotic dynamics is Shilnikov-type on either side, which can be therefore aptly named as ‘two-sided Shilnikov chaos’. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Indrava Roy: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Resources, Data curation, Writing - original draft, Writing - review & editing, Visualization, Supervision, Project administration, Funding acquisition. Mahashweta Patra: Methodology, Methodology, Software, Validation, Formal analysis, Investigation, Resources, Data curation, Writing - original draft, Writing - review & editing, Visualization. Soumitro Banerjee: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Resources, Writing - review & editing, Supervision, Project administration, Funding acquisition.

Acknowledgments IR thanks the Indian Science and Engineering Research Board (SERB) for support via MATRICS project MTR/2017/0 0 0835. SB acknowledges financial support in the form of J C Bose Fellowship by the SERB, Govt. of India, no. SB/S2/JCB-023/2015. IR thanks A. Prasad for reading a preliminary draft of the paper. References [1] Arecchi F, Meucci R, Gadomski W. Laser dynamics with competing instabilities. Phys Rev Lett 1987;58(21):2205. [2] Avrutin V, Zhusubaliyev ZT, Saha A, Banerjee S, Sushko I, Gardini L. Dangerous bifurcations revisited. Int J Bifurcation Chaos 2016;26(14):1630040. [3] Braun T, Lisboa JA, Gallas JA. Evidence of homoclinic chaos in the plasma of a glow discharge. Phys Rev Lett 1992;68(18):2770. [4] De S, Dutta PS, Banerjee S, Roy AR. Local and global bifurcations in three-dimensioanl, continuous, piecewise smooth maps. Int J Bifurcation Chaos 2011;21:1617–36. [5] de Tomasi F, Hennequin D, Zambon B, Arimondo E. Instabilities and chaos in an infrared laser with saturable absorber: experiments and vibrorotational model. JOSA B 1989;6(1):45–57. [6] di Bernardo M. Normal forms of border collisions in high-dimensional nonsmooth maps. In: Circuits and systems, 2003. ISCAS’03. Proceedings of the 2003 international symposium on, volume 3, pages III–III. IEEE; 2003. [7] Hsieh P.-F., Sibuya Y.. Fundamental theorems of ordinary differential equations. 1999. [8] Jameson GJ. Counting zeros of generalised polynomials: Descartes’ rule of signs and Laguerre’s extensions. Math Gaz 2006;90(518):223–34. [9] Patra M, Banerjee S. Bifurcation of quasiperiodic orbit in a 3d piecewise linear map. Int J Bifurcation Chaos 2017;27(10):1730033. [10] Pisarchik A, Meucci R, Arecchi F. Discrete homoclinic orbits in a laser with feedback. Phys Rev E 20 0 0;62(6):8823. [11] Pisarchik A, Meucci R, Arecchi F. Theoretical and experimental study of discrete behavior of Shilnikov chaos in a Co2 laser. Eur Phys J D 2001;13(3):385–91. [12] Roy I, Roy AR. Border collision bifurcations in three-dimensional piecewise smooth systems. Int J Bifurcation Chaos 2008;18:577–86. [13] Saha A., Banerjee S. Existence and stability of periodic orbits in n-dimensional piecewise linear continuous maps. 2015. arXiv:1504.01899. [14] Shilnikov LP. A case of the existence of a denumerable set of periodic motions. Dokl Akad Nauk SSSR 1965;160:558–61. [15] Silva CP. Shil’nikov’s theorem-a tutorial. IEEE Trans Circuits Syst I 1993;40(10):675–82. [16] Tossavainen T. The lost cousin of the fundamental theorem of algebra. Math Mag 2007;80(4):290–4. [17] Tresser C. About some theorems by lp šil’nikov. Ann Inst H Poincaré Phys Théor 1984;40(4):441–61. [18] Viktorov EA, Klemer DR, Karim MA. Shil’nikov case of antiphase dynamics in a multimode laser. Opt Commun 1995;113(4–6):441–8. [19] Zhou T, Chen G, Yang Q. Constructing a new chaotic system based on the Silnikov criterion. Chaos Solitons Fractals 2004;19(4):985–93. [20] Gonchenko A, Gonchenko S, Kazakov A, Turaev D. Simple scenarios of onset of chaos in three-dimensional maps. Int J Bifurcation Chaos 2014;24(8):1440 0 05.