A high-precision numerical method to simulating the fractional-order EI Ni ñ o chaotic systems with Riemann–Liouville fractional derivative

Abstract

Chaotic systems arise everywhere in control theory and nonlinear vibration. This paper uses a high-precision numerical approach for capturing chaotic attractors of the fractional EI Niño chaotic systems. The numerical results indicate that some of the chaotic attractors of fractional-order systems are topologically equivalent to those discovered in the integer-order systems, and some unusual attractors of fractional EI Niño chaotic systems are found by using the present method.

Introduction

The vallis model^1–2 for El Niño consists of a set of three autonomous first-order nonlinear ordinary differential equations. The El Niño chaotic dynamical system is given by Borghezan and Rech ³

$$\begin{array}{l}\{\begin{array}{l}\frac{dx\left(t\right)}{dt}=\beta y\left(t\right)-\gamma x\left(t\right)+\gamma \rho ,\hfill \\ \frac{dy\left(t\right)}{dt}=x\left(t\right)z\left(t\right)-y\left(t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t\in \left[0,T\right],\hfill \\ \frac{dz\left(t\right)}{dt}=-x\left(t\right)y\left(t\right)-z\left(t\right)+1,\hfill \end{array}\hfill \end{array}$$

(1)

with the initial conditions

$$x\left(0\right)=c_1,y\left(0\right)=c_2,z\left(0\right)=c_3,$$

(2)

which has been used to describe an anomalous ^4-6 of the Peru and Ecuador coastal waters that occurs roughly about every 3–6 years, which has a great impact on global climate. Although the model has been applied with significant degrees of success, irregular and anomalous phenomena cannot be adequately described with the local or classical derivatives, which have been suggested in many research papers, books, and monographs.^7–9 Fractional calculus (FC) can accurately describe physical process. It is a very interesting topic of investigation in the field of nonlinear research. Hence, for a more accurate representation, some authors present the fractional EI Niño chaotic systems (1) is the following form

$$\begin{array}{l}\{\begin{array}{l}{D}_t^{{\alpha }_1}x\left(t\right)=\beta y\left(t\right)-\gamma x\left(t\right)+\gamma \rho ,\hfill \\ D_t^{{\alpha }_2}y\left(t\right)=x\left(t\right)z\left(t\right)-y\left(t\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t\in \left[0,T\right],\hfill \\ D_t^{{\alpha }_3}z\left(t\right)=-x\left(t\right)y\left(t\right)-z\left(t\right)+1,\hfill \end{array}\hfill \end{array}$$

(3)

Some numerical methods have been studied for this equation. In Ref., ⁴ authors carried out numerical simulations by using estimates of the largest Lyapunov exponent to characterize the dynamical behavior of the systems (1). In Refs.,^5–6 authors researched the integer-order systems (1). Although some numerical methods of the fractional differential equations have been announced, such as spectral method,^10–13 HPM,^14–18 RKM,^19–23 and so on.^24–25 However, due to the nonlocal property of fractional derivative, these methods involved in solving the fractional chaotic dynamical system are time consuming and tedious. High-precision numerical methods have always been the direction that scholars strive to pursue. Based on this problem, this paper uses a high-precision numerical approach for the fractional-order EI Ni $\tilde{n}$ o chaotic systems. We observe many novel dynamic behaviors of fractional-order systems in numerical experiments which are unlike any that have been previously discovered in numerical experiments or theoretical studies.

Fractional vibration systems^26–32 and fractional equation^33–36 have been widely concerned. There are many definitions^37–38 of fractional derivatives (FDs). The most famous of these FD that have been widely popularized in the world of fractional calculus are the Riemann–Liouville (RL) FD, Gr $\ddot{u}$ nwald–Letnikov (GL) FD, and Caputo FD.

Definition 1.1. RLFD of order α is defined as

$${\hspace{0.17em}}_{t_0}^{RL}{D}_t^{\alpha }g\left(t\right)=\{\begin{array}{l}\frac{1}{\mathrm{\Gamma }\left(n-\alpha \right)}\frac{d^n}{d{t}^n}{\displaystyle {\int }_{t_0}^t\frac{g\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha -n+1}}}d\tau ,0\le n-1<\alpha <n,\hfill \\ \frac{d^{n}f\left(t\right)}{d{t}^n},\alpha =n\in N.\hfill \end{array}$$

(4)

Definition 1.2. GLFD of order α is defined as^37–38

$${\hspace{0.17em}}_{t_0}^{GL}\hspace{0.17em}{D}_t^{\alpha }g\left(t\right)=\underset{\tau \to 0}{\text{lim }}\frac{1}{{\tau }^{\alpha }}{\displaystyle \sum _{j=0}^{\left[\left(t-t_0\right)/\tau \right]}{\left(-1\right)}^j\left({\hspace{0.17em}}_j^{\alpha }\right)g\left(t-j\tau \right)},$$

(5)

where ${\left(-1\right)}^j\left({\hspace{0.17em}}_j^{\alpha }\right)=\frac{{\left(-1\right)}^j\mathrm{\Gamma }\left(\alpha +1\right)}{\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(\alpha -j+1\right)}$ . Note

$$\begin{array}{l}{\rho }_j^{\left(\alpha \right)}={\left(-1\right)}^j\left({\hspace{0.17em}}_j^{\alpha }\right)=\frac{{\left(-1\right)}^j\mathrm{\Gamma }\left(\alpha +1\right)}{\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(\alpha -j+1\right)},j=\mathrm{0,1},\dots ,n.\hfill \end{array}$$

(6)

Therefore, GLFD of α-order is transformed into the following form

$${\hspace{0.17em}}_{t_0}^{GL}{D}_t^{\alpha }g\left(t\right)=\underset{\tau \to 0}{\text{lim }}\frac{1}{{\tau }^{\alpha }}{\displaystyle \sum _{j=0}^{\left[\left(t-t_0\right)/\tau \right]}{\rho }_j^{\left(\alpha \right)}}g\left(t-j\tau \right).$$

(7)

Using Newton’s binomial theorem, we can know

$$\begin{array}{l}{\left(1-z\right)}^n={\displaystyle \sum _{j=0}^n{\left(-1\right)}^j}\left({\hspace{0.17em}}_j^n\right)z^j={\displaystyle \sum _{j=0}^n{\rho }_j^{\left(n\right)}}{z}^j.\hfill \end{array}$$

(8)

Note

$$\begin{array}{l}{\displaystyle \sum _{j=0}^{\infty }{\rho }_j^{\left(\alpha \right)}}{z}^j={\left(1-z\right)}^{\alpha }.\hfill \end{array}$$

(9)

First-order generating function (GF) of GLFD is (1 − z)^α.

$${\hspace{0.17em}}_{t_0}^{GL}{D}_t^{\alpha }y\left(t\right)=\begin{array}{l}\frac{1}{{\tau }^{\alpha }}{\displaystyle \sum _{i=0}^{\left[\left(t-t_0\right)/\tau \right]}{\rho }_j^{\left(\alpha \right)}}y\left(t-i\tau \right)+o\left(\tau \right).\hfill \end{array}$$

(10)

GLFD and RLFD have the following relationship

$${\hspace{0.17em}}_{t_0}^{GL}{D}_t^{\alpha }y\left(t\right)={\hspace{0.17em}}_{t_0}^{RL}{D}_t^{\alpha }y\left(t\right).$$

(11)

For convenience of expression, we denote

$${\hspace{0.17em}}_{t_0}^{GL}{D}_t^{\alpha }g\left(t\right){\hspace{0.17em}=\hspace{0.17em}}_{t_0}^{RL}{D}_t^{\alpha }g\left(t\right)=D_t^{\alpha }g\left(t\right).$$

In order to obtain higher precision, author has given the construction method of GF.^37,39-41

Numerical method

Definition 2.1. A polynomial function of n-order f _n (t) is defined as $f_n\left(t\right)={\displaystyle \sum _{k=1}^{n}1/k{\left(1-t\right)}^k}$ .

Theorem 2.2. The polynomial function f _n (t) could be written as $f_n\left(t\right)={\displaystyle \sum _{k=0}^n{f}_k}{t}^k$ , where f _k is the solution of the following equations

$$\begin{array}{l}\{\begin{array}{l}\begin{array}{ll}{f}_0\hfill & +f_1+f_2+\cdots +f_n=0,\hfill \\ f_0\hfill & +2f_1+3f_2+\cdots +\left(n+1\right)f_n=-1,\hfill \\ f_0\hfill & +2^2{f}_1+3^2{f}_2+\cdots +{\left(n+1\right)}^2{f}_n=-2,\hfill \\ \hfill & ⋮\hfill \\ f_0\hfill & +2^n{f}_1+3^n{f}_2+\cdots +{\left(f+1\right)}^n{f}_n=-n.\hfill \end{array}\hfill \end{array}\hfill \end{array}$$

Proof. For the details of the proof, one may refer to Refs.^37‐41

Definition 2.3. A GF $f_n^{\alpha }\left(t\right)$ of GLFD is defined as

$$\begin{array}{c}{f}_n^{\alpha }\left(t\right)={\left(f_0+f_{1}t+\cdots +f_n{t}^n\right)}^{\alpha }.\hfill \end{array}$$

(12)

Theorem 2.4. The GF $f_n^{\alpha }\left(t\right)$ can be written as

$$\begin{array}{l}{f}_n^{\alpha }\left(t\right)={\displaystyle \sum _{k=0}^{\infty }{\rho }_k^{(\alpha ,n)}}{t}^k.\hfill \end{array}$$

(13)

where ρ _k ^{( α,n )} can be calculated by the following formula

$$\{\begin{array}{l}{\rho }_k^{\left(\alpha ,n\right)}=f_0^{\alpha },k=0,\hfill \\ {\rho }_k^{\left(\alpha ,n\right)}=-\hspace{0.17em}\frac{1}{p_0}{\displaystyle \sum _{i=1}^n{f}_i}\left(1-i\frac{1+\alpha }{k}\right){\rho }_{k-i}^{\left(\alpha ,n\right)},k=\mathrm{1,2},\dots ,\hfill \\ {\rho }_k^{\left(\alpha ,n\right)}=0,k<0.\hfill \end{array}$$

(14)

Proof. In equation (12), if t = 0, then it has ${\rho }_0^{\left(\alpha ,n\right)}=f_0^{\alpha }$ . Rewriting equation (12), it can be found that

$$\begin{array}{l}{(f_0+f_{1}t+\cdots +f_n{t}^n)}^{\alpha }={\displaystyle \sum _{k=-\infty }^{\infty }{\rho }_k^{\left(\alpha ,n\right)}}{t}^k,\hfill \end{array}$$

(15)

where k < 0 , then ${\rho }_k^{(\alpha ,n)}=0$ .

Calculating the first-order derivative on both sides of equation (15), it can be obtained that

$$\begin{array}{l}\alpha \left(f_1+2f_{2}t+\cdots +n{f}_n{t}^{n-1}\right){\left(f_0+f_{1}t+\cdots +f_n{t}^n\right)}^{\alpha -1}={\displaystyle \sum _{k=-\infty }^{\infty }k}{\rho }_k^{(\alpha ,n)}{t}^{k-1},\hfill \end{array}$$

(16)

Multiplying ( f ₀ + f ₂ z + ⋯ + f _n t ⁿ ) on both sides of equation (16), it can be obtained that

$$\begin{array}{l}\alpha \left(f_1+\cdots +n{f}_n{t}^{n-1}\right){\left(f_0+\cdots +f_n{t}^n\right)}^{\alpha }=\left(f_0+\cdots +f_n{t}^n\right){\displaystyle \sum _{k=-\infty }^{\infty }k}{\psi }_k^{(\alpha ,n)}{t}^{k-1}.\hfill \end{array}$$

(17)

Substituting equation (15) into equation (17), we can know

$$\begin{array}{l}\alpha \left(f_1+\cdots +n{f}_n{t}^{n-1}\right){\displaystyle \sum _{k=-\infty }^{\infty }{\rho }_k^{(\alpha ,n)}}{t}^k=\left(f_0+\cdots +f_n{t}^n\right){\displaystyle \sum _{k=-\infty }^{\infty }k}{\rho }_k^{(\alpha ,n)}{t}^{k-1}.\hfill \end{array}$$

(18)

Using the $t^d{\rho }_k^{\left(\alpha ,n\right)}={\rho }_{k-d}^{\left(\alpha ,n\right)}$ property, the left side of equation (18) can be written as

$$\begin{array}{l}{\displaystyle \sum _{k=-\infty }^{\infty }\alpha }\left(f_1{\rho }_k^{(\alpha ,n)}+2f_2{\rho }_{k-1}^{(\alpha ,n)}+\cdots +n{f}_n{\rho }_{k-n+1}^{(\alpha ,n)}\right)t^k,\hfill \end{array}$$

(19)

then its right side can be written as

$$\begin{array}{l}{\displaystyle \sum _{k=-\infty }^{\infty }\left[\left(k+1\right)f_0{\rho }_{k+1}^{(\alpha ,n)}+k{f}_1{\rho }_k^{(\alpha ,n)}+\cdots +\left(k-n+1\right)f_n{\rho }_{k-n+1}^{(\alpha ,n)}\right]}{t}^k.\hfill \end{array}$$

(20)

Comparing the coefficient of t in equation (19), we can get

$$\begin{array}{l}\alpha \left(f_1{\rho }_k^{(\alpha ,n)}+\cdots +n{f}_n{\rho }_{k-n+1}^{(\alpha ,n)}\right)=\left(k+1\right)f_0{\rho }_{k+1}^{(\alpha ,n)}+\cdots +\left(k-n+1\right)f_n{\rho }_{k-n+1}^{(\alpha ,n)}.\hfill \end{array}$$

(21)

Moving equation (21) back one step, then the equation can become

$$\begin{array}{l}{f}_{0}k{\rho }_k^{(\alpha ,n)}+f_1\left(k-1-\alpha \right){\rho }_{k-1}^{(\alpha ,n)}+\cdots +f_n\left(k-n-n\alpha \right){\rho }_{k-n}^{(\alpha ,n)}=0.\hfill \end{array}$$

(22)

If k ≠ 0, it is thus clear that

$$\begin{array}{l}{\rho }_k^{(\alpha ,n)}=-\frac{1}{f_0}\left[f_1\left(1-\frac{1+\alpha }{k}\right){\rho }_{k-1}^{(\alpha ,n)}+f_2\left(1-2\frac{1+\alpha }{k}\right){\rho }_{k-2}^{(\alpha ,n)}+\cdots +f_n\left(1-n\frac{1+\alpha }{k}\right){\rho }_{k-n}^{(\alpha ,n)}\right].\hfill \end{array}$$

(23)

where, when k < 0 , ${\rho }_k^{(\alpha ,n)}=0$ . So the theorem is proved. For the details of the proof, one may refer to Refs.^37‐41

Corollary 2.5. The GF $f_n^{\alpha }\left(t\right)$ of GLFD could be written as

$$\begin{array}{l}{f}_n^{\alpha }\left(t\right)={\displaystyle \sum _{j=0}^{\infty }{\rho }_j^{(\alpha ,n)}}{t}^j.\hfill \end{array}$$

(24)

where

$$\{\begin{array}{l}{\rho }_k^{\left(\alpha ,n\right)}=f_0,k=0,\hfill \\ {\rho }_k^{\left(\alpha ,n\right)}=-\hspace{0.17em}\frac{1}{f_0}{\displaystyle \sum _{i=0}^{k-1}\left(1-i\frac{1+\alpha }{k}\right)}{f}_i{\rho }_{k-i}^{\left(\alpha ,n\right)},k=\mathrm{1,2},\dots .,n-1,\hfill \\ {\rho }_k^{\left(\alpha ,n\right)}=-\hspace{0.17em}\frac{1}{f_0}{\displaystyle \sum _{i=0}^n\left(1-i\frac{1+\alpha }{k}\right)}{f}_i{\rho }_{k-i}^{\left(\alpha ,n\right)},k=n,n+1,n+2,\dots .,\hfill \\ {\rho }_k^{\left(\alpha ,n\right)}=0,k<0,\hfill \end{array}$$

(25)

GLFD with n-order GF $f_n^{\alpha }\left(t\right)$ is given as

$${\hspace{0.17em}}_0{D}_t^{\alpha }g\left(t\right)=\underset{\tau \to 0}{\text{lim}}\frac{1}{{\tau }^{\alpha }}{\displaystyle \sum _{i=0}^{[t/\tau ]}{\rho }_i^{(\alpha ,n)}}g\left(t-i\tau \right).$$

(26)

Applying (26), an approximate computation scheme of GLFD with n-order GF $f_{\left(t\right)}^{\alpha }$ is given as

$$D_t^{\alpha }g\left(t\right)\simeq y_m=\begin{array}{l}\frac{1}{{\tau }^{\alpha }}{\displaystyle \sum _{j=0}^m{\rho }_j^{(\alpha ,n)}}g\left(t-j\tau \right)=\frac{1}{{\tau }^{\alpha }}g\left(t_k\right)-\frac{1}{{\tau }^{\alpha }}{\displaystyle \sum _{i=1}^m{\rho }_j^{(\alpha ,n)}}g\left(t_{k-j}\right).\hfill \end{array}$$

(27)

Next, we consider the following nonlinear fractional ordinary differential equations

$$\{\begin{array}{l}{D}_t^{{\alpha }_1}{z}_1\left(t\right)=f_1\left(t,z_1\left(t\right),z_2\left(t\right)\dots ,z_n\left(t\right)\right),\hfill \\ D_t^{{\alpha }_2}{z}_2\left(t\right)=f_2\left(t,z_1\left(t\right),z_2\left(t\right)\dots ,z_n\left(t\right)\right),\hfill \\ ⋮\hfill \\ D_t^{{\alpha }_n}{z}_n\left(t\right)=f_n\left(t,z_1\left(t\right),z_2\left(t\right)\dots ,z_n\left(t\right)\right),\hfill \end{array}$$

(28)

subject to the initial condition given by ${[z_1\left(0\right),z_2\left(0\right),\dots ,z_n\left(0\right)]}^T={[c_1,c_2,\dots ,c_n]}^T$ .

Constructing the following linearized iterative method

$$\{\begin{array}{l}{z}_1\left(t_k\right)={\tau }^{{\alpha }_i}{f}_1\left(t_k,z_1\left(t_{k-1}\right),z_2\left(t_{k-1}\right),\dots ,z_n\left(t_{k-1}\right)\right)+z_1\left(0\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{\left({\alpha }_1,p\right)}\left(z_1\left(t_{k-j}\right)-z_1\left(0\right)\right)},\hfill \\ z_2\left(t_k\right)={\tau }^{{\alpha }_i}{f}_2\left(t_k,z_1\left(t_k\right),z_2\left(t_{k-1}\right),\dots ,z_n\left(t_{k-1}\right)\right)+z_2\left(0\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{\left({\alpha }_2,p\right)}\left(z_2\left(t_{k-j}\right)-z_2\left(0\right)\right)},\hfill \\ ⋮\hfill \\ z_n\left(t_k\right)={\tau }^{{\alpha }_i}{f}_n\left(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_{k-1}\right)\right)+z_n\left(0\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{\left({\alpha }_n,p\right)}\left(z_n\left(t_{k-j}\right)-z_n\left(0\right)\right).}\hfill \end{array}$$

(29)

Theorem 2.6. If f₁ (t, z₁(t), z₂(t), …, z_n(t)) is a continuous function on 0 ≤ t ≤ T, − ∞ ≤ z_j ≤ ∞ and satisfies the local Lipschitz condition, then the solution of equation (35) is convergent.

Proof. In this section, for convenience of expression, we assume that ${[z_1\left(0\right),z_2\left(0\right),\dots ,z_n\left(0\right)]}^T={[\mathrm{0,0},\dots ,0]}^T$ .

In equation (28), let t = t_k, we get

$$\{\begin{array}{l}{D}_t^{{\alpha }_1}{z}_1\left(t_k\right)=f_1\left(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_k\right)\right),\hfill \\ D_t^{{\alpha }_2}{z}_2\left(t_k\right)=f_2\left(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_k\right)\right),\hfill \\ ⋮\hfill \\ D_t^{{\alpha }_n}{z}_n\left(t_k\right)=f_n\left(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_k\right)\right).\hfill \end{array}$$

(30)

Substituting (26) into (30), we get

$$\{\begin{array}{c}\underset{\tau \to 0}{\text{lim}}\frac{1}{{\tau }^{{\alpha }_1}}{\displaystyle \sum _{j=0}^{[(t-t_0)/h]}{\rho }_i^{({\alpha }_1,p)}}{z}_1\left(t_{k-i}\right)=f_1(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_k\right)),\\ \underset{\tau \to 0}{\text{lim}}\frac{1}{{\tau }^{{\alpha }_2}}{\displaystyle \sum _{j=0}^{[(t-t_0)/h]}{\rho }_i^{({\alpha }_2,p)}}{z}_2\left(t_{k-i}\right)=f_2(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_k\right)), \\ ⋮\\ \underset{\tau \to 0}{\text{lim}}\frac{1}{{\tau }^{{\alpha }_n}}{\displaystyle \sum _{j=0}^{[(t-t_0)/h]}{\rho }_i^{({\alpha }_2,p)}}{z}_n\left(t_{k-i}\right)=f_n(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_k\right)).\end{array}$$

(31)

Using the approximate computation scheme (27), we obtain

$$\{\begin{array}{l}\frac{1}{{\tau }^{{\alpha }_1}}{\displaystyle \sum _{i=0}^m{\rho }_i^{\left({\alpha }_1,p\right)}}{\tilde{z}}_1\left(t_{k-i}\right)=f_1\left(t_k,{\tilde{z}}_1\left(t_k\right){\tilde{z}}_2\left(t_k\right),\dots ,{\tilde{z}}_n\left(t_k\right)\right),\hfill \\ \frac{1}{{\tau }^{{\alpha }_2}}{\displaystyle \sum _{i=0}^m{\rho }_i^{\left({\alpha }_2,p\right)}}{\tilde{z}}_2\left(t_{k-i}\right)=f_2(t_k,{\tilde{z}}_1\left(t_k\right){\tilde{z}}_2\left(t_k\right),\dots ,{\tilde{z}}_n\left(t_k\right)),\hfill \\ ⋮\hfill \\ \frac{1}{{\tau }^{{\alpha }_n}}{\displaystyle \sum _{i=0}^m{\rho }_i^{{\hspace{0.17em}}_i^{\left({\alpha }_2,p\right)}}}{\tilde{z}}_n\left(t_{k-i}\right)=f_n(t_k,{\tilde{z}}_1\left(t_k\right){\tilde{z}}_2\left(t_k\right),\dots ,{\tilde{z}}_n\left(t_k\right)).\hfill \end{array}$$

(32)

It is easy to see by induction that

$$\{\begin{array}{l}{\tilde{z}}_1\left(t_k\right)={\tau }^{{\alpha }_1}{f}_1\left(t_k,{\tilde{z}}_1\left(t_k\right),{\tilde{z}}_2\left(t_k\right),\dots ,{\tilde{z}}_n\left(t_k\right)\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_1,p)}}{\tilde{z}}_1\left(t_{k-j}\right),\hfill \\ {\tilde{z}}_2\left(t_k\right)={\tau }^{{\alpha }_2}{f}_2\left(t_k,{\tilde{z}}_1\left(t_k\right),{\tilde{z}}_2\left(t_k\right),\dots ,{\tilde{z}}_n\left(t_k\right)\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_2,p)}}{\tilde{z}}_2\left(t_{k-j}\right),\hfill \\ ⋮\hfill \\ {\tilde{z}}_n\left(t_k\right)={\tau }^{{\alpha }_n}{f}_n\left(t_k,{\tilde{z}}_1\left(t_k\right),{\tilde{z}}_2\left(t_k\right),\dots ,{\tilde{z}}_n\left(t_k\right)\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_n,p)}}{\tilde{z}}_n\left(t_{k-j}\right).\hfill \end{array}$$

