# Dynamical Behaviors of Coupled Memristor-Based Oscillators with Identical and Different Nonlinearities

Authors: A. Elsonbaty; A. Abdelkhalek; A. Elsaid
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In this work, the interesting dynamics of coupled nonlinear memristor-based oscillators in ring configuration are explored. The mathematical models are derived to describe the possible cases of employing identical or different nonlinearities. Analytical and numerical techniques, involving perturbation methods, normal forms, phase portraits, and Lyapunov exponents are used to investigate various types of dynamical behaviors along with their stability regions in parameters space. The effects of time-delayed coupling on the proposed system are numerically studied. It is demonstrated that the coupled oscillators show rich dynamics including periodic orbits, quasiperiodicity, two-dimensional, and three-dimensional tori.


---

## Body

## 1. Introduction

Many applications in engineering, physics, chemistry, biology, and economy are modeled using nonlinear differential equations (DEs) where their exact forms of solutions are very difficult to obtain in general or they do not provide tools to understand the underlying dynamics. However, dynamical systems methods provide powerful techniques to study the qualitative behaviors of mathematical model solution rather than quantitative ones [1].One of the interesting problems in nonlinear theory is to study the dynamics of systems involving coupled nonlinear oscillators. This problem appears in many applications in biology [2], neuroscience [3–6], chemistry [7, 8], and physics [9–11]. Researchers in the field of dynamical systems are interested in studying the stability of coupled oscillators systems with nonlinear properties since 1950s till now [12–15].In 1970, Chua introduced the fourth electronic element called the memristor [16]. It is a two-terminal element that relates the magnetic flux ϕ between the terminals with the electric charge q that passes through it [16, 17]. Studies on the dynamical properties, i.e., stability and bifurcation of the memristor-based oscillator circuit, have received increasing interest [18–21]. The main reason is that coupled memristor-based oscillatory circuits appear in many applications such as filtering, energy storage, energy transmission, isolation impedance transformation, learning networks that require a synapse-like function, secure communications, and image stabilization [22–31].In [32], the stability of three nonlinear identical memristor-based oscillators coupled by inductors in ladder circuit configuration was studied. In this paper, we examine a more generalized case by considering coupled three memristor-based oscillators in ring configuration where nonidentical memristors are allowed. Two cases are studied where both identical and different nonlinearities are considered. The memristor used is of a flux-controlled type that has memductance function w(ϕ) representing the flux-dependent rate of change of charge (1)wϕ=dqϕdϕ.The two cases studied are according to the definition of q(ϕ). In the first case the memristors are described by three third-order polynomials, whereas in the second case one memristor is described by third-order polynomial while the remaining memristors are characterized by fourth-order polynomials.Time delays play an important role in mathematical modeling of real systems and studying synchronization such that their effects cannot be neglected in real life; see, for example, [33–35]. Therefore, numerical investigations of time-delayed coupling influences on coupled oscillators dynamics are carried out in this work. The rest of the paper is organized as follows: In Section 2, mathematical models of all cases studied in this work are presented. In Section 3, the multiple scales method is applied to acquire the normal form for each system. In Section 4, various types of dynamical behaviors that can be extracted from normal forms in addition to their stability analysis are obtained. In Section 5, numerical simulations are carried out to verify theoretical results. The general discussion and main conclusion of this work are presented in Section 6.


## 2. Ring Model Equations

In this section, we derive the fundamental equations describing the three ring coupled oscillators circuit. Figure1 shows the memristor-based oscillator circuit we adopt in this work. Applying Kirchhoff’s voltage law (KVL) and Kirchhoff’s current law (KCL) to the circuit in Figure 2, the model can be described as follows:(2)I1-I3+C1V1˙+g1V1+IL1+Iw1V1=0,I2-I1+C2V2˙+g2V2+IL2+Iw2V2=0,I3-I2+C3V3˙+g3V3+IL3+Iw3V3=0,I1˙-V1-V2L0=0,I2˙-V2-V3L0=0,I3˙-V3-V1L0=0,IL1˙-V1L1=0,IL2˙-V2L2=0,IL3˙-V3L3=0.We consider two cases of system (2) according to the equations describing memristors.Figure 1
Theith memristor-based oscillator circuit in coupled network.Figure 2
Ring type schematic of three coupled memristor-based oscillators where each oscillator is referred to as osci, i=1,2,3.


### 2.1. Case 1

In this case, the nonlinearities of the three memristors are described as follows:(3)qwiϕi=miϕi3-niϕi,mi,ni>0,i=1,2,3.By integrating both sides of each equation of system (2) between times t0 and t while using (3) and the relations(4)qt=∫-∞tIudu,(5)ϕt=∫-∞tVudu,we get(6)
ϕ
1
˙
+
g
1
ϕ
1
C
1
+
q
L
1
C
1
+
m
1
ϕ
1
3
C
1
-
n
1
ϕ
1
C
1
+
q
1
C
1
-
q
3
C
1
=
0
,
(7)
ϕ
2
˙
+
g
2
ϕ
2
C
2
+
q
L
2
C
2
+
m
2
ϕ
2
3
C
2
-
n
2
ϕ
2
C
2
-
q
1
C
2
+
q
2
C
2
=
0
,
(8)
ϕ
3
˙
+
g
3
ϕ
3
C
3
+
q
L
3
C
3
+
m
3
ϕ
3
3
C
3
-
n
3
ϕ
3
C
3
+
q
3
C
3
-
q
2
C
3
=
0
,
(9)
q
1
˙
-
ϕ
1
-
ϕ
2
L
0
=
0
,
(10)
q
2
˙
-
ϕ
2
-
ϕ
3
L
0
=
0
,
(11)
q
3
˙
-
ϕ
3
-
ϕ
1
L
0
=
0
,
(12)
q
L
1
˙
-
ϕ
1
L
1
=
0
,
(13)
q
L
2
˙
-
ϕ
2
L
2
=
0
,
(14)
q
L
3
˙
-
ϕ
3
L
3
=
0
.By differentiating (6)–(8) and using (9)–(14), we obtain the following equations that depict our intended case:(15)ϕj¨-ϕj˙nj-gjCj-3mjϕj2Cj+ϕj1CjLj+2CjL0-∑k=1,k≠j3ϕkCjL0=0,j=1,2,3.


### 2.2. Case 2

