Tracking the chaotic behaviour of fractional-order Chua’s system by Mexican hat wavelet-based artificial neural network

Abstract

In this paper, a process is devised systematically to scrutinize the scrolling chaotic behaviour of fractional-order Chua's system. The process is composed of fractional Laplace transformation, artificial neural network with Mexican hat wavelet as an activation function and simulated annealing. Sequentially, the parametric expansion of fractional Laplace transform is employed to convert the governing fractional system into an ordinary differential system. Next, artificial neural network and simulated annealing approximate and optimize the attained system and produce accurate solutions. The predictability and elaboration of double scrolling chaotic structures of fractional-order Chua's system are also studied using the Lyapunov exponent and fifth–fourth Runge–Kutta method. Moreover, the mean absolute error and root mean square error are measured for the convergence analysis of the proposed scheme. On the whole, the accurate approximate solutions, the phase plots of Lyapunov exponent spectrum and bifurcation maps of the dynamical evolution of fractional Chua's system are a triumph of this endeavour.

Keywords

Fractional calculus Chua’s system chaotic dynamics bifurcation analysis

Introduction

The chaos theory seems to be everywhere; wild swirls of cigarette smoke, back-forth snaps of the flag in the wind, the behaviour of cars on an expressway, etc. are real-life examples of this theory.^1,2 Its prototype took initiative from the study of three-body motion, simulation of weather predictions, fluctuations in economic systems, to the recent analysis of nonlinear systems, existing in all the multidisciplinary areas of science. The exploration of chaotic behaviour and its patterns, using several methodologies and computer-based algorithms has been assessed extensively.³ Similarly, the applications of the chaotic systems are also widely found in electrical engineering, where Chua’s system is one of the most expedient example of electric circuits that exhibit classical chaotic behaviour.^4,5 It is a system of nonlinear ordinary differential equations, which defines the voltage across the capacitors and the electric current in the inductor of the system. Chua's system produces a double scrolling chaotic pattern, which brought about an evolution in the study of chaos theory. Due to many other interesting hidden chaotic attractors of Chua's system, it has undergone several generalizations.^6,7

Among many other methods for measuring the existence of chaos in a dynamical system, Lyapunov exponent (LE) is considered to be the most appropriate tool for this purpose. It is used to a great extent to measure the exponential divergence or convergence rate of the trajectories, starting from nearby initial points. Similarly, LE can also be used to check the stability of limit sets and sensitivity to the initial conditions of a nonlinear system. The signs of the LEs provide a qualitative picture of a system's dynamics, as found in Singh and Roy,⁸ where the (+, +, −, −) nature of LE are discussed for the hyperchaotic behaviour of a four-dimensional hyperchaotic system. A single LE characterizes a one-dimensional map into chaos, a marginally stable orbit and a periodic orbit, by retaining positive, zero and negative signs, respectively. The existence of at least one positive LE identifies the presence of chaos, with the magnitude of the exponent reflecting the phase at which the system becomes unpredictable.⁹

Modelling different real-world phenomena, using fractional definitions, has become the most highly appreciated areas of realistic sciences. This is because the nonlocal properties of fractional operator enable these differential models to condense the information, about recent and historical situations. For the last few decades, many advancements have been made in this regard to enhance the definitions and properties of fractional calculus to overcome the inadequacies, e.g. He’s fractional derivative,¹⁰ Atangana–Baleanu fractional derivative,¹¹ conformable derivative,¹² etc. Consequently, these novel aspects enrich the capabilities of fractional differential models by bringing diverse physical significances to light.^13–15 Hence, by means of different theories of fractional derivatives, the behaviours of many fractional differential equations have been studied and various techniques have been developed,^16–19 but still, there are many things that can be done in this area.

In this endeavour, we provide the study of chaotic characteristics of Chua’s system with fractional order derivative. Our main aim is to acquaint an approximating and optimizing algorithm, which enables to accurately simplify the fractional Chua's system and corroborate LE to review the chaotic attractor of the system. In this regard, a parametric expansion of fractional Laplace transform of Caputo derivative has been utilized to deal with the fractional operators of the system. This new expanded form of fractional Laplace transform has been greatly used to solve fractional order systems nowadays, as it converts the problem in the ordinary differential equation, remaining in a fractional environment.²⁰ Here, after applying the fractional Laplace transformation, artificial neural network (ANN)²¹ with simulated annealing (SA)²² and fifth–fourth Runge–Kutta (FFRK) method²³ are employed to further attain the numerical solutions of the voltages across the capacitors and electric current of the Chua's circuit. In this study, the Mexican hat wavelet (MHW)²⁴ is taken into account for the activation of neurons of ANN. MHW is classified into a continuous wavelet transforms, which divides a given continuous-time signal into different scale components. It has significant applications in the representation of continuous-time signals of fractional bands and for resolving component mixtures.²⁵ Besides, the amalgamation of ANN and SA has also been significantly utilized these days to numerically optimize a number of integers, non-integer ordinary and partial differential equations.²¹

Moreover, the remaining paper contains: the basic definitions of the Caputo fractional derivative, the parametric expansion of fractional Laplace transform and the demonstrative structure of MHW in the following section. Next, the formulation of fractional Chua's model and implementation of proposed schemes are described in detail. Furthermore, a comprehensive discussion of bifurcation investigations of the fractional Chua's system, graphical solutions of the parameters along with its phase diagrams, is imparted. Conclusively, effective outcomes of the whole attempt are transliterated.

Prologue

This section comprises the illustrative preliminaries of the paper, which include some major definitions and properties of Caputo derivative, Laplace transformation and MHWs.

Caputo fractional operator

Assume $ν > 0, m = ⌈ ν ⌉$ and $f (t) \in C^{m} [0, 1]$ , then the Caputo fractional derivative of $f (t)$ with respect to $t$ is defined as

