Chaotic Dynamics and Complexity in Real and Physical Systems

Abstract

Emergence of chaos and complex behavior in real and physical systems has been discussed within the framework of nonlinear dynamical systems. The problems investigated include complexity of Child’s swing dynamics , chaotic neuronal dynamics (FHN model), complex Food-web dynamics, Financial model (involving interest rate, investment demand and price index) etc. Proper numerical simulations have been carried out to unravel the complex dynamics of these systems and significant results obtained are displayed through tables and various plots like bifurcations, attractors, Lyapunov exponents, topological entropies, correlation dimensions, recurrence plots etc. The significance of artificial neural network (ANN) framework for time series generation of some dynamical system is suggested.

1. Introduction

In this chapter, we investigate the dynamical complexity of several real physical systems. We present our analysis of various problems considered here and present results graphically based on actual numerical simulation for various system. We revisit the analysis of complexity of nonlinear pendulum dynamics and its application to unravel the complex oscillations observed in a swing pumped by a child. For the analysis, we use various tools e.g., phase plot, bifurcation diagram, Poincare surface of section and maps, Lyapunov exponent (LCE) etc., of theory of nonlinear dynamical system. Next we consider the problem of prey-predator system with Allee effect and introduce correlation dimension and topological entropy to characterize the fractal structure and the associated complexity in its dynamics. Further, beside the normal analysis used to understand the complex neuronal dynamics, say using Fitzhug-Nagumo model (FHN), recurrence plots (RPs) have been used along with the phase plot analysis and bifurcation diagram to picturise the transition of spike occurrence from periodic to quasi-periodic and chaotic oscillations in the presence of external periodic stimulation. Significance of multi-scale permutation entropy analysis to characterize nonlinear dynamical complexity of real system is also suggested while analyzing a financial system involving interest rate, investment demand and interest rates. Finally, we describe the utility of time series generation of dynamical variables of chaotic system, such as Lorentz system, using artificial neural network.

2. Pendulum motion and dynamics behind a swing

2.1 Pendulum motion

Motion of nonlinear driven pendulum with friction are widely discussed through numerous literature in Physics and applied Mathematics, (e.g. here, [1, 2, 3]). The nonlinear analysis of driven nonlinear pendulum provides a basis for understanding the complexity of various nonlinear dynamical systems. Regular and chaotic motions are observed in such pendulums depending on the numerical values assigned to the parameters associated in their equations of motion. A swing dynamics is very similar to that of nonlinear driven pendulum, [4, 5, 6, 7]. In the present text regular and chaotic motion of a pendulum and that of the child’s swing is discussed mathematically. Numerical results are presented in various forms of graphics. The equation of motion of a driven pendulum having angular displacement, $\theta$, from vertical with linear damping expressed as

$$\frac{d^2\theta }{{\mathit{dt}}^2}+k\frac{\mathit{d\theta }}{\mathit{dt}}+{\omega }_0^{2}sin\theta =Fcos\left(\mathit{\omega t}\right)$$E1

where $F$ and $\omega$ are respectively the amplitude and frequency of the driving force and $k$ is the damping coefficient and ${\omega }_0=\sqrt{\frac{g}{L}}$ is the natural frequency for free small-amplitude oscillations. Here $g$ refers to acceleration due to gravity and $L$ the length of the pendulum. Often it is convenient to express frequency in units of ${\omega }_0$ by setting ${\omega }_0\to 1$ and rescaling the time unit accordingly. The periodic force $Fcos\mathit{\omega t}$ is active and influence the motion of the pendulum.

The Eq. (1) can easily be replaced by equation

$$\frac{d^2\theta }{{\mathit{dt}}^2}+2\beta \frac{\mathit{d\theta }}{\mathit{dt}}+{\omega }_0^{2}sin\theta =f{\omega }_0^{2}cos\left(\mathit{\omega t}\right)$$E2

Here $k=2\beta$, $\beta =\frac{\gamma }{m}$, i.e., ratio of damping coefficient per unit mass $m$, ${\omega }_0^2=\frac{g}{L}$ and $f=\frac{F}{{\omega }_0^2}$. One obtains a bifurcation diagram for Eq. (2), shown in Figure 1. Bifurcation scenario indicates a period doubling phenomena followed by chaos. This implies the pendulum oscillations may be regular or chaotic depending on the magnitudes of external forcing.

Figure 1.

Bifurcation scenario of damped and driven pendulum for values $\beta =0.75,\omega =2\pi ,{\omega }_0=3\pi$ and $1.06\le f\le 1.09$.

System (1) or (2) are very common structure with very few degrees of freedom. The simple forced pendulum is periodic in $\theta$ when the driving force applied is smaller than a certain threshold $f_C$ and chaotic when the force $f$ is greater than this threshold. An interesting case would also be when the driving force becomes comparable to the weight.

Keeping fixed values for $k=0.125,{\omega }_0=1,\omega =2/3$ and then varying forcing amplitude $F$, following regular (periodic) and chaotic motion of the pendulum is observed:

Case (a): for a value $F=0.2$, a time-series plot, phase plot, surface of section and Poincaré map, are drawn as shown in Figure 2.

Figure 2.

Time-series plot, phase plot, surface of section and Poincaré map for periodic motion for $F=0.2$.

Case (b): for a value $F=0.8$, corresponding figures of case (a) are obtained as a time-series plot, phase plot, surface of section and Poincaré map, are drawn as shown in Figure 3.

Figure 3.

Time-series plot, phase plot, surface of section and Poincaré map for periodic motion for $F=0.8$.

The plots shown in Figure 3 indicate at a value $F=0.8$ the pendulum oscillation is chaotic and this leads to unpredictability. Pendulum may be whirling irregularly or overturn or show very irregular oscillations.

As an application of the foregoing analysis, in the following section, we extend the formalism to discuss the problem of swing oscillation where the length of the pendulum varies periodically.