In this case, the nonlinearities of the three memristors are different and described as follows:(16)qw1ϕ1=m1ϕ13-n1ϕ1,qw2ϕ2=h2ϕ24+m2ϕ23-n2ϕ2,qw3ϕ3=h3ϕ34-n3ϕ3.mi,hi+1,nj>0,i=1,2,j=1,2,3.Using the same procedure as in Section 2.1, we obtain(17)ϕ1¨-ϕ1˙n1-g1C1-3m1ϕ12C1+ϕ11C1L1+2C1L0-ϕ2+ϕ3C1L0=0,ϕ2¨-ϕ2˙n2-g2C2-3m2ϕ22C2-4h2ϕ23C2+ϕ21C2L2+2C2L0-ϕ1+ϕ3C2L0=0,ϕ3¨-ϕ3˙n3-g3C3-4h3ϕ33C3+ϕ31C3L3+2C3L0-ϕ2+ϕ1C3L0=0.


## 2.1. Case 1

In this case, the nonlinearities of the three memristors are described as follows:(3)qwiϕi=miϕi3-niϕi,mi,ni>0,i=1,2,3.By integrating both sides of each equation of system (2) between times t0 and t while using (3) and the relations(4)qt=∫-∞tIudu,(5)ϕt=∫-∞tVudu,we get(6)
ϕ
1
˙
+
g
1
ϕ
1
C
1
+
q
L
1
C
1
+
m
1
ϕ
1
3
C
1
-
n
1
ϕ
1
C
1
+
q
1
C
1
-
q
3
C
1
=
0
,
(7)
ϕ
2
˙
+
g
2
ϕ
2
C
2
+
q
L
2
C
2
+
m
2
ϕ
2
3
C
2
-
n
2
ϕ
2
C
2
-
q
1
C
2
+
q
2
C
2
=
0
,
(8)
ϕ
3
˙
+
g
3
ϕ
3
C
3
+
q
L
3
C
3
+
m
3
ϕ
3
3
C
3
-
n
3
ϕ
3
C
3
+
q
3
C
3
-
q
2
C
3
=
0
,
(9)
q
1
˙
-
ϕ
1
-
ϕ
2
L
0
=
0
,
(10)
q
2
˙
-
ϕ
2
-
ϕ
3
L
0
=
0
,
(11)
q
3
˙
-
ϕ
3
-
ϕ
1
L
0
=
0
,
(12)
q
L
1
˙
-
ϕ
1
L
1
=
0
,
(13)
q
L
2
˙
-
ϕ
2
L
2
=
0
,
(14)
q
L
3
˙
-
ϕ
3
L
3
=
0
.By differentiating (6)–(8) and using (9)–(14), we obtain the following equations that depict our intended case:(15)ϕj¨-ϕj˙nj-gjCj-3mjϕj2Cj+ϕj1CjLj+2CjL0-∑k=1,k≠j3ϕkCjL0=0,j=1,2,3.


## 2.2. Case 2

In this case, the nonlinearities of the three memristors are different and described as follows:(16)qw1ϕ1=m1ϕ13-n1ϕ1,qw2ϕ2=h2ϕ24+m2ϕ23-n2ϕ2,qw3ϕ3=h3ϕ34-n3ϕ3.mi,hi+1,nj>0,i=1,2,j=1,2,3.Using the same procedure as in Section 2.1, we obtain(17)ϕ1¨-ϕ1˙n1-g1C1-3m1ϕ12C1+ϕ11C1L1+2C1L0-ϕ2+ϕ3C1L0=0,ϕ2¨-ϕ2˙n2-g2C2-3m2ϕ22C2-4h2ϕ23C2+ϕ21C2L2+2C2L0-ϕ1+ϕ3C2L0=0,ϕ3¨-ϕ3˙n3-g3C3-4h3ϕ33C3+ϕ31C3L3+2C3L0-ϕ2+ϕ1C3L0=0.


## 3. Normal Forms

The purpose of this section is to acquire the normal forms corresponding to each case of ring model derived in previous section. The powerful perturbation techniques are used to achieve this goal. First, systems (15) and (17) are written in perturbed form; then each system is reduced to an equivalent first-order system [36].System (15) can be written in perturbed form as follows:(18)xj¨-εxj˙aj-αjxj2+βjxj-bjε∑k=1,k≠j3xk=0,j=1,2,3.And system (17) can be written in perturbed form as follows:(19)x1¨-εx1˙a1-α1x12+β1x1-b1εx2+x3=0,x2¨-εx2˙a2-α3x23-α2x22+β2x2-b2εx1+x3=0,x3¨-εx3˙a3-α4x33+β3x3-b3εx1+x2=0.By using the substitutions(20)xjt=yj¯t+yjt,xj˙t=iβjyjt-yjt¯,j=1,2,3,where yj¯(t) is the complex conjugate of yj(t), systems (18) and (19) are reduced to the more simplified first-order systems of DEs. More specifically, system (18) is reduced to the following system:(21)yj˙-12εyj-yj¯aj-αjyj¯+yj2+ibjε∑k=1,k≠j3yk0¯+yk02βj-iβjyj=0,j=1,2,3.And system (19) is reduced to the following system:(22)y1˙-12εy1-y1¯a1-α1y1¯+y12+ib1εy2¯+y3¯+y2+y32β1-iβ1y1=0,y2˙-12εy2-y2¯-α3y2¯+y23-α2y2¯+y22+a2+ib2εy1¯+y3¯+y1+y32β2-iβ2y2=0,y3˙-12εy3-y3¯a3-α4y3¯+y33+ib3εy1¯+y2¯+y1+y22β3-iβ3y3=0.The normal forms corresponding to systems (21) and (22) are obtained by applying the method of multiple scales [36].


### 3.1. The Normal Form of Case 1 of Ring Model

