Maysoon M. Aziz¹, Arqam S. Mahmood^2
1Department of Mathematics, College of Sciences Computer and Mathematics, University of Mosul, Mosul, Iraq.
²Directorate General of Education Nineveh, Ministry of Education, Mosul, Iraq.
DOI: 10.4236/oalib.1115716   PDF    HTML   XML   33 Downloads   102 Views   Citations

Abstract

The study and analysis of dynamical is currently important, especially in applied fields. The research presents a numerical and analytical study of the classical model for under standing the dynamic behavior of the suggested a three dimensional biomathematical autonomous system, which made up of three ordinary equations, that represent the model variables. The dynamic properties of the suggested system are verified by equilibrium points and its stability such as the roots of the characteristic equation, Routh-Hurwitz criteria, and Lyapunov function, dissipativity, bifurcation and Kaplan-York dimension. Simultaneously, the system parameters were estimated using genetic algorithm based on real EEG data to improve the models fit to the target data. All simulations were performed using MATLAB and R-K4 numerical method. The 0 - 1 test was applied to simulated system data and real EEG data, and the results were compared, the comparison showed similar characteristic in both cases confirming and supporting the system’s ability to represented complex biological data. The results indicate the system’s efficiency as an effective tools. Showed EEG system is unstable and chaotic biological system.

Keywords

Rossler System, Stability Analysis, Chaos, Lyapunov Exponents, EEG Signals Modeling

Share and Cite:

1. Introduction

The analysis of nonlinear and chaotic systems is a major topic in applied mathematics and modern physics due to their complex behavior and high sensitivity to initial conditions [1] [2]. Among well-known continuous models, the three-dimension Rössler system stands out as a simple-structured model capable of generating rich dynamics, including periodic orbits and chaotic attractors [3] [4]. This system has been widely used in the study of local stability, bifurcations, Lyapunov exponents, and chaos detection tests. It also has applications in bio signal processing and time-based modeling [5]-[8]. Furthermore, the simplicity of its equations compared to more complex systems makes it a flexible model suitable for theoretical analysis and numerical verification in many recent studies. In this research, the ordinary system is re-examined through Roth-Heuertz equilibrium point analysis, the study of the roots of the characteristic equation, the calculation of Lyapunov exponents, and the application of the (0 - 1) chaos test. With the adoption of parametric optimization techniques to increase the model’s conformity with real-world data by using genetic algorithm to estimate system parameters [9]-[16], this convergence and conformity indicate the system’s importance in linking abstract mathematical models with complex applied phenomena.

2. Model Formulation

Recently Rossler constructed the 3-D system [1]. the proposed system described by:

$\begin{array}{l}\dot{X}=-y-z\\ \dot{Y}=x+ay\\ \dot{Z}=b+z\left(x-c\right)\end{array}$ (1)

Where x, y, z are dynamical variables of the proposed system where (x = F_z: frontal midline, y = C_z: vertex, z = P_z: parietal midline) in EEG signal, and (a,b and c) are parameters of the system are estimating from real EEG data for epilepsy people taken from https://zenodo.org/records/10259996?utm_source=chatgpt.com by using Genetic Algorithm [3] as the following:

Form of Residual error as: $r_i=x_{real\left(i\right)}-x_{sim\left(i\right)}$

Fitness function = RMSE = $\sqrt{\frac{1}{N}}{\displaystyle {\sum }_{i=1}^N{\left(\text{Residual error}\right)}^2}$

General formula of fitness function as the following:

RMSE $\sqrt{\frac{1}{N}}{\displaystyle \sum }_{i=1}^N{\left(x_{EEGi}-x_{simi}\right)}^2$ , where: f_z = x, C_z = y, P_z = z.

And $F=w_x\ast F_x+w_y\ast F_y+w_z\ast F_z$ , if the weights $w_x=w_y=w_z=1$