\frac{d^{ν}}{d t^{ν}} f (t) = {\begin{array}{l} I_{t}^{m - ν} \frac{d^{m}}{d t^{m}} f (t), & m - 1 < ν \leq m, \\ \frac{d^{m}}{d t^{m}} f (t), & ν = m \in N \end{array}

(1)

where

I_{t}^{ν}

is the Riemann–Liouville integral, given as

\begin{array}{l} I_{t}^{ν} f (t) = \frac{1}{Γ (ν)} \int_{0}^{t} {(t - ψ)}^{ν - 1} f (ψ) d ψ, \\ I_{t}^{0} f (t) = f (t) \end{array}

(2)

we use the notation

D_{t}^{ν}

in replacement of

\frac{d^{ν}}{d t^{ν}}

for the Caputo fractional derivative throughout the study.

Fractional Laplace transform

The Laplace transform of Caputo fractional derivative of order $ν$ is outlined as

L [D_{t}^{ν} f (t)] = s^{ν} \hat{f} (s) - \sum_{k = 0}^{n - 1} s^{ν - k - 1} f^{(k)} (0), n - 1 \leq ν \leq n, n \in N

(3)

where

\hat{f} (s)

is the Laplace transform of

f (t)

, i.e.

L [f (t)] = \hat{f} (s) = \int_{0}^{\infty} e^{- s t} f (t) d t

. On considering

n = 1

, equation (3) gives

L [D^{ν} f (t)] = s^{ν} \hat{f} (s) - s^{ν - 1} f (0), 0 \leq ν \leq 1

(4)

MHW function

Let $ψ \in L^{2} (ℜ)$ is a continuous nonorthogonal wavelet transform, described as^24,25

ψ_{a, h} (x) = {| a |}^{- 1 / 2} ψ (\frac{x - h}{a})

(5)

for

a, h \in ℜ

, such that

a \neq 0

, then the so-called MHW is defined by

ψ (τ) = \frac{- 2}{\sqrt[4]{π} \sqrt{3 σ}} (\frac{τ^{2}}{σ} - 1) e^{(- \frac{τ^{2}}{σ^{2}})}

(6)

Here $σ$ represents the width and smoothness of MHW function, as shown in Figure 1.

Figure 1.

MHW at different values of $σ$ .

Fractional-order Chua’s model

Chua’s system is a simple electronic system which contains a nonlinear resistor (Chua’s diode) with a piecewise-linear characteristic.³ As mentioned previously, Chua's system is a paradigm for chaos and has received considerable attention for more than a decade. Its mathematical structure has been widely studied with the definitions of fractional derivatives by several well-known authors.^5,7 In general, the normalized dimensionless dynamical equation of Chua's circuit, with Caputo fractional derivative can be written as

\begin{array}{l} D_{t}^{α_{1}} y_{1} (t) = c_{1} (y_{2} (t) - y_{1} (t) - g (y_{1} (t))) \\ D_{t}^{α_{2}} y_{2} (t) = c_{2} (y_{1} (t) - y_{2} (t) + y_{3} (t)), 0 \leq t \leq p, p \in N \\ D_{t}^{α_{3}} y_{3} (t) = c_{3} y_{2} (t) \end{array}

(7)

with initial conditions

y_{1} (0) = ρ_{1}, y_{2} (0) = ρ_{2}, y_{3} (0) = ρ_{3}

(8)

where

α_{1}

α_{2}

and

α_{3}

are the fractional order of Caputo derivatives, such that

0 < α_{1}, α_{2}, α_{3} \leq 1

c_{1}

c_{2}

c_{3}

are normalized dimensionless parameters, and

g (y_{1} (t))

is a piecewise linear function. The Chua's system considered in the present study differs from the usual one, i.e. in system (7) the piecewise linear function is replaced by an appropriate cubic nonlinear function

\begin{array}{l} D_{t}^{α_{1}} y_{1} (t) = c_{1} (y_{2} (t) + \frac{y_{1} (t) - 2 y_{1}^{3} (t)}{7}) \\ D_{t}^{α_{2}} y_{2} (t) = c_{2} (y_{1} (t) - y_{2} (t) + y_{3} (t)), 0 \leq t \leq p, p \in N \\ D_{t}^{α_{3}} y_{3} (t) = c_{3} y_{2} (t) \end{array}

(9)

This generalization is found in Cafagna and Grassi,⁷ where the substitution of the continuous piecewise-linear function is made with a smooth function, such as a cubic polynomial etc. Here, we obtain different representations of solutions of the fractional-order Chua’s model (9), by using fractional Laplace transformation along with ANN and an optimization technique SA, elaborated in the next section.

Methodical scheme

In this section, we systematically illustrate the strategy that is proposed to evaluate the governing Chua's system. Firstly, a parametric expansion of the fractional-order Chua’s system into a new integer order system is carried out by employing fractional Laplace transform. Then, an ANN technique with SA is exerted to achieve an approximate analytical solution for the system. In fact, the parametric reduction of the fractional model fulfilled at the first step can further enhance the computational productivity.

Implementation of fractional Laplace transform

Consider the Caputo fractional operator $D_{t}^{α_{i}} y_{i} (t)$ for $i = 1, 2, 3$ in system (9), then its fractional Laplace transformation can be expanded as

\begin{matrix} L [D_{t}^{α_{1}} y_{1} (t)] = s^{α_{1}} {\hat{y}}_{1} (t) - s^{α_{1} - 1} y_{1} (0) \\ L [D_{t}^{α_{2}} y_{2} (t)] = s^{α_{2}} {\hat{y}}_{2} (t) - s^{α_{2} - 1} y_{2} (0) \\ L [D_{t}^{α_{3}} y_{3} (t)] = s^{α_{3}} {\hat{y}}_{3} (t) - s^{α_{3} - 1} y_{3} (0) \end{matrix}

(10)

where

{\hat{y}}_{i} (t)

represents the fractional Laplace transform of

y_{i} (t)

. Here, we assume that

0 < α_{1}, α_{2}, α_{3} \leq 1

, so the terms

α_{i}

maybe linearized by following the method defined for interval arithmetic in Ren et al.,²⁰ i.e.

s^{α_{i}} \approx α_{i} s + (1 - α_{i}) s^{α_{i}} = α_{i} s + (1 - α_{i})

(11)

On replacing equation (11) into equation (10) and using the inverse Laplace transform we attain the parametric expansions as

\begin{matrix} D_{t}^{α_{1}} y_{1} (t) = α_{1} y_{1}' (t) - (1 - α_{1}) (y_{1} (t) - y_{1} (0)) \\ D_{t}^{α_{2}} y_{2} (t) = α_{2} y_{2}' (t) - (1 - α_{2}) (y_{2} (t) - y_{2} (0)) \\ D_{t}^{α_{3}} y_{3} (t) = α_{3} y_{3}' (t) - (1 - α_{3}) (y_{3} (t) - y_{3} (0)) \end{matrix}

(12)

Taking into account equation (12), the fractional-order Chua’s system (9) is simplified to the following differential system

\begin{array}{l} α_{1} y_{1}' (t) - (1 - α_{1}) (y_{1} (t) - y_{1} (0)) = c_{1} (y_{2} (t) + \frac{y_{1} (t) - 2 y_{1}^{3} (t)}{7}) \\ α_{2} y_{2}' (t) - (1 - α_{2}) (y_{2} (t) - y_{2} (0)) = c_{2} (y_{1} (t) - y_{2} (t) + y_{3} (t)), 0 \leq t \leq p, 0 \leq α_{i} \leq 1, p \in N \\ α_{3} y_{3}' (t) - (1 - α_{3}) (y_{3} (t) - y_{3} (0)) = c_{3} y_{2} (t) \end{array}

(13)

This process of linearization is found to be significantly helpful in reducing computational costs that occur in different complicated approximations. Thenceforth, we exercise the proposed scheme to get the analytical solutions of system (13) with conditions (8), which will consequently lead to an approximate analytical solution of the fractional-order Chua’s system (9).

Approximations of MHW-based ANN

ANN predominantly approximates the functions by introducing a trial solution, which features two important parts: the first part contains an ANN with a vector, encompassing all its corresponding weights, whereas the second part must satisfy the initial/boundary conditions. Therefore, assume $\bar{y_{i}} (t; Ω_{i})$ , to be the trial solution for the problem under consideration, with $Ω_{i}$ as a weight vector and second component satisfy the conditions defined in equation (8). Hence, the mathematical formulation of the trial solution $\bar{y_{i}} (t; Ω_{i})$ for equation (13) is articulated as

\begin{matrix} \bar{y_{1}} (t; Ω_{1, i}) = y_{1} (0) + t \sum_{i = 1}^{n} λ_{i} N (℘_{1}) \\ \bar{y_{2}} (t; Ω_{2, i}) = y_{2} (0) + t \sum_{i = 1}^{n} μ_{i} N (℘_{2}) \\ \bar{y_{3}} (t; Ω_{3, i}) = y_{3} (0) + t \sum_{i = 1}^{n} η_{i} N (℘_{3}) \end{matrix}

(14)

where weight vectors

Ω_{1, i}

Ω_{2, i}

and

Ω_{3, i}

contain all the unknown weights that are to be determined and

N (℘)

is called the activation function. Here,

N (℘)

is taken as a non-orthogonal MHW function defined in equation (6), where

℘_{1}

℘_{2}

and

℘_{3}

are the linear sums that can be written as

\begin{matrix} ℘_{1} = \sum_{i = 2}^{n} u_{i} t + k_{i} \\ ℘_{2} = \sum_{i = 2}^{n} v_{i} t + l_{i} \\ ℘_{3} = \sum_{i = 2}^{n} w_{i} t + q_{i} \end{matrix}

(15)

Such that $λ_{i}, μ_{i}, η_{i}, u_{i}, v_{i}$ and $w_{i}$ denote network adaptive coefficients (or weights) and $k_{i}, l_{i}$ , and $q_{i}$ are the biases of the network, which are absorbed in weight vectors $Ω_{1, i}$ , $Ω_{2, i}$ , $Ω_{3, i}$ . Substituting equation (15) into equation (14), we get the perceptron of ANN that provides a universal approximator of continuous functions of equation (13).

Fitness function

The construction of fitness function is based on the residual error of equation (13), which, after substitutions of trial solutions, can be written as

\begin{array}{l} R_{i 1} (t, Ω_{1, i}) = \sum_{i = 1}^{n} (α_{1} \bar{y_{1}'} (t; Ω_{1, i}) - (1 - α_{1}) (\bar{y_{1}'} (t; Ω_{1, i}) - ρ_{1}) \\ - c_{1} (\bar{y_{2}} (t, Ω_{2, i}) + \frac{\bar{y_{1}} (t, Ω_{1, i}) - 2 {\bar{y_{1}}}^{3} (t, Ω_{1, i})}{7})) \\ R_{i 2} (t, Ω_{2, i}) = \sum_{i = 1}^{n} (α_{2} \bar{y_{2}'} (t; Ω_{2, i}) - (1 - α_{2}) (\bar{y_{2}'} (t; Ω_{2, i}) - ρ_{2}) \\ - c_{2} (\bar{y_{1}} (t; Ω_{1, i}) - \bar{y_{2}} (t; Ω_{2, i}) - \bar{y_{3}} (t; Ω_{3, i}))) \\ R_{i 3} (t, Ω_{3, i}) = \sum_{i = 1}^{n} (α_{3} \bar{y_{3}'} (t; Ω_{3, i}) - (1 - α_{3}) (\bar{y_{3}'} (t; Ω_{3, i}) - ρ_{3}) - c_{3} \bar{y_{3}} (t; Ω_{3, i})) \end{array}

