Stability and slow-fast oscillation in fractional-order Belousov-Zhabotinsky reaction with two time scales

1. Introduction

One of the best-studied chemical oscillation systems is the Belousov-Zhabotinsky (BZ) reaction, which was elucidated by 20 chemical equations to explain the reaction mechanism and was simplified to three-variable differential equations [1, 2]. Subsequently, many works about physical and chemical mechanism, numerical simulation and experimental research on BZ reaction appeared abundantly [3-5]. After the 1990s, the slow-fast oscillation was found in many chemical reactions [6, 7], the reason was that the catalyst could make the reaction process involve in different time scales with large gap. One of the classical slow-fast oscillation, i.e. bursting oscillation was observed by Strizhak in the experiment of BZ reaction [8]. However, most of the researches on the slow-fast oscillation in chemical reaction were limited to numerical simulation and experimental investigation.

The better method to qualitatively reveal the bifurcation mechanism of the slow-fast phenomenon should be the slow-fast dynamical analysis method proposed by Rinzel [9]. By use of the method, multiple time-scale systems had been developed in overwhelming growth in last three decades [10-12]. For example, Izhikevich [13] provided all co-dimensional one bifurcation about slow-fast oscillation. Shilnikov [14] had summarized the qualitative methods on Hindmarsh-Rose model, and presented the different bursting oscillations under Hopf bifurcation. Chumakov [15] established a kinetic model of catalytic hydrogen oxidation, and studied the generation mechanism of slow-fast oscillation. Simpson [16] analyzed bursting behavior in a stochastic piecewise system, and discussed the influence of noise on slow-fast oscillation. The bifurcation patterns about neuron system with three time scales were studied by Lu [17]. Li and Bi [18] proposed enveloping slow-fast analysis method to reveal the bursting oscillation mechanism of the system with periodic excitation. Zhang [19] used differential inclusion theory to analyze the bifurcation mechanism of non-smooth systems with multiple time scales.

Actually slow-fast oscillation could be involved not only in the integer-order system but also in the fractional-order system. The reason lies in that fractional calculus would make the dynamical system more complicated, and those slow-fast oscillations have been found in various fields such as chemical, physical, mechanical, electrical, biological, economical and control engineering [20].

In the last three decades, fractional-order dynamical systems had been developed rapidly. Ahmad [21] analyzed the fractional-order Wien-bridge oscillator associated with the periodic oscillations. Ahmed et al. [22] proved the existence and uniqueness of solutions in fractional-order predator-prey system. Wang [23] gave the similarities and differences of the feedbacks between fractional-order and integer-order SDOF linear damped oscillator. Shen et al. [24-29] studied analytically and numerically the primary resonances of van der Pol (VDP) and Duffing oscillators with fractional-order derivative, and investigated the dynamical of linear single degree-of-freedom oscillator with fractional-order derivative. Liu and Duan [30] studied a fractional-order oscillator by Laplace transform and its complex inversion integral formula. Elouahab [31] extended the nonlinear feedback control in fractional-order financial system to eliminate the chaotic behaviors. Many works on stochastic dynamical system with fractional-order derivative had also been done [32-34]. Furthermore, the classical Brusselator with fractional-order derivative was investigated by Gafiychuk [35] and Li [36], in which the stability conditions and limit cycle were discussed.

However, to our best knowledge, the fractional-order BZ reaction is rarely studied, especially when the different time scales is involved. Here we will focus on the fractional-order BZ reaction with two time scales. The paper is organized as follow. In Section 2, the mathematical model of fractional-order BZ reaction is given and the bifurcation condition is discussed. The stabilities of integer-order and fractional-order systems are analyzed in details in Section 3. Then the slow-fast oscillation phenomenon and the corresponding generation mechanism are discussed by use of slow-fast dynamical analysis method in Section 4. At last, the main conclusions of this paper are made.

2. Mathematical model and bifurcation analysis

The photosensitive version of BZ reaction, i.e. Oregonator, was proposed by Seliguchi et al [37], and it was completed by the well-known reaction steps from the FKN mechanism [38], described as:

where ${}_{}{}^{*}Ru(bpy{)}_3^{2+}$ is the light-excited molecule of $Ru(bpy{)}_3^{2+}\text{.}$ The corresponding mathematical model was given in the form:

where $x$, $y$ and $z$ represent the concentrations of $HBr{O}_2$, $B{r}^{-}$, and $Ru\left(bpy\right)$ respectively, and $s$, $q$, $f$ and $k$ are the dimensionless constants related to the reaction rates. Because these parameters are closely related to the reaction condition such as temperature, pressure, feed rate, etc., the dynamical behaviors of reaction process caused by the parameter variation are very important.

The fractional-order version of BZ reaction could be established as:

in which $D^{\alpha }$ is the operator of fractional derivative with the order $\alpha \in (0,1]$. Obviously, Eq. (1) is a special case of Eq. (2). Here, we adopt Caputo’s definition as:

The equilibria $(x_0,y_0,z_0)$ of fractional-order BZ reaction can be obtained as:

where $m=\sqrt{f^2+6fq-2f+q^2+2q+1}$.

The corresponding Jacobian matrix at the equilibrium is:

The stabilities of $E_1$, $E_2$ and $E_3$ can be determined by the associated characteristic equation:

where:

For convenience, some expressions are defined as follow:

Based on Shengjin’s formulas [39], a method to solve the univariate cubic equation, three real roots can exist in Eq. (3) for $\mathrm{\Delta }\le 0$. While for $\mathrm{\Delta }>0$, there are a single real root and a pair of conjugate complex roots, shown as:

in which:

The absolute values of arguments of the complex roots are expressed as:

Based on the stability theory of fractional-order differential equation [40], the stability condition for the equilibria of Eq. (2) should be:

In the range of $\alpha \in \left(\mathrm{0,1}\right]$, the stability of Eq. (2) would keep unchanged for $\mathrm{\Delta }\le 0$, whereas it may be related to the fractional order for $\mathrm{\Delta }>0$. The critical condition for losing stability of Eq. (2) can be expressed as:

where:

Because of $\pi \alpha /2\in (0,\pi /2]$, the critical condition Eq. (4) can be ${\alpha }_0={\alpha }_{01}$.

If we fix the parameters as $s=\text{2}$, $q=\text{0.05}$ and $f=\text{1.61}$ in Eq. (2), three equilibria are ca$$lculated and denoted as follow respectively:

$E_2$ is not appropriate in practice and $E_3$ is always unstable, so that only the equilibrium point $E_1$ will be discussed for practical significance. The coefficients of Eq. (3) can be obtained:

and:

The double-parameter bifurcation diagram of Eq. (2) with respect to the parameters $k$ and $\alpha$ is plotted in Fig. 1(a), where the equilibria are stable in region (I). When the parameters pass across the critical curve into region (II), the equilibria will become unstable, and a stable limit cycle will appear in the system. It is obviously that the parameter range of stable equilibria is bigger than stable limit cycles. However, if $0<k\ll 1$, the parameter interval of stable limit cycle is larger than stable equilibria. The phase diagram of the stable limit cycle for $k=\text{0.15}$ and $\alpha =\text{0.995}$ is presented in Fig. 1(b).

Fig. 1The bifurcation diagram and the phase diagram of limit cycle: a) the bifurcation diagram with respect to parameters k and α; b) the phase diagram for the parameters k=0.15 and α=0.995

a)

b)

3. Stability analysis of the integer-order and fractional-order system

The parameter $k$ is closely related to the stability of the system (2) although it is irrelevant to the equilibrium point of this system. Therefore, the variation of the parameter $k$ will result into the change of the dynamical behaviors. The stability of the system for $\alpha =\text{1}$ and $\alpha =\text{0.95}$ when $\text{0.0005}<k<\text{1.1}$ will be analyzed in details in the following parts.

Considering $\alpha =1$, the values of other parameters are the same as those in Section 2. The real and imaginary parts of the eigenvalues of Eq. (2) are plotted in Fig. 2, denoted by the solid line and stars respectively. For $k\in \text{(0.0005,}\text{0.00375)}$, there are two positive and one negative real eigenvalues, so that $E_1$ is unstable. A negative real eigenvalue and two complex conjugate roots are found for $k\in \text{(0.00375,}\text{1.1)}$. Here we would like to point that the real part of complex conjugate roots undergoes variation with the increase of the parameter $k$. For $k<\text{0.1782}$, the real part of conjugate complex roots is larger than zero, while it is less than zero for $k>\text{0.1782}$. This means that the equilibrium is unstable for $k\in \text{(0.00375,}\text{0.1782)}$ and stable for $k\in \text{(0.1782,}\text{1.1)}$. Therefore, Hopf bifurcation happens at the critical parameter value $k=\text{0.1782}$ denoted by point A in Fig. 2(a). On the other hand, the necessary condition of Hopf bifurcation is that a pair of pure imaginary roots appear in the system, which can be obtained by $a_1>0$, $a_2>0$ and $a_1{a}_2-a_0=0$. Because there is stable limit cycle when $k<\text{0.1782}$ and stable equilibrium for $k>\text{0.1782}$, the Hopf bifurcation at critical point A is subcritical. This means that stable periodic reaction could exist when the reaction rate $k$ is small enough, while large $k$ will make the stable periodic reaction disappear, and the concentrations $x$, $y$ and $z$ will approach constants.