(33)

So, constructing the following linearized iterative method

$$\{\begin{array}{l}{z}_1\left(t_k\right)={\tau }^{{\alpha }_1}{f}_1(t_k,z_1\left(t_{k-1}\right),z_2\left(t_{k-1}\right),\dots ,z_n\left(t_{k-1}\right))-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_1,p)}}{z}_1\left(t_{k-j}\right),\hfill \\ z_2\left(t_k\right)={\tau }^{{\alpha }_2}{f}_2(t_k,z_1\left(t_k\right),z_2\left(t_{k-1}\right),\dots ,z_n\left(t_{k-1}\right))-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_2,p)}}{z}_2\left(t_{k-j}\right),\hfill \\ ⋮\hfill \\ z_n\left(t_k\right)={\tau }^{{\alpha }_n}{f}_n(t_k,z_1\left(t_k\right),z_2\left(t_k\right),\dots ,z_n\left(t_{k-1}\right))-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_n,p)}}{z}_n\left(t_{k-j}\right).\hfill \end{array}$$

(34)

From (2) and (33), we have the Theorem (2.6).

So, numerical approach of the fractional-order EI Ni $\tilde{n}$ o chaotic systems (3) is given by

$$\{\begin{array}{l}x\left(t_k\right)=\beta {\tau }^{{\alpha }_1}\left(\left(y\left(t_{k-1}\right)\right)-\gamma x\left(t_k\right)\right)+\gamma \rho x\left(0\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_1,p)}}\left(x\left(t_{k-1}\right)-x\left(0\right)\right),\hfill \\ y\left(t_k\right)={\tau }^{{\alpha }_2}\left(\left(-y\left(t_{k-1}\right)\right)+x\left(t_k\right)z\left(t_{k-1}\right)\right)+y\left(0\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_2,p)}\left(y\left(t_{k-1}\right)-y\left(0\right)\right)},\hfill \\ z\left(t_k\right)={\tau }^{{\alpha }_3}\left(-x\left(t_k\right)y\left(t_k\right)-\beta z\left(t_{k-1}\right)\right)+z\left(0\right)-{\displaystyle \sum _{j=1}^m{\rho }_j^{({\alpha }_3,p)}}\left(z\left(t_{k-1}\right)-z\left(0\right)\right).\hfill \end{array}$$