2.2 Problem of Swing oscillation

Oscillations of a swing pumped by a child is very familiar to us. Every time the swing passes through its lowest point the child pumps it over and again. The dynamics of weightless rod with a point mass sliding along the length mimics like a pendulum swing whose length varies periodically with time. The motion of the swing governed by the dynamical system written as [5]:

$$\frac{d}{\mathit{dt}}\left[{\mathit{ml}}^2\frac{\mathit{d\theta }}{\mathit{dt}}\right]+\gamma l^2\frac{\mathit{d\theta }}{\mathit{dt}}+\mathit{mgl}sin\theta =0$$E3

where $m$ is the mass, $l$ is the length, $\theta$ is the angle made by the swing from the vertical position and $g$ is the acceleration due to gravity. As the length of the swing varies periodically with time, one assumes

$$l=l_0+\mathit{a\varphi }\left(\mathrm{\Omega }t\right)$$E4

where $l_0$ is the mean length of the swing and is constant, $a$ and $\mathrm{\Omega }$ are, respectively, the amplitude and frequency of excitation. The function $\varphi \left(\mathrm{\Omega }t\right)$ should be a periodic function of time. Then, by introducing the following dimensionless parameters and variables

$$\tau =\mathrm{\Omega }t,\epsilon =\frac{a}{l_0},{\mathrm{\Omega }}_0=\sqrt{\frac{g}{l_0}},\omega =\frac{{\mathrm{\Omega }}_0}{\mathrm{\Omega }},\beta =\frac{\gamma }{m{\mathrm{\Omega }}_0},$$

equation of motion of the swing in dimensionless form written as:

$$\ddot{\theta }+\left[\frac{2\epsilon \dot{\varphi }\left(\tau \right)}{1+\mathit{\epsilon \varphi }\left(\tau \right)}+\mathit{\beta \omega }\right]\dot{\theta }+\frac{{\omega }^2}{1+\mathit{\epsilon \varphi }\left(\tau \right)}sin\theta =0$$E5

where $\left\{.\right\}$ in Eq. (5) corresponds to differentiation with respect $\tau$.

Since $\varphi \left(\tau \right)$ is a periodic function, we may take $\varphi \left(\tau \right)=Asin\left(\mathit{\lambda \tau }\right)$ and thus the foregoing equation may be rewritten as:

$$\ddot{\theta }+\left[\frac{2\mathit{k\lambda }cos\left(\mathit{\lambda \tau }\right)}{1+ksin\left(\mathit{\lambda \tau }\right)}+\mathit{\beta \omega }\right]\dot{\theta }+\frac{{\omega }^2}{1+ksin\left(\mathit{\lambda \tau }\right)}sin\theta =0$$E6

where $k=\mathit{\epsilon A}$.

For stability of motion of the swing a linear stability analysis is applied. We may write Eq. (6) as the following two first order equations:

$$\begin{array}{ll}\dot{\theta }& =u\equiv f\left(\theta ,u\right)\\ \dot{u}& =-\left[\frac{2\mathit{k\lambda }cos\left(\mathit{\lambda \tau }\right)}{1+ksin\left(\mathit{\lambda \tau }\right)}+\mathit{\beta \omega }\right]u-\frac{{\omega }^2}{1+ksin\left(\mathit{\lambda \tau }\right)}sin\theta \equiv g\left(\theta ,u\right)\end{array}$$E7

The Jacobian for the above system may be written as,

$$J=\left(\begin{array}{cc}0& 1\\ -\frac{{\omega }^2}{1+ksin\left(\mathit{\lambda \tau }\right)}cos\theta & -\left[\frac{2\mathit{k\lambda }cos\left(\mathit{\lambda \tau }\right)}{1+ksin\left(\mathit{\lambda \tau }\right)}+\mathit{\beta \omega }\right]\end{array}\right)$$

When an external periodic force $Fcos\left(\mathit{\vartheta \tau }\right)$ is applied to pump the swing, final form of the equation of motion stands as

$$\ddot{\theta }+\left[\frac{2\mathit{k\lambda }cos\left(\mathit{\lambda \tau }\right)}{1+ksin\left(\mathit{\lambda \tau }\right)}+\mathit{\beta \omega }\right]\dot{\theta }+\frac{{\omega }^2}{1+ksin\left(\mathit{\lambda \tau }\right)}sin\theta =Fcos\left(\mathit{\vartheta \tau }\right)$$E8

2.3 Regular and Chaotic motion of the swing

Figure 4.

A time-series and phase plots and plots of surface of section and Poincaré map for regular motion of the swing for $F=0.8,\beta =0.5,k=0.1,\lambda =0.05,\omega =1,\vartheta =2/3.$

Figure 5.

Regular two periodic motion of the swing when $F=0.9,\beta =0.5,k=0.3,\lambda =1,\omega =1,\vartheta =2/3$.

The swing, Eq. (8), oscillates in regular motion for significant contribution of friction, (i.e. when the frictional coefficient $\beta$ has sufficiently higher value) and it is in chaotic motion in case of small friction and higher values of driving force. Figures 4 and 5 showing the case of regular motion.

When the frictional contribution is insignificant, swing oscillations are chaotic and unpredictable. Figure 6 stands for such chaotic motion of the swing when $\beta =0$.

Figure 6.

Chaotic oscillation of the swing when $F=0.2,\beta =0,k=0.1,\lambda =0.05,\omega =1,\vartheta =2/3$.

Figure 7 show chaotic oscillation when $\beta$ is not zero but small. Surface of section and Poincare map shown in this figure are interesting showing typical chaotic behavior.

Figure 7.

Chaotic oscillation of the swing when $F=0.2,\beta =0.01,k=0.1,\lambda =0.05,\omega =1,\vartheta =2/3$.