For $\alpha =\text{0.95}$, there are a negative real eigenvalue ${\lambda }_1$ and a pair of complex conjugate roots ${\lambda }_{\mathrm{2,3}}$, and all the arguments are shown in Fig. 3. The argument of real eigenvalue ${\lambda }_1$ is always $\pi$, and the arguments of the complex conjugate roots may vary with the increase of parameter $k$. The equilibria are unstable for $k\in \text{(0.0005,}\text{0.1049)}$ because of $\left|arg\left({\lambda }_{\mathrm{2,3}}\right)\right|<\alpha \pi /2$. And they are locally asymptotical stable for $k\in \text{(0.1049,}\text{1.1)}$ because the stable condition $\left|arg\lambda \right|>\alpha \pi /2$ is met. The critical point of losing stability takes place at $k=\text{0.1049}$, which can also be verified by Fig. 1.

Comparing the abovementioned two situations, it could be found that the stabilities of the fractional-order and integer-order systems are uniform in the most cases. However, there is a parameter interval, $k\in$(0.1049, 0.1782), in which the stabilities of the two systems are totally different. In this interval, $E_1$ is unstable for $\alpha =\text{1}$, while it is stable for $\alpha =\text{0.95}$. The numerical simulations of time history for $k=\text{0.15}$ are plotted in Fig. 4, which coincide with the abovementioned theoretical analysis. Furthermore, with the decrease of fractional order $\alpha$, the interval length of different stability may become larger and larger, which can be found from Fig. 1.

Fig. 2The stability of the integer-order system: a) The real and imaginary parts of the eigenvalues of Eq. (2) for α=1; b) The enlarge figure near k= 3.75×10-3

a)

b)

Fig. 3The critical value and the absolute values of eigenvalue’s argument for α= 0.95

Fig. 4The stability of Eq. (2): a) the time history of the fractional-order system; b) the time history of the integer-order system

a)

b)

4. Fold/Fold slow-fast oscillation and bifurcation mechanism

Considering $k=\text{0.001}$, Eq. (2) may involve in two time scales, so that the whole system may behave in the typical slow-fast phenomenon. For $\alpha =\text{0.95}$, the periodic slow-fast oscillation appears in the whole system. The corresponding phase diagram and time history are plotted in Fig. 5. In the periodic process, there are twice instantaneous jumping behaviors, denoted by the arrows in Fig. 5(a). Furthermore, the fast and slow variables present different dynamical features. For the fast variable $x$, the instantaneous jumping behaviors form the spiking state, and the quiet state takes up the most time in the periodic oscillation, which can be found in Fig. 5(b). While the instantaneous jumping phenomenon doesn’t happen in the slow variable $z$, and it will change uniformly, which is presented in Fig. 5(c). Here we would like to point that the periodic oscillation is stable, which is produced by subcritical Hopf bifurcation. The time histories under different initial values can illustrate the stability of the slow-fast periodic oscillation, as shown in Fig. 6.

Fig. 5The periodic slow-fast oscillation for α=0.95

a) The phase diagram

b) The time history for the fast variable $x$

c) The time history for the slow variable $z$

Fig. 6The limit cycle for α=0.95 with different initial values (x*,y*,z*)

a)$(x^{*},y^{*},z^{*})=\text{(8,}\text{1.2254,}\text{3.5)}$

b)$(x^{*},y^{*},z^{*})=\text{(16,}\text{1.2254,}\text{2)}$

In order to reveal the generation mechanism, we will analyze the system by use of the bifurcation theory. Obviously, the whole system Eq. (2) can be divided into a fast subsystem (FS) and a slow subsystem (SS). The FS is given by the fast variables $x$ and $y$, while the SS is modeled by the slow variable $z$. Slow variable is taken as the bifurcation parameter of the fast variables. The FS can be written as:

The equilibria of FS can be determined by:

or

Letting:

one could obtain the extreme points. The extreme points of the equilibrium line is used to analyze the bifurcation of FS. For $q=\text{0.05}$, $f=\text{1.61}$ and $s=\text{2}$, we can obtain two extreme points, denoted by LP1(2.7425, 3.1569) and LP2(8.6777, 3.8461) respectively. Then the equilibrium line of FS is divided into three branches and plotted in Fig. 7, where the points on the branch (I) and (III) are stable nodes, and the ones on the branch (II) are unstable saddle points. The details can be verified by Fig. 8, where the eigenvalues of the equilibrium curve are presented. Therefore, LP1 and LP2 are the critical points of the Fold bifurcation of FS. Here we would like to point out that the equilibria are either nodes or saddles under the taken parameter condition. The imaginary parts of eigenvalues are always zero, so that the stabilities of FS may keep unchanged with the variation of the fractional order.

Fig. 7Bifurcation diagram of FS

Fig. 8The eigenvalues of FS in Fig. 7

a) The eigenvalues for branch (I)

b) The eigenvalues for branch (II)

c) The eigenvalues for branch (III)

Because the system possesses fast and slow subsystems, the slow-fast dynamical analysis method is used here. By overlapping the bifurcation diagram of FS with the phase diagram of the whole system, Fig. 9 is obtained, which can be used to explain the generation mechanism of the periodic slow-fast oscillation.

Fig. 9Overlapping of phase diagram of the fractional-order system with bifurcation diagram of FS for α=0.95

Now we describe the periodic oscillation in details. The trajectory beginning at point D may keep in quiescent state (QS), because the trajectory may move slowly along the stable equilibrium line branch (I) of FS. The QS will be interrupted at point A, where the trajectory moves to the minimal value, i.e. the critical point LP1 of Fold bifurcation of FS. At the same time, the system may be attracted by the stable attractors on branch (III), which results into the jumping phenomenon called as spiking state (SP). When the trajectory jumps to point B, the system response may enter QS again, shown as the slow movement along the stable equilibrium line branch (III). When the system response arrives at the critical point LP2 of FS, the QS terminates and is followed by the SP characterized by jumping to point D because of the attraction from the stable equilibrium line branch (I). The whole procedure forms one period of oscillation. In a word, the twice Fold bifurcations result into twice transitions between the QS and FS, so that the periodic oscillation should belong to Fold/Fold type slow-fast oscillation.

Fig. 10The generation mechanisms of periodic slow-fast oscillation for different fractional orders

a)$\alpha =\text{1}$

b)$\alpha =\text{0.85}$

Furthermore, the effect of fractional order on the periodic slow-fast oscillation is compared here. The generation mechanisms for $\alpha =\text{1}$ and $\alpha =\text{0.85}$ are shown in Fig. 10. From the figure it could be observed that the generation mechanisms of the periodic oscillation are almost the same with the variation of the fractional order. The primary reason should be that the type of the equilibrium point of the subsystem keeps unchanged with the variation of the fractional order. However, the oscillation period may become much longer with the decrease of fractional order $\alpha$. The differences between the periods of integer-order and fractional-order systems lie in the fact, power law stability $t^{-\alpha }$ [41] is used to define the asymptotical stability of fractional-order system instead of traditional exponential stability. The details associated with the time histories can be found in Fig. 11.

Fig. 11The time histories of periodic slow-fast oscillation for different fractional orders

a)$\alpha =\text{1}$

b)$\alpha =\text{0.85}$

5. Conclusions

The fractional-order BZ reaction is investigated by use of theoretical analysis and numerical simulation. Based on the stability condition of fractional-order system, the double-parameter bifurcation diagram with respect to fractional order is firstly given, and subcritical Hopf bifurcation is found in fractional-order BZ reaction. By comparing factional-order and integer-order systems, we present the parameter interval about different stability of the two systems, and find that the different features may become more obvious with the increase of fractional order.

Furthermore, the slow-fast oscillation phenomenon is firstly discussed in fractional-order BZ reaction with two time scales coupled. Based on the slow-fast dynamical analysis method, it is found that the fast subsystem possesses twice Fold bifurcations, which leads the system to jump transiently. The switch between different stable equilibrium line branches results into the transition between QS and SP. Accordingly, the slow-fast oscillation belongs to Fold/Fold type. It is also found that the effect of fractional order on trajectory shape and generation mechanism is small, because the type of the equilibrium points of FS keeps unchanged with the variation of the fractional order. However, power law stability $t^{-\alpha }$ in fractional-order system will make the period become longer.