(35)

Numerical experiment

We consider the model (3) with the initial conditions x(0) = y(0) = z(0) = 0. We set h = 0.01, T = 100, p = 20, the numerical results with different parameters and different fractional derivatives are shown in Figure 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10. Figure 1 shows numerical results at α = [1, 1, 1], [β, γ, ρ] = [140, 3, 0.45, 1], it is the integer-order systems. The simulation results are consistent with those of other methods, and confirm the validity of the present method. Next, we set different parameters and different fractional derivatives to capture chaotic attractors of the fractional EI Niño chaotic systems. We set parameters α = [1.1, 1.1, 0.95], [β, γ, ρ] = [102, 3, 0.45, 1], projected on the (x, y), (x, z), (y, z)-plane and Phase diagram are shown in Figure 3. Figure 2 shows projected on the (x, y), (x, z), (y, z)-plane and Phase diagram in parameters α = [1.1, 1.1, 0.95], [β, γ, ρ] = [102, 3, 0.45, 1]. The three-dimensional space graphs of the systems (3) with different parameters are shown in Figure 4. Time series plots and three-dimensional space graphs are shown Figure 9. We set parameters α = [0.95, 0.95, 1.05], [β, γ, ρ] = [102, 3, 0.45, 1], projected on the (x, y), (x, z), (y, z)-plane are shown in Figure 9. We set parameters α = [0.9, 1.1, 0.95], [β, γ, ρ] = [120, 5, 0.45, 1], projected on the (x, y), (x, z), (y, z)-plane are shown in Figure 10.