We may thus conclude that the swing oscillates smoothly when the frictions are higher but for no friction or insignificant friction, swing oscillations would be chaotic or unpredictable. In such a case whirling, overturn or any unpredictable situation may happen.

Beside the application of bifurcation diagram, phase plot and Poincare surface of section technique, we introduce the idea of Lyapunov characteristic exponents (LCE), correlation dimension and topological entropy which provide further insight of a complex dynamical system. In the following section, we analyze the complexity of Prey-predator system using such tools.

3. Complexity in prey-predator system with Allee effect

In recent years many type of predator-prey problems, originated in Biological sciences, investigated which depend on various environmental and social conditions, [8, 9, 10, 11]. Some problems solved by the application of Allee effect, which is an interesting phenomenon, to some predator-prey systems appear to be very interesting, [12, 13, 14, 15, 16]. The Allee effect on prey-predator system is a phenomenon in biology which characterizes certain correlation between population size or density and the mean individual fitness of a population or species. In the following study we investigate the complexity in a predator – prey problem with the Allee effect.

3.1 Discrete prey-predator model

A model for the prey-predator problem with Allee effect can represented as

$$\begin{array}{l}{X}_{n+1}=X_n+{\mathit{rX}}_n\left(1-X_n\right)\left(1-exp\left(-\epsilon X_n\right)\right)-{\mathit{aX}}_n{Y}_n\\ Y_{n+1}=Y_n+{\mathit{aY}}_n\left(X_n-Y_n\right)\left(\frac{Y_n}{\mu +Y_n}\right)\end{array}$$E9

where $X_n$ and $Y_n$ refers to the density of prey and predators. Further, $r$ correspond to the growth rate parameter of the prey population and $a$ the predation prameter. Here,

• $1-exp\left(-\epsilon X_n\right)$ stands for mate finding Allee effect on prey population, here $\epsilon$ is defined as the Allee effect constant and the term

• $\frac{Y_n}{\mu +y_n}$ stands for the Allee effect on predator and here, $\mu$ is the Allee effect constant. Bigger $\mu$ means the stronger the Allee effect on predator population.

For assumed values of parameters $a=2.0,r=2.4$, fixed points of system (9) are obtained, approximately, as $P_1^{\ast }\left(0,0\right),P_2^{\ast }\left(1,0\right),P_3^{\ast }\left(\mathrm{0.545455,0.545455}\right)$ and by using stability analysis, we find all are unstable.

3.2 Bifurcation diagrams

The phenomena of bifurcation provide a qualitative change in the behavior of a system during evolution. Such a change occurs when a particular parameter is varied while keeping other parameters constant. Bifurcation diagram shows the splitting of stable solutions within a certain range of values of the parameter. During the processes of bifurcation, one observes different cycles of evolution which leading to the chaotic situation. Phenomena like bistability, periodic windows within chaos etc. may also be observed for some systems. A bifurcation can be taken as a tool to analyze the regular, chaotic as well as complexity within the system. For $a=2.0$ and $1.8\le r\le 2.4$ , Figure 8 shows bifurcation of system (9), where some interesting phenomena observed that the system is not producing a period doubling bifurcation scenario which is very common for many nonlinear systems.

Figure 8.

Bifurcation diagram of system (9), (a) Prey densities, (b) Predator densities for $a=2.0$ and $1.8\le r\le 2.4$.

3.3 Numerical simulations

Simulation of the foregoing models provides,

3.3.1 Attractors

Keeping parameters $a$ and $r$ fixed viz., $a=2.0,r=2.4$, attractors for different cases are obtained through numerical technique [16], and shown in Figure 9.

Figure 9.

Plots of regular and chaotic attractors for $a=2.0$ and $r=2.4$; (i) plot (a) without Allee effect, (ii) plot (b) with Allee effect on prey only, $\epsilon =4.5$ , (iii) plot (c) Allee effect on predator only $\mathit{mu}=0.1$, and (iv) plot (d) Allee effect on prey as well as on predator, $\epsilon =4.5,\mu =0.1$.

Looking plots of attractor of Figure 9, one finds a chaotic attractor, figure (a), when Allee effect is not in consideration, for $a=2.0,r=2.4$ . But, the application of Allee effect to either of the population or to both population, system returned to regularity, e.g. figures (b), (c) and (d) are no more chaotic. This also follow from the plots of LCEs given below.

3.3.2 Lyapunov exponents (LCEs)

The phase space dynamics of a nonlinear chaotic physical system is very complex in general. One of the important feature of such a system is its sensitivity to initial conditions i.e., two very nearby trajectory in phase space show divergence exponentially. Such divergence are characterized by LCEs. To indicate chaotic and regular evolution, an appropriate measure is to find Lyapunov exponents (LCEs) which are obtained for different cases by using appropriate procedure. Plots of LCEs are shown in Figure 10.

Figure 10.

Plots of Lyapunov exponents for $a=2.0,r=2.4$ and (i) figure (a) without Allee effect, (ii) figure (b) with $\epsilon =4.5,\mu =0$ , (iii) figure (c) Allee effect on predator only with $\mu =0.1$ , (iv) Allee effect on both populations $\epsilon =4.5,\mu =0.1$ .

3.3.3 Correlation dimension

Lorenz attractor provides an example of a fractal object with noninteger dimension. The correlation dimension permits us to quantify the space filling property and provides the measure of dimensionality of the chaotic attractor. It is expressed as

$$D^c=\frac{dlogC\left(R\right)}{dlog\left(R\right)}$$

where $C\left(R\right)$ is defined as

$$C\left(R\right)=\frac{1}{n\left(n-1\right)}\sum _{i=1}^n\sum _{j=1}^n\left[\mathrm{\Theta }\left(R-\parallel x\left[i\right]-x\left[j\right]\parallel \right)\right]$$