(16)

Now, using the capability of approximation theory in mean square error (MSE) sense, the fitness/objective function is developed as

E (Ω) = \frac{1}{3 m} \sum_{j = 1}^{m} \sum_{i = 1}^{n} (R_{i j} {(t_{j}, Ω_{1, i})}^{2} + R_{i j} {(t_{j}, Ω_{2, i})}^{2} + R_{i j} {(t_{j}, Ω_{3, i})}^{2})

(17)

where

t_{j} = (t_{1}, t_{2}, …, t_{n})

are the input data. As the fitness function

E

, given in equation (17), approaches zero, the approximate solution,

\bar{y_{i}} (t; Ω_{i, j})

approaches the exact solution

y_{i} (t)

of equation (13).

Next, for the required minimum value of fitness function, we apply an SA algorithm,²² which is a popular heuristic optimizer.

SA algorithm

SA algorithm imitates the annealing process, which is followed in material sciences, where a metal is cooled and frozen into a crystalline state with the minimum energy and larger crystal sizes so as to refine the metallic structures. Here, we use SA to compute the unknown terms of the series and optimize the performance index. This process encompasses the following features:

The unknown parameters whose values are to be determined.

A conditional equation or an objective function. Its values play a crucial part in measuring the probability.

Boltzmann probability distribution, which is an exponential function that assigns the probability to each value of the conditional equation.