Then $\text{RMSE}\left(a,b,c\right)=\sqrt{\frac{1}{3}+{\displaystyle {\sum }_i^3{\left(x_{EEG}-x_{sim}\right)}^2+{\left(y_{EEG}-y_{sim}\right)}^2+{\left(z_{EEG}-z_{EEG}\right)}^2}}$

Obtaining on the parameters (a = 0.25030, b = 0.25017, c = 6.5).

3. Numerical and Grapgical Analysis

3.1. Runge-Kutta (4th Order) Numerical Method

The Runge-Kutta 4^th order method can be considered one of the most widely used numerical methods for solving ordinary differential equations, due to its good accuracy, numerical stability, and ease of programmatic application. This method relies on estimating the value of the solution in the next step through a weighted average of four calculated steps within the specified time period, making it more efficient than other simple methods.

Using the time step (h), which allows for the systematic division of the time period, the four slope formulas for the method are as follows:

$K_1=f\left(t_n,y_n\right)$ . $K_3=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_2\right)$ .

$K_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1\right)$ . $K_4=f\left(t_n+h,y_n+h{k}_3\right)$ .

$Y_{n+1}=y_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)$ (2)

where (n) is number of iterations, (h) is step size.

This method was used in this study to obtain accurate numerical solutions for the trajectories of the three-part Rossler system and to clearly represent its dynamic behavior, as well as to generate data used to simulate and compare real biological systems.

3.2. Graphical Analysis

Figure 1 shows the phase space of variables over time, where the variable x representing the frontal midline in the EEG, the variable y representing the vertex, and the variable z representing the preitial midline at the proposed system (1).

Figure 2 shows the two-dimensional phase-space diagrams of the proposed system(1), through the projections (x, y), (x, z), (y, z). Traditional projections reveal the basic dynamical behavior,

3.3. Attractor of Dynamical System (1)

The chaotic attractor is the most obvious geometric representation of complex behavior in continuous nonlinear systems. It describes the evolution of time-time paths within phase space in a finite, non-periodic structure. In the three-stage Rössler system, the attractor appears as an extended spiral twist, reflecting the high sensitivity to initial conditions and the irregularity of the motion. Figure 3 illustrates the resulting attractor for the system, highlighting the characteristic geometric structure and chaotic dynamics of the paths.

4. Equilibrium Points

Equilibrium points are obtained by the algebraic system’s equations as follows: $-y-z=0$ , $x+ay=0$ , $b+z\left(x-c\right)=0$ from the first equation we get $x=-ay$

and $z=-y$ substituting in to the third equation yields getting the quadratic:

$a{z}^2-cz+b=0$ (3)

substitute values of Parameters of system(1) by: (a = 0.25, b = 0.25, c = 6.5) in

Figure 1. Show the state space of (x vs t) , (y vs t) and (z vs t) for dynamical system (1).

Figure 2. Show the phase space (x vs y), (x vs z) and (y vs z) for dynamical system (1).

Figure 3. 3D graphic shows the attractor of the dynamical system (1).

quadratic Equation (3): $0.25z^2-6.5z+0.25=0$ and solved to find the values of z, and from it we obtain the values of x and y.

$z_{1,2}=\frac{-b\pm \sqrt{c^2-4ab}}{2a}$ $\to$ $z=\frac{-b\pm \sqrt{\Delta }}{2a}$ , where $\Delta =c^2-4ab>0$ ,

$x_{1,2}=a{z}_{1,2}$ , $y_{1,2}=-z_{1,2}$

$\Delta ={\left(-6.5\right)}^2-4\times 0.25\times 0.25=42.25-0.25\approx 42$

$\sqrt{\Delta }=\sqrt{42}\approx 6.48$ $\to$ $z_1=\frac{6.5+6.48}{2\times 0.25}=25.96$ , and $z_2=\frac{6.5-6.48}{2\times 0.25}=0.04$