corresponds to the correlation sum and is a measure of total number of points contained within a hypersphere of radius $R$ as a function of $R$ normalized to the total number of points squared. Using the algorithm [17, 18], the correlation dimension can be determined from the scaling region found in the plot of $logC\left(R\right)$ as a function of $log\left(R\right)$. For the Lorenz system with parameters $\sigma =10,\rho =28$ and $b=8/3$, the correlation dimension $D^c$ is found to be $2.069$. The correlation dimension of the chaotic attractor Figure 9(a) is found to be $D^c=0.571$ (Figure 11).

Figure 11.

Plot of correlation integral data.

3.3.4 Topological entropies

As explained in the beginning, topological entropy measures the complexity of the system. More topological entropy implies system is more complex. Presence of complexity does not mean the system is chaotic and vice versa. In Figure 12, we have plots of topological entropy for different cases. In figure (a), topological entropy increases for $r>2$ but bifurcation diagrams and calculations of LCEs indicate the system is regular within $2.0\le r\le 2.2$. Similar observation can be made looking at figures (b) and (c). In figure (d) one finds no fluctuations of topological entropy, it establishes a steady state situation.

Figure 12.

Plots of topological entropies for $a=2.0$ and $1.8\le r\le 2.6$ : (i) figure (a) with no Allee effect, (ii) figure (b) when $\epsilon =4.5,\mu =0$, figure (c) with Allee effect on predator only $\mu =0.1$ , (iv) when $\epsilon =4.5,\mu =01$.

The results obtained through bifurcation plots, Figure 8, and those of LCEs plots, Figure 10, show that the Allee effect stabilize the motion from chaos to regularity. The correlation dimension of the chaotic attractor is obtained as $D^c\cong 0.571$. Through this study we find the existence of complexity within the system, even when system behavior is regular, we find significant amount of increase in topological entropy. This implies the fact that the system may be regular but may exhibit complexity.

4. Recurrence plot

Natural system exhibits periodicities and also irregular cyclicities. Usually measures such as Lyapunov characteristic exponent (LCE), correlation dimension, Kolmogorov- Sinai (KS) entropy etc., have been used to characterize the complexity of observed nonlinear dynamical behavior of a system. But the analysis based on application of the foregoing tools inherently assumes the system to be noise free and stationary. An alternative framework based on the idea of recurrence plot was introduced in [19] for visualization of the dynamical behavior of a system in phase space and subsequently the formalism has been extended to quantify the recurrence plots to unravel the observed complexities i.e., regular, quasi-periodic, chaotic transition etc. For a discrete time series P with N data points such that

$$P:\left\{x_1,x_2,x_3,\cdots ,x_N\right\},$$E10

where $x_i,i=1,2,..\dots ,N$ refers to observed values at time $t_1,t_1+\mathrm{\Delta }t,.\dots ,t_1+n\mathrm{\Delta }t$. If the system has true dimension $m$ , a sequence of vectors may be constructed from the time series as:

$${\mathbf{X}}_{\mathbf{i}}=\left(x_i,x_{i+\tau },x_{i+2\tau },\cdots ,x_{i+\left(m-1\right)\tau }\right);i=1,2,\cdots ,n;n=N-\left(m-1\right)\tau$$E11

where $\tau$ corresponds to time lag or delay and $m$ - the embedding dimension of the phase space. By considering the distances in $m$- dimensional reconstructed points, we construct a recurrence plot (RP). In fact RP is an $n\times n$ symmetrical array where a dot is marked at a point $\left(i,j\right)$ if ${\mathbf{X}}_{\mathbf{i}}$ is close to another point ${\mathbf{X}}_{\mathbf{j}}$. We may write

$$R_{i,j}\left(r\right)=\mathrm{\Theta }\left(r-\parallel {\mathbf{X}}_{\mathbf{i}}-{\mathbf{X}}_{\mathbf{j}}\parallel \right);i,j=1,2,\cdots ,n$$E12

where $\mathrm{\Theta }$ is a Heavyside function, $r$ is a small threshold distance between neighboring points and $\parallel .\parallel$ is an Euclidean norm.

A characteristic pattern emerges in RPs which characterizes a dynamical system. The method of RP is suitable for both stationary and non-stationary dynamical system. Since a trajectory may return to a point or close to it in phase space, the deterministic dynamical system shows recurrent behavior and RP therefore exhibits both horizontal and vertical lines. For a stochastic dynamical system, such lines in RP are of very small size and in fact appear by chance. Therefore the distribution of such points appear to be homogeneous. In case of periodic system the RP is filled with longer diagonal lines. Various measures that quantify RPs are mainly, RR, DET, ENT, DIV, LAM, TT which refers to density of recurrence points, determinism, divergence, entropy, laminarity and trapping time respectively. For a periodic system diagonal lines are longer which for chaotic system RP shows broken short lines. Recently [20] has provided a very useful description of applying RPs and recurrence quantification analysis to unravel the complex dynamics of general problem of three species interaction in ecology. In the following section, we extend the analysis of complexity to the problem associated with neuronal dynamics, an area of current interest in neuro-bio-science [21]. We however restrict ourselves to only RPs to supplement the analysis of complexity using phase portrait, bifurcation diagram etc.

5. Regular and Chaotic neuronal dynamics

The neuronal communication is known to be mediated by electrical pulses called spikes. Studies of various spiking patterns reveal nonlinear characteristics of slow-fast neuronal dynamics. A considerable amount of information regarding neuronal activity has been obtained by studying the dynamics of spiking pattern [22]. The phenomenon of tonic firing, mixed mode (bursting and spiking) etc., are typical responses exhibited by an excitable neuron [23]. Cortical neurons have been reported to show tonic bursting wherein the neurons periodically switches between firing state and resting state. The mixed mode firing is observed in mammalian neocortex [24]. Spike generation in fact depends on the firing threshold and the stimulus intensity. In recent years, the perception regarding constancy of neuron’s firing threshold has changed to dynamic [21]. In this work, we first briefly introduce the Fitzhugh-Nagumo model (FHN) that have been proposed for spike generation like well known Hodgkin-Huxley model. It is however to be noted that FHN model reproduce the experimental results less accurately. Our interest in FHN model emanates mainly due to its showing complex spiking pattern even though it is mathematically simple. The basic FHN model assumes the threshold to be constant. We also study the changes caused in the spiking pattern as a result of time varying threshold. This study assumes significance as such a model may throw insight into the model the dynamics of cortisol secretion from hypothalamus [25]. It is to be noted that the neuronal firings may take place at regular interval or randomly due to inherent mechanism or may be due to its interaction with the neighborhood neurons or result of exogenous stimulus [26, 27, 28].