Using the substitutions(23)yjt;ε=∑k=02εkyjkT0,T1,T2,Tj-1=εj-1t,ddt=∑k=02Dkεk,Dj-1=∂Tj-1,j=1,2,3,in system (21) and equating coefficients of equal powers of ε, we obtain(Orderε0)(24)D0yj0-iβjyj0=0,j=1,2,3.(Orderε1)(25)D0yj1-iβjyj1=-D1yj0+-ibj∑k=1,k≠j3yk0¯+yk02βj+12yj0¯-yj0αjyj0¯+yj02-aj,j=1,2,3.(Orderε2)(26)
D
0
y
12
-
i
β
1
y
12
=
-
i
b
1
y
21
¯
+
y
31
¯
+
y
21
+
y
31
2
β
1
-
D
2
y
10
-
D
1
y
11
+
1
2
y
11
α
1
y
10
¯
2
-
2
y
10
y
10
¯
-
3
y
10
2
+
a
1
-
y
11
¯
α
1
-
3
y
10
¯
2
-
2
y
10
y
10
¯
+
y
10
2
+
a
1
,
D
0
y
22
-
i
β
2
y
22
=
-
i
b
2
y
11
¯
+
y
31
¯
+
y
11
+
y
31
2
β
2
-
D
2
y
20
-
D
1
y
21
+
1
2
y
21
α
2
y
20
¯
2
-
2
y
20
y
20
¯
-
3
y
20
2
+
a
2
-
y
21
¯
α
2
-
3
y
20
¯
2
-
2
y
20
y
20
¯
+
y
20
2
+
a
2
,
D
0
y
32
-
i
β
3
y
32
=
-
i
b
3
y
11
¯
+
y
21
¯
+
y
11
+
y
21
2
β
3
-
D
2
y
30
-
D
1
y
31
+
1
2
y
31
α
3
y
30
¯
2
-
2
y
30
y
30
¯
-
3
y
30
2
+
a
3
-
y
31
¯
α
3
-
3
y
30
¯
2
-
2
y
30
y
30
¯
+
y
30
2
+
a
3
.Solving system (24) of DEs we get(27)yj0=eiβjT0Aj0T1,T2,j=1,2,3.By substituting these solutions in (25) we get the conditions needed to remove the secular terms in the following forms:(28)D1Aj0=ajAj02-12αjAj02Aj0¯,j=1,2,3.Using (28) to solve system (25) and substituting from the results in (26), we get the following conditions to remove the secular terms(29)D2A10=-iA10β1-β2β1-β3-12α1a1A10A10¯+11α12A102A10¯2+2a12-8b1b3β1-β2+b2β1-β316β1β1-β2β1-β3,D2A20=iA20β1-β2β2-β3-12α2a2A20A20¯+11α22A202A20¯2+2a22-8b2b3β1-β2+b1β3-β216β2β2-β1β2-β3,D2A30=-iA30β2-β3β1-β3-12α3a3A30A30¯+11α32A302A30¯2+2a32+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.Substituting from (28) and (29) in Aj0˙=D2ε2Aj0+D1εAj0, where j=1,2 and 3, we obtain the normal form of system (21) which can be expressed in the following form:(30)rj0˙=12εrj0aj-αjrj02,j=1,2,3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.


### 3.2. The Normal Form of Case 2 of Ring Model

Using the same procedure applied in Section3.1, the following conditions are obtained to remove the secular terms of order ε1(31)D1A10=a1A102-12α1A102A10¯,D1A20=a2A202-12α2A202A20¯,D1A30=a3A302,and the conditions needed to remove the secular terms of order ε2 are given by(32)D2A10=-iA10β1-β2-12α1a1A10A10¯+11α12A102A10¯2+2a12-8b1b216β1β1-β2,D2A20=iA20β1-β2-60α2a2A20A20¯+A202A20¯256α32A20A20¯+55α22+10a22+40b1b280β2β2-β1,D2A30=iA30β2-β3-a3β1-β3a3-α3A30A30¯-4b1b3+4b2b3β3-β18β3β3-β1β3-β2.So, we get the following normal form:(33)r10˙=12εr10a1-α1r102,r20˙=12εr20a2-α2r202,r30˙=12a3εr30,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,θ30˙=-ε2β2-β3a3β1-β3a3-α3r302+4b1b3+4b2b3β1-β38β3β3-β1β3-β2.


## 3.1. The Normal Form of Case 1 of Ring Model

Using the substitutions(23)yjt;ε=∑k=02εkyjkT0,T1,T2,Tj-1=εj-1t,ddt=∑k=02Dkεk,Dj-1=∂Tj-1,j=1,2,3,in system (21) and equating coefficients of equal powers of ε, we obtain(Orderε0)(24)D0yj0-iβjyj0=0,j=1,2,3.(Orderε1)(25)D0yj1-iβjyj1=-D1yj0+-ibj∑k=1,k≠j3yk0¯+yk02βj+12yj0¯-yj0αjyj0¯+yj02-aj,j=1,2,3.(Orderε2)(26)
D
0
y
12
-
i
β
1
y
12
=
-
i
b
1
y
21
¯
+
y
31
¯
+
y
21
+
y
31
2
β
1
-
D
2
y
10
-
D
1
y
11
+
1
2
y
11
α
1
y
10
¯
2
-
2
y
10
y
10
¯
-
3
y
10
2
+
a
1
-
y
11
¯
α
1
-
3
y
10
¯
2
-
2
y
10
y
10
¯
+
y
10
2
+
a
1
,
D
0
y
22
-
i
β
2
y
22
=
-
i
b
2
y
11
¯
+
y
31
¯
+
y
11
+
y
31
2
β
2
-
D
2
y
20
-
D
1
y
21
+
1
2
y
21
α
2
y
20
¯
2
-
2
y
20
y
20
¯
-
3
y
20
2
+
a
2
-
y
21
¯
α
2
-
3
y
20
¯
2
-
2
y
20
y
20
¯
+
y
20
2
+
a
2
,
D
0
y
32
-
i
β
3
y
32
=
-
i
b
3
y
11
¯
+
y
21
¯
+
y
11
+
y
21
2
β
3
-
D
2
y
30
-
D
1
y
31
+
1
2
y
31
α
3
y
30
¯
2
-
2
y
30
y
30
¯
-
3
y
30
2
+
a
3
-
y
31
¯
α
3
-
3
y
30
¯
2
-
2
y
30
y
30
¯
+
y
30
2
+
a
3
.Solving system (24) of DEs we get(27)yj0=eiβjT0Aj0T1,T2,j=1,2,3.By substituting these solutions in (25) we get the conditions needed to remove the secular terms in the following forms:(28)D1Aj0=ajAj02-12αjAj02Aj0¯,j=1,2,3.Using (28) to solve system (25) and substituting from the results in (26), we get the following conditions to remove the secular terms(29)D2A10=-iA10β1-β2β1-β3-12α1a1A10A10¯+11α12A102A10¯2+2a12-8b1b3β1-β2+b2β1-β316β1β1-β2β1-β3,D2A20=iA20β1-β2β2-β3-12α2a2A20A20¯+11α22A202A20¯2+2a22-8b2b3β1-β2+b1β3-β216β2β2-β1β2-β3,D2A30=-iA30β2-β3β1-β3-12α3a3A30A30¯+11α32A302A30¯2+2a32+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.Substituting from (28) and (29) in Aj0˙=D2ε2Aj0+D1εAj0, where j=1,2 and 3, we obtain the normal form of system (21) which can be expressed in the following form:(30)rj0˙=12εrj0aj-αjrj02,j=1,2,3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.