$x_1=a{z}_1=0.25\times 25.96=6.49$ , $x_2=a{z}_2=0.25\times 0.04=0.01$

$y_1=-z_1=–25.96$ , $y_2=-z_2=-0.04$ ,

then: E₀ = (0, 0, 0), E₁ = (6.49, −25.96, 25.96), E₂ = (0.01, −0.04, 0.04), are three equilibrium points.

5. Stability Analysis

5.1. Characteristic Equation

Linearizing system (1) around equation to determine the Jacobian matrix J, its corresponding eigenvalue ${\lambda }_i$ , $i=1,2,3$ , are found by solving the characteristic equation: $\left|J-\lambda I\right|=0$

To obtain the characteristic equation, the differential equations of the dynamical system (1) are partially derived to obtain the Jacobi matrix, definition of system (1) equations is:

$\dot{X}=s_1\left(x,y,z\right)$ , $\dot{Y}=s_2\left(x,y,z\right)$ , $\dot{Z}=s_3\left(x,y,z\right)$ and $S=\left|\begin{array}{c}{s}_1\left(x,y,z\right)\\ s_2\left(x,y,z\right)\\ s_3\left(x,y,z\right)\end{array}\right|$

$J\left(x,y,z\right)=\left(\begin{array}{ccc}\frac{\partial s_1}{\partial x}& \frac{\partial s_1}{\partial y}& \frac{\partial s_1}{\partial z}\\ \frac{\partial s_2}{\partial x}& \frac{\partial s_2}{\partial y}& \frac{\partial s_2}{\partial z}\\ \frac{\partial s_3}{\partial x}& \frac{\partial s_3}{\partial y}& \frac{\partial s_3}{\partial z}\end{array}\right)=\left(\begin{array}{ccc}0& -1& -1\\ 1& a& 0\\ z& 0& x-c\end{array}\right)$

$\to \left|J-\lambda I\right|=0$

$\to \left|\left(\begin{array}{ccc}0& -1& -1\\ 1& a& 0\\ z& 0& x-c\end{array}\right)-\left(\begin{array}{ccc}\lambda & 0& 0\\ 0& \lambda & 0\\ 0& 0& \lambda \end{array}\right)\right|=\left|\left(\begin{array}{ccc}-\lambda & -1& -1\\ 1& a-\lambda & 0\\ z& 0& \left(x-c\right)-\lambda \end{array}\right)\right|=0$

${\lambda }^3-\left(x+a-c\right){\lambda }^2+\left(ax-ac+z+1\right)\lambda -\left(az+x-c\right)=0$ (characteristic eq.)(4)

Substitute values of the parameters and variables of E₀ for dynamical system (1) in eq,(4)and getting on:

${\lambda }^3+6.25{\lambda }^2-0.625\lambda +6.5=0$ (5)

solve it to getting on: ${\lambda }_1=-6.5$ , ${\lambda }_{2,3}=0.125\mp i0.992157$ (Two imaginary roots conjugated with a positive real part and one negative real root), then E₀ of system (1) is (hyperbolic saddle-focus unstable).

Substitute values of the parameters and variables of E₁ for dynamical system (1) in eq,(4)and getting on:

${\lambda }^3-0.24037{\lambda }^2+26.96\lambda -6.480741=0$ (6)

solve it to getting on: ${\lambda }_1\approx 0.240392>0$ , ${\lambda }_{2,3}=-0.000107\mp i5.19221$ (Two imaginary roots conjugated with a negative real part and one positive real root). then E₁ of system (1) is (hyperbolic saddle-focus unstable).

Substitute values of the parameters and variables of E₂ for dynamical system (1) in (4) and getting:

${\lambda }^3+6.24{\lambda }^2-0.5825\lambda +6.48=0$ (7)

solve it to getting on: ${\lambda }_1\approx -6.49<0$ and ${\lambda }_{2,3}\approx 0.122096\mp i0.992221$ (Two imaginary roots conjugated with a positive real part and one negative real root), then E₂ of system (1) is (hyperbolic saddle-focus unstable).