The detailed algorithm of SA scheme is defined as follows:

Objective function $E (Ω) = E (Ω_{1, i}, Ω_{2, i}, Ω_{3, i})$

Initialize the initial temperature $T_{0}$ and initial guess $Ω_{0, i}$

Set the final temperature $T_{f}$ and the max number of iterations $j$

Define the cooling schedule $T \to η T, 0 < η \leq 1$ while ( $T > T_{f}$ and $d < M$ )

Drawn from a Boltzmann probability distribution

Move randomly to a new location: $Ω_{d + 1} = Ω_{d} + δ$ (random walk)

Calculate $Δ E = E_{d + 1} (Ω_{d + 1}) - E_{d} (Ω_{d})$

Accept the new solution if betterif not improved

Generate a random number $r$

Accept if $P = \exp (\frac{- Δ E}{T}) > r$ end if

Update the best $Ω_{*}$ and $E_{*}$ $d = d + 1$ end while

Suppose that $Ω$ is the optimal solution of the optimization problem (17), then by substituting it into equation (14), the approximate numerical solutions of equation (13) is achieved.

Performance index

The accuracy of the algorithm is assessed in terms of the two different performance measures, namely, mean absolute error (MAE) and the root mean square error (RMSE), i.e.

L_{i, a b s} = | y_{i}_{, R K} - {\bar{y}}_{i, t r i a l} | for i = 1, 2, 3

(18)

L_{i, \infty} = \max L_{i, a b s}

(19)

L_{i, r m s} = \sqrt{\frac{1}{m} {\sum_{j = 1}^{m} | {\hat{y}}_{i, R K} - {\bar{y}}_{i, j} |}^{2}}

(20)

where

y_{i, R K}

are the solutions obtained by using FFRK method and

{\bar{y}}_{i,_{t r i a l}}

is the approximate solution obtained by the proposed method. We choose the computational domain

[0, 1]

, for the governing model; thus, the values of performance indices, based on MAE and RMSE, should be equal to zero.

Application and discussion

Now, we investigate the dynamical behaviour of system (13) by the method described in the section Methodical scheme and the classical FFRK method. Scientifically, all the numerical manipulations and constructions of figures and tables are carried out using Mathematica 11. Following the aforementioned procedure, the trial solution of equation (13) for $n = 5$ , given in equation (14), is expressed as

\begin{array}{l} \bar{y_{1}} (t; Ω_{1}) = 0.6 + t \sum_{i = 1}^{5} \frac{- 2 λ_{i}}{\sqrt[4]{π} \sqrt{3 p}} (\frac{{(u_{i} t + k_{i})}^{2}}{p} - 1) e^{(- \frac{{(u_{i} t + k_{i})}^{2}}{p^{2}})} \\ \bar{y_{2}} (t; Ω_{2}) = 0.1 + t \sum_{i = 1}^{5} \frac{- 2 μ_{i}}{\sqrt[4]{π} \sqrt{3 p}} (\frac{{(v_{i} t + l_{i})}^{2}}{p} - 1) e^{(- \frac{{(v_{i} t + l_{i})}^{2}}{p^{2}})} \\ \bar{y_{3}} (t; Ω_{3}) = - 0.6 + t \sum_{i = 1}^{5} \frac{- 2 η_{i}}{\sqrt[4]{π} \sqrt{3 p}} (\frac{{(w_{i} t + q_{i})}^{2}}{p} - 1) e^{(- \frac{{(w_{i} t + q_{i})}^{2}}{p^{2}})} \end{array}

(21)

Here, $ρ_{1} = 0.6$ , $ρ_{2} = 0.1$ and $ρ_{3} = - 0.6$ are taken to be sensitive chaotic initial conditions. In addition, the residual error function is attained in the following form