## 3.2. The Normal Form of Case 2 of Ring Model

Using the same procedure applied in Section3.1, the following conditions are obtained to remove the secular terms of order ε1(31)D1A10=a1A102-12α1A102A10¯,D1A20=a2A202-12α2A202A20¯,D1A30=a3A302,and the conditions needed to remove the secular terms of order ε2 are given by(32)D2A10=-iA10β1-β2-12α1a1A10A10¯+11α12A102A10¯2+2a12-8b1b216β1β1-β2,D2A20=iA20β1-β2-60α2a2A20A20¯+A202A20¯256α32A20A20¯+55α22+10a22+40b1b280β2β2-β1,D2A30=iA30β2-β3-a3β1-β3a3-α3A30A30¯-4b1b3+4b2b3β3-β18β3β3-β1β3-β2.So, we get the following normal form:(33)r10˙=12εr10a1-α1r102,r20˙=12εr20a2-α2r202,r30˙=12a3εr30,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,θ30˙=-ε2β2-β3a3β1-β3a3-α3r302+4b1b3+4b2b3β1-β38β3β3-β1β3-β2.


## 4. Stability Analysis of Equilibrium Solutions

In this section, the stability analysis of the equilibrium solutions is carried out using the normal forms obtained in (30) and (33).


### 4.1. Stability Analysis of Case 1 of Ring Model

The Jacobian matrix of system (30) is expressed as(34)J=12εa1-r102α1-εr102α100012εa2-r202α2-εr202α200012εa3-r302α3-εr302α3.System (30) has the following equilibrium solutions.


#### 4.1.1. Equilibrium Point

The following equilibrium solution,(35)ri0=0,i=1,2,3,is corresponding to an equilibrium point at origin. Evaluating the eigenvalues of Jacobian matrix (34) at equilibrium point (35), we get the stability region of this point in parameters space as follows:(36)a1<0,a2<0,a3<0.


#### 4.1.2. Periodic Solution I

The second equilibrium solution,(37)r10=0,r20=0,r30=a3α3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a periodic solution in the phase space of case 1 of ladder model. The second order approximate solution of this periodic orbit has an amplitude given by a3/α3 and an angular frequency computed as(38)β3-ε2β2-β3a32β1-β3+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution I (37), we get its stability region in parameters space as follows:(39)a1<0,a2<0,a3>0.


#### 4.1.3. Periodic Solution II

The third equilibrium solution,(40)r10=0,r20=a2α2,r30=0,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,represents a periodic solution in the phase space of case 1 of ring model. The second order approximate solution of this periodic orbit has an amplitude given by a2/α2 and an angular frequency computed as(41)β2-ε28b2b3β1-β2+b1β3-β2-a22β1-β2β2-β316β2β2-β1β2-β3.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution II (40), we get its stability region in parameters space as follows:(42)a1<0,a2>0,a3<0.


#### 4.1.4. Periodic Solution III

The fourth equilibrium solution,(43)r10=a1α1,r20=0,r30=0,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,corresponds to a periodic solution in the phase space of case 1 of ring model. The second order approximate solution of this periodic orbit has an amplitude given by a1/α1 and an angular frequency computed as(44)β1+ε28b1b3β1-β2+b2β1-β3-a12β1-β2β1-β316β1β1-β2β1-β3.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution III (43), we get its stability region in parameters space as follows:(45)a1>0,a2<0,a3<0.


#### 4.1.5. 2D Torus I

The fifth equilibrium solution,(46)r10=a1α1,r20=a2α2,r30=0,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (44) and (41). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus I (46), we get the stability region in parameters space as follows:(47)a1>0,a2>0,a3<0.


#### 4.1.6. 2D Torus II

The sixth equilibrium solution,(48)r10=a1α1,r20=0,r30=a3α3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (44) and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus II (48), we get the stability region in parameters space as follows:(49)a1>0,a2<0,a3>0.


#### 4.1.7. 2D Torus III

The seventh equilibrium solution,(50)r10=0,r20=a2α2,r30=a3α3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (41) and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus III (50), we get the stability region in parameters space as follows:(51)a1<0,a2>0,a3>0.


#### 4.1.8. 3D Torus

The eighth equilibrium solution,(52)r10=a1α1,r20=a2α2,r30=a3α3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 3D torus in the phase space of case 1 of ring model. The second order approximate solution of this 3D torus has three fundamental angular frequencies given by (44), (41), and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 3D torus (52), we get the stability region in parameters space as follows:(53)a1>0,a2>0,a3>0.


### 4.2. Stability Analysis of Case 2 of Ring Model

The Jacobian matrix of system (33) is given by(54)12a1-r102α1ε-r102α1ε00012a2-r202α2ε-r202α2ε000a3ε2.System (33) has the following equilibrium solutions.


#### 4.2.1. Equilibrium Point

The equilibrium solution(55)ri0=0,i=1,2,3,constitutes an equilibrium point at origin for case 2 of ring model. Evaluating the eigenvalues of Jacobian matrix (54) at equilibrium point (55), we get the following region of stability in parameters space:(56)a1<0,a2<0,a3<0.


#### 4.2.2. Periodic Solution I

The equilibrium solution,(57)r10=0,r20=a2α2,r30=0,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,constitutes a periodic solution in the phase space of case 2 of ring model. The second order approximate solution of this periodic orbit has an amplitude a2/α2 and an angular frequency given as(58)β2+ε2a22β1-β256a2α32+5α23+40α23b1b280α23β2β2-β1.Evaluating the eigenvalues of Jacobian matrix (54) at periodic solution I (57), we get its stability region in parameters space as follows:(59)a1<0,a2>0,a3<0.


#### 4.2.3. Periodic Solution II

The equilibrium solution,(60)r10=a1α1,r20=0,r30=0,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,constitutes a periodic solution in case 2 of ring model phase space. The second order approximate solution of this periodic orbit has an amplitude a1/α1 and an angular frequency given as(61)β1+ε2a12β2-β1+8b1b216β1β1-β2.Evaluating the eigenvalues of Jacobian matrix (54) at periodic solution II (60), we get its stability region in parametric space as follows:(62)a1>0,a2<0,a3<0.


#### 4.2.4. 2D Torus

The equilibrium solution,(63)r10=a1α1,r20=a2α2,r30=0,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,corresponds to a 2D torus in the phase space of case 2 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (61) and (58). Evaluating the eigenvalues of Jacobian matrix (54) at 2D torus (63), we get the following stability region in parametric space:(64)a1>0,a2>0,a3<0.


