A new multistage numerical method based on blending a Gauss-Siedel relaxation method and Chebyshev pseudospectral method, for solving complex dynamical systems exhibiting hyperchaotic behavior, is presented. The proposed method, called the multistage spectral relaxation method (MSRM), is applied for the numerical solution of three hyperchaotic systems, namely, the Chua, Chen, and Rabinovich-Fabrikant systems. To demonstrate the performance of the method, results are presented in tables and diagrams and compared to results obtained using a Runge-Kutta-(4,5)-based MATLAB solver, ode45, and other previously published results.

1. Introduction

Systems exhibiting chaotic behaviour have received increasing attention from various researchers in recent years due to the challenge associated with computing their solutions and their applications in a number of areas in science such as in electrical circuits, lasers, fluid dynamics, mechanical devices, population growth, and many other areas. Chaotic behavior was first observed by Lorenz [1] in 1963 in a system of ordinary differential equations modelling weather phenomena.

Chaotic systems are complex dynamical systems which are characterised by rapidly changing solutions and high sensitivity to small perturbations of the initial data and hence computing their solutions has proven to be challenging. A chaotic system consists of only one positive Lyapunov exponent. Chaotic systems with at least two positive Lyapunov exponents are said to be hyperchaotic. Hyperchaotic systems generally have more complex dynamical behaviours than ordinary chaotic systems. The concept of hyperchaos was first introduced by Rössler [2] in the ordinary differential equations for modeling chemical reactions. The quest for finding accurate and efficient methods for solving chaotic systems is what drives the research of scientists who specialise in the development of new solution methods or seek to optimize existing methods.

For rapidly oscillating dynamical systems like the chaotic and hyperchaotic systems, standard analytical iterative methods are not guaranteed to give solutions valid globally in time. Recent research efforts have been able to overcome this lack of global convergence by modifying the standard analytical iterative schemes. This is achieved by implementing the schemes on sequences of subintervals whose union make up the domain of the underlying problem. These modified methods are known as multi-stage (or piece-wise) methods. Examples which have recently been applied to chaotic systems include the multi-stage Adomian decomposition method [36], multistage homotopy analysis method [7, 8], multi-stage differential transformation method [911], multi-stage variational iteration method [1214], multistage homotopy perturbation methods [1518]. The analytical multi-stage approaches are limited in their applications because they seek to obtain explicit analytical solutions at each of the multiple intervals. This process involves time-consuming and tedious computational operations and if too many small intervals are considered, as may be the case when dealing with highly oscillatory systems, the analytical integration process will be computationally too costly and impossible to resolve even with the use of symbolic scientific software. To overcome the inefficiency of analytical approaches in solving chaotic systems, attempts have recently been made to use the multi-stage implementation idea on some numerical methods of solution. Examples of the multi-stage numerical methods that have recently been reported for solving chaotic systems include the piecewise-spectral parametric iteration method [19] and the piecewise successive linearization method (PSLM) [2022].

In this work we present a new multi-stage iterative scheme which is based on blending a Gauss-Seidel type relaxation method with spectral collocation integration. The new method is called the multi-stage spectral relaxation method (MSRM). The MSRM is based on simple decoupling and rearrangement of the governing IVPs and numerically integrating the resulting equations in multiple intervals. We examine the applicability of the proposed MSRM in hyperchaotic chaotic systems of IVPs which include the Chua, Chen, and Rabinovich-Fabrikant systems. The computed numerical simulations are presented and it is noted that the method is accurate, efficient, and very easy to implement because its algorithm is easy to derive as it does not require any linearization of the original equations. The numerical results are also compared with some Matlab built-in Runge-Kutta solvers and good agreement is observed.

2. Multistage Spectral Relaxation Method

In this section, we give a brief description of how the multi-stage spectral relaxation method (MSRM) algorithm is developed for the solution of common chaotic systems governed by nonlinear systems of first order IVPs. Consider a chaotic system defined by nonlinear first order differential equations in the form subject to the initial conditions where are the unknown variables and are the corresponding initial conditions, are known constant input parameters and is the nonlinear component of the th equation and the dot denotes differentiation with respect to time .