The equilibrium points of system (1) are classified as hyperbolic saddle-focus, according to eigenvalues, this is the most important structural feature of system (1), we observe that all equilibrium points of the dynamic system are unstable, and initially this explains the emergence of a strange attraction and the presence of chaos.

5.2. Routh Stability

The Routh array was originally formulated by the British mathematician Edward Routh in the year 1876 [5].

Theorem

If the signs in the first column are different, then the system is unstable.

Proof: by characteristic Equation (5): a₃ = 1, a₂ = −0.24, a₁ = 26.96, a₀ = −6.48

$b_1=\frac{a_2\times a_1-a_3\times a_0}{a_2}$ , $b_2=\frac{a_1\times a_4-a_5\times a_0}{a_1}=0$ ,

$c_1=\frac{b_1\times a_0-a_2\times b_2}{b_1}$ , $c_2=\frac{b_1{a}_5-a_1{b}_3}{b_1}=0$

Since the initial column of the constructed Routh array exhibits coefficients of mixed sign, system (1) is therefore deemed unstable at E₁.

Proceeding in a similar manner, the analysis revealed that the other equilibria, E₁ and E₂, also exhibit instability, as detailed in Table 1. Therefore, the system (1) under consideration fails to be stable.

Table 1. Classification detected equilibrium points and its stability for suggested system (1).

Equilibrium points	Characteristic equation Roots 𝜆3 – (x + a – c) 𝜆2 + (ax – ac + z + 1) 𝜆 – (az + x – c) = 0	Routh stability	Hurwitz stability > 0	Lyapunov function	Type of equilibrium Point
E0(0,0,0,0)	𝜆1 ≈ –6.5 𝜆2 , 3 = 0.125 ∓ i 0.99215	a0 = 1 a1 = 6.25 a2 = –0.625 a3 = 6.5 b1 = –1.3384 b2 = 0 c1 = 6.5, c2 = 0	∆1 = 6.25 ∆2 = –10.406 ∆3 = –67.640	0	Hyperbolic saddle-focus unstable
	Unstable	Unstable	Unstable	Unstable
E1 = (6.49, –25.96, 25.96)	𝜆1 ≈ 0.2403 > 0 𝜆2 , 3 = –0.000107 ∓ i 5.19221	a0 = 1 a1 = –0.24037 a2 = 26.96 a3 = –6.48 b1 = –0.00014 b2 = 0 c1 = 1, c2 = 0	∆1 = –0.24037 ∆2 = –0.00375 ∆30.0024313	V ≈ –0.249216	Hyperbolic saddle-focus unstable
	Unstable	Unstable	Unstable	Unstable	Unstable
E2 = (0.01, –0.04, 0.04)	𝜆1 ≈ –6.49 < 0 𝜆2,3 ≈ 0.122096 ∓ i 0. 992221	a0 = 1 a1 = 6.24 a2 = –0.5825 a3 = 6.48 b1 = 17.3644 b2 = 0 c1 = 1, c2 = 0	∆1 = 6.24 ∆2 –10.115 ∆3 = –65.5439	V ≈ 0.000384	hyperbolic saddle-focus unstable
	Unstable	Unstable	Unstable	Unstable	Unstable

5.3. Hurwitz Stability

The procedure was contributed by Adolf Hurwitz, a German mathematician, in 1895 [14].

Considering the real polynomial p evaluated at the equilibrium point E₂, the Hurwitz matrix (H) given in (8) is employed in order to compute its leading principal minors.:

$H=\left[\begin{array}{ccc}{a}_1& a_3& a_5\\ a_0& a_2& a_4\\ 0& a_1& a_3\end{array}\right]$ (8)

In order for the system to attain stability, all of the successive principal minors of the Hurwitz matrix must take positive values.

${\Delta }_1\left(p\right)=\left|a_1\right|=-0.24037<0$ ,