5.1 Basic dynamics of FHN model

The FHN model describes the interaction between the voltage $v$ across the axon membrane driven by input current $I$ and the recovery variable $w$. The recovery variable $w$ is the result of mainly the reflecting outward potassium current ($K^{+}$) that results in hyperpolarization of the axon after each spike occurrence. We may write the FHN model equation as [28, 29, 30]:

$$\begin{array}{cc}\frac{\mathit{dv}}{\mathit{dt}}& =\mathit{\alpha v}\left(\beta -v\right)\left(v-b_0\right)-\mathit{\sigma w}+I\\ \frac{\mathit{dw}}{\mathit{dt}}& =\epsilon \left(v-\mathit{\delta w}\right)\end{array}$$E13

where $\delta >0$ and the parameter $\alpha >0$ scales the amplitude of the membrane potential $v$, and $\epsilon$ is used here to control the recovery variable $w$ with respect to action potential $v$. The parameter $b_0$ i.e., $\left(0<b_0<1\right)$, corresponds to the threshold value that controls the excitable behavior of the neuron. Also $\beta$ and $\sigma$ are constants for the system.

In our analysis of Eq. (13), we take $\beta =\sigma =1$ and for the case of no external input current, $I=0$ the dynamical system (13) has three equilibrium points or fixed points $\left(v_e,w_e\right)$ as:

$$\begin{array}{c}{E}_1=\left(0,0\right),\\ E_{2,3}=\frac{\left(1+b_0\right)}{2}\pm \frac{\sqrt{{\left(1-b_0\right)}^2-\frac{4}{\mathit{\alpha \delta }}}}{2}.\end{array}$$

It may be noted that if

$${\left(1-b_0\right)}^2-\frac{4}{\mathit{\alpha \delta }}<0,$$E14

then the system possess $E_1$ as the only equilibrium point.

Further defining

$$h\left(v\right)=v\left(1-v\right)\left(v-b_0\right)$$E15

the Jacobian matrix $J$ of the system may be written as:

$$J=\left|\begin{array}{cc}\alpha h^{\prime }\left(v\right)& -1\\ \epsilon & -\mathit{\epsilon \delta }\end{array}\right|.$$E16

For the equilibrium or fixed points $\left(v_e,w_e\right)$, the eigenvalues ${\lambda }_{1,2}$ of the Jacobian matrix are given by

$${\lambda }_{1,2}=\frac{-\left(\mathit{\epsilon \delta }-\alpha b_1\right)\pm \sqrt{{\left(\mathit{\epsilon \delta }-\alpha b_1\right)}^2-4\epsilon \left(1-\alpha b_1\delta \right)}}{2}$$E17

where $b_1=h^{\prime }\left(v_e\right)$.

Therefore (1) if $\alpha b_1\delta <1$, the equilibrium point $\left(v_e,w_e\right)$ is asymptotically stable if $\alpha b_1<\mathit{\epsilon \delta }$, a repellor if $\alpha b_1>\epsilon$, (2) if $\alpha b_1\delta >1$ the equilibrium point is a saddle point and (3) if $\alpha b_1\delta =1$ then the equilibrium point is stable (unstable) if $\alpha b_1<\mathit{\delta \epsilon }\left(\alpha b_1>\mathit{\delta \epsilon }\right)$.

In case the parameters of the system are such that condition Eq. (14) holds then using Eq. (17) we find the equilibrium point $E_1$ i.e., origin, to be asymptotically stable if

$$\mathit{Re}\left[-\left(\alpha b_0+\mathit{\epsilon \delta }\right)\pm \sqrt{{\left(\alpha b_0+\mathit{\epsilon \delta }\right)}^2-4\epsilon \left(\mathit{\alpha \delta }{b}_0+1\right)}\right]<0.$$E18

Based on Routh’s criteria we may write equivalently

$$\alpha b_0+\mathit{\epsilon \delta }>0,\mathit{\alpha \delta }{b}_0+1>0.$$E19

Therefore if origin is the only fixed point of the system, then following [30], it is observed that the system has no limit cycle if ${\left(b_0-\frac{1}{2}\right)}^2+\frac{3}{4}-\frac{3\mathit{\delta \epsilon }}{\alpha }<0$. This result however assumes $\alpha b_0+\mathit{\epsilon \delta }>0$. In case $\alpha b_0+\mathit{\epsilon \delta }<0$, the origin becomes unstable and the system can be shown to exhibit one stable limit cycle. It is noted here that on varying the threshold parameter $b_0$, the system may exhibit Andronov-Hopf bifurcation when $\alpha b_0+\mathit{\epsilon \delta }=0$ as per Eq. (18) and at this point the origin of the system becomes unstable causing a bifurcation to at least one stable limit cycle. It is to be noted that gaps exist in parameter space when origin is the only asymptotically fixed point and where limit cycles may exists. Interestingly, we numerically show the existence of bistable behavior [18, 30] in terms of occurrence of double cycle bifurcation by taking $\epsilon =0.015,\delta =3.5,\alpha =1.0$ and allowing the threshold value $b_0<0$ i.e., $-0.044$(Figure 13).

Figure 13.

Phase portrait of FHN system showing bistability between limit cycle and stable fixed point.