\begin{array}{l} R_{i 1} (t, Ω_{1}) = \sum_{i = 1}^{5} α_{1} {\bar{y_{1}}}^{'} (t, Ω_{1, i}) - (1 - α_{1}) (\bar{y_{1}} (t, Ω_{1, i}) - ρ_{1}) - c_{1} (\bar{y_{2}} (t, Ω_{2, i}) + \frac{\bar{y_{1}} (t, Ω_{1, i}) - 2 {\bar{y_{1}}}^{3} (t, Ω_{1, i})}{7}) \\ R_{i 2} (t, Ω_{2}) = \sum_{i = 1}^{5} α_{2} {\bar{y_{2}}}^{'} (t, Ω_{2, i}) - (1 - α_{2}) (\bar{y_{2}} (t, Ω_{2, i}) - ρ_{2}) - c_{2} (\bar{y_{1}} (t, Ω_{1, i}) - \bar{y_{2}} (t, Ω_{2, i}) + \bar{y_{3}} (t, Ω_{3, i})) \\ R_{i 3} (t, Ω_{3}) = \sum_{i = 1}^{5} α_{3} {\bar{y_{3}}}^{'} (t, Ω_{3, i}) - (1 - α_{3}) (\bar{y_{3}} (t, Ω_{3, i}) - ρ_{3}) - c_{3} \bar{y_{2}} (t, Ω_{2, i}) \end{array}

(22)

Substituting equation (21) into equation (22) and on discretizing $t$ in 20 equally spaced points in the interval $[0, 1]$ , we get the fitness function (17). Initially, set $α_{1} = α_{2} = α_{3} = α$ , $T_{0}$ and $Ω_{0}$ , and then process the SA algorithm, defined in the section SA algorithm. The set of weight vector $Ω = {Ω_{i}}$ for $i = 1, 2, 3, 4, 5$ , in which $Ω_{i} = (λ_{i}, μ_{i}, η_{i}, u_{i}, v_{i}, w_{i}, k_{i}, l_{i}, q_{i})$ , are shown graphically in Figure 2(a) to (h). The minimum values of fitness function $E (Ω)$ are calculated by equation (17) and tabulated in Table 1 for different fractional values and number of neurons around 10⁻³ to 10⁻⁷. The smaller the value of fitness function, the better the performance of the algorithm, which is calculated using equations (18) to (20) and depicted in Tables 2 to 4. It is clear from the results of Tables 5 to 7 for MAE that the accuracy and convergence of the solutions acquired from the proposed method are in a good agreement. Moreover, the simulated path of fitness function is also presented graphically in Figure 3(a) to (b). In order to determine dynamical evolution of the Chua’s model, the optimized values of neurons are substituted in the approximate solutions $\bar{y_{i}} (t; Ω_{i})$ , for $i = 1, 2, 3$ , given in equation (21).

Figure 2.

Set of weights of ANN for different values of $α$ . (a) $α = 0.91$ ; (b) $α = 0.92$ ; (c) $α = 0.925$ ; (d) $α = 0.93$ ; (e) $α = 0.97$ ; (f) $α = 0.98$ ; (g) $α = 0.99$ ; (h) $α = 1$ .

Figure 3.

Fitness function $E$ for different values of $α$ , for n = 5. (a) α = 0.92,0.93,0.95,0.97,0.99; (b) α = 0.91,0.925,0.94,0.96,0.98.

Table 1.

The optimized values of the fitness function $E$ with different number of neurons.

α	n=2	n=3	n=4	n=5
0.91	6.03135 × 10⁻⁴	3.89984 × 10⁻³	1.09027 × 10⁻⁵	2.31533 × 10⁻⁶
0.92	5.75683 × 10⁻⁴	4.01213 × 10⁻³	1.52583 × 10⁻⁵	4.40363 × 10⁻⁷
0.925	1.4246 × 10⁻³	3.96938 × 10⁻³	2.19853 × 10⁻⁵	2.2456 × 10⁻⁶
0.93	1.34297 × 10⁻³	4.01196 × 10⁻³	9.9448 × 10⁻⁷	1.26264 × 10⁻⁵
0.94	1.3472 × 10⁻³	1.7422 × 10⁻⁵	5.2273 × 10⁻⁶	1.68169 × 10⁻⁷
0.95	5.45476 × 10⁻⁴	6.6479 × 10⁻⁵	4.28122 × 10⁻⁶	2.42878 × 10⁻⁷
0.96	4.45721 × 10⁻⁴	2.26685 × 10⁻⁴	4.52505 × 10⁻⁶	2.17893 × 10⁻⁷
0.97	5.77965 × 10⁻⁴	2.41879 × 10⁻⁶	5.63894 × 10⁻⁶	1.7401 × 10⁻⁷
0.98	5.51202 × 10⁻⁴	2.58464 × 10⁻⁶	2.2085 × 10⁻⁶	2.24173 × 10⁻⁶
0.99	5.78318 × 10⁻⁴	2.08444 × 10⁻⁴	4.12595 × 10⁻⁷	2.23617 × 10⁻⁶
1.00	1.36125 × 10⁻³	1.97773 × 10⁻⁴	1.75325 × 10⁻⁶	4.98837 × 10⁻⁶

Table 2.

The RMSE values of $y_{1} (t)$ for $α_{1} = α_{2} = α_{3} = α$ with different number of neurons.

α	n=2	n=3	n=4	n=5
0.910	0.	0.	0.	0.
0.920	1.79527 × 10⁻²	1.25237 × 10⁻¹	6.76873 × 10⁻⁴	5.47821 × 10⁻⁴
0.925	3.27488 × 10⁻²	2.12629 × 10⁻¹	4.5492 × 10⁻⁴	4.32552 × 10⁻⁴
0.930	3.1245 × 10⁻²	3.12747 × 10⁻¹	1.27574 × 10⁻³	3.41395 × 10⁻⁴
0.940	2.15593 × 10⁻²	4.98516 × 10⁻¹	1.46017 × 10⁻³	5.53936 × 10⁻⁴
0.950	2.35032 × 10⁻²	7.37568 × 10⁻¹	5.84684 × 10⁻⁴	3.23442 × 10⁻⁴
0.960	4.57716 × 10⁻²	9.86811 × 10⁻¹	1.63991 × 10⁻³	5.56365 × 10⁻⁴
0.970	7.2416 × 10⁻²	1.2252 × 10⁻²	1.26801 × 10⁻³	5.59064 × 10⁻⁴
0.980	1.24246 × 10⁻¹	1.47731 × 10⁻²	1.12602 × 10⁻³	1.05325 × 10⁻³
0.990	1.80935 × 10⁻¹	1.82389 × 10⁻²	1.92778 × 10⁻³	8.32138 × 10⁻⁴
1.00	2.2321 × 10⁻¹	2.11697 × 10⁻²	3.55936 × 10⁻³	1.27568 × 10⁻³

Table 3.

The RMSE values of $y_{2} (t)$ for $α_{1} = α_{2} = α_{3} = α$ with different number of neurons.

$α$	$n = 2$	$n = 3$	$n = 4$	$n = 5$
0.91	0.	0.	0.	0.
0.92	1.33458 × 10⁻²	1.13999 × 10⁻¹	2.7943 × 10⁻⁴	2.70204 × 10⁻⁴
0.925	1.99349 × 10⁻²	2.14987 × 10⁻¹	2.52746 × 10⁻³	2.16135 × 10⁻⁴
0.93	1.58829 × 10⁻²	3.03364 × 10⁻¹	1.86508 × 10⁻³	2.23496 × 10⁻⁴
0.94	1.32971 × 10⁻²	3.78038 × 10⁻¹	8.50141 × 10⁻⁴	1.93201 × 10⁻⁴
0.95	1.88583 × 10⁻²	4.37521 × 10⁻¹	2.60616 × 10⁻³	1.89552 × 10⁻⁴
0.96	3.22436 × 10⁻²	4.81326 × 10⁻¹	2.10578 × 10⁻³	1.78957 × 10⁻⁴
0.97	3.98189 × 10⁻²	5.10128 × 10⁻¹	5.59056 × 10⁻⁴	3.09481 × 10⁻⁴
0.98	3.76075 × 10⁻²	5.25076 × 10⁻¹	1.28543 × 10⁻⁴	1.64271 × 10⁻⁴
0.99	2.8184 × 10⁻²	5.2849 × 10⁻¹	1.40648 × 10⁻³	3.20518 × 10⁻⁴
1.00	2.50841 × 10⁻²	5.54005 × 10⁻¹	1.40011 × 10⁻³	2.95743 × 10⁻⁴

Table 4.

The RMSE values of $y_{3} (t)$ for $α_{1} = α_{2} = α_{3} = α$ with different number of neurons.

$α$	$n = 2$	$n = 3$	$n = 4$	$n = 5$
0.91	0.	0.	0.	0.
0.92	2.08247 × 10⁻²	4.71525 × 10⁻¹	2.43805 × 10⁻³	5.7026 × 10⁻⁴
0.925	6.00806 × 10⁻²	9.81781 × 10⁻¹	1.41115 × 10⁻³	8.81843 × 10⁻⁴
0.93	1.11236 × 10⁻¹	1.46856 × 10⁻²	1.00846 × 10⁻²	1.07394 × 10⁻³
0.94	1.41984 × 10⁻¹	1.86864 × 10⁻²	1.2046 × 10⁻²	9.32497 × 10⁻⁴
0.95	1.25074 × 10⁻¹	2.13658 × 10⁻²	5.31119 × 10⁻³	5.73164 × 10⁻⁴
0.96	8.08439 × 10⁻²	2.25274 × 10⁻²	2.7826 × 10⁻³	7.69002 × 10⁻⁴
0.97	1.038 × 10⁻¹	2.22233 × 10⁻²	5.70426 × 10⁻³	6.62432 × 10⁻⁴
0.98	1.51721 × 10⁻¹	2.07042 × 10⁻²	4.96893 × 10⁻³	1.04998 × 10⁻³
0.99	2.36533 × 10⁻¹	2.58315 × 10⁻²	7.38548 × 10⁻³	1.24823 × 10⁻³
1.00	2.93588 × 10⁻¹	3.14597 × 10⁻²	1.29591 × 10⁻²	1.51511 × 10⁻³

Table 5.

The MAE values of $y_{1} (t)$ for $α_{1} = α_{2} = α_{3} = α$ with different number of neurons.

$α$	$n = 2$	$n = 3$	$n = 4$	$n = 5$
0.91	1.0192 × 10⁻³	8.45526 × 10⁻³	8.34352 × 10⁻⁶	1.19794 × 10⁻⁷
0.92	9.67988 × 10⁻⁴	1.01238 × 10⁻³	8.61912 × 10⁻⁷	2.41177 × 10⁻⁷
0.925	9.27058 × 10⁻³	7.92797 × 10⁻⁴	8.72194 × 10⁻⁶	3.6872 × 10⁻⁷
0.93	8.62086 × 10⁻³	2.40082 × 10⁻⁴	1.21634 × 10⁻⁷	4.47545 × 10⁻⁷
0.94	8.85064 × 10⁻³	1.27115 × 10⁻⁴	1.7323 × 10⁻⁶	9.69204 × 10⁻⁸
0.95	8.05044 × 10⁻⁴	1.3348 × 10⁻⁵	4.5164 × 10⁻⁷	2.63846 × 10⁻⁷
0.96	9.46187 × 10⁻⁴	1.40401 × 10⁻⁵	8.28316 × 10⁻⁷	6.04626 × 10⁻⁸
0.97	1.63153 × 10⁻³	1.19406 × 10⁻⁴	1.55015 × 10⁻⁶	7.24064 × 10⁻⁸
0.98	1.37203 × 10⁻³	1.42808 × 10⁻⁴	6.19608 × 10⁻⁸	2.16951 × 10⁻⁷
0.99	1.34295 × 10⁻³	4.09688 × 10⁻⁴	1.3274 × 10⁻⁷	2.32376 × 10⁻⁷
1.00	5.13437 × 10⁻²	2.07393 × 10⁻⁴	1.03094 × 10⁻⁷	5.06675 × 10⁻⁷

Table 6.

The MAE values of $y_{2} (t)$ for $α_{1} = α_{2} = α_{3} = α$ with different number of neurons.

α	n=3	n=2	n=4	n=5
0.91	2.81948 × 10⁻²	3.62599 × 10⁻⁵	6.17461 × 10⁻⁷	3.33582 × 10⁻⁸
0.92	3.93335 × 10⁻²	3.38324 × 10⁻⁵	1.04763 × 10⁻⁷	2.07015 × 10⁻⁸
0.925	4.34069 × 10⁻²	6.63781 × 10⁻⁵	5.01667 × 10⁻⁷	1.11443 × 10⁻⁷
0.93	6.49941 × 10⁻³	6.17928 × 10⁻⁴	3.45252 × 10⁻⁸	7.51071 × 10⁻⁸
0.94	1.1406 × 10⁻²	6.05695 × 10⁻⁴	1.73936 × 10⁻⁷	3.53147 × 10⁻⁹
0.95	1.48353 × 10⁻¹	3.3616 × 10⁻⁵	1.51261 × 10⁻⁷	1.30539 × 10⁻⁸
0.96	1.07002 × 10⁻¹	6.72586 × 10⁻⁵	7.59788 × 10⁻⁸	5.55576 × 10⁻⁹
0.97	9.40824 × 10⁻²	6.81374 × 10⁻⁵	1.01052 × 10⁻⁷	6.4805 × 10⁻⁹
0.98	2.56723 × 10⁻²	5.75292 × 10⁻⁵	3.71695 × 10⁻⁸	1.42432 × 10⁻⁸
0.99	6.64585 × 10⁻²	5.81643 × 10⁻⁵	1.08081 × 10⁻⁸	1.54103 × 10⁻⁸
1.00	4.39674 × 10⁻⁴	1.63473 × 10⁻³	2.98973 × 10⁻⁸	4.48061 × 10⁻⁸

Table 7.

The values of MAE $y_{3} (t)$ for $α_{1} = α_{2} = α_{3} = α$ with different number of neurons.

α	n=2	n=3	n=4	n=5
0.91	1.47302 × 10⁻³	2.72487 × 10⁻¹	1.52671 × 10⁻⁵	2.19408 × 10⁻⁷
0.92	1.35661 × 10⁻³	3.22595 × 10⁻³	5.60809 × 10⁻⁷	2.78145 × 10⁻⁷
0.925	2.86158 × 10⁻²	1.32671 × 10⁻³	1.08043 × 10⁻⁵	1.41903 × 10⁻⁶
0.93	2.58598 × 10⁻²	9.63953 × 10⁻³	3.57978 × 10⁻⁷	6.68145 × 10⁻⁷
0.94	2.56086 × 10⁻²	6.3173 × 10⁻¹	3.37428 × 10⁻⁶	8.48716 × 10⁻⁸
0.95	1.05601 × 10⁻³	5.561 × 10⁻¹	1.9535 × 10⁻⁶	3.22159 × 10⁻⁷
0.96	1.2242 × 10⁻³	5.10432 × 10⁻¹	1.48546 × 10⁻⁶	7.61834 × 10⁻⁸
0.97	2.22051 × 10⁻³	5.48337 × 10⁻¹	2.28513 × 10⁻⁶	1.04404 × 10⁻⁷
0.98	1.85189 × 10⁻³	1.31341 × 10⁻³	4.39618 × 10⁻⁷	2.32253 × 10⁻⁷
0.99	1.81678 × 10⁻³	3.98598 × 10⁻³	1.53801 × 10⁻⁷	2.91462 × 10⁻⁷
1.00	5.04879 × 10⁻²	6.37526 × 10⁻³	2.96405 × 10⁻⁷	7.086 × 10⁻⁷

The dynamical behaviour of the time series of equation (13), attained by using FFRK method, are plotted in Figures 4 to 6 for $y_{1} (t)$ , $y_{2} (t)$ and $y_{3} (t)$ , respectively, at different fractional values of $α$ . These portrayals illustrate that with the increase in the time interval, i.e. the value of $p$ , Chua’s system exhibits chaotic behaviour. It can be visually demonstrated by the phase portraits of peculiar attractors, shown in Figure 7(a) to (f) for $y_{1} (t)$ versus $y_{2} (t)$ , Figure 8(a) to (f) for $y_{3} (t)$ versus $y_{2} (t)$ , Figure 9(a) to (f) for $y_{3} (t)$ versus $y_{1} (t)$ , for $p = 20$ . Consequently, with respect to the fractional value of $α$ , system (13) has different periodic attractors for the different fractional values of $α$ . Moreover, Figure 10(a) to (f) elaborates the three-dimensional phase portraits of the Chua's system. Hence it is clearly seen that, the system (13) exhibits different chaotic and periodic behaviour with the increase of parameter $p$ and $0.925 \leq α \leq 1$ .

Figure 4.

Time series of $y_{1} (t)$ for the Chua’s attractor obtained by FFRK method. (a) α = 0.96; (b) α = 0.97; (c) α = 0.99; (d): α = 1.

Figure 5.

Time series of $y_{2} (t)$ for the Chua’s attractor obtained by FFRK method. (a) α = 0.96; (b) α = 0.97; (c): α = 0.99; (d) α = 1.

Figure 6.

Time series of $y_{3} (t)$ for the Chua’s attractor obtained by FFRK method. (a) α = 0.96; (b) α = 0.97; (c) α = 0.99; (d) α = 1.

Figure 7.

Phase plane projections, for $y_{1} (t)$ versus $y_{2} (t)$ , for the state space configuration of Chua's system. (a) $α = 0.95$ ; (b) $α = 0.96$ ; (c) $α = 0.97$ ; (d) $α = 0.98$ ; (e) $α = 0.99$ ; (f) $α = 1.00$ .

Figure 8.

Phase plane projections, for $y_{3} (t)$ versus $y_{2} (t)$ , for the state space configuration of Chua's system. (a) $α = 0.95$ ; (b) $α = 0.96$ ; (c) $α = 0.97$ ; (d) $α = 0.98$ ; (e) $α = 0.99$ ; (f) $α = 1.00$ .

Figure 9.

Figure 10.

3D Phase portrait of Chua’s system at different order $α$ . (a) $α = 0.95$ ; (b) $α = 0.96$ ; (c) $α = 0.97$ ; (d) $α = 0.98$ ; (e) $α = 0.99$ ; (f) $α = 1.00$ .

Furthermore, the bifurcation diagram of the system (13) is also plotted in Figure 11, with respect to the parameter $c_{1}$ and visited asymptotically from initial conditions equation (8). It can be observed that when $c_{1} \in (0.04, 0.30)$ the system performs the chaotic behaviour. Additionally, to study the effect of fractional derivatives on the dynamics of the system (13), the state space configuration is considered as well, to verify chaos by computing LE.^8,9 The spectacle of this configuration is shown in Figure 12, which is helpful to further make a comprehensive bifurcation study.

Figure 11.

Bifurcation diagram of the fractional-order Chua’s system.

Figure 12.

LEs of the Chua’s system with cubic nonlinearity.

Conclusion

In the present work, we carried out a systematic stratagem to scrutinize the fractional Chua's system and also comprehensively discussed its chaotic behaviour. The parametric linearization of fractional Laplace transform and two well-known methodologies, MHW-based ANN with SA and FFRK method, were proposed to attain the solutions of fractional Chua's model. Comprehensively, from the entire study, following lucrative outcomes were achieved:

The advantage of fractional Laplace transform expansion, which converted the fractional Chua's model into integer order differential system that remains in the fractional environment, made its simulation much easier.

The MHW-based ANN with SA algorithm efficiently solved the cubic nonlinearity and produced optimized values of the dynamical variables of the Chua's circuit for different fractional values that ascertains its constructive applicability to solve nonlinear fractional systems.

Through LE entropy we analysed the stability and sensitivity of the initial conditions and presented chaotic predictions of fractional Chua's model. The successful plots of LE for the cubic nonlinearity of Chua's system with the fractional factor verify the ability of LE to study the stability of fractional order systems as well.

The FFRK method made it possible to attain phase portraits of chaotic attractors in different dimensions and for different fractional values of the Chua's circuit, which illustrates a notable capability of this method to describe the chaotic behaviour of the systems.

Additionally, by means of visualization of the bifurcation plots and chaotic attractor of fractional Chua's circuit, the key objective of this attempt is achieved.

In the future, we will study dynamical behaviours of many other nonlinear systems that appear in different fields of science. Technically, we will consider some other amalgamation of neural networks with optimization algorithms to solve fractional partial differential equations as well.

Footnotes

Acknowledgements

The corresponding author acknowledges the referees and Editor Prof JH He for his comments and suggestions.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Hirsch

Smale

Devaney

RL.

Differential equations, dynamical systems, and an introduction to chaos. 2nd ed. San Diego: Elsevier, 2004.

Aguilar

JFG

López

MGL

Martínez

VMA

, et al. Chaos in a cancer model via fractional derivatives with exponential decay and Mittag-Leffler Law. Entropy 2017; 19: 681.

Pivka

Huang

Chua’s oscillator: a compendium of chaotic phenomena. J Franklin Inst 1994; 331: 705–741.

Matsumoto

A chaotic attractor from Chua’s circuit. IEEE Trans Circuits Syst 1984; 31: 1055–1058.

Alkahtani

BST.

Chua’s circuit model with Atangana–Baleanu derivative with fractional order. Chaos Solitons Fractals 2016; 89: 547–551.

Menacer

Lozi

Chua

LO.

Hidden bifurcations in the multispiral Chua attractor. Int J Bifurcation Chaos 2016; 26: 1630039.

Cafagna

Grassi

Fractional-order Chua’s circuit: time-domain analysis, bifurcation, chaotic behavior and test for chaos. Int J Bifurcation Chaos 2008; 18: 615–639.

Singh

Roy

BK.

The nature of Lyapunov exponents is (+, +, −, −). Is it a hyperchaotic system?

Chaos Solitons Fractals 2016; 92: 73–85.

Baleanu

Luo

WH.

Lyapunov functions for Riemann–Liouville-like fractional difference equations. Appl Math Comput 2017; 314: 228–236.

10.

JH.

A new fractal derivation. Therm Sci 2011; 15: S145–S147.

11.

Atangana

Baleanu

New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Therm Sci 2016; 20: 763–769.

12.

Khalil

Horani

Yousef

, et al. A new definition of fractional derivative. J Comput Appl Math 2014; 264: 65–70.

13.

Wang

Liu

SY.

He's fractional derivative for non-linear fractional heat transfer equation. Therm Sci 2016; 20: 793–796.

14.

Sayevand

Pichaghchi

Analysis of nonlinear fractional KdV equation based on He’s fractional derivative. Nonlinear Sci Lett A 2016; 7: 77–85.

15.

Khan

Hameed

Khan

, et al. A heuristic optimization method of fractional convection reaction: an application to diffusion process. Therm Sci 2018; 22: S243–S252.

16.

Khan

Razzaq

Ayaz

Some properties and applications of conformable fractional Laplace transform (CFLT.). J Frac Calc Appl 2018; 9: 72–81.

17.

Wang

Zhang

Liu

, et al. A short review on analytical methods for fractional equations with He's fractional derivative. Therm Sci 2017; 21: DOI: 10.2298/TSCI160513036W

18.

JH.

Fractal calculus and its geometrical explanation. Results Phys 2018; 10: 272–276.

19.

Khan

Razzaq

Ayaz

Notes on fuzzy fractional Sumudu transform. J Math Computer Sci 2018; 18: 63–73.

20.

Ren

Sun

Dai

New approximations for solving the Caputo-type fractional partial differential equations. Appl Math Model 2016; 40: 2625–2636.

21.

Khan

Shaikh

Sultan

, et al. Numerical simulation using artificial neural network on fractional differential equations. Numerical Simulation InTech, http://dx.doi.org/10.5772/64151 (accessed on 1 November 2018).

22.

Yang

XS.

Simulated annealing. In: XS

Yang

(ed.) Nature-inspired optimization algorithms. Elsevier: Luniver Press, 2014, pp. 67–75.

23.

Butcher

JC.

Numerical methods for ordinary differential equations in the 20th century. J Comput Appl Math 2000; 125: 1–29.

24.

Masooda

Majeeda

Samara

, et al. Design of Mexican hat wavelet neural networks for solving Bratu type nonlinear systems. Neurocomputing 2017; 221: 1–14.

25.

Dinç

Duarte

Machado

JAT

, et al. Application of continuous wavelet transform to the analysis of the modulus of the fractional Fourier transform bands for resolving two component mixture. Signal Image Video Process 2015; 9: 801–807.