Following [2022], we begin by decomposing the interval of integration into nonoverlapping intervals where and . Let be the solution of (2.1) in the first subinterval and be the solutions in the subsequent sub-intervals (). Equation (2.2) is used as the initial condition for obtaining the solution in the first subinterval . After that we use the continuity condition between neighbouring sub-intervals to obtain the initial conditions for solving (2.1) in the rest of the sub-intervals. Thus, in each interval we must solve subject to where is the Kronecker delta. We remark that, for the convenience of developing the MSRM scheme, it is necessary to express the governing equation in the form (2.3).

The MSRM algorithm uses the ideas of the Gauss-Siedel method to decouple systems of equations then apply the Chebyshev spectral collocation method to discretize and solve the resulting decoupled subsystems. We illustrate the development of the MSRM for a four equation IVP system (with ).

In the interval , the proposed MSRM iteration scheme is given as subject to the initial conditions where is the estimate of the solution after iterations. We remark that the iteration scheme (2.5)–(2.8) is only suitable for IVP systems in which the nonlinear term does not contain the variable in the th equation. This is usually the case for most of the well-known chaotic and hyperchaotic systems that have been reported in literature.

A suitable initial guess to start the iteration scheme (2.5)–(2.8) is one that satisfies the initial condition (2.9). A convenient choice of initial guess that was found to work in the numerical experiments considered in this work is

The system of (2.5)–(2.8) form a decoupled set of linear ordinary differential equations whose solutions can easily be found analytically using standard techniques of solving differential equations. However, since the solutions have to be computed in intervals which, to ensure accuracy and convergence, must be chosen to be very small, it may not be practical to seek analytical solutions. For this reason, we use a numerical approach based on pseudo-spectral methods to solve the equations. Spectral methods are powerful tools for solving ODEs if the physical domain is simple and the solution is smooth. They have the distinct advantage of having spectral accuracy and thus are significantly more accurate than related numerical methods such as finite differences and finite elements. Details of properties of spectral methods can be found in books by Canuto et al. [23], Fornberg [24], and Trefethen [25].

We use the Chebyshev spectral method to solve (2.5)–(2.8) on each interval . We first transform the region to the interval on which the spectral method is defined by using the linear transformation in each interval for . After the transformation, the interval is discretized using the Chebyshev-Gauss-Lobatto collocation points [23, 25] which are the extrema of the th order Chebyshev polynomial

The Chebyshev spectral collocation method is based on the idea of introducing a differentiation matrix which is used to approximate the derivatives of the unknown variables at the collocation points as the matrix vector product where and are the vector functions at the collocation points . The entries of the Chebyshev differentiation matrix are defined as with

Applying the Chebyshev spectral collocation method in (2.5)–(2.7) gives with where is an identity matrix of order . Thus, starting from the initial approximation (2.10), the recurrence formula can be used to obtain the solution in the interval . The solution approximating in the entire interval is given by

3. Numerical Examples

In this section, we apply the proposed MSRM to systems of IVPs with chaotic behavior to illustrate its effectiveness. In particular, we consider the hyperchaotic Chua, Chen, Rabinovich-Fabrikant systems. The results obtained are compared to results obtained by the built-in Matlab solvers, ode45.

3.1. Chua System

The Chua system was originally created as an electic circuit by Chua in 1983 [26]. The system is a set of equations with a smooth nonlinearity given by

Based on the Chua oscillator, Rech and Albuquerque [27] constructed a new four-dimensional system by introducing a fourth variable which is an adequate feedback controller to the third equation in system (3.1) to obtain

When , the system (3.2) has two positive Lyapunov exponents and hence exhibits a hyperchaotic behavior [28]. The hyperchaotic system (3.2) was solved for the initial conditions .