For the case of $I\ne 0$, the equilibrium points may be one, two or three and their stability may be analyzed following the foregoing analysis. Taking the parameter values: $a=0.06,b_0=0.50,\epsilon =14$, as in [29], the phase portrait were obtained using numerical integration of the system, Eq. (13), for different $I$ values (Figure 14). The appearance of limit cycle behavior is due to supercritical Hopf bifurcation and as a consequence of loss of stability of the unique equilibrium point that exist for $I<4.2$ [29]. Figure 14 also suggest that the amplitude of limit cycles first increases and subsequently decreases with increase in values of $I$. At around $I\sim 12.45$ the second bifurcation occurs and system is led to a stable equilibrium. [29] has provided a detailed discussion on the richness of various bifurcation event as $I$ is varied.

Figure 14.

Phase portrait of FHN system showing limit cycle for $4.2\le I\le 12.45$. $\left(a\right)I=4.23,$$\left(b\right)I=4.8,$$\left(c\right)I=5.5,$$\left(d\right)I=6.8,$$\left(e\right)I=7.5,$$\left(f\right)I=9.5,$$\left(g\right)I=11.75,$$\left(h\right)I=12.42,$ and $\left(i\right)I=12.45$.

5.2 FHN neuron in the presence of external periodic electrical stimulation

Chaotic systems exhibit complexity and are sensitively dependent on initial condition of the system under investigation and also unpredictable. Chaos as a nonlinear phenomenon has attracted researchers from different disciplines e.g., physics, biology, ecology, neurobioscience etc. In this section, we investigate the effect of periodic electrical stimulation on the dynamics of an FHN system, Eq. (13).

The basic equation that governs the dynamics of FHN system in the presence of external periodic stimulation, $I\left(t\right)$, may be written as:

$$\begin{array}{l}\frac{\mathit{dv}}{\mathit{dt}}=\mathit{\alpha v}\left(\beta -v\right)\left(v-b_0\right)-\mathit{\sigma w}+I\left(t\right),\\ \frac{\mathit{dw}}{\mathit{dt}}=\epsilon \left(v-\mathit{\delta w}\right).\end{array}$$E20

Here, we take the external periodic stimulation as given by $I\left(t\right)=\left[\frac{I_0}{2\mathit{\pi \nu }}\right]cos\left(2\mathit{\pi \nu t}\right)$, where $I_0,\nu$ refers to the amplitude and frequency of the input stimulus. Further, we present the simulation results of the system, Eq. (20), by taking $\alpha =10,\beta =1,b_0=0.10,\delta =0.25,\epsilon =1$, and $\sigma =1$ and varying both $I_0$ and $\nu$. The variation of both amplitude and frequency of the external periodic stimulus is found to result in the membrane potential $v$ exhibiting regular or chaotic temporal behavior. The regular or periodic neuron spiking could be classified as $p:q$ phase-locking, where $p$ and $q$ corresponds to the number of spikes and number of periodic stimuli per unit response period. For instance Figures 15 and 16 illustrates respectively the response of the neuronal spiking corresponding to $1:1$ and $1:2$ phase locked rhythm.

Figure 15.

$1:1$ phase locking rhythm of spiking neuron. (a) Time series of membrane potential with $I_0=0.1,\nu =0.05$. (b) $v-w$ phase portrait with same parameters as in (a).

Figure 16.

$1:1$ phase locking rhythm of spiking neuron. (a) Time series of membrane potential with $I_0=0.1,\nu =0.1015$. (b) $v-w$ phase portrait with same parameters as in (a).

The response of the single FHN neuron to external periodic response could also be chaotic for certain values of the amplitude $I_0$ and frequency $\nu$ of driving stimulus i.e., i.e., $I_0=0.183,\nu =0.1931$, as shown in Figure 17.

Figure 17.

$1:1$ phase locking rhythm of spiking neuron. (a) Time series of membrane potential with $I_0=0.183,\nu =0.1931$. (b) $v-w$ phase portrait with same parameters as in (a).

The observed dynamical transition from regular/periodic to chaotic of membrane potential $v$ with increase in amplitude and frequency of external could be further seen by constructing the RPs. The construction of RP, discussed earlier in section 4, involves the reconstruction of phase space using a time series of a dynamical variable, say the membrane potential $v$, based on the information regarding the delay parameter $\tau$ and the embedding dimension $m$. The delay parameter $\tau$ for the time series of $v$ could be obtained using the method of mutual information (MI) [31] and the embedding dimension $m$ may be determined using the algorithm of [32, 33]. The time series of Figures 15–17 for the membrane potential $v$ have been used to construct the RP shown in Figure 18. The change in spiking patterns caused by external periodic stimulation from regular to chaotic is well indicated in RP of almost equally spaced diagonal lines to irregularly occurring broken diagonal lines of varying length.

Figure 18.

RP of of the membrane potential with (a) $1:1$ phase locking rhythm, (b) $1:2$ phase locking and (c) chaotic rhythm.

5.3 FHN neuron with time varying threshold

The dynamics of cortisol secretion from hypothalamus could be modeled using FHN system with time varying threshold [25]. Complexities of spike dynamics of FHN neuron has been earlier investigated in [34] incorporating the time varying threshold. Models of observed tonic firing and spike bursting were simulated by considering both periodic and noisy form of the threshold variation. The effect of a mixed mode threshold on the spiking FHN system was also investigated. Here we present results of bifurcation analysis of different states of neuronal firing of FHN neuron by considering a discrete form of the system [34, 35]. Following [34], the discrete form of the FHN system may be written as:

$$v_{n+1}=v_n+\mathrm{\Delta }\alpha \left[-v_n\left(v_n-1\right)\left(v_n-b_n\right)-w+I\right]$$E21

$$w_{n+1}=w_n+\mathrm{\Delta }\left(v_n-\delta w_n\right)$$E22