${\Delta }_2\left(p\right)=\left|\begin{array}{cc}{a}_1& a_3\\ a_0& a_2\end{array}\right|= \left|\begin{array}{cc}-0.24037& -6.48\\ 1& 26.96\end{array}\right|=-0.00375<0$

${\Delta }_3\left(p\right)=\left|\begin{array}{ccc}-0.24037& -6.48& 0\\ 1& 26.96& 0\\ 0& -0.24073& -6.48\end{array}\right|=41.9928313-41.9904=0.0024313>0$

Accordingly, the dynamical model described by (1) proves to be unstable at E₂.

By repeating the same procedure for the remaining equilibrium points E₀ and E₁, system (1) proves to be unstable, the outcomes being tabulated in Table 1.

5.4. Lyapunov Stability Function

$V\left(x,y,z\right)=\frac{1}{2}\left(x^2+y^2+z^2\right)$ , $\dot{V}=x{s}_1+y{s}_2+z{s}_3$ , $\dot{V}>0$ .

$\dot{V}=x\left(-y–z\right)+y\left(x+ay\right)+z\left(b+z\left(x-c\right)\right)\to \dot{V}=a{y}^2+\left(x-c\right)z^2-xz+bz$ (9)

Since the equilibrium points are not locally stable, and the system’s trajectories do not converge towards these equilibrium points but rather move towards a strange attractor, we conclude that the system does not achieve asymptomatic stability (locally or globally) with respect to the equilibrium points. Therefore, a dynamic system does not stabilize at an equilibrium point but rather stabilizes at a limited strange attractor, which may be chaotic.

6. Dissipativity

Dispersion is an important property of dynamical systems because it indicates whether the magnitudes of paths in phase space are time-dependent. When a system is dispersed, the motion remains confined within a limited region of space, often leading to complex attractors and chaotic behavior. To study dispersion, we calculate the vector field divergence ∇.S as follows:

$\nabla \cdot S=\frac{\partial s_1}{\partial x}+\frac{\partial s_2}{\partial y}+\frac{\partial s_3}{\partial z}=Tr\left(J\right)$ , $\nabla \cdot S=0+a+\left(x-c\right)=a+x-c$

the system is dissipative iff $\nabla \cdot S<0$ (dissipation condition). $a+x-c<0$ $x<c-a$ , when $a=0.25$ , $x<c=6.5$ then $x<6.5-0.25$ $\to$ $x<6.25$ in E_0: $x=0<6.25$ $\to$ the system is dissipative in E₀ in E₁: $x=6.49>6.25$ $\to$ the system is unbounded (undissipated) in E₁ in E₂: $x=0.01<6.25$ $\to$ the system is dissipative in E₂ ($\nabla \cdot S<0$ ).

7. The Bifurcation

This section deals with the bifurcation analysis of system (1), bifurcation diagram is the tool generally used to classify the dynamical systems. In order to observe the bifurcation in parietal midline in EEG signal, the bifurcation analysis conducted through numerical simulation using R-K4 numerical method over parameter $c\in \left[3,9\right]$ , the numerical analysis started with initial condition (0.1, 0.1, 0.1), The numerical integration is carried out with the initial time t₀ = 0, a fixed step size of Δt = 0.01, and a final time tend = 300 s, system (1) exhibits rich behaviors as: period-doubling routes to chaos, intermittency, quasi-periodic oscillations and full-scale attractors. In this section explores those behaviors through numerical simulation using R-K4 numerical method, Figure 4 shown bifurcation behavior for proposed system (1) with respected to $c\in \left[3,9\right]$ , (t step = h = 0.01) [11] [12].

Figure 4. Shows bifurcation of dynamical system (1) with respected to c ∈ [3,9] and (h=0.01).

8. Chaos Analysis

8.1. Lyapunov Exponent