In this example, the parameters used in the MSRM iteration (2.5)–(2.7) are

3.2. Chen System

Qi et al. [29] presented a four-dimensional Chen system with cubic nonlinearity in each equation. The system is given by

The system exhibits hyperchaotic behavior since it contains more than one Lyapunov exponent. Qi et al. [29] analysed the basic properties of (3.4) using Lyapunov exponents and bifurcation diagrams and found that the system can generate cyclic, periodic, and chaotic attractors depending on the choice of the parameters.

Following Qi et al. [29], in implementing the MSRM on (3.4) we use the parameters . The parameter is varied between and for the cyclic, periodic and chaotic attractors respectively. The initial conditions used are .

In this example, the parameters used in the MSRM iteration (2.5)–(2.7) are

3.3. Rabinovich-Fabrikant System

The Rabinovich-Fabrikant system models the dynamical behavior arising from the modulation instability in a non-equilibrium dissipative medium [30, 31]. The system was introduced by Rabinovich and Fabrikant [32]. The Rabinovich-Fabrikant equation possesses multiple chaotic attractors. The system is described by the following set of equations: with parameters . Luo et al. [31] reported that different chaotic behaviors are observed for different values of and .

In this example, the parameters used in the MSRM iteration (2.5)–(2.7) are

Two cases are considered for this system:

Case 1. with initial conditions .

Case 2. with initial conditions .

4. Results and Discussion

We now present and discuss numerical results from the hyperchaotic systems introduced in the last section. The MSRM results are compared with those obtained from MATLAB solver ode45. The results are presented in table and graphical form. In generating all the results presented in this work, it was found that collocation points in each interval was sufficient to give good accuracy. The MSRM algorithm was run repeatedly in each interval until the norm of the difference between successive iterations was less than .

4.1. Chua System

Comparison of results and computation times for the Chua system is depicted in Table 1 and Figure 1. From observation of these, it is clear that the MSRM and ode45 yield comparable results. However, the MSRM is faster than ode45, as can be observed from the computation times. Figures 2 and 3 shows two and three dimensional phase portraits of the Chua system respectively.

4.2. Chen System

For the Chen system, computed MSRM results are compared with the piece-wise successive linearisation method (PSLM) [22] results. The results are depicted in Tables 2, 3, and 4 for the cyclic, periodic, and chaotic respectively. Comparable results are obtained in all the cases. The performance comparison for the PSLM and ode45 was made by Motsa and Sibanda [22] and a good agreement of the results was observed. The phase diagrams for the cyclic, periodic and chaotic of the Chen systems are shown by Figures 4, 5, and 6, respectively. The phase plots are similar to those of Motsa and Sibanda [22].

4.3. Rabinovich-Fabrikant System

For the Rabinovich-Fabrikant, two cases were considered depending on the value of , which gives rise to different chaotic behaviors. Results for Case 1 () are shown in Table 5 and results for Case 2 () are shown in Table 6. Furthermore, Figures 7 and 9 show graphical solutions for Case 1 and Case 2, respectively. As seen from the figures and tables the MSRM and ode45 results are comparable. The different chaotic behaviors can clearly be seen from Figures 8 and 10 which depicts the phase plots for Case 1 and Case 2, respectively. The phase diagrams are similar to those of Luo et al. [31].

5. Conclusion

In this work we have successfully computed solutions of three hyperchaotic systems namely Chua, Chen, and Rabinovich-Fabrikant systems, using a new method which is based on blending Gauss-Siedel relaxation method and the Chebyshev pseudo-spectral method. The method, called the multi-stage spectral relaxation method (MSRM) is a multi-stage method which is adapted to solve complex dynamical systems like the hyperchaotic systems. The results presented in table and graphical form are comparable to results obtained using the Runge-Kutta-(4,5)-based MATLAB built-in solver, ode45, as well as other previously published results. The results also show that the proposed MSRM is accurate, computationally efficient, and a reliable method for solving complex dynamical systems with hyperchaotic behavior.