where $\mathrm{\Delta }$ refers to the integral step size and is treated here as a bifurcation parameter.

In case of mixed mode threshold variation, the membrane potential $v$ exhibits a complex behavior as shown in bifurcation diagram (Figure 19a). It is readily observed that the temporal behavior is chaotic in the region $0.42\le \mathrm{\Delta }\le 0.68$. Thereafter windows of regular and chaotic regimes are observe till $\mathrm{\Delta }=0.8$ for $I=0$. A slight increase in $I$ changes the dynamics to a more complex behavior as shown in Figure 19b wherein windows of quasi-periodic behavior that sets sets in say at $\delta \sim 0.68$ for $I=0$ (Figure 20a) make a transition to chaos as shown in Figure 20b.

Figure 19.

Bifurcation of membrane potential. (a) I = 0; (b) I = 1.

Figure 20.

Phase portrait showing quasiperiodic behavior for (a) I = 0 and chaotic behavior for (b) I = 1.

In the following section, we study the complex dynamics in economics, an area actively pursued by researchers, by introducing another measure called ’multi-scale permutation entropy’ by considering a nonlinear financial model.

6. Chaotic dynamics in finance model

In mainstream economics, economic dynamics has assume great importance in recent years in view of the availability of market and other data. Economic dynamics has therefore influenced both micro- and macroeconomics. Therefore lot of research output has poured in explaining irregular micro-economic fluctuation, erratic business cycles, irregular growth and aperiodic behavior of economic data etc. Nonlinear systems provides an alternative simple and deterministic framework that easily can explain aperiodic or chaotic behaviors of various financial systems. One of the important features of nonlinear system is that the irregular/chaotic behavior supports an endogenous mechanism for the observed complexity in economic time series. As a result nonlinear dynamic framework has been applied to economic modeling and several examples are available in [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46].

In the present work we revisit the synthetic chaotic financial model discussed in [45, 46] which is based on interest rate, investment demand and price index as dynamical variables. We numerically explored and analyze the complexity of the model using the multiscale entropy (MPE) frame work. In this section, we briefly describe the chaotic financial model and its basic characteristics. We also outlines the procedure of MPE for analyzing the complexity of the finance model.

6.1 Chaotic financial model

We consider a dynamic finance model composed of three coupled first order differential equation. This model describes the temporal evolution of the state variables viz. the interest rate $X$, the investment demand $Y$ and the price index $Z$. The model is described as [39]:

$$\begin{array}{cc}\frac{\mathit{dX}}{\mathit{dt}}& =Z+\left(Y-a\right)X,\\ \frac{\mathit{dY}}{\mathit{dt}}& =1-\mathit{bY}-X^2,\\ \frac{\mathit{dZ}}{\mathit{dt}}& =-X-\mathit{cZ}.\end{array}$$E23

Here $a$, $b$ and $c$ are positive constants and represent the saving amount, cost per investment and elasticity of demand of the commercial markets. First equation appears, representing the changes in $X$, as a result of contradiction in the investment market and structural adjustment from the goods prices. Second equation representing the changes in $Y$ appears due to proportionality to the rate of investment and also to an inversion of the cost of investment and interest rate. The third equation emerges due to contradiction between supply and demand in commercial markets which is influenced by interest rates.

Eq. (23) has been numerically integrated using fourth-order Runge-Kutta method to obtain the time series of the dynamic variables $X,Y$ and $Z$, shown in Figure 21 with $a=3.0,b=0.1,c=1.0$ and initial condition $\left(X_0,Y_0,Z_0\right)=\left(2,3,2\right)$. Similar choice of the parameter were made in [45].

Figure 21.

Temporal evolution of Finance model.

The representation of two dimensional phase portrait $\left(X,Y\right)$ and the attractor are shown in Figure 22. Obviously they represent chaotic dynamics of the temporal behavior shown in Figure 22.

Figure 22.

(a) Phase portrait of Finance model, (b) 3-D attractor of the Finance model.

6.2 Complexity analysis using multiscale permutation entropy (${MPE}$) method

The multiscale permutation method involves two steps. A “coarse graining” is applied first to a time series $X_i,i=1,N$ to construct a consecutive coarse-grained time series. The coarse-grained process involves averaging a successively increasing number of data points in non-over lapping windows. The elements of each of the coarse grained time series $y_j^{\left(s\right)}$ is computed as,

$$y_j^{\left(s\right)}=\frac{1}{s}\sum _{i=\left(j-1\right)\tau +1}^{\mathit{js}}{X}_i$$E24

where $1\le i\le N/s$ and $s$ defines the scale factor. Each time series length is of size that is an integral multiple of $N/s$. For $s=1$, the coarse-grained time series is just the original time series.

The second step involves the computation of permutation entropy [47] for each of the coarse-grained time series. For a coarse-grained time series $y_j$ we first consider the series of vector of length $m$, and obtain $S^m\left(n\right)=\left[y_n,y_{n+1},\cdots ,y_{n+m-1}\right],1\le n\le \left(N/s\right)-m+1$. Subsequently, $S^m\left(n\right)$ is arranged in an increasing order viz., $\left[y_{n+j_1+1}\le y_{n+j_2+1}\le \cdots \le y_{n+j_n+1}\right]$. For $m$ different numbers, there will be $m!$ possible order patterns/structures $\mathrm{\Pi }$ which are termed as permutations. If $f\left(\mathrm{\Pi }\right)$ denotes the frequency of order pattern $\mathrm{\Pi }$, then the relative frequency and hence the probability $p=f\left(\mathrm{\Pi }\right)/\left(N/s-m+1\right)$. The permutation entropy $H\left(m\right)$ therefore is given by

$$H\left(m\right)=-\sum _{\mathrm{\Pi }=1}^{m!}p\left(\mathrm{\Pi }\right)\mathit{ln}\left(\mathrm{\Pi }\right).$$E25

The maximum value of $H\left(m\right)$ is $\mathit{log}\left(m!\right)$ thus showing all permutations to have equal probability. Also, the time series is termed as regular if minimum value of $H\left(m\right)$ is zero. Therefore $H\left(m\right)$ the permutation entropy provide a quantitative measure of dynamical complexity of a time series as it refers to its local structures. It may be noted that the permutation entropy depends on the chosen value of $m$. For $m<3$ , there will be very few distinct states and the foregoing scheme does not work satisfactorily. In the present analysis we have considered sufficiently large time series and chosen $m=6$ to estimate the complexity measure $\mathit{MPE}$.

For the financial model, Eq. (23), we have simulated the permutation entropy as a function of the scale $s$ for $m=6$. The simulated results have been shown in Figure 23a where we observe a saturation behavior with increased value of the scale factor $s$. It is also observed that the permutation entropy at any scale $s$ for the interest rate $X$ time series is higher than the investment demand time series $Y$ which is further higher than that of the price index $Z$ time series. Therefore we may conclude that the complexity measure relation for the considered financial model time series can be expressed as ${\mathit{MPE}}_X>{\mathit{MPE}}_Y>{\mathit{MPE}}_Z$. The behavior of the complexity measure of the considered finance model has been found to be quite similar to that of chaotic Rossler attractor [with parameters $a=0.15;b=0.20;c=10.0$] (Figure 23b).

Figure 23.

Multiscale permutation of (a) Financial model time series, (b) Rossler model.

The simulation results for the multi-scale permutation entropy, $\mathit{MPE}$, presented for the financial model and the Rossler chaotic model exhibit long term correlation of the respective time series of a dynamical variable. Such inference is made in view of the increasing trend of $\mathit{MPE}$ with scale factor $s$ for a given $m$. In case of a standard financial model, the efficacy of such model could be made on comparing the $\mathit{MPE}$ trend of resulting simulated time series for interest rate ($X$), investment demand (Y) and that of price index (Z) with the availability of the real time series data for the corresponding dynamical variables. Finally, we introduce the idea of generation of time series of a nonlinear chaotic dynamical system, say a Lorenz system, using artificial neural network.

7. Time series generation using artificial neural network (ANN)

A large class of different architecture have been used in neural network for various application. Among these application an issue relates the approximation of a nonlinear mapping $f\left(\mathbf{x}\right)$ with the network $f_{\mathit{ANN}}\left(\mathbf{x}\right),\mathbf{x}\in {\mathbf{R}}^{\mathbf{K}}$ where $K$ corresponds to the size of te input. Besides the Radial Basis Function (RBF), a Multi Layer Perceptron (MLP) has been used extensively in function approximation. A MLP neural network comprises an input layer, several hidden layers and an output layer as shown in Figure 24.

Figure 24.

A multilayer perceptron network.

An $\mathit{MLP}$ comprises inputs $x_i,i=1,2,\cdots ,K$ to the neurons gets multiplied with weights $w_{\mathit{ki}}$ and summed up along with the bias ${\theta }_i$. The resulting $n_i$ is then acts as an input to the activation function $g$ which could be chosen as a sigmoid function or a $tanh$ function. The output at node $i$ is given by $y_i=g\left[{\sum }_{j=1}^K{w}_{{\mathit{jix}}_j+{\theta }_i}\right]$. Figure 24 illustrates a typical $\mathit{MLP}$ network where the output is given by

$$y_i=g\left[\sum _{j=1}^3{w}_{\mathit{ji}}^{2}g\left(\sum _{k=1}^K{w}_{\mathit{kj}}^1{x}_k+{\theta }_j^1\right)+{\theta }_j^2\right].$$E26

Several algorithms are available to determine the network parameters e.g., weights ($w_{\mathit{ji}}^k)$) and biases(${\theta }_j^k$). Such algorithms are termed as teaching or learning algorithms. The basic procedure involving the learning algorithm of an $\mathit{MLP}$ network are: (a) Define the network structure, selecting the activation function and initializing the weights and biases, (b) providing the error estimates and number of epochs for training algorithm before running the training algorithm, (c) the output is simulated using input data to the network and compared with the given output, and (d) finally validating the result with independent data.

In this work, using the inputs as $x,y$ and $z$ time series from Lorenz system exhibiting chaotic dynamics and using the $\mathit{\text{newff}},\mathit{\text{train}}$ and $\mathit{sim}$$\mathit{\text{MATLAB}}$ commands [31], we simulated each of these time series. In our simulation we take the learning parameters viz., $\mathit{net}.\mathit{\text{trainParam}}.\mathit{\text{show}}=50;$$\mathit{net}.\mathit{\text{trainParam}}.\mathit{lr}=0.05;$$\mathit{net}.\mathit{\text{trainParam}}.\mathit{\text{epochs}}=1000;$$\mathit{net}.\mathit{\text{trainParam}}.\mathit{\text{goal}}=1e^{-3}$, and use $100$ neurons and $3$ output layers [48]. Figure 25 shows the $\mathit{MLP}$ network generated time series of Lorenz variables and the corresponding deviations from the input time series.

Figure 25.

Neural network generated Lorenz time series ($\mathit{\text{black}},\mathit{red},\mathit{\text{brown}}$) and respective deviation from input time series ($\mathit{\text{blue}}$).

8. Conclusion

In this chapter, we have applied phase portrait, bifurcation diagram, Poincare surface of section, LCEs, correlation dimension, topological entropy and multi-scale permutation entropy method to unravel the complexity of various physical systems e.g., nonlinear forced pendulum, child’s swing problem, prey-predator system, periodically stimulated FHN neuron model and nonlinear financial model. Important characterization of transition from regular to chaotic dynamics have been made using the foregoing methods. Finally artificial neural network based on multi-layer perceptron network have been shown to satisfactorily generate the time series of dynamical variable of chaotic system such as Lorenz system.