Figure 1.

Numerical results for the systems (3) at α = [1, 1, 1], [β, γ, ρ] = [140, 3, 0.45, 1].

Figure 2.

Numerical results for the systems (3) at α = [1.1, 1.1, 0.95], [β, γ, ρ] = [102, 3, 0.45, 1].

Figure 3.

Numerical results for the systems (3) at α = [1.1, 1.1, 0.96], [β, γ, ρ] = [102, 3, 0.45].

Figure 4.

The three-dimensional space graphs of the systems (3) with different parameters.

Figure 5.

Numerical results for the systems (3) at α = [0.95, 0.95, 1.05], [β, γ, ρ] = [102, 3, 0.45].

Figure 6.

Numerical results for the systems (3) at α = [1.1, 1.1, 0.95], [β, γ, ρ] = [100, 3, 0.45].

Figure 7.

Numerical results for the systems (3) at α = [0.79, 0.9, 1.1], [β, γ, ρ] = [140, 15, 0.45].

Figure 8.

Numerical results for the systems (3) at α = [0.9, 1.1, 0.95], [β, γ, ρ] = [120, 5, 0.45].

Figure 9.

Numerical results for the systems (3) at α = [0.95, 0.95, 1.05], [β, γ, ρ] = [102, 3, 0.45, 1].

Figure 10.

Numerical results for the systems (3) at α = [0.9, 1.1, 0.95], [β, γ, ρ] = [120, 5, 0.45, 1].

Conclusions and remarks

This paper effectively captures chaotic attractors of the fractional EI Ni $\tilde{n}$ o chaotic systems in different parameters β, γ, ρ and different fractional derivatives α₁, α₂, α₃ by using a high-precision numerical approach. The numerical results indicate that some of the chaotic attractors of fractional-order EI Niño systems are topologically equivalent to those discovered in the integer-order systems, and some unusual attractors of fractional-order chaotic systems are found by using the present method.