## 4.1. Stability Analysis of Case 1 of Ring Model

The Jacobian matrix of system (30) is expressed as(34)J=12εa1-r102α1-εr102α100012εa2-r202α2-εr202α200012εa3-r302α3-εr302α3.System (30) has the following equilibrium solutions.


### 4.1.1. Equilibrium Point

The following equilibrium solution,(35)ri0=0,i=1,2,3,is corresponding to an equilibrium point at origin. Evaluating the eigenvalues of Jacobian matrix (34) at equilibrium point (35), we get the stability region of this point in parameters space as follows:(36)a1<0,a2<0,a3<0.


### 4.1.2. Periodic Solution I

The second equilibrium solution,(37)r10=0,r20=0,r30=a3α3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a periodic solution in the phase space of case 1 of ladder model. The second order approximate solution of this periodic orbit has an amplitude given by a3/α3 and an angular frequency computed as(38)β3-ε2β2-β3a32β1-β3+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution I (37), we get its stability region in parameters space as follows:(39)a1<0,a2<0,a3>0.


### 4.1.3. Periodic Solution II

The third equilibrium solution,(40)r10=0,r20=a2α2,r30=0,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,represents a periodic solution in the phase space of case 1 of ring model. The second order approximate solution of this periodic orbit has an amplitude given by a2/α2 and an angular frequency computed as(41)β2-ε28b2b3β1-β2+b1β3-β2-a22β1-β2β2-β316β2β2-β1β2-β3.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution II (40), we get its stability region in parameters space as follows:(42)a1<0,a2>0,a3<0.


### 4.1.4. Periodic Solution III

The fourth equilibrium solution,(43)r10=a1α1,r20=0,r30=0,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,corresponds to a periodic solution in the phase space of case 1 of ring model. The second order approximate solution of this periodic orbit has an amplitude given by a1/α1 and an angular frequency computed as(44)β1+ε28b1b3β1-β2+b2β1-β3-a12β1-β2β1-β316β1β1-β2β1-β3.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution III (43), we get its stability region in parameters space as follows:(45)a1>0,a2<0,a3<0.


### 4.1.5. 2D Torus I

The fifth equilibrium solution,(46)r10=a1α1,r20=a2α2,r30=0,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (44) and (41). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus I (46), we get the stability region in parameters space as follows:(47)a1>0,a2>0,a3<0.


### 4.1.6. 2D Torus II

The sixth equilibrium solution,(48)r10=a1α1,r20=0,r30=a3α3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (44) and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus II (48), we get the stability region in parameters space as follows:(49)a1>0,a2<0,a3>0.


### 4.1.7. 2D Torus III

The seventh equilibrium solution,(50)r10=0,r20=a2α2,r30=a3α3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (41) and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus III (50), we get the stability region in parameters space as follows:(51)a1<0,a2>0,a3>0.


### 4.1.8. 3D Torus

The eighth equilibrium solution,(52)r10=a1α1,r20=a2α2,r30=a3α3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 3D torus in the phase space of case 1 of ring model. The second order approximate solution of this 3D torus has three fundamental angular frequencies given by (44), (41), and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 3D torus (52), we get the stability region in parameters space as follows:(53)a1>0,a2>0,a3>0.


## 4.1.1. Equilibrium Point

The following equilibrium solution,(35)ri0=0,i=1,2,3,is corresponding to an equilibrium point at origin. Evaluating the eigenvalues of Jacobian matrix (34) at equilibrium point (35), we get the stability region of this point in parameters space as follows:(36)a1<0,a2<0,a3<0.


## 4.1.2. Periodic Solution I

The second equilibrium solution,(37)r10=0,r20=0,r30=a3α3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a periodic solution in the phase space of case 1 of ladder model. The second order approximate solution of this periodic orbit has an amplitude given by a3/α3 and an angular frequency computed as(38)β3-ε2β2-β3a32β1-β3+8b1b3+8b2b3β1-β316β3β3-β1β3-β2.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution I (37), we get its stability region in parameters space as follows:(39)a1<0,a2<0,a3>0.


## 4.1.3. Periodic Solution II

The third equilibrium solution,(40)r10=0,r20=a2α2,r30=0,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,represents a periodic solution in the phase space of case 1 of ring model. The second order approximate solution of this periodic orbit has an amplitude given by a2/α2 and an angular frequency computed as(41)β2-ε28b2b3β1-β2+b1β3-β2-a22β1-β2β2-β316β2β2-β1β2-β3.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution II (40), we get its stability region in parameters space as follows:(42)a1<0,a2>0,a3<0.


## 4.1.4. Periodic Solution III

The fourth equilibrium solution,(43)r10=a1α1,r20=0,r30=0,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,corresponds to a periodic solution in the phase space of case 1 of ring model. The second order approximate solution of this periodic orbit has an amplitude given by a1/α1 and an angular frequency computed as(44)β1+ε28b1b3β1-β2+b2β1-β3-a12β1-β2β1-β316β1β1-β2β1-β3.Evaluating the eigenvalues of Jacobian matrix (34) at periodic solution III (43), we get its stability region in parameters space as follows:(45)a1>0,a2<0,a3<0.


## 4.1.5. 2D Torus I

The fifth equilibrium solution,(46)r10=a1α1,r20=a2α2,r30=0,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (44) and (41). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus I (46), we get the stability region in parameters space as follows:(47)a1>0,a2>0,a3<0.


## 4.1.6. 2D Torus II

The sixth equilibrium solution,(48)r10=a1α1,r20=0,r30=a3α3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (44) and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus II (48), we get the stability region in parameters space as follows:(49)a1>0,a2<0,a3>0.


## 4.1.7. 2D Torus III

The seventh equilibrium solution,(50)r10=0,r20=a2α2,r30=a3α3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 2D torus in the phase space of case 1 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (41) and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 2D torus III (50), we get the stability region in parameters space as follows:(51)a1<0,a2>0,a3>0.


## 4.1.8. 3D Torus

The eighth equilibrium solution,(52)r10=a1α1,r20=a2α2,r30=a3α3,θ10˙=ε28b1b3β1-β2+b2β1-β3-β1-β2β1-β3-12a1α1r102+2a12+11α12r10416β1β1-β2β1-β3,θ20˙=-ε28b2b3β1-β2+b1β3-β2-β1-β2β2-β3-12a2α2r202+2a22+11α22r20416β2β2-β1β2-β3,θ30˙=-ε2β2-β3β1-β3-12a3α3r302+2a32+11α32r304+8b1b3+8b2b3β1-β316β3β3-β1β3-β2,corresponds to a 3D torus in the phase space of case 1 of ring model. The second order approximate solution of this 3D torus has three fundamental angular frequencies given by (44), (41), and (38). Evaluating the eigenvalues of Jacobian matrix (34) at 3D torus (52), we get the stability region in parameters space as follows:(53)a1>0,a2>0,a3>0.


## 4.2. Stability Analysis of Case 2 of Ring Model

The Jacobian matrix of system (33) is given by(54)12a1-r102α1ε-r102α1ε00012a2-r202α2ε-r202α2ε000a3ε2.System (33) has the following equilibrium solutions.


### 4.2.1. Equilibrium Point

The equilibrium solution(55)ri0=0,i=1,2,3,constitutes an equilibrium point at origin for case 2 of ring model. Evaluating the eigenvalues of Jacobian matrix (54) at equilibrium point (55), we get the following region of stability in parameters space:(56)a1<0,a2<0,a3<0.


### 4.2.2. Periodic Solution I

The equilibrium solution,(57)r10=0,r20=a2α2,r30=0,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,constitutes a periodic solution in the phase space of case 2 of ring model. The second order approximate solution of this periodic orbit has an amplitude a2/α2 and an angular frequency given as(58)β2+ε2a22β1-β256a2α32+5α23+40α23b1b280α23β2β2-β1.Evaluating the eigenvalues of Jacobian matrix (54) at periodic solution I (57), we get its stability region in parameters space as follows:(59)a1<0,a2>0,a3<0.


### 4.2.3. Periodic Solution II

The equilibrium solution,(60)r10=a1α1,r20=0,r30=0,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,constitutes a periodic solution in case 2 of ring model phase space. The second order approximate solution of this periodic orbit has an amplitude a1/α1 and an angular frequency given as(61)β1+ε2a12β2-β1+8b1b216β1β1-β2.Evaluating the eigenvalues of Jacobian matrix (54) at periodic solution II (60), we get its stability region in parametric space as follows:(62)a1>0,a2<0,a3<0.


### 4.2.4. 2D Torus

The equilibrium solution,(63)r10=a1α1,r20=a2α2,r30=0,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,corresponds to a 2D torus in the phase space of case 2 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (61) and (58). Evaluating the eigenvalues of Jacobian matrix (54) at 2D torus (63), we get the following stability region in parametric space:(64)a1>0,a2>0,a3<0.


## 4.2.1. Equilibrium Point

The equilibrium solution(55)ri0=0,i=1,2,3,constitutes an equilibrium point at origin for case 2 of ring model. Evaluating the eigenvalues of Jacobian matrix (54) at equilibrium point (55), we get the following region of stability in parameters space:(56)a1<0,a2<0,a3<0.


## 4.2.2. Periodic Solution I

The equilibrium solution,(57)r10=0,r20=a2α2,r30=0,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,constitutes a periodic solution in the phase space of case 2 of ring model. The second order approximate solution of this periodic orbit has an amplitude a2/α2 and an angular frequency given as(58)β2+ε2a22β1-β256a2α32+5α23+40α23b1b280α23β2β2-β1.Evaluating the eigenvalues of Jacobian matrix (54) at periodic solution I (57), we get its stability region in parameters space as follows:(59)a1<0,a2>0,a3<0.


## 4.2.3. Periodic Solution II

The equilibrium solution,(60)r10=a1α1,r20=0,r30=0,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,constitutes a periodic solution in case 2 of ring model phase space. The second order approximate solution of this periodic orbit has an amplitude a1/α1 and an angular frequency given as(61)β1+ε2a12β2-β1+8b1b216β1β1-β2.Evaluating the eigenvalues of Jacobian matrix (54) at periodic solution II (60), we get its stability region in parametric space as follows:(62)a1>0,a2<0,a3<0.


## 4.2.4. 2D Torus

The equilibrium solution,(63)r10=a1α1,r20=a2α2,r30=0,θ10˙=ε28b1b2-β1-β2-12a1α1r102+2a12+11α12r10416β1β1-β2,θ20˙=ε2β1-β2-60a2α2r202+10a22+56α32r206+55α22r204+40b1b280β2β2-β1,corresponds to a 2D torus in the phase space of case 2 of ring model. The second order approximate solution of this 2D torus has two fundamental angular frequencies given by (61) and (58). Evaluating the eigenvalues of Jacobian matrix (54) at 2D torus (63), we get the following stability region in parametric space:(64)a1>0,a2>0,a3<0.


## 5. Numerical Simulations

In this section, numerical simulations are carried out to verify the theoretical results obtained in Section4.


### 5.1. Numerical Simulation of Case 1 of Ring Model

Some numerical examples of system (18) are given to address the following scenarios. First, when parameters a1, a2, and a3 have negative values, it is found that a stable equilibrium point at the origin exists as shown in Figure 3. This scenario always occurs for values of parameters satisfying conditions (36). Second, a stable periodic orbit appears in system dynamics when parameters a1 and a2 have negative values while a3 is positive. The equilibrium point at the origin loses its stability and this scenario arises for values of parameters satisfying conditions in (39). Periodic orbits are still observed when parameters a1, a2, and a3 satisfy conditions given in (42) and (45). Figure 4 shows the stable periodic orbit that emerges at a1=-0.4, a2=-0.5, and a3=0.7. In third scenario, parameters a1 and a2 are positive while a3 is negative. In this case, a 2D torus appears as shown in Figure 5. This agrees with previously obtained stability conditions in parameter space given by (47). The Lyapunov exponents spectrum plot associated with the 2D torus case is shown in Figure 6. Note that 2D tori can also be spotted in system dynamics when conditions (49) and (51) are satisfied. Finally, the parameters a1, a2, and a3 are taken positive in numerical simulations. In particular, Figure 7 clarifies the occurrence of 3D torus in phase space of the system at a1=0.4, a2=0.4, and a3=0.5. This case is corresponding to the conditions given in (53). The associated Lyapunov exponents spectrum is plotted and shown in Figure 8.Figure 3
A stable equilibrium point of case1 of ring model obtained when ε=0.1, a1=-1, a2=-0.5, a3=-1, β1=1, β2=2, β3=2, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 4
A stable periodic orbit of case1 of ring model obtained when ε=0.1, a1=-0.4, a2=-0.5, a3=0.7, β1=1, β2=2, β3=3, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 5
A stable 2D torus of case1 of ring model obtained when ε=0.1, a1=0.42, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 6
Lyapunov exponents of the equilibrium solution in Figure5.Figure 7
A stable 3D torus of case1 of ring model obtained when ε=0.1, a1=0.4, a2=0.4, a3=0.5, β1=1, β2=3, β3=2, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 8
Lyapunov exponents of the equilibrium solution in Figure7.


### 5.2. Numerical Simulation of Case 2 of Ring Model

The numerical examples of system (19) are shown in Figures 9–12 to illustrate the following scenarios. First, a stable equilibrium point at origin occurs when parameters a1, a2, and a3 are all negative as shown in Figure 9. Setting the values of parameters to satisfy conditions (56) leads to this scenario. Second, when parameters a1 and a3 are negative while a2 is positive, a stable periodic orbit is obtained. The stable equilibrium point at origin becomes unstable and this scenario arises for the values of parameters conditions in (59). A stable periodic orbit is shown in Figure 10 when a1=-0.7, a2=0.4, and a3=-0.5. Periodic orbits also occur when parameters a1, a2, and a3 satisfy conditions in (62). Finally, a 2D torus appears in system dynamics when parameters a1 and a2 are positive while a3 is negative. This is clarified in Figure 11 where a1=0.45, a2=0.3, and a3=-0.2. The values of parameters a1, a2, and a3 in this example satisfy the conditions in (64). The associated Lyapunov exponents spectrum is plotted and shown in Figure 12.Figure 9
A stable equilibrium point of case2 of ring model obtained when ε=0.1, a1=-1, a2=-1, a3=-1, β1=3, β2=1, β3=2, and b1=b2=b3=α1=α2=α3=α4=1 in system (19).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 10
A stable periodic orbit of case2 of ring model obtained when ε=0.1, a1=-0.7, a2=0.4, a3=-0.5, β1=1, β2=2, β3=1, and b1=b2=b3=α1=α2=α3=α4=1 in system (19).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 11
A stable 2D torus of case2 of ring model obtained when ε=0.1, a1=0.45, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1, and b1=b2=b3=α1=α2=α3=α4=1 in system (19).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 12
Lyapunov exponents of the equilibrium solution in Figure11.


### 5.3. Effects of Time-Delayed Coupling

In this part, the influences of time-delayed coupling on the aforementioned numerical examples are numerically investigated. For case1 of ring model the mathematical model can be rewritten in the following form:(65)xj¨-εxj˙aj-αjxj2+βjxj-bjε∑k=1,k≠j3xkτ=0,j=1,2,3.xkτ=xk(t-τ) and τ denotes constant time delay value. It is observed that for the values of parameters used in Figure 13 the equilibrium point loses its stability when time delay value increases over a certain value (approximately 0.4) and a periodic orbit is induced.Figure 13
Phase portraits of case1 of ring model with time-delayed coupling obtained for three values of τ, ε=0.1, a1=-1, a2=-0.5, a3=-1, β1=1, β2=2, β3=2, and b1=b2=b3=α1=α2=α3=1 in system (65). Time delay values are for (a) τ=0.3, (b) τ=0.4, and (c) τ=0.6.
(a) (b) (c)However, the effects of time delay coupling on other dynamics of system are almost negligible as shown in Figure14.Figure 14
Phase portraits of case1 of ring model with time-delayed coupling obtained for (a) τ=2, ε=0.1, a1=-0.4, a2=-0.5, a3=0.7, β1=1, β2=2, β3=3 and b1=b2=b3=α1=α2=α3=1, (b) τ=2, ε=0.1, a1=0.42, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1 and b1=b2=b3=α1=α2=α3=1, and (c) τ=1, ε=0.1, a1=0.4, a2=0.4, a3=0.5, β1=1, β2=3, β3=2 and b1=b2=b3=α1=α2=α3=1.
(a) (b) (c)For case2 of ring model the rate equations can be reformulated as follows:(66)x1¨-εx1˙a1-α1x12+β1x1-b1εx2τ+x3τ=0,x2¨-εx2˙a2-α3x23-α2x22+β2x2-b2εx1τ+x3τ=0,x3¨-εx3˙a3-α4x33+β3x3-b3εx1τ+x2τ=0.Although the effect of time-delayed coupling on the stable equilibrium point at the origin and torus is measly as shown in Figure 15, it is observed that it has an obvious impact on periodic solution that exists when conditions are satisfied as shown in Figure 16.Figure 15
Phase portraits of case2 of ring model with time-delayed coupling obtained for (a) τ=2, ε=0.1, a1=-1, a2=-1, a3=-1, β1=3, β2=1, β3=2 and b1=b2=b3=α1=α2=α3=α4=1 and (b) τ=0.5, ε=0.1, a1=0.45, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1 and b1=b2=b3=α1=α2=α3=α4=1.
(a) (b)Figure 16
Phase portraits of case2 of ring model with time-delayed coupling obtained for two values of τ, ε=0.1, a1=-1, a2=-0.5, a3=-1, β1=1, β2=2, β3=2 and b1=b2=b3=α1=α2=α3=1 in system (66). Time delay values are for (a) τ=0.7 and (b) τ=0.1.
(a) (b)


## 5.1. Numerical Simulation of Case 1 of Ring Model

Some numerical examples of system (18) are given to address the following scenarios. First, when parameters a1, a2, and a3 have negative values, it is found that a stable equilibrium point at the origin exists as shown in Figure 3. This scenario always occurs for values of parameters satisfying conditions (36). Second, a stable periodic orbit appears in system dynamics when parameters a1 and a2 have negative values while a3 is positive. The equilibrium point at the origin loses its stability and this scenario arises for values of parameters satisfying conditions in (39). Periodic orbits are still observed when parameters a1, a2, and a3 satisfy conditions given in (42) and (45). Figure 4 shows the stable periodic orbit that emerges at a1=-0.4, a2=-0.5, and a3=0.7. In third scenario, parameters a1 and a2 are positive while a3 is negative. In this case, a 2D torus appears as shown in Figure 5. This agrees with previously obtained stability conditions in parameter space given by (47). The Lyapunov exponents spectrum plot associated with the 2D torus case is shown in Figure 6. Note that 2D tori can also be spotted in system dynamics when conditions (49) and (51) are satisfied. Finally, the parameters a1, a2, and a3 are taken positive in numerical simulations. In particular, Figure 7 clarifies the occurrence of 3D torus in phase space of the system at a1=0.4, a2=0.4, and a3=0.5. This case is corresponding to the conditions given in (53). The associated Lyapunov exponents spectrum is plotted and shown in Figure 8.Figure 3
A stable equilibrium point of case1 of ring model obtained when ε=0.1, a1=-1, a2=-0.5, a3=-1, β1=1, β2=2, β3=2, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 4
A stable periodic orbit of case1 of ring model obtained when ε=0.1, a1=-0.4, a2=-0.5, a3=0.7, β1=1, β2=2, β3=3, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 5
A stable 2D torus of case1 of ring model obtained when ε=0.1, a1=0.42, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 6
Lyapunov exponents of the equilibrium solution in Figure5.Figure 7
A stable 3D torus of case1 of ring model obtained when ε=0.1, a1=0.4, a2=0.4, a3=0.5, β1=1, β2=3, β3=2, and b1=b2=b3=α1=α2=α3=1 in system (18).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 8
Lyapunov exponents of the equilibrium solution in Figure7.


## 5.2. Numerical Simulation of Case 2 of Ring Model

The numerical examples of system (19) are shown in Figures 9–12 to illustrate the following scenarios. First, a stable equilibrium point at origin occurs when parameters a1, a2, and a3 are all negative as shown in Figure 9. Setting the values of parameters to satisfy conditions (56) leads to this scenario. Second, when parameters a1 and a3 are negative while a2 is positive, a stable periodic orbit is obtained. The stable equilibrium point at origin becomes unstable and this scenario arises for the values of parameters conditions in (59). A stable periodic orbit is shown in Figure 10 when a1=-0.7, a2=0.4, and a3=-0.5. Periodic orbits also occur when parameters a1, a2, and a3 satisfy conditions in (62). Finally, a 2D torus appears in system dynamics when parameters a1 and a2 are positive while a3 is negative. This is clarified in Figure 11 where a1=0.45, a2=0.3, and a3=-0.2. The values of parameters a1, a2, and a3 in this example satisfy the conditions in (64). The associated Lyapunov exponents spectrum is plotted and shown in Figure 12.Figure 9
A stable equilibrium point of case2 of ring model obtained when ε=0.1, a1=-1, a2=-1, a3=-1, β1=3, β2=1, β3=2, and b1=b2=b3=α1=α2=α3=α4=1 in system (19).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 10
A stable periodic orbit of case2 of ring model obtained when ε=0.1, a1=-0.7, a2=0.4, a3=-0.5, β1=1, β2=2, β3=1, and b1=b2=b3=α1=α2=α3=α4=1 in system (19).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 11
A stable 2D torus of case2 of ring model obtained when ε=0.1, a1=0.45, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1, and b1=b2=b3=α1=α2=α3=α4=1 in system (19).
(a)
x 1-x2-x3 (b)
x 1-x2 (c)
x 1-x3 (d)
x 2-x3Figure 12
Lyapunov exponents of the equilibrium solution in Figure11.


## 5.3. Effects of Time-Delayed Coupling

In this part, the influences of time-delayed coupling on the aforementioned numerical examples are numerically investigated. For case1 of ring model the mathematical model can be rewritten in the following form:(65)xj¨-εxj˙aj-αjxj2+βjxj-bjε∑k=1,k≠j3xkτ=0,j=1,2,3.xkτ=xk(t-τ) and τ denotes constant time delay value. It is observed that for the values of parameters used in Figure 13 the equilibrium point loses its stability when time delay value increases over a certain value (approximately 0.4) and a periodic orbit is induced.Figure 13
Phase portraits of case1 of ring model with time-delayed coupling obtained for three values of τ, ε=0.1, a1=-1, a2=-0.5, a3=-1, β1=1, β2=2, β3=2, and b1=b2=b3=α1=α2=α3=1 in system (65). Time delay values are for (a) τ=0.3, (b) τ=0.4, and (c) τ=0.6.
(a) (b) (c)However, the effects of time delay coupling on other dynamics of system are almost negligible as shown in Figure14.Figure 14
Phase portraits of case1 of ring model with time-delayed coupling obtained for (a) τ=2, ε=0.1, a1=-0.4, a2=-0.5, a3=0.7, β1=1, β2=2, β3=3 and b1=b2=b3=α1=α2=α3=1, (b) τ=2, ε=0.1, a1=0.42, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1 and b1=b2=b3=α1=α2=α3=1, and (c) τ=1, ε=0.1, a1=0.4, a2=0.4, a3=0.5, β1=1, β2=3, β3=2 and b1=b2=b3=α1=α2=α3=1.
(a) (b) (c)For case2 of ring model the rate equations can be reformulated as follows:(66)x1¨-εx1˙a1-α1x12+β1x1-b1εx2τ+x3τ=0,x2¨-εx2˙a2-α3x23-α2x22+β2x2-b2εx1τ+x3τ=0,x3¨-εx3˙a3-α4x33+β3x3-b3εx1τ+x2τ=0.Although the effect of time-delayed coupling on the stable equilibrium point at the origin and torus is measly as shown in Figure 15, it is observed that it has an obvious impact on periodic solution that exists when conditions are satisfied as shown in Figure 16.Figure 15
Phase portraits of case2 of ring model with time-delayed coupling obtained for (a) τ=2, ε=0.1, a1=-1, a2=-1, a3=-1, β1=3, β2=1, β3=2 and b1=b2=b3=α1=α2=α3=α4=1 and (b) τ=0.5, ε=0.1, a1=0.45, a2=0.3, a3=-0.2, β1=1, β2=2, β3=1 and b1=b2=b3=α1=α2=α3=α4=1.
(a) (b)Figure 16
Phase portraits of case2 of ring model with time-delayed coupling obtained for two values of τ, ε=0.1, a1=-1, a2=-0.5, a3=-1, β1=1, β2=2, β3=2 and b1=b2=b3=α1=α2=α3=1 in system (66). Time delay values are for (a) τ=0.7 and (b) τ=0.1.
(a) (b)


## 6. Conclusion

In this work, we have investigated new analytical results concerning the dynamical behaviors of three coupled memristor-based oscillators with identical and different nonlinearities in ring configuration. The stability regions in parameters space are attained for each type of induced dynamics. The analytical perturbation techniques were employed to achieve this goal analytically. Furthermore, numerical simulations are performed and they agreed with analytical results obtained. It was demonstrated that the proposed models exhibit complex dynamics including periodic orbits, 2D tori, and 3D tori in ring coupling configuration. A detailed analytical investigations of the effects of time-delayed couplings, with different time delays, can be considered in future work. Also, the mathematical modeling of coupled memristor-based oscillators using fractional order derivatives, in order to accurately incorporate memory effects, can be examined.



---

*Source: 4394058-2018-10-21.xml*