The Lyapunov exponent is one of the most important quantitative tools used in diagnosing chaotic behavior and measuring the high sensitivity exhibited by the system upon the initial conditions. The presence of a maximum positive Lyapunov exponent is a clear indicator of chaos. These exponents also illustrate the rate of divergence or convergence of adjacent paths in phase space over time. Figure 5(a) shows the convergence of the three Lyapunov exponents of the three-part Rossler system to reach stable values. Figure 5(b) also shows the change of each exponent individually as a function of the parameter (c = 6.5), revealing the regions of stability and transition to chaotic behavior within the studied domain, the values of three Lyapunov exponents are:

LE₁ = 0.0915, LE₂ = −0.0018, LE₃ = −6.2007

8.2. Zero-One Test (0-1)

The binary test (0 1) for chaos is a simple and effective numerical tool that helps to so as to discriminate between (regular dynamic, whether periodic or quasi-periodic) dynamics and chaotic dynamics, including the dynamic system (1), it works directly with the time series and does not require any phase space reconstruction:

Step(1)- compute translation variables: Given a time series $\varphi \left(j\right)$ , compute the translation variables p(n) and q(n) using $p\left(n\right)={\displaystyle {\sum }_{i=1}^n\varphi \left(j\right)\cos \left(jc\right)}$ , and $q\left(n\right)={\displaystyle {\sum }_{i=1}^n\varphi \left(j\right)\sin \left(jc\right)}$ , where c is a randomly chosen constant in $\left(0,\pi \right)$ , $n=1,2,\cdots ,L$ .

(a) (b)

Figure 5. (a) Shows the convergence of the three Lyapunov exponents of proposed system (1), (b) Shows the change of each exponent.

Step(2)- compute Mean Square Displacement (MSD): calculate M(n) the Mean Square Displacement of the (p,q) trajectory.

$M\left(n\right)=\underset{L\to \infty }{\lim }\left(\frac{1}{L}\right){\displaystyle {\sum }_{j=1}^L\left[{\left(P\left(j+n\right)-P\left(j\right)\right)}^2+{\left(q_n\left(j+n\right)-q_n\left(j\right)\right)}^2\right]}$

Wherein n takes integer values from $1,2,\cdots ,\left[L/10\right]$ .

$Vosc\left(n\right)={\left[E\left(\varphi \right)\right]}^2\times \frac{1-\cos \left(nc\right)}{1-\cos \left(c\right)}$ , where $E\left(\varphi \right)=\underset{L\to \infty }{\lim }\left(\frac{1}{L}\right){\displaystyle {\sum }_{j=1}^n\varphi \left(j\right)}$

The quantity D(n) is defined through the relation $D\left(n\right)=M\left(n\right)-Vosc\left(n\right)$

$Kcorr=Kc=\underset{n\to \infty }{\lim }\frac{\log Mc\left(n\right)}{\log \left(n\right)}$

Step(3)- Estimate Growth Rate K: Determine the asymptotic growth rate K of M(n),

Kc states: If $Kc\cong 0$ the dynamic is regular (periodic or quasi-periodic).

and if $Kc\cong 1$ then the dynamic is chaotic. with Kc takes value in the interval $\in \left[0,1\right]$ .

Through the use of the MATLAB programming environment, the following results are obtained for simulation data: $Kc=0.92695\cong 1$ Then the system is a chaotic. and we get for real data: $Kc=0.9586\cong 1$ , it is chaotic too.

The 0 - 1 test results obtained from the Rössler simulations exhibit a strong agreement with those taken from the real EEG data, demonstrating a high level of similarity in their chaotic signatures. This close correspondence confirms the capability of the Rossler system to reliably reproduce the underlying dynamical characteristics of the measured brain signals. Figure 6(a)-(b) shown the compared between real data of (EEG) signal and simulation data of Dynamics System(1).

(a) For real data (b) For simulation data

Figure 6. Shows the chaotic of real data and simulation data for proposed system.

9. Conclusion

This study demonstrated that the typical three-dimension proposed system exhibits rich and complex dynamic behavior despite its simple mathematical formulation. Analysis of equilibrium points which appeared as hyperbolic saddle-focus by means of the roots of the characteristic equation and the Routh-Hurwitz criteria showed that the system loses stability within the selected parameters, which explains its transition towards more complex motion patterns. Furthermore, the Lyapunov exponents confirmed the presence of a clear chaotic state. On the other hand, the (0 - 1) chaotic test was used, which supported the exponents results, yielding a mean value approaching one for both the real data (EEG) and the simulated data generated from the system using the numerical R-K4 method. The median value reached approximately 0.9586 for the real data and 0.92695 for the simulated data, indicating a confirmed chaotic nature in both cases. Additionally, the use of the intelligent genetic algorithm demonstrated good efficiency in estimating the system’s parameters and improving the model’s fit to the target data. In general, the results confirm that the 3-D proposed system is an effective model for studying chaos and representing nonlinear dynamical systems, with the possibility of employing it in the future in studying control and synchronization applications, biological modeling, and fractional-order or multidimensional expansions.

Acknowledgements

The authors are very grateful to the University of Mosul/College of Computers Sciences and Mathematics for their provided facilities, which helped to improve the quality of this work.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Rössler, O.E. (1976) An Equation for Continuous Chaos. Physics Letters A, 57, 397-398.[CrossRef]
[2]	Strogatz, S.H. (2018) Nonlinear Dynamics and Chaos. 2nd Edition, CRC Press.
[3]	Ott, E. (2002) Chaos in Dynamical Systems. 2nd Edition, Cambridge University Press.[CrossRef]
[4]	Sparrow, C. (1982) The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors. Springer.
[5]	(2007) ECE 680 Modern Automatic Control. Routh’s Stability Criterion, 1-6.
[6]	Wolf, A., Swift, J.B., Swinney, H.L. and Vastano, J.A. (1985) Determining Lyapunov Exponents from a Time Series. Physica D: Nonlinear Phenomena, 16, 285-317.[CrossRef]
[7]	Rosenstein, M.T., Collins, J.J. and De Luca, C.J. (1993) A Practical Method for Calculating Largest Lyapunov Exponents from Small Data Sets. Physica D: Nonlinear Phenomena, 65, 117-134.[CrossRef]
[8]	Gottwald, G.A. and Melbourne, I. (2004) A New Test for Chaos in Deterministic Systems. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460, 603-611.[CrossRef]
[9]	Gottwald, G.A. and Melbourne, I. (2009) On the Implementation of the 0-1 Test for Chaos. SIAM Journal on Applied Dynamical Systems, 8, 129-145.[CrossRef]
[10]	Ogata, K. (2010) Modern Control Engineering. 5th Edition, Prentice Hall.
[11]	MathWorks (2024) MATLAB Documentation. The MathWorks, Inc.
[12]	Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P. (2007) Numerical Recipes: The Art of Scientific Computing. 3rd Edition, Cambridge University Press.
[13]	Aziz, M.M. and Mahmood, A.S. (2023) Mathematical Model of Epidemic Disease Covid-19. AIP Conference Proceedings, 2414, Article ID: 040072.[CrossRef]
[14]	Aziz, M.M. and Mahmood, A.S. (2021) Analysis of Dynamical Behavior for Epidemic Disease COVID-19 with Application. Turkish Journal of Computer and Mathematics Education (TURCOMAT), 12, 568-577.[CrossRef]
[15]	Haupt, R.L. and Haupt, S.E. (2004) Practical Genetic Algorithms. 2nd Edition, Wiley.[CrossRef]
[16]	Aziz, M. (2023) Mathematical Model for the Effect of Buoyancy Forces on the Stability of a Fluid Flow. European Journal of Pure and Applied Mathematics, 16, 983-996.[CrossRef]