# Synchronization in Dynamically Coupled Fractional-Order Chaotic Systems: Studying the Effects of Fractional Derivatives

Authors: J. L. Echenausía-Monroy; C. A. Rodríguez-Martíne; L. J. Ontañón-García; J. Alvarez; J. Pena Ramirez
Journal: Complexity
 (2021)

---

## Abstract

This study presents the effectiveness of dynamic coupling as a synchronization strategy for fractional chaotic systems. Using an auxiliary system as a link between the oscillators, we investigate the onset of synchronization in the coupled systems and we analytically determine the regions where both systems achieve complete synchronization. In the analysis, the integration order is considered as a key parameter affecting the onset of full synchronization, considering the stability conditions for fractional systems. The local stability of the synchronous solution is studied using the linearized error dynamics. Moreover, some statistical metrics such as the average synchronization error and Pearson’s correlation are used to numerically identify the synchronous behavior. Two particular examples are considered, namely, the fractional-order Rössler and Chua systems. By using bifurcation diagrams, it is also shown that the integration order has a strong influence not only on the onset of full synchronization but also on the individual dynamic behavior of the uncoupled systems.


---

## Body

## 1. Introduction

Synchronization is an emergent physical phenomenon caused by the interaction of two or more dynamic entities that pervade the natural world [1–4]. In the case of oscillating units, the synchronization phenomenon can be defined as the adjustment of temporal evolution to a common rhythm.For the case of integer-order systems, there exists a vast and mature literature where we can find different interconnection schemes for synchronizing dynamic systems, like, for example, master-slave synchronization scheme, adaptive synchronization, and synchronization based on state observers, to name a few [5–10]. Although each of these strategies is effective, there are limitations in their applications, e.g., there are cases where these schemes have marginal ranges for which the synchronous response is achieved or have poor robustness to maintain a stable synchronous state under the influence of external disturbances. This is one of the reasons why dynamic interconnections have emerged as an alternative to the classical static schemes. In this case, the interaction between agents is indirectly achieved through a suitably designed dynamic coupling [11–13]. This type of synchronization strategy has shown better performance than static couplings. In particular, dynamic coupling increases the intervals of coupling strength values for which it is possible to achieve synchronized behavior, and it may also be possible to synchronize systems that cannot be synchronized with static coupling [11].On the other hand, the use of fractional calculus has been extensively studied in nonlinear systems (see, e.g., [14–18]) and also, there exist notable contributions related to the study of synchronization in fractional-order systems (see, e.g., [16, 19–23]). For example, there is work based on applying sliding modes to fractional-order models to achieve synchronization [24–29]. The modeling and analytical study of fractional-order systems is also a fruitful field, e.g., the use of the Razumikhin approximation for fractional-order systems with delay [30, 31], the extrapolation of Lyapunov theory to fractional systems [32, 33], and the existence and uniqueness of equilibrium points of the Mittag–Leffler criteria [34, 35]. However, the use of dynamic couplings in the context of fractional-order systems seems to be unexplored so far.Consequently, in this study we present a synchronization scheme of fractional order based on dynamic coupling. In particular, a master-slave interconnection is considered, in which the systems have an indirect interaction via a fractional-order linear system. In the analysis, the Rössler equation [36] and the Chua double-scroll oscillator [37] are considered. Among the questions to be addressed is whether a dynamic coupling designed for integer-order systems is also effective in its fractional-order version? If so, how does the derivative order influence the occurrence of synchronization in the coupled systems? The local stability of the synchronous solution in the coupled systems is investigated by analyzing the error dynamics, and furthermore, the onset of synchronization is also numerically investigated by computing some statistical metrics like Pearson’s correlation between time series. Additionally, using bifurcation diagrams, we have investigated the dynamic behavior of the uncoupled systems. The obtained results show that the integration order has a strong influence on the stability of the synchronous solution, and interestingly, it also produces a period-doubling cascade route to chaos in the uncoupled systems.The rest of the study is organized as follows: Section2 describes the basics of fractional calculus and gives a brief introduction to fractional-order systems. Then, Section 3 describes the proposed synchronization scheme and the local stability of the synchronous solution in the coupled systems is discussed. Subsequently, in Section 4, the performance of the dynamic coupling is investigated using the Rössler and Chua systems as application examples. Finally, Sections 6 and 7 are dedicated to the discussion and conclusions, respectively.


## 2. Preliminaries

This section presents a brief overview of some basic concepts about fractional-order systems. In particular, the Caputo derivative, the general representation of a fractional-order system, and the stability of linear time-invariant fractional-order systems are revisited.


### 2.1. Fractional Caputo Derivative

In the literature, there are various definitions of fractional-order derivatives, the most common being the Riemann–Liouville and Caputo operators [38, 39]. The fractional Caputo derivative for a time-invariant system described by the vector field fx is defined as(1)aD0qfx=1Γn−q∫axfntx−tq−n+1dt,with n=⌈q⌉ for the integration order 0<q<1, being Γ the gamma function defined as follows:(2)Γz=∫0∞tz−1e−tdt.


### 2.2. Fractional-Order Dynamic System

A commensurate fractional-order time-invariant system can be described, in general, as follows:(3)D0nkxt=ft,xt,D0n1xt,D0n2xt,…,D0nk−1xt,subject to initial conditions(4)xj0=x0j,withj=0,1,…,nk−1,where n1,n2,…,nk are rational numbers, such that nk>nk−1>⋯>n1>0, nj−nj−1≤1 for all j=2,3,…,k, and 0<n1≤1. The least common multiple of the denominator of n1,n2,…,nk is defined by M and set q=1/M and N=Mnk. Then, equation (3) can be expressed as follows [38]:(5)D0qx0t=x1t,D0qx1t=x2t,⋮D0qxN−2t=xN−1t,D0qxN−1t=ft,x0t,xn1/qt,…,xnk−1/qt.


### 2.3. Stability of Time-Invariant Fractional-Order System

A linear time-invariant fractional-order system is described by(6)dqxtdtq=Ax,where x∈ℝn is the state vector, A∈ℝn×n is a constant matrix, and 0<q<1 is the fractional commensurate derivative order. Then, the stability of the system described by equation (6) is determined as follows [38]:(i)
System (6) is stable, if and only if argλj≥qπ/2,∀j=1,2,…,n(ii)
System (6) is asymptotically stable, if and only if argλj>qπ/2,∀j=1,2,…,n(iii)
System (6) is unstable, if and only if argλj<qπ/2, for at least one λj, j=1,2,…,nFrom these conditions, it is clear that the local stability of fractional-order systems depends on the integration orderq, so that the stability of an equilibrium point can be modified by the fractional order, and therefore, the stability region at the complex plane is as in Figure 1 [40].Figure 1
Stability region in a fractional-order linear time-invariant system for0<q<1.


## 2.1. Fractional Caputo Derivative

In the literature, there are various definitions of fractional-order derivatives, the most common being the Riemann–Liouville and Caputo operators [38, 39]. The fractional Caputo derivative for a time-invariant system described by the vector field fx is defined as(1)aD0qfx=1Γn−q∫axfntx−tq−n+1dt,with n=⌈q⌉ for the integration order 0<q<1, being Γ the gamma function defined as follows:(2)Γz=∫0∞tz−1e−tdt.


## 2.2. Fractional-Order Dynamic System

A commensurate fractional-order time-invariant system can be described, in general, as follows:(3)D0nkxt=ft,xt,D0n1xt,D0n2xt,…,D0nk−1xt,subject to initial conditions(4)xj0=x0j,withj=0,1,…,nk−1,where n1,n2,…,nk are rational numbers, such that nk>nk−1>⋯>n1>0, nj−nj−1≤1 for all j=2,3,…,k, and 0<n1≤1. The least common multiple of the denominator of n1,n2,…,nk is defined by M and set q=1/M and N=Mnk. Then, equation (3) can be expressed as follows [38]:(5)D0qx0t=x1t,D0qx1t=x2t,⋮D0qxN−2t=xN−1t,D0qxN−1t=ft,x0t,xn1/qt,…,xnk−1/qt.


## 2.3. Stability of Time-Invariant Fractional-Order System

A linear time-invariant fractional-order system is described by(6)dqxtdtq=Ax,where x∈ℝn is the state vector, A∈ℝn×n is a constant matrix, and 0<q<1 is the fractional commensurate derivative order. Then, the stability of the system described by equation (6) is determined as follows [38]:(i)
System (6) is stable, if and only if argλj≥qπ/2,∀j=1,2,…,n(ii)
System (6) is asymptotically stable, if and only if argλj>qπ/2,∀j=1,2,…,n(iii)
System (6) is unstable, if and only if argλj<qπ/2, for at least one λj, j=1,2,…,nFrom these conditions, it is clear that the local stability of fractional-order systems depends on the integration orderq, so that the stability of an equilibrium point can be modified by the fractional order, and therefore, the stability region at the complex plane is as in Figure 1 [40].Figure 1
Stability region in a fractional-order linear time-invariant system for0<q<1.


## 3. Proposed Synchronization Scheme Based on Dynamic Coupling

The dynamic interconnection considered in this study has been presented in [11] for integer-order systems and is adapted here to the fractional-order case. The scheme, where the interaction between the systems is indirect via a dynamical system (Figure 2) is described by the following set of equations:(7)Dqxm=Fxm,Dqxs=Fxs−B1h,Dqh=Gh−kB2xm−xs,where xm,xs∈ℝn represents the state vectors of both the master and slave systems and h=h1,h2T for hi∈ℝ,i=1,2 is the state variables of the dynamic coupling. It is assumed that the vector field F is smooth enough, which can be either linear or nonlinear, and the coupling force between the systems is denoted by k.Figure 2
Schematic representation of unidirectionally coupled fractional-order systems interacting via fractional-order dynamic coupling.On the other hand, the design of a dynamic coupling involves two coupling matrices, denotedB1∈ℝn×2 and B2∈ℝ2×n. These matrices are generated under the premise that only one of the elements of each of these matrices is equal to 1, and the other entries are zero, which means that the coupling is applied only in one state variable of the slave system and that the coupling considers only one measured variable.Finally, the matrixG from equation (7) is given by(8)G=−αc1−γ1−γ2,where γ1,γ2, and αc are design parameters of the dynamic coupling. The construction of the coupling system, for the integer-order case, is inspired by the so-called Huygen’s coupling [11, 41], which in its simplest form can be interpreted as a damped oscillator.Since in this study the focus is on synchronization, it is necessary to give the following definition.Definition 1.
The coupled systems, equation (7), are said to be asymptotically synchronized if(9)limt⟶∞xm−xs=0,limt⟶∞h=0.


### 3.1. Local Stability Analysis

In order to investigate the stability of the synchronous solution defined in equation (9), we proceed as follows. First, it is assumed that function F in equation (7) can be written as the sum of linear and nonlinear components, i.e., it is assumed that(10)Fxi=Pxi+Exi,where P∈ℝn×n is a constant matrix and Exi∈ℝn is a vector containing nonlinear terms.Next, the synchronization error is defined asep:=xm−xs,hT. Note that in the definition of the error, we have included the state h of the dynamic coupling. The reason is because the parameters in the dynamic coupling should be chosen such that, when the systems synchronize, the coupling vanishes. Then, by replacing equations (10) into (7), and computing the corresponding synchronization error dynamics, we obtain(11)Dqep=A˜ep+gpt,ep,where(12)A˜=PB1−KB2G,gpt,ep=fxm−fxsO,where O=0,0T. Furthermore, note that the term gpt,ep is a vanishing perturbation [42] because gpt,0=0. Then, the stability properties of system equation (11) are fully determined by the eigenvalues of the matrix A˜. In particular, following the results presented in Section 2.3, we have that the synchronization error dynamics (equation (11)) is locally asymptotically stable if(13)argλj>qπ2,∀j=1,2,…,n.Thus, if it is possible to find values ofk,γ1,γ2, and αc such that the above condition is satisfied, then the coupled systems described by equation (7) will achieve complete synchronization, according to Definition 1.


### 3.2. Statistical Metrics Used for Detecting Complete Synchronization

In this study, the onset of synchronization in the coupled systems (equation (7)) is also numerically studied by computing the following synchronization index.(14)S=∑i=1nei+∑i=1nCi,where n is the dimension of the systems to be synchronized, ei is the absolute value of the synchronization error between the i-th state variable of the master and slave systems, i.e.,(15)ei=xmi−xsi,and Ci is the Pearson correlation [43], computed from the i-th state variable of the master and slave systems described by(16)Ci=σNmiNsiσNmiσNsi,where σNmiNsi is the covariance between the data obtained from the time series of the state variables of the master and slave systems, and σNmiσNsi is the standard deviation obtained from the i-th state variable of the master (slave) oscillator. Finally, when S=n, the systems are synchronized.In the next section, the onset of synchronization with dynamical coupling is studied for two particular fractional-order chaotic systems, namely the Rössler and the Chua oscillators.


## 3.1. Local Stability Analysis

In order to investigate the stability of the synchronous solution defined in equation (9), we proceed as follows. First, it is assumed that function F in equation (7) can be written as the sum of linear and nonlinear components, i.e., it is assumed that(10)Fxi=Pxi+Exi,where P∈ℝn×n is a constant matrix and Exi∈ℝn is a vector containing nonlinear terms.Next, the synchronization error is defined asep:=xm−xs,hT. Note that in the definition of the error, we have included the state h of the dynamic coupling. The reason is because the parameters in the dynamic coupling should be chosen such that, when the systems synchronize, the coupling vanishes. Then, by replacing equations (10) into (7), and computing the corresponding synchronization error dynamics, we obtain(11)Dqep=A˜ep+gpt,ep,where(12)A˜=PB1−KB2G,gpt,ep=fxm−fxsO,where O=0,0T. Furthermore, note that the term gpt,ep is a vanishing perturbation [42] because gpt,0=0. Then, the stability properties of system equation (11) are fully determined by the eigenvalues of the matrix A˜. In particular, following the results presented in Section 2.3, we have that the synchronization error dynamics (equation (11)) is locally asymptotically stable if(13)argλj>qπ2,∀j=1,2,…,n.Thus, if it is possible to find values ofk,γ1,γ2, and αc such that the above condition is satisfied, then the coupled systems described by equation (7) will achieve complete synchronization, according to Definition 1.


## 3.2. Statistical Metrics Used for Detecting Complete Synchronization

In this study, the onset of synchronization in the coupled systems (equation (7)) is also numerically studied by computing the following synchronization index.(14)S=∑i=1nei+∑i=1nCi,where n is the dimension of the systems to be synchronized, ei is the absolute value of the synchronization error between the i-th state variable of the master and slave systems, i.e.,(15)ei=xmi−xsi,and Ci is the Pearson correlation [43], computed from the i-th state variable of the master and slave systems described by(16)Ci=σNmiNsiσNmiσNsi,where σNmiNsi is the covariance between the data obtained from the time series of the state variables of the master and slave systems, and σNmiσNsi is the standard deviation obtained from the i-th state variable of the master (slave) oscillator. Finally, when S=n, the systems are synchronized.In the next section, the onset of synchronization with dynamical coupling is studied for two particular fractional-order chaotic systems, namely the Rössler and the Chua oscillators.


## 4. Application Example 1: Rössler System

The fractional-order version of the well-known Rössler system [44] is given by(17)Dqx1=−x2−x3,Dqx2=x1+ax2,Dqx3=b+x1x3−cx3,where xi,i=1,2,3 denotes the state variables and a,b,c are constants.It has been reported that every system has a limit of integration order for which it is possible to use a fractional-order derivative without stabilizing its dynamics [16, 17, 45]. This could be interpreted to mean that the dynamics described in equation (17) are restricted to have at least one eigenvalue in the unstable region, and this being true only if argλ<qπ/2 for at least one of its eigenvalues λ. With a=0.2,b=0.2,c=5.7, the eigenvalues obtained by the Jacobian matrix evaluated at the equilibrium point E1=0.007,−0.0351,0.0351 are λ=0.0970±0.9952i;−5.6870. Since minargλ=1.4736, the critical order of integration is defined as qc=minargλ2/π (Section 2), and by substitution it is obtained as qc=0.9381. This result is confirmed by the bifurcation diagram as shown in Figure 3, where the local maxima in x1 are plotted as a function of the variation of the integration order. Note that the derivative order is the same for all state variables of the system. The inset shows the attractors and the corresponding integration order considered for the developed analysis.Figure 3
Bifurcation diagram of the Rössler system by modifying the integration order, equation (16). The abscissa axis indicates that the three oscillator state variables have the same integration order.It is worth noting that the bifurcation diagram shown in Figure3 was numerically calculated using the Adams–Bashforth–Moulton (ABM) method [46] and following the guides for plotting a proper bifurcation diagram [47].


### 4.1. Dynamically Coupled Fractional-Order Rössler Systems

Now, we consider a pair of identical Rössler systems [44] described by(18)master systemDqmx1m=−x2m−x3m,Dqmx2m=x1m+ax2m,Dqmx3m=b+x1mx3m−cx3m,(19)slave systemDqsx1s=−x2s−x3s−h2,Dqsx2s=x1s+ax2s,Dqsx3s=b+x1sx3s−cx3s,(20)dynamic couplingDqdh1=−αch1+h2,Dqdh2=−γ1h1−γ2h2−kx2m−x2s,where xm,s denotes the state variables of the master and slave systems, respectively, h1,2 denotes the states of the dynamic coupling, k indicates the coupling force between the oscillators, a,b,c are constants of the Rössler model, and qj, j=m,s,d denotes the integration order of the master, slave, and dynamic coupling, respectively.Here, we consider the case in which all orders of integration are equal, i.e.,qm=qs=qd. Then, the dynamic behavior of system equations (18)–(20) is numerically studied as a function of the coupling strength k and the integration order q. For this purpose, equations (18)–(20) are numerically integrated with the following parameter values cf. [11]: αc=a,γ1=k,γ2=k,a=0,2,b=0,2,c=5,7 and initial conditions xim,xis,hi=0.1,0.1,0.1,−0.2,−0.1,−0.1,0,0. On the other hand, the coupling strength is varied in the interval 0≤k≤20 with a test size Δk=0.2 and the integration order is varied in the interval 0.96≤q≤1 at Δq=2e−4.The obtained results are shown in Figure4(a), where the colors indicate the value of the synchronization index S (see equation (14)). Synchronous behavior is indicated by the dark blue color (S=3), while the remaining colors indicate unsynchronized dynamics.Figure 4
Numerical results for equations (18)–(20). (a) Limit behavior of the system as a function of coupling strength k and integration order q. (b) Synchronization index equation (14) as a function of integration order, for a fixed k=10 and k=15. A value of S=3 indicates complete synchronization.
(a)(b)Moreover, it is evident from Figure4(a) that there are regions where, for a fixed coupling strength, the limit behavior is completely determined by the integration order. For example, for a fixed coupling strength of k=10, the synchronization index S abruptly changes when the integration order q is within the interval 0.988<q<0.995, as can be seen in Figure 4(b), but if the coupling strength is increased up to k=15, the systems always achieve complete synchronization.To validate the synchronization regions obtained by the time series analysis, we conduct a stability analysis following the results presented in Section3.1. First, let the following synchronization errors be defined as follows: ej=xjm−xjs, for j=1,2,3, and e4=h1,e5=h2. Then, the resulting error dynamics are given by(21)Dqe1=−e2−e3+e5,Dqe2=e1+ae2,Dqe3=x1mx3m−x1sx3s−ce3,Dqe4=−αce4+e5,Dqe5=−γ1e4−γ2e5−ke2.Note that the error dynamics in equation (21) can be written in the form of equation (11) with(22)A˜=0−1−1011a00000−c00000−αc10−k0−γ1−γ2,gpt,ep=00x1mx3m−x1sx3s00.Note that the termgpt,ep in equation (21) is indeed a vanishing perturbation since when the systems synchronize it follows that x1m=x1s,x2m=x2s, and x3m=x3s and therefore gpt,0=0. Consequently, the local stability of the synchronous solution in the coupled systems described by equations (18)–(20) can be determined from the condition in equation (13). In particular, we compute condition equation (13) as a function of the integration order q and the coupling strength k. The obtained results are shown in Figure 5(a), where the blue region corresponds to values of k and q for which condition (13) is not satisfied and thus the synchronous solution is unstable, whereas on the white region, condition (13) is satisfied and then the synchronous solution is expected to be stable. For the sake of comparison, Figure 5(b) shows the overlap of Figure 4 with Figure 5(a). It can be seen that there is a good agreement between the numerical and the analytical results.Figure 5
(a) Stability analysis from linear error dynamics for the coupled Rössler system shown in equation (21). (b) Overlap of the analytical result in the synchronization map shown in Figure 4(a).
(a)(b)Remark 1.
In the previous analysis, we have considered the case where the integration order of the systems and the dynamic coupling are the same. However, we also have conducted a numerical study in which the integration orders are different. In particular, we have numerically integrated equations (18)–(20) using the parameter values considered before, as a function of the integration orders of the master and slave systems, while the integration order of the dynamic coupling remains fixed. The integration orders qm,qs are varied in the interval 0.98≤qj<1, j=m,s, and considering the integration order qd=0.985. The obtained results are shown in Figure 6(a) where the blue areas correspond to synchronization (S=3). From the obtained results, it is clear to see that the integration order of the master and slave systems should be almost the same to observe a synchronized behavior and that larger differences are tolerated as long as the integration order of both systems approaches to one.
On the other hand, Figure6(b) shows the obtained results for the case that only the integration order of the dynamic coupling is varied, while the oscillators are assumed to have integer order, i.e., qm=qs=1. In this case, equations (18)–(20) are numerically integrated by varying the integration order qd of the dynamic coupling in the interval 0.8<qd<1 and the corresponding synchronization index, given by equation (14), is calculated. The obtained results are shown in Figure 6(b) for two different coupling forces, where the choice of these k values corresponds to those reported in [11] and those obtained in Figure 4(b).Figure 6
(a) Synchronization map as a function of integration ordersqm and qs of master and slave Rössler systems, equations (18)–(20), and considering a fixed integration order of qd=0.985 in the dynamic coupling. (b) Synchronization map for system equations (18)–(20) as a function of the integration order qd of the dynamic coupling, while considering that the master and slave oscillators have fixed integer order, i.e., q=1.
(a)(b)


## 4.1. Dynamically Coupled Fractional-Order Rössler Systems

Now, we consider a pair of identical Rössler systems [44] described by(18)master systemDqmx1m=−x2m−x3m,Dqmx2m=x1m+ax2m,Dqmx3m=b+x1mx3m−cx3m,(19)slave systemDqsx1s=−x2s−x3s−h2,Dqsx2s=x1s+ax2s,Dqsx3s=b+x1sx3s−cx3s,(20)dynamic couplingDqdh1=−αch1+h2,Dqdh2=−γ1h1−γ2h2−kx2m−x2s,where xm,s denotes the state variables of the master and slave systems, respectively, h1,2 denotes the states of the dynamic coupling, k indicates the coupling force between the oscillators, a,b,c are constants of the Rössler model, and qj, j=m,s,d denotes the integration order of the master, slave, and dynamic coupling, respectively.Here, we consider the case in which all orders of integration are equal, i.e.,qm=qs=qd. Then, the dynamic behavior of system equations (18)–(20) is numerically studied as a function of the coupling strength k and the integration order q. For this purpose, equations (18)–(20) are numerically integrated with the following parameter values cf. [11]: αc=a,γ1=k,γ2=k,a=0,2,b=0,2,c=5,7 and initial conditions xim,xis,hi=0.1,0.1,0.1,−0.2,−0.1,−0.1,0,0. On the other hand, the coupling strength is varied in the interval 0≤k≤20 with a test size Δk=0.2 and the integration order is varied in the interval 0.96≤q≤1 at Δq=2e−4.The obtained results are shown in Figure4(a), where the colors indicate the value of the synchronization index S (see equation (14)). Synchronous behavior is indicated by the dark blue color (S=3), while the remaining colors indicate unsynchronized dynamics.Figure 4
Numerical results for equations (18)–(20). (a) Limit behavior of the system as a function of coupling strength k and integration order q. (b) Synchronization index equation (14) as a function of integration order, for a fixed k=10 and k=15. A value of S=3 indicates complete synchronization.
(a)(b)Moreover, it is evident from Figure4(a) that there are regions where, for a fixed coupling strength, the limit behavior is completely determined by the integration order. For example, for a fixed coupling strength of k=10, the synchronization index S abruptly changes when the integration order q is within the interval 0.988<q<0.995, as can be seen in Figure 4(b), but if the coupling strength is increased up to k=15, the systems always achieve complete synchronization.To validate the synchronization regions obtained by the time series analysis, we conduct a stability analysis following the results presented in Section3.1. First, let the following synchronization errors be defined as follows: ej=xjm−xjs, for j=1,2,3, and e4=h1,e5=h2. Then, the resulting error dynamics are given by(21)Dqe1=−e2−e3+e5,Dqe2=e1+ae2,Dqe3=x1mx3m−x1sx3s−ce3,Dqe4=−αce4+e5,Dqe5=−γ1e4−γ2e5−ke2.Note that the error dynamics in equation (21) can be written in the form of equation (11) with(22)A˜=0−1−1011a00000−c00000−αc10−k0−γ1−γ2,gpt,ep=00x1mx3m−x1sx3s00.Note that the termgpt,ep in equation (21) is indeed a vanishing perturbation since when the systems synchronize it follows that x1m=x1s,x2m=x2s, and x3m=x3s and therefore gpt,0=0. Consequently, the local stability of the synchronous solution in the coupled systems described by equations (18)–(20) can be determined from the condition in equation (13). In particular, we compute condition equation (13) as a function of the integration order q and the coupling strength k. The obtained results are shown in Figure 5(a), where the blue region corresponds to values of k and q for which condition (13) is not satisfied and thus the synchronous solution is unstable, whereas on the white region, condition (13) is satisfied and then the synchronous solution is expected to be stable. For the sake of comparison, Figure 5(b) shows the overlap of Figure 4 with Figure 5(a). It can be seen that there is a good agreement between the numerical and the analytical results.Figure 5
(a) Stability analysis from linear error dynamics for the coupled Rössler system shown in equation (21). (b) Overlap of the analytical result in the synchronization map shown in Figure 4(a).
(a)(b)Remark 1.
In the previous analysis, we have considered the case where the integration order of the systems and the dynamic coupling are the same. However, we also have conducted a numerical study in which the integration orders are different. In particular, we have numerically integrated equations (18)–(20) using the parameter values considered before, as a function of the integration orders of the master and slave systems, while the integration order of the dynamic coupling remains fixed. The integration orders qm,qs are varied in the interval 0.98≤qj<1, j=m,s, and considering the integration order qd=0.985. The obtained results are shown in Figure 6(a) where the blue areas correspond to synchronization (S=3). From the obtained results, it is clear to see that the integration order of the master and slave systems should be almost the same to observe a synchronized behavior and that larger differences are tolerated as long as the integration order of both systems approaches to one.
On the other hand, Figure6(b) shows the obtained results for the case that only the integration order of the dynamic coupling is varied, while the oscillators are assumed to have integer order, i.e., qm=qs=1. In this case, equations (18)–(20) are numerically integrated by varying the integration order qd of the dynamic coupling in the interval 0.8<qd<1 and the corresponding synchronization index, given by equation (14), is calculated. The obtained results are shown in Figure 6(b) for two different coupling forces, where the choice of these k values corresponds to those reported in [11] and those obtained in Figure 4(b).Figure 6
(a) Synchronization map as a function of integration ordersqm and qs of master and slave Rössler systems, equations (18)–(20), and considering a fixed integration order of qd=0.985 in the dynamic coupling. (b) Synchronization map for system equations (18)–(20) as a function of the integration order qd of the dynamic coupling, while considering that the master and slave oscillators have fixed integer order, i.e., q=1.
(a)(b)


## 5. Application Example 2: Chua System

If the Chua system described in [37] is modeled with derivatives of fractional order, then the system described in equation (23) is analyzed.(23)Dqx1=σx2−x1−ϕx1,Dqx2=x1−x2+x3,Dqx3=−βx2,where Dq is the fractional-order derivative by Caputo’s definition, xi,i=1,2,3 denotes the system state variables, σ,β are constants of the Chua circuit, and ϕx1 is a nonlinear function defined in equation (24), with constant values a,b.(24)ϕx1=bx1+12a−bx1+1−x1−1.In the same way as for the Rössler fractional-order model, the system described in (23) is analyzed to identify the minimum fractional order that can be modeled without stabilizing the dynamics, namely, qc=0.9541, since σ=10,β=14.87,a=−1.27,b=0.68. This result is confirmed by the bifurcation diagram as shown in Figure 7, where the local maxima in x1 are plotted as a function of the integration order variation. The inset shows the attractors and the corresponding integration order.Figure 7
Bifurcation diagram of the Chua system by modifying the integration order, equation (23). The abscissa axis indicates that the three oscillator state variables have the same integration order.


### 5.1. Dynamically Coupled Fractional-Order Chua Systems

In the same way as for the Rössler system described above, a pair of Chua oscillators [48] coupled by dynamical connections and defined by the system of equations (25)–(27) is considered.(25)master systemDqmx1m=σx2m−x1m−ϕx1m,Dqmx2m=x1m−x2m+x3m,Dqmx3m=−βx2m,(26)slave systemDqsx1s=σx2s−x1s−ϕx1s,Dqsx2s=x1s−x2s+x3s,Dqsx3s=−βx2s−h1,(27)dynamic couplingDqdh1=−αch1+h2,Dqdh2=−γ1h1−γ2h2−kx3m−x3s,where xm,s denotes the state variables of the master and slave systems, respectively, h1,2 is the states of the dynamic coupling, ϕx1m,s is the nonlinear function defined in equation (28), and k is the coupling force between the systems.(28)ϕx1i=bx1i+12a−bx1i+1−x1i−1,fori=m,s.As reported in [11], the following values are used in this study for all analyzes developed; σ=10,β=14.87,a=−1.27,b=−0.68 for αc=1,γ1=k,γ2=k/5, and initial conditions xim,xis,hi=0.1,0.1,0.1,−0.2,−0.1,−0.1,0,0. On the other hand, the coupling strength is varied in the interval 0≤k≤60 with Δk=0.5941 and the integration order is varied in the interval 0.965≤q≤1 for step size.To identify synchronization regions in the coupled Chua systems of fractional order, an analysis of the coupling force as a function of the integration order is developed. The obtained results are shown in Figure8. Note that the color map represents the value of the metric S and that the darkest shade of blue represents S=3, which means that the systems have reached full synchronization.Figure 8
Numerical result equations (25)–(27). (a) Limit behavior of the system as a function of coupling strength k and integration order q. (b) Synchronization index equation (14) as a function of integration order, for a fixed k=13 and k=20. A value of S=3 indicates complete synchronization.
(a)(b)Analogous to the stability analysis performed for the coupled Rössler system, the dynamic error model of the coupled pair of Chua oscillators described in equations (25)–(27) is described by(29)Dqe1=σe2−e1−ϕx1m−ϕx1s,Dqe2=e1−e2+e3,Dqe3=−βe2−h1,Dqe4=−αe4+e5,Dqe5=−γ1e4−γ2e5−ke3.Due to the nature of the nonlinearity of the Chua circuit, it is not possible to perform the same analysis as in the Rössler system; instead, it is necessary to use the Jacobian of the error model, equation (30), evaluated in one of the equilibrium points of the system. Since the Chua circuit has symmetric equilibrium points located at E1=−1.841,0.0004474,2.179 for the previously defined values, the choice of one of these points does not affect the analysis.(30)J˜=−σ−σpxσ0001−11000−β010000−αc1000−γ1−γ2,where px=1/2bsignx−1+b−1/2asignx−1−a.After defining the system shown in equation (30), it is possible to perform the stability system analysis, where the stability of the dynamic model of the Chua coupling error is sought by modifying the integration order and the coupling force. The analytical result is shown in Figure 9(a) and then compared with the map obtained from the time series analysis in Figure 8(b). As with the Rössler case, the analytical result is able to describe the boundary at which the system is unstable.Figure 9
(a) Stability analysis from linear error dynamics for the coupled Chua system as shown in equation (29). (b) Overlay of the analytical result in the synchronization map as shown in Figure 8(a).
(a)(b)Remark 2.
Similarly, to the case of the Rössler systems, we have also investigated the onset of synchronization in the fractional-order Chua oscillators as a function of the derivative orders in the oscillators while keeping the coupling system at fixedqd, resulting in the map as shown in Figure 10(a). The results from the analysis of the behavior of the Chua systems under the dynamic coupling integration order variation are shown in Figure 10(b), where the integration orders of the oscillators remain fixed at qm,s=1.Figure 10
Synchronization map as a function of integration orders for the Chua systems, equations (25)–(27). (a) The integration order of the master and slave is varied, whereas the integration order of the dynamic coupling remains fixed. (b) The order of the dynamic coupling is varied, while the master and slave systems have fixed integer order q=1.
(a)(b)


## 5.1. Dynamically Coupled Fractional-Order Chua Systems

In the same way as for the Rössler system described above, a pair of Chua oscillators [48] coupled by dynamical connections and defined by the system of equations (25)–(27) is considered.(25)master systemDqmx1m=σx2m−x1m−ϕx1m,Dqmx2m=x1m−x2m+x3m,Dqmx3m=−βx2m,(26)slave systemDqsx1s=σx2s−x1s−ϕx1s,Dqsx2s=x1s−x2s+x3s,Dqsx3s=−βx2s−h1,(27)dynamic couplingDqdh1=−αch1+h2,Dqdh2=−γ1h1−γ2h2−kx3m−x3s,where xm,s denotes the state variables of the master and slave systems, respectively, h1,2 is the states of the dynamic coupling, ϕx1m,s is the nonlinear function defined in equation (28), and k is the coupling force between the systems.(28)ϕx1i=bx1i+12a−bx1i+1−x1i−1,fori=m,s.As reported in [11], the following values are used in this study for all analyzes developed; σ=10,β=14.87,a=−1.27,b=−0.68 for αc=1,γ1=k,γ2=k/5, and initial conditions xim,xis,hi=0.1,0.1,0.1,−0.2,−0.1,−0.1,0,0. On the other hand, the coupling strength is varied in the interval 0≤k≤60 with Δk=0.5941 and the integration order is varied in the interval 0.965≤q≤1 for step size.To identify synchronization regions in the coupled Chua systems of fractional order, an analysis of the coupling force as a function of the integration order is developed. The obtained results are shown in Figure8. Note that the color map represents the value of the metric S and that the darkest shade of blue represents S=3, which means that the systems have reached full synchronization.Figure 8
Numerical result equations (25)–(27). (a) Limit behavior of the system as a function of coupling strength k and integration order q. (b) Synchronization index equation (14) as a function of integration order, for a fixed k=13 and k=20. A value of S=3 indicates complete synchronization.
(a)(b)Analogous to the stability analysis performed for the coupled Rössler system, the dynamic error model of the coupled pair of Chua oscillators described in equations (25)–(27) is described by(29)Dqe1=σe2−e1−ϕx1m−ϕx1s,Dqe2=e1−e2+e3,Dqe3=−βe2−h1,Dqe4=−αe4+e5,Dqe5=−γ1e4−γ2e5−ke3.Due to the nature of the nonlinearity of the Chua circuit, it is not possible to perform the same analysis as in the Rössler system; instead, it is necessary to use the Jacobian of the error model, equation (30), evaluated in one of the equilibrium points of the system. Since the Chua circuit has symmetric equilibrium points located at E1=−1.841,0.0004474,2.179 for the previously defined values, the choice of one of these points does not affect the analysis.(30)J˜=−σ−σpxσ0001−11000−β010000−αc1000−γ1−γ2,where px=1/2bsignx−1+b−1/2asignx−1−a.After defining the system shown in equation (30), it is possible to perform the stability system analysis, where the stability of the dynamic model of the Chua coupling error is sought by modifying the integration order and the coupling force. The analytical result is shown in Figure 9(a) and then compared with the map obtained from the time series analysis in Figure 8(b). As with the Rössler case, the analytical result is able to describe the boundary at which the system is unstable.Figure 9
(a) Stability analysis from linear error dynamics for the coupled Chua system as shown in equation (29). (b) Overlay of the analytical result in the synchronization map as shown in Figure 8(a).
(a)(b)Remark 2.
Similarly, to the case of the Rössler systems, we have also investigated the onset of synchronization in the fractional-order Chua oscillators as a function of the derivative orders in the oscillators while keeping the coupling system at fixedqd, resulting in the map as shown in Figure 10(a). The results from the analysis of the behavior of the Chua systems under the dynamic coupling integration order variation are shown in Figure 10(b), where the integration orders of the oscillators remain fixed at qm,s=1.Figure 10
Synchronization map as a function of integration orders for the Chua systems, equations (25)–(27). (a) The integration order of the master and slave is varied, whereas the integration order of the dynamic coupling remains fixed. (b) The order of the dynamic coupling is varied, while the master and slave systems have fixed integer order q=1.
(a)(b)


## 6. Discussion

From the results presented in this study, it is possible to confirm the research question formulated in the introduction of this study, according to which the use of dynamic couplings in chaotic systems of fractional order is able to induce complete synchronization, as it has been reported in their counterparts of integer order. Likewise, it is noteworthy to mention that the transition from unsynchronized behavior to synchronization is not abrupt, since there exists a region where the coupled systems may exhibit some sort of intermittency phenomenon. These areas are indicated by the blurred areas in Figures4 and 8.If the oscillators are modeled with fractional derivatives but a fixed integration order is maintained in the coupling system, the desired synchronous behavior is achieved only under the condition that both oscillators have the same integration order. Small variations in the integration order in some of the models cause the systems to lose their synchrony, as shown in Figures6(a) and 10(a). In contrast, the synchronization seems to have some robustness against variations in the integration order of the dynamic coupling, provided that the oscillators have the same integration order, as shown in the numerical results presented in Figures 6(b) and 10(b).It is also noteworthy that in theq,k plane there exist regions where the systems are easier to synchronize. This is explained by the bifurcation diagrams of the isolated oscillators, given in Figures 3 and 7, where the modification in the integration order causes important qualitative changes in the dynamics of the system, where both models are able to present chaotic or periodic behavior, for a set of parameters where the integer-order dynamics is always chaotic, only due to the modification of the derivative order. This is not only an indication that the synchronization between the systems requires lower coupling forces for periodic and quasiperiodic behaviors but also a clear indication that the change in the integration order can be associated with a modification of the vector field, which can also be achieved in the integer-order system by modifying the system parameters [17].The local stability analysis is in good agreement with the numerical analysis, as shown in Figures5 and 9. It should be noted, however, that the stability conditions are only necessary conditions. This result is similar when talking about the stability of the equilibrium points in a chaotic system, where obtaining unstable saddle points with index 2 favors the occurrence of chaotic behavior but does not guarantee the occurrence of a strange attractor in the system [49].It is worth mentioning that in the cases of analysis where the equations to solve do not have the same derivative order, the algorithm proposed by Petráš is implemented [50]; otherwise, the Adams–Bashforth–Moulton (ABM) method is used [46], which is a generalization of the classical ABM integrator that is well known in the resolution of first-order switching system problems [40, 51].Notice that the results presented here have been obtained under the assumption of identical oscillators. It is still necessary to determine the robustness of the dynamic coupling against parameter mismatches or external disturbances in the oscillators.


## 7. Conclusions

We have analyzed the onset of synchronization in fractional-order chaotic systems interacting via a linear time-invariant dynamic coupling, which also is described by fractional derivatives. The obtained results have shown the ability of the dynamic coupling to induce synchronization in the systems, and the strong influence of the integration order on the onset of synchronized behavior has been demonstrated.Among the observed limitations is that dynamic coupling is sensitive to variations in the integration order of the systems. Therefore, the master and slave systems should have the same integration order. Furthermore, it has been shown that the linearization approach used here to study the local stability of the synchronous solution only provides necessary conditions. Further investigation is needed to derive stronger stability conditions. Perhaps the use of transverse Lyapunov exponents can solve this problem. It remains as future work to extend these results to the bidirectional case and also to the case of networks. Also, it would be interesting to investigate whether any emergent behavior or other types of synchronous behaviors can occur in the coupled systems depending on the integration order of both the systems and the coupling.Finally, we would like to point out that the results presented here apply to the Caputo definition of fractional derivative. Moreover, these results were obtained using two different numerical approaches, namely, the ABM and the Petráš integrator method. As future work, we plan to compare these results with different fractional operators, such as the Riemann–Liouville operator or the Atangana–Baleanu operator, to name a few.



---

*Source: 7242253-2021-12-24.xml*