La Niña-Southern oscillation model

Journal Pre-proof Multiple attractors and periodicity on the Vallis model for El Niño/La Niña-Southern oscillation model J.F. Gómez-Aguilar PII:

S1364-6826(19)30440-7

DOI:

https://doi.org/10.1016/j.jastp.2019.105172

Reference:

ATP 105172

To appear in:

Journal of Atmospheric and Solar-Terrestrial Physics

Received Date: 7 August 2019 Revised Date:

11 November 2019

Accepted Date: 12 November 2019

Please cite this article as: Gómez-Aguilar, J.F., Multiple attractors and periodicity on the Vallis model for El Niño/La Niña-Southern oscillation model, Journal of Atmospheric and Solar-Terrestrial Physics, https://doi.org/10.1016/j.jastp.2019.105172. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.

Multiple attractors and periodicity on the Vallis model for El Ni˜no/La Ni˜na-Southern oscillation model J.F. G´omez-Aguilar CONACyT-Tecnol´ogico Nacional de M´exico/CENIDET. Interior Internado Palmira S/N, Col. Palmira, C.P. 62490, Cuernavaca, Morelos, M´exico. [email protected] In this paper, we obtain multiple attractors and periodicity using differential and integral operators with power-law and Mittag-Leffler law for the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model and the continuous-time Vallis model for El Ni˜ no. Also, we consider the extension of these models considering a stochastic approach where the given parameters are converted to normal distributions. Additionally, we consider for both models novel differential and integral operators with fractional order and fractal dimension. These novel operators predict chaotic behaviors involving the fractal derivative in convolution with power-law and the Mittag-Leffler function, also, these operators can capture self-similarities for both chaotic attractors. We have presented the conditions of existence of uniquely exact solutions of the system using the fixed-point theorem approach. Each model is solved numerical via the Adams-Bashforth-Moulton, Adams-Moulton and the Atangana-Toufik schemes. We presented numerical simulations for different values of fractional order to show the applicability and computational efficiency of these methods. The results obtained presents more information that were not revealed in the models with local derivative. Keywords Fractional calculus; El Ni˜ no/La Ni˜ na-Southern oscillation model; Mathematical modelling; Liouville-Caputo fractional derivative; Atangana-Baleanu fractional derivative; Fractal-fractional operators.

1

Introduction

Fractional calculus (FC) can accurately describe physical process since it has historical memory and spatial global correlations. It is a very interesting topic of investigation in the field of nonlinear research and has extensive applications in electromagnetic oscillations, system control, physics, chaos, economics, engineering, biochemistry, microbiology, epidemiology and other applications [1]-[20]. Three different type of fractional operators involving power law, exponential decay law and the generalized Mittag-Leffler function law had been constructed. The Riemann-Liouville-Caputo fractional operator consider the power law. Caputo and Fabrizio consider the exponential decay law to developed a new operator non-singular but with local kernel. Atangana and Baleanu generalized the exponential decay law and developed a new fractional operator characterized by non-singular and non-local kernel based in the Mittag-Leffler law [21]-[23]. Recently, Atangana introduced the fractal-fractional calculus, this kind of calculus with power law, exponential decay law and the Mittag-Leffler kernel generalized the classical differentiation considering the concept of fractal differentiation. These operators involving naturally two fractional parameters, one for the fractional order (describing power-law, exponential decay and Mittag-leffler memories) and other describe the fractal dimension (involving the fractal derivative) [24]-[25]. El Ni˜ no/La Ni˜ na-Southern oscillation is a quasiperiodic climate pattern that arises across the tropical pacific ocean on the coast of Peru and Ecuador every five years. It is combined with two phases, El Ni˜ no describe the warm oceanic phase, and La Ni˜ na describe the cold phase. This model permits analyze the global ocean climate and atmospheric physics [26]-[33]. Another mathematical model that permits describe El Ni˜ no phenomena is the continuous-time Vallis model, which consists of a set of three autonomous first-order nonlinear ordinary differential equations. This model is described in [34] for Borghezan and Rech, the authors studied chaos and periodicity for this model, also reported the existence of periodic structures embedded in a chaotic region. Alkahtani and

1

Atangana consider this model and comparing the classical and fractional behaviors involving the Liouville-Caputo and Caputo-Fabrizio derivatives. The authors presented the existence and unique solutions of the model [35]. Considering the stochastic approach, the Liouville-Caputo and Atangana-Baleanu derivatives, the aim of the present work is to obtain novel chaotic behaviors for El Ni˜ no/La Ni˜ na-Southern oscillation model and the continuoustime vallis model for el Ni˜ no. We studied in detail the exactness and uniqueness of the solutions of these nonlinear models by using the fixed-point theorem. Also, considering fractal-fractional operators with power law and MittagLeffler kernel we obtain novel chaotical behaviors. These models based on fractional derivatives will be more descriptive than the one with classical derivative.

2 2.1

El Ni˜ no/La Ni˜ na-Southern oscillation models Classical and random simulations.

Classical model. The coupled dynamical El Ni˜ no/La Ni˜ na-Southern Oscillation (ENSO) model was observed to recount the oscillating physical mechanism as given by [26]-[30]. This model is defined as [31] x(t) ˙ = Cx(t) + ηy(t) − ξx3 (t), y(t) ˙ = −ςx(t) − γy(t),

(1)

where C, η, γ and ς are physical constants,  represents a perturbation coefficient which is defined in the interval 0 <  < 1. x(t) represents the temperature of the eastern equatorial pacific sea surface and y(t) represents the thermo-cline depth anomaly. For Eq. (1), we consider the following dimensionless parameters set: C = 1, η = 1; 1.5; 2, ς = 1; 2; 3, γ = 1; 0.85; 0.70 and ξ = 0.1; 0.3; 0.6. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1 and y(0) = 1. The system (1) was solved numerically using the fourth-order Runge-Kutta integration scheme. The numerical results given in Figs. 1a-1i shows the dynamics and phase portraits for the classical coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model. Random coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model. Now, we consider that the coefficients considered in the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model are computed via statistical analysis of data got from different countries that make use of this model. In the real life situations, the coefficients of this model are not constant and therefore change in every different situation. These coefficients are random numbers in nature, therefore we present the random coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model transforming the coefficients of system (1) to random variables [36]. Hence, by using normal distribution, we will obtain the new set of random coefficients A ∼ N (a1 , s21 );

B ∼ N (a2 , s22 );

C ∼ N (a3 , s23 );

D ∼ N (a4 , s24 );

E ∼ N (a5 , s25 ),

where (ai , si ), i ∈ (1, 5) corresponds to the means and the standard deviations of the normal distributions, respectively. Random variables that are expected to be a sum of independent quantities often have a normal distribution. The mean values of these distributions will be chosen according to the numerical values of the coefficients. The distributed random variable Ξ = N (ai , s2i ), i ∈ (1, 5) can be written as Ξ = ai + si χi , i ∈ (1, 5), where χ ∼ N (0, 1) is the standard normally distributed random variable. Now, these random variables can be rewritten in the following manner A = a1 + s1 χ1 ;

B = a2 + s2 χ2 ;

C = a3 + s3 χ3 ;

D = a4 + s4 χ4 ;

E = a5 + s5 χ5 ,

where the initial conditions for this model are x(0) = 1 and y(0) = 1. For the independent random variables χi , i ∈ (1, 5), the distribution is N (0, 1). Setting the appropriate values of (ai , si ), yields A = 1 + 0.03χ1 ;

B1 = 1 + 0.038χ2 ;

B2 = 1.5 + 0.04χ2 ; 2

B3 = 2 + 0.043χ2 ,

2.5

1.5

1.5

2

1

1

0.5

0.5

0

0

1.5

0.5

y(t)

y(t)

x(t)

1

-0.5

-0.5

-1

-1

-1.5

-1.5

0 -0.5 -1 -1.5

-2 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

120

140

160

180

-2 -1.5

200

1

1

1

0.5

0.5

0.5

0

y(t)

1.5

0

-0.5

-0.5

-1

-1

-1

-1.5

-1.5 60

80

100

120

0.5

140

160

180

200

0

20

40

60

80

t (days)

100

120

140

160

180

200

t (days)

(d)

1

1

0.5

0.5

y(t)

1.5

0

2

-1.5 -1.5

-1

-0.5

0

0.5

1

(f)

0

-0.5

-0.5

-1

-1

-1.5

1.5

x(t)

(e)

1.5

1

-1.5 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

100

120

140

160

180

200

t (days)

(g)

(h)

(i)

Figure 1: Numerical simulation of classical coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model given by Eq. (1). In (a)-(c) C = 1, η = 1, ς = 1, γ = 1, and ξ = 0.1. In (d)-(f) C = 1, η = 1.5, ς = 2, γ = 0.85, and ξ = 0.3. In (g)-(i) C = 1, η = 2, ς = 3, γ = 0.70, and ξ = 0.6.

3

2.5

0

-0.5

40

0

(c)

1.5

20

-0.5

x(t)

1.5

0

-1

(b)

y(t)

x(t)

(a)

x(t)

100

t (days)

1.5

C1 = 1 + 0.016χ3 ; D1 = 1 + 0.034χ4 ; ;

C2 = 2 + 0.039χ3 ;

C3 = 3 + 0.017χ3 ,

D2 = 0.85 + 0.0235χ4 ;

E1 = 0.1 + 0.0363χ5 ;

(2)

D3 = 0.70 + 0.0363χ4 ,

E2 = 0.3 + 0.0235χ5 ;

E3 = 0.6 + 0.0175χ5 .

Substituting the parameters described in (2) into the model (1), we get dx(t) = Ax(t) + Bj y(t) − Ej x3 (t), , dt dy(t) = −Cj x(t) + Dj y(t). (3) dt The solution of the model (3) can be obtained applying the Adams-Bashforth method [37]. The Numerical scheme is given by o 1 n o 3 n xn+1 = x0 + h Axn (tn ) + Bj yn (tn ) − Ej x3n (tn ) − h Axn−1 (tn−1 ) + Bj yn−1 (tn−1 ) − Ej x3n−1 (tn−1 ) , 2 2 n o o 3 1 n (4) yn+1 = y0 + h − Cj xn (tn ) + Dj yn (tn ) − h − Cj xn−1 (tn−1 ) + Dj yn−1 (tn−1 ) , 2 2 where j ∈ (1, 3). Now, we consider the dimensionless parameters given by Eq. (2), the time simulation t = 200 [days], step size h = 1 × 10−3 , and initial conditions: x(0) = 1 and y(0) = 1. The system (3) was solved numerically using the Adams-Bashforth method (4). The numerical results given in Figs. 2a-2i shows the dynamics and phase portraits for the random coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model. The continuous-time vallis model for El Ni˜ no consists of a set of three autonomous first-order nonlinear ordinary differential equations. Therefore, the vallis model for El Ni˜ no phenomenon is given by [34] x(t) ˙ = βy(t) − C(x(t) + p), y(t) ˙ = x(t)z(t) − y(t),

(5)

z(t) ˙ = −x(t)y(t) − z(t) + 1, where x(t) represents the current generated by the temperature gradient, y(t) represents the difference of temperatures of the east and west parts of the ocean, and z(t) is the sum of these temperatures. For Eq. (5), we consider the following dimensionless parameters set: β = 102, C = 3 and p = 0.5. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1, y(0) = 1 and z(0) = 1. The system (5) was solved numerically using the fourth-order Runge-Kutta integration scheme. The numerical results given in Figs. 3a-3i shows the dynamics and phase portraits for the classical continuous-time vallis model for El Ni˜ no. Random continuous-time vallis model for El Ni˜ no. Now, for this model we also consider that the coefficients involved in the vallis model for El Ni˜ no are computed via statistical analysis of data got from different countries that make use of this model. Hence, by using normal distribution, we will obtain the new set of random coefficients A ∼ N (a1 , s21 );

B ∼ N (a2 , s22 );

C ∼ N (a3 , s23 ),

where (ai , si ), i ∈ (1, 3) corresponds to the means and the standard deviations of the normal distributions, respectively. The distributed random variable Ξ = N (ai , s2i ), i ∈ (1, 3) can be written as Ξ = ai + si χi , i ∈ (1, 3), where χ ∼ N (0, 1) is the standard normally distributed random variable. Now, these random variables can be rewritten in the following manner A = a1 + s1 χ1 ; B = a2 + s2 χ2 ; C = a3 + s3 χ3 , where the initial conditions for this model are x(0) = 1; y(0) = 1 and z(0) = 1.

4

2.5

1.5

1.5

2

1

1

0.5

0.5

0

0

1.5

0.5

y(t)

y(t)

x(t)

1

-0.5

-0.5

0 -0.5 -1 -1.5

-1

-1

-1.5

-1.5

-2 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

100

120

140

160

180

-2 -1.5

200

-1

-0.5

0

t (days)

(a)

(b)

1.5

1

0.5

1

1.5

2

2.5

x(t)

(c)

1.5

1.5

1

1

0.5

0.5

0

0

0

y(t)

y(t)

x(t)

0.5

-0.5

-0.5

-1

-1

-1.5

-1.5

-0.5

-1

-1.5

-2 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

(d)

120

140

160

180

-2 -1.5

200

1

1

1

0.5

0.5

0.5

y(t)

1.5

0

0

-0.5

-0.5

-1

-1

-1

-1.5 40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

100

t (days)

(g)

(h)

0.5

1

1.5

0.5

1

1.5

0

-0.5

-1.5

0

(f)

1.5

20

-0.5

x(t)

1.5

0

-1

(e)

y(t)

x(t)

100

t (days)

120

140

160

180

200

-1.5 -1.5

-1

-0.5

0

x(t)

(i)

Figure 2: Numerical simulation of random coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model given by Eq. (3).

5

20

1.5

1

1

0.5

0.5

0

15

y(t)

x(t)

5

0

z(t)

10

0

-0.5

-0.5

-1

-5

-10

-15

-1 0

20

40

60

80

100

120

140

160

180

200

-1.5 0

20

40

60

80

t (days)

100

120

140

160

180

200

0

20

40

60

80

t (days)

(a)

100

120

140

160

180

200

t (days)

(b)

(c)

1

0.5

z(t)

0

-0.5

-1

-1.5 -1.5

-1

-0.5

0

0.5

1

1.5

y(t)

(d)

(e)

(f)

1.5

20

1 10

x(t)

y(t)

0.5 0

0

-0.5 -10 -1 -1.5 1

-20 1 20

0

2

0

10 0

-1

1 0

-1

-10

z(t)

-2

-20

-1

z(t)

x(t)

(g)

(h)

-2

-2

y(t)

(i)

Figure 3: Numerical simulation of classical continuous-time vallis model for El Ni˜ no given by Eq. (5) for β = 102, C = 3 and p = 0.5.

6

For the independent random variables χi , i ∈ (1, 3), the distribution is N (0, 1). Setting the appropriate values of (ai , si ), yields A = 102 + 0.65χ1 ; B = 3 + 0.45χ2 ; C = 0.5 + 0.03χ3 . (6) Substituting the parameters described in (6) into the model (5), we get dx(t) = Ax(t) + Bj y(t) − Ej x3 (t), , dt dy(t) = −Cj x(t) + Dj y(t). (7) dt The solution of the model (7) can be obtained applying the Adams-Bashforth method [37]. The Numerical scheme is given by o 1 n o 3 n xn+1 = x0 + h Axn (tn ) + Bj yn (tn ) − Ej x3n (tn ) − h Axn−1 (tn−1 ) + Bj yn−1 (tn−1 ) − Ej x3n−1 (tn−1 ) , 2 2 n o o 1 n 3 (8) yn+1 = y0 + h − Cj xn (tn ) + Dj yn (tn ) − h − Cj xn−1 (tn−1 ) + Dj yn−1 (tn−1 ) . 2 2 Now, we consider the following dimensionless parameters given by Eq. (6), the time simulation t = 200 [days], step size h = 1 × 10−3 , and initial conditions: x(0) = 1, y(0) = 1 and z(0) = 1. The system (7) was solved numerically using the the Adams-Bashforth method (8). The numerical results given in Figs. 4a-4i shows the dynamics and phase portraits for the random continuous-time vallis model for El Ni˜ no.

2.2

Coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model.

Liouville-Caputo fractional-order derivative. The fractional-order coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model corresponding to the system (1) is as follows C α 3 0 Dt x(t) = Cx(t) + ηy(t) − ξx (t), C α 0 Dt y(t)

= −ςx(t) − γy(t).

(9)

The model is subject to initial conditions x(0) (t) = x(0);

y(0) (t) = y(0),

(10)

where α ∈ (0, 1]. Fractional derivative used in the model (9) are all in Liouville-Caputo sense, C α 0 Dt {f (t)}

1 = Γ(n − α)

Z

t

0

dn f (σ)(t − σ)n−α−1 dσ, dtn

n − 1 < α ≤ n.

C α 0 Dt

[21] (11)

Existence and uniqueness of the solutions. Assume a convex, bounded and closed subset Ω of a Banach space Φ and ξ : Ω → Ω a condensing map, where Φ have a fixed point in Ω. We consider the initial value problem on the cylinder ∆ = {(t, p) ∈ R × Φ : t ∈ [0, T ], x ∈ Ω(0, φ)} for some + fixed T > 0, φ > 0 and suppose T 1 that ∃ ∆ ∈ (0, α), x, y L1 ∈ L1/∆ ([0, T ], R ). Also, The following functions x10 , y10 , x11 , y11 ∈⊂ (R, Φ) Lloc (R, Φ) such that x = x10 + x11 , y = y10 + x11 and the following assumptions hold 1. x10 and y10 are bounded and Lipszchitz. 2. x11 and y11 are compact and bounded. 3. |R(t, p) − R(t, z)| ≤ L1 (t)||p − z|| ∀ (t, p), (t, z) ∈ R.

7

20

1.5

15

1

1

0.5

10 0.5 0

z(t)

y(t)

x(t)

5 0

0

-0.5 -0.5

-5 -1

-1

-10

-15

-1.5 0

20

40

60

80

100

120

140

160

180

200

-1.5 0

20

40

60

80

100

t (days)

120

140

160

180

200

0

20

40

60

t (days)

(a)

1

100

120

140

160

180

200

t (days)

(b)

1.5

80

(c)

1

1

0.5

0.5

0

0

0

z(t)

z(t)

y(t)

0.5

-0.5

-0.5

-1

-1

-0.5

-1

-1.5 -15

-10

-5

0

5

10

15

-1.5 -15

20

-10

-5

0

x(t)

5

10

15

-1.5 -1.5

20

x(t)

(d)

-0.5

0

0.5

y(t)

(e)

1.5

-1

(f)

20

1 10

x(t)

y(t)

0.5 0

0

-0.5 -10 -1 -1.5 1

-20 2 20

0

10 0

-1

1

1 0.5

0

0 -1

-10

z(t)

-2

-20

x(t)

-0.5

y(t)

-2

(g)

-1 -1.5

(h)

z(t)

(i)

Figure 4: Numerical simulation of random continuous-time vallis model for El Ni˜ no given by Eq. (7).

8

1

1.5

Considering the Riemann-Liouville integral [38] to both sides of Eq. (9), we get Z t Z t 1 1 α−1 (t − σ) x10 (σ, x(σ))dσ + (t − σ)α−1 x11 (σ, x(σ))dσ, x(t) = x(0) + Γ(α1 ) 0 Γ(α) 0 Z t Z t 1 1 y(t) = y(0) + (t − σ)α−1 y10 (σ, y(σ))dσ + (t − σ)α−1 y11 (σ, y(σ))dσ. Γ(α) 0 Γ(α) 0

(12)

Theorem 1. Based on the assumptions 1 and 2 the intial value problem possesses at least one solution in the interval [0, T ] subject to the condition ω||L||1/O T M < 1, (13) ψ= Γ(α)  1−O 1−O where M = α − O and Ω = α−O . 1 ω(||H1 ||1/O + ||H2 ||1/O )T M ≤ µ and supposed Ωµ = {p : ||p|| ≤ µ}, the Proof. Considering µ such that ξ(0)| + Γ(α) closed ball in a Banach space ([0, T ], Φ) with sup || · ||. Now, we consider p : Ωµ → a Banach space ([0, T ], Φ), p → x10 p + x11 p with Z t 1 (t − σ)α−1 x10 (σ, p(σ))dσ, x10 (t) = x(0) + Γ(α) 0 Z t 1 x11 (t) = x(0) + (t − σ)α−1 x11 (σ, p(σ))dσ, Γ(α) 0 Z t 1 y10 (t) = y(0) + (t − σ)α−1 y10 (σ, p(σ))dσ, Γ(α) 0 Z t 1 (t − σ)α−1 y11 (σ, p(σ))dσ. (14) y11 (t) = y(0) + Γ(α) 0

Now, we proved that x(t) and y(t) are condensing, we can prove the existence of a fixed point of x(t) and y(t). A) We need prove that x(Ωµ ) ⊂ Ωµ . For p ∈ Ωµ , yields Z t 1 ||x(t)|| ≤ |x(0)| + (t − σ)α−1 x10 (σ, p(σ))dσ, Γ(α) 0 !1−O Z !O !1−O Z t Z t t α−1 α−1 1 1 1/O (t−σ) 1−O dσ J1 (σ)dσ + (t−σ) 1−O dσ ≤ |x(0)|+ Γ(α) Γ(α) 0 0 0 ≤ |x(0)| +

Z

t

!O 1/O J2 (σ)dσ

, (15)

0

ω1 (||J1 ||1/O + ||J2 ||1/O ) M1 T ≤ µ1 , Γ(α)

following the same procedure, we obtain ||y(t)|| ≤ |y(0)| +

ω2 (||J3 ||1/O + ||J4 ||1/O ) M2 T ≤ µ2 , Γ(α)

and therefore xi (Ωµ ) ⊂ Ωµ for i = 1, 2. B) We need prove that x and y are constractions. For p, z ∈ Ωµ , we obtain Z t 1 ||xp(t) − xz(t)|| ≤ (t − σ)α−1 L(σ)|p(σ) − z(σ)|dσ, Γ(α) 0 !1−O Z t α−1 1 1−O ≤ (t − σ) dσ (L1/O (σ)dσ)O ||p − z|| ≤ Ψi ||p − z||, Γ(α) 0 9

(16)

||yp(t) − yz(t)|| ≤ ≤

1 Γ(α)

Z

t

Z t 1 (t − σ)α−1 L(σ)|p(σ) − z(σ)|dσ, Γ(α) 0 !1−O

α−1

(t − σ) 1−O dσ

(L1/O (σ)dσ)O ||p − z|| ≤ Ψi ||p − z||,

(17)

0

where, ωi ||L||1/O T µi < 1, for i = 1, 2. (18) Γ(α) The above equation prove that x and y are contractions and satisfying, particularly ||x(p) − x11 (z)|| ≤ Ψ||p − z||. Ψi =

C) We need prove that x11 and y11 are compact. For 0 ≤ p1 ≤ p2 ≤ T , yields Z Z p1 1 p2 α−1 α1 −1 (p2 − σ) x11 (σ, p(σ))dσ − (p1 − σ) x11 (σ, p(σ))dσ , ||x11 p(p1 ) − x11 z(p2 )|| ≤ Γ(α) 0 0 i α−O α−O 1−O ωi h α−O ωi p11−O − p2 1−O + (p2 − p1 ) 1−O (p2 − p1 )α−O ||J11 ||1/O , ||J11 ||1/O + (19) Γ(α) Γ(α) i 2ωi ||J11 ||1/O α−O 1−O ωi ωi h ||J11 ||1/O + (p2 − p1 )α−O ||J11 ||1/O ≤ (p2 − p1 )α−O , (p2 − p1 ) 1−O ≤ Γ(α) Γ(α) Γ(α) following the same procedure, we obtain i 2ωi ||J21 ||1/O α−O 1−O ωi ωi h ||J21 ||1/O + ||y11 p(p1 ) − y11 z(p2 )|| ≤ (p2 − p1 ) 1−O (p2 − p1 )α−O ||J21 ||1/O ≤ (p2 − p1 )α−O , Γ(α) Γ(α) Γ(α) (20) for ωi ∈ [1, 2]. From Arzela-Ascoli principle [39], can be concluded that x11 (Ωµ ) and y11 (Ωµ ) are relatively compact which implies that x11 and y11 are compact. Since x10 and y10 are contractions and x11 and y11 are compact and hence completely continuos [40], the maps x = x10 + x11 and y = y10 + y11 are condensing on Ωµ , we get the existence of fixed points of x and y, respectively. ≤

D) We need prove that the given initial value problem possesses the solution on the real interval [0, T ]. For this prove, we consider the assumption 3, the condition (18) and the map ϑ given by Z t 1 (t − σ)α−1 x(σ, x(σ))dσ, ϑ[x(t)] = x(0) + Γ(α) 0 Z t 1 (t − σ)α−1 y(σ, y(σ))dσ. (21) ϑ[y(t)] = y(0) + Γ(α) 0 For x10 (t), x11 (t), y10 (t), y11 (t) ∈ Ωµ , we obtain Z t 1 |ϑ[x10 (t)] − ϑ[x11 (t)]| ≤ (t − σ)α−1 L1 (σ)|x10 (σ) − x11 (σ)|dσ, Γ(α) 0 #O "Z #1−O " Z t t α−1 1 1/O L1 (σ)dσ , (22) ≤ (t − σ) 1−O dσ Γ(α) 0 0 ≤

ω1 ||L1 ||1/O T M1 ||x10 − x11 ||, Γ(α)

following the same procedure, we obtain |ϑ[y10 (t)] − ϑ[y11 (t)]| ≤

ω2 ||L2 ||1/O T M2 ||y10 − y11 ||. Γ(α)

(23)

For the above cases, the condition (18) is guaranteed and the existence of the unique solution is verified for the coupled dynamical El Ni˜ no/La Ni˜ na-southern oscillation model involving the Liouville-Caputo fractional-order derivative. 10

Predictor-corrector Adams-Bashforth-Moulton method. The Adams-Bashforth-Moulton (ABM) method is a linear multi-step integration method. fractional-order problem of the form C γ 0 Dt f (t)

f k (0) = f0k ,

= Υ(t, f (t)),

We consider a

k = 0, 1, ..., n − 1.

(24)

Considering that Eq. (24) has a unique solution defined on the interval t ∈ [0, T ], Eq. (24) satisfying the following Volterra integral equation f (t) =

n−1 X

k

1 + k! Γ(γ)

(k) t f0

k=0

Zt

(t − σ)γ−1 Ξ(σ, f (σ))dσ,

t < T,

(25)

0

γ where γ > 0 and C 0 Dt is the Liouville-Caputo fractional derivative given by Eq. (11). The solution scheme used to solve Eq. (24) is given by [37]

P fk+1

fk+1

j n−1 P tk+1

k 1 P bj,k+1 Ξ(tj , fj ), Γ(γ) j=0 j=0 j! ! j n−1 k P tk+1 (j) P 1 P = f0 + aj,k+1 Ξ(tj , fj ) + ak+1,k+1 Ξ(tk+1 , fk+1 ) , Γ(γ) j=0 j=0 j!

=

(j)

f0 +

(26)

where, aj,k+1

 γ+1 − (k − γ)(k + 1)γ ) j = 0,  (k hγ ((k − j + 2)γ+1 + (k − j)γ+1 − 2(k − j + 1)γ+1 ) 1 ≤ j ≤ k, · = γ(γ + 1)  1 j = k + 1, hγ ((k + 1 − j)γ − (k − j)γ ), γ

bj,k+1 =

j = 0, 1, 2, ..., k.

(27)

Following this procedure, we can propose a solution for system given by Eq. (9) using the Adams-BashforthMoulton method given by Eq. (26) we get

x(t) =

n−1 X

k (k) t

x(0)

k=0

y(t) =

1 + k! Γ(α)

Zt

h i (t − σ)α−1 Cx(σ) + ηy(σ) − ξx3 (σ) dσ,

0 n−1 X k=0

y(0)(k)

1 tk + k! Γ(α)

Zt

h i (t − σ)α−1 − ςx(σ) − γy(σ) dσ.

(28)

0

We shown the chaos mechanism in the coupled dynamical El Ni˜ no/La Ni˜ na-southern oscillation model. We consider the following dimensionless parameters set: C = 1, η = 1; 1.5; 2, ς = 1; 2; 3, γ = 1; 0.85; 0.70 and ξ = 0.1; 0.3; 0.6. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1 and y(0) = 1. The numerical results given in Figs. 5a-5i shows the phase portraits for the fractional-order coupled dynamical El Ni˜ no/La Ni˜ na-southern oscillation model in the Liouville-Caputo sense applying the solution scheme (28) for several values of α ∈ (0, 1], arbitrarily chosen. Atangana-Baleanu-Caputo fractional-order derivative. The fractional-order coupled dynamical El Ni˜ no/La Ni˜ na-southern oscillation model corresponding to the system (1) is as follows ABC α Dt x(t) = Cx(t) + ηy(t) − ξx3 (t), 0 ABC α Dt y(t) 0

= −ςx(t) − γy(t). 11

(29)

2.5

1.5 =1 =0.95 =0.90 =0.85

2

1.5 =1 =0.95 =0.90 =0.85

1

1.5

=1 =0.95 =0.90 =0.85

1

0.5

0.5

0

0

0.5

y(t)

y(t)

x(t)

1

-0.5

-0.5

-1

-1

-1.5

-1.5

0 -0.5 -1 -1.5

-2 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

100

120

140

160

180

-2 -1.5

200

(a)

1

1

0.5

0.5

0.5

0

0

y(t)

x(t) 0

-0.5

-0.5

-1

-1

100

120

140

160

180

2

2.5

=1 =0.95 =0.90 =0.85

-1 =1 =0.95 =0.90 =0.85

-1.5

-2 80

1.5

-0.5

-1.5

-1.5

1

y(t)

1

1.5

60

0.5

(c)

1.5

40

0

x(t)

1.5 =1 =0.95 =0.90 =0.85

20

-0.5

(b)

2

0

-1

t (days)

200

0

20

40

60

80

t (days)

100

120

140

160

180

-2 -1.5

200

-1

-0.5

t (days)

(d)

0

0.5

1

1.5

2

x(t)

(e)

(f)

1.5 =1 =0.95 =0.90 =0.85

1.5

1

=1 =0.95 =0.90 =0.85

1.5

1

=1 =0.95 =0.90 =0.85

1 0.5

0

y(t)

0.5

y(t)

x(t)

0.5

0

0 -0.5

-0.5

-0.5 -1

-1

-1

-1.5 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

t (days)

80

100

t (days)

(g)

(h)

120

140

160

180

200

-1.5 -1.5

-1

-0.5

0

0.5

1

x(t)

(i)

Figure 5: Numerical simulation for the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model given by Eq. (9). In (a)-(c) C = 1, η = 1, ς = 1, γ = 1, and ξ = 0.1. In (d)-(f) C = 1, η = 1.5, ς = 2, γ = 0.85, and ξ = 0.3. In (g)-(i) C = 1, η = 2, ς = 3, γ = 0.70, and ξ = 0.6, we consider different values for α ∈ (0, 1], arbitrarily chosen.

12

1.5

The model is subject to initial conditions x(0) (t) = x(0);

y(0) (t) = y(0),

(30)

where αi ∈ (0, 1] for i = 1, 2 are the orders of the fractional derivatives. Fractional derivative used in the model given by Eq. (29) are all in Atangana-Baleanu-Caputo sense, [23] Z h (t − σ)α i B(α) t dn ABC α f (σ)Eα − α dσ, n − 1 < α ≤ n, D {f (t)} = 0 t n n − α 0 dt n−α

ABC α Dt 0

(31)

where B(α) is a normalization function, B(0) = B(1) = 1. This fractional operator uses the Mittag-Leffler law as nonsingular and nonlocal kernel. Existence and uniqueness of the solution. The Atangana-Baleanu fractional integral of order α of a function f (t) is defined as [23] AB α 0 It {f (t)}

= f (t) − f (0) =

1−α α f (t) + B(α) B(α)Γ(α)

Z

t

f (σ)(t − σ)α−1 dσ,

(32)

0

where B(α) is a normalization function as in the case commented above. Consider Eq. (32), the system (29) is equivalent to the following Z t n o o 1 − αn α 3 x(t) − g1 (t) = (t − σ)α−1 Cx(σ) + ηy(σ) − ξx3 (σ) dσ, Cx(t) + ηy(t) − ξx (t) + B(α) B(α)Γ(α) 0 y(t) − g2 (t) =

Z t o n o 1 − αn α − ςx(t) − γy(t) + (t − σ)α−1 − ςx(σ) − γy(σ) dσ. B(α) B(α)Γ(α) 0

(33)

The above system is convert to iterative routine considering x(0) (t) = g1 (t);

x(0) (t) = g2 (t),

(34)

obtaining the following iterative formula x(n+1) (t) =

o α 1 − αn Cxn (t) + ηyn (t) − ξx3n (t) + B(α) B(α)Γ(α)

y(n+1) (t) =

1 − αn B(α)

o − ςxn (t) − γyn (t) +

α B(α)Γ(α)

t

Z

n o (t − σ)α−1 Cxn (σ) + ηyn (σ) − ξx3n (σ) dσ,

0

Z

t

n o (t − σ)α−1 − ςxn (σ) − γyn (σ) dσ.

(35)

0

Taking the limit for a large value of n we expect to obtain the exact solution. Theorem 2. The initial value problem given by Eq. (29) possesses at least one solution in the interval [0, T ]. Proof. Considering the Picard-Lindelof approach, we consider the following operator Ξ1 (t, x) = Cx(t) + ηy(t) − ξx3 (t), Ξ2 (t, y) = −ςx(t) − γy(t), where Ξ1 (t, x) and Ξ2 (t, y) are contractions respect to x and y, respectively. Let N1 = sup||Ψa,b1 Ξ1 (t, x)||; N2 = sup||Ψa,b2 Ξ2 (t, y)||,

(36)

(37)

where, Ψa,b1 = |t − a, t + a| × [x − b1 , x + b1 ] = A1 × B1 ,

Ψa,b2 = |t − a, t + a| × [y − b2 , y + b2 ] = A1 × B2 .

13

(38)

Considering the Picard’s operator, we have Θ : Ψ(A1 , B1 , B2 , B3 ) → Ψ(A1 , B1 , B2 , B3 ),

(39)

defined as follows ΘΩ(t) = Ω0 (t) + ∆(t, Ω(t))

1−α α + B(α) B(α)Γ(α)

Z

t

(t − σ)α−1 ∆(σ, Ω(σ))dσ,

(40)

0

where Ω(t) = {x(t), y(t)}, Ω0 (t) = {g1 , g2 } and ∆(t, Ω(t)) = {Ξ1 (t, x(t)), Ξ2 (t, y(t))}. Now we assume that all the solutions are bounded within a period of time ||Ω(t)||∞ ≤ max{B1 , B2 }, Z t 1−α α (t − σ)α−1 ∆(σ, Ω(σ))dσ ||Ω(t) − Ω0 (t)|| = ∆(t, Ω(t)) + B(α) B(α)Γ(α) 0 Z t 1−α α ≤ (t − σ)α−1 ||∆(σ, Ω(σ))||dσ ||∆(t, Ω(t))|| + B(α) B(α)Γ(α) 0 ≤

(41)

1−α α Z = max{B1 , B2 } + Z$α ≤ $Z ≤ B = max{B1 , B2 }. B(α) B(α)

Here we request that B . Z Using the fixed point theorem of Banach space together with the metric, we have $<

||ΘΩ1 − ΘΩ2 ||∞ = sup||t∈A |Ω1 − Ω2 |,

(42)

(43)

Z t 1−α α ||ΘΩ1 − ΘΩ2 || = {∆(t, Ω1 (t)) − ∆(t, Ω2 (t))} + (t − σ)α−1 {∆(σ, Ω1 (t)) − ∆(σ, Ω2 (t))}dσ , B(α) B(α)Γ(α) 0 Z t α% 1−α q||Ω1 (t) − Ω2 (t)|| + (t − σ)α−1 ||Ω1 (t)) − Ω2 (t))||dσ, (44) ≤ B(α) B(α)Γ(α) 0 ≤ $%||Ω1 (t)) − Ω2 (t))||, with % < 1. Since Ω is a contraction we have that $% < 1, thus the defined operator Θ is a contraction too. We conclude that the system (29) as a unique set of solution. Adams-Moulton method. The numerical approximation of Atangana-Baleanu fractional integral (32) using the Adams-Moulton rule is given by [41] ∞ α X h f (tν+1 ) − f (tν ) i α 1 − α h f (tn+1 ) − f (tn ) i AB α + ξν , (45) 0 It [f (tn+1 )] = B(α) 2 Γ(α) ν=0 2 where ξνα = (ν + 1)1−α − (ν)1−α . Using the above numerical scheme, we have ( ! ! !3 )) x (t) − x (t) y (t) − y (t) x (t) − x (t) 1 − α (n+1) (n) (n+1) (n) (n+1) (n) x(n+1) (t)−x(n) (t) = x(0)n (t)+ C +η −ξ + B(α) 2 2 2 (

! ! !3 ) ( ∞ x(ν+1) (t) − x(ν) (t) y(ν+1) (t) − y(ν) (t) x(ν+1) (t) − x(ν) (t) α X α + ξ C +η −ξ , B(α) ν=0 ν 2 2 2

14

(

! !)) x (t) − x (t) y (t) − y (t) (n+1) (n) (n+1) (n) y(n+1) (t) − y(n) (t) = y(0)n (t) + −ς −γ + 2 2 ( ! !) ∞ x(ν+1) (t) − x(ν) (t) y(ν+1) (t) − y(ν) (t) α X α ξ −ς −γ . + B(α) ν=0 ν 2 2 1−α B(α)

(

(46)

We shown the chaos mechanism in the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model. We consider the following dimensionless parameters set: C = 1, η = 1; 1.5; 2, ς = 1; 2; 3, γ = 1; 0.85; 0.70 and ξ = 0.1; 0.3; 0.6. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1 and y(0) = 1. The numerical results given in Figs. 6a-6i shows the phase portraits for the fractional-order coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model in the Atangana-Baleanu-Caputo sense applying the solution scheme (46) for several values of α ∈ (0, 1], arbitrarily chosen. Fractal-Fractional power-law derivative. The fractional-fractional coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model involving the power-law derivative corresponding to the system (1) is given by F F −P α,τ 0 Dt x(t)

= Ξ1 (x(t), y(t), t),

F F −P α,τ 0 Dt y(t)

= Ξ2 (x(t), y(t), t),

(47)

3

where Ξ1 (x(t), y(t), t) = Cx(t) + ηy(t) − ξx (t), Ξ2 (x(t), y(t), t) = −ςx(t) − γy(t) and α, τ ∈ (0, 1]. The model is subject to initial conditions x(0) (t) = x(0);

y(0) (t) = y(0).

The power-law fractal-fractional derivative is given by [25] Z t 1 d F F −P α,τ D {f (t)} = f (σ)(t − σ)n−α−1 dσ, 0 t Γ(n − α) dtτ 0

(48)

n − 1 < α, τ ≤ n ∈ N,

(49)

where, d f (t) − f (σ) f (σ) = lim . (50) τ t→σ dt tτ − σ τ In order to consider integer-order initial conditions, we replace the power-law Riemann-Liouville derivative to power-law Liouville-Caputo derivative. Eq. (47) can be converted to the Volterra integrals and the numerical scheme at tn+1 is given by n

xn+1 (t) = x( 0) +

τ X Γ(α) j=0 n

τ X yn+1 (t) = y(0) + Γ(α) j=0

Z

tj+1

σ τ −1 (tn+1 − σ)α−1 Ξ1 (x, y, σ)dσ,

tj

Z

tj+1

σ τ −1 (tn+1 − σ)α−1 Ξ2 (x, y, σ)dσ.

(51)

tj

Within the finite interval [tj , tj+1 ], we approximate the functions σ τ −1 Ξi (x, y, σ) ∀ i = 1, 2 using the Lagrangian piece-wise interpolation such that Θj (σ) =

σ − tj−1 τ −1 σ − tj τ −1 tj Ξi (xj , y j , tj ) − t Ξi (xj−1 , y j−1 , tj−1 ). tj − tj−1 tj − tj−1 j−1

Thus, we obtain n

τ X xn+1 (t) = x(0) + Γ(α) j=0

Z

tj+1

tj

15

σ τ −1 (tn+1 − σ)α−1 Θj (σ)dσ,

(52)

2.5

1.5 =1 =0.95 =0.90 =0.85

2

1.5 =1 =0.95 =0.90 =0.85

1

1.5

=1 =0.95 =0.90 =0.85

1

0.5

0.5

0

0

0.5

y(t)

y(t)

x(t)

1

-0.5

-0.5

-1

-1

-1.5

-1.5

0 -0.5 -1 -1.5

-2 0

20

40

60

80

100

120

140

160

180

200

0

20

40

60

80

t (days)

(a)

120

140

160

180

-2 -1.5

200

1

1

0.5

0.5

0.5

0

0

y(t)

x(t) 0

-0.5

-0.5

-1

-1

-1.5 100

120

140

160

180

200

0

20

40

60

80

t (days)

(d)

100

120

140

160

180

-2 -1.5

200

0

-0.5

-1

-1

-1

-1.5 60

80

100

120

140

160

180

200

1.5

2

0

-0.5

40

1

0.5

-0.5

-1.5

0.5

1

y(t)

y(t)

0

0

1.5 =1 =0.95 =0.90 =0.85

0.5

20

-0.5

(f)

1

0.5

0

-1

x(t)

1.5

1

2.5

-1.5

(e)

=1 =0.95 =0.90 =0.85

2

=1 =0.95 =0.90 =0.85

t (days)

1.5

1.5

-1 =1 =0.95 =0.90 =0.85

-2 80

1

-0.5

-1.5

60

0.5

y(t)

1

1.5

40

0

(c)

1.5

20

-0.5

x(t)

1.5 =1 =0.95 =0.90 =0.85

0

-1

(b)

2

x(t)

100

t (days)

0

20

40

60

t (days)

80

100

t (days)

(g)

(h)

120

140

160

180

200

-1.5 -1.5

=1 =0.95 =0.90 =0.85

-1

-0.5

0

0.5

1

x(t)

(i)

Figure 6: Numerical simulation for the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model given by Eq. (29). In (a)-(c) C = 1, η = 1, ς = 1, γ = 1, and ξ = 0.1. In (d)-(f) C = 1, η = 1.5, ς = 2, γ = 0.85, and ξ = 0.3. In (g)-(i) C = 1, η = 2, ς = 3, γ = 0.70, and ξ = 0.6, we consider different values for α ∈ (0, 1], arbitrarily chosen.

16

1.5

n

τ X yn+1 (t) = y(0) + Γ(α) j=0

Z

tj+1

σ τ −1 (tn+1 − σ)α−1 Θj (σ)dσ,

(53)

tj

Solving the integrals of the right hand side, we obtain the following numerical scheme n   τ (∆t)α X h τ −1 tj Ξ1 (xj , y j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − Γ(α + 2) j=0  i −1 tτj−1 Ξ1 ((xj−1 , y j−1 , tj−1 ) (n + 1 − j)α+1 − (n − j)α (n − j + 1 + α) ,

xn+1 (t) = x(0) +

n   τ (∆t)α X h τ −1 tj Ξ2 (xj , y j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − Γ(α + 2) j=0 (54)  i τ −1 j−1 j−1 α+1 α tj−1 Ξ2 (x , y , tj−1 ) (n + 1 − j) − (n − j) (n − j + 1 + α) .

yn+1 (t) = y(0) +

We shown the chaos mechanism in the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model. We consider the following dimensionless parameters set: C = 1, η = 1; 1.5; 2, ς = 1; 2; 3, γ = 1; 0.85; 0.70 and ξ = 0.1; 0.3; 0.6. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1 and y(0) = 1. The numerical results given in Figs. 7a-7i shows the phase portraits for the fractional-order coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model involving the fractal-fractional Liouville-Caputo derivative considering the solution scheme (54) for several values of α − τ ∈ (0, 1], arbitrarily chosen. Fractal-Fractional Mittag-Leffler derivative. The fractional-fractional coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model involving the MittagLeffler memory corresponding to the system (1) is given by F F −M α,τ 0 Dt x(t)

= Ξ1 (x(t), y(t), t),

F F −M α,τ 0 Dt y(t)

= Ξ2 (x(t), y(t), t),

(55)

3

where Ξ1 (x(t), y(t), t) = Cx(t) + ηy(t) − ξx (t), Ξ2 (x(t), y(t), t) = −ςx(t) − γy(t) and α, τ ∈ (0, 1]. The model is subject to initial conditions x(0) (t) = x(0);

y(0) (t) = y(0).

(56)

The Atangana-Baleanu-Caputo fractal-fractional derivative is given by [25] F F −M

where

d dtγ f (σ)

α,γ 0 Dt {f (t)} =

B(α) d n − α dtγ

Z 0

t

h (t − σ)α i dσ, f (σ)Eα − α n−α

n − 1 < α, γ ≤ n ∈ N,

(57)

is given by Eq. (50).

Applying the Atangana-Baleanu integral, we have τ tτ −1 (1 − α) ατ Ξ1 (x, y, t) + AB(α) AB(α)Γ(α)

Z

x(t) = x(0) +

y(t) = y(0) +

ατ τ tτ −1 (1 − α) Ξ2 (x, y, t) + AB(α) AB(α)Γ(α)

Z

t

Ψτ −1 (t − σ)α−1 Ξ1 (x, y, σ)dσ,

0 t

Ψτ −1 (t − σ)α−1 Ξ2 (x, y, σ)dσ,

0

at tn+1 , we have n

X τ tτ −1 (1 − α) ατ xn+1 (t) = x(0) + n Ξ1 (xn , y n , tn ) + AB(α) AB(α)Γ(α) j=0 17

Z

tj+1

tj

σ τ −1 (tn+1 − σ)α−1 Ξ1 (x, y, σ)dσ,

(58)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 7: Numerical simulation for the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model given by Eq. (47). In (a)-(c) C = 1, η = 1, ς = 1, γ = 1, and ξ = 0.1. In (d)-(f) C = 1, η = 1.5, ς = 2, γ = 0.85, and ξ = 0.3. In (g)-(i) C = 1, η = 2, ς = 3, γ = 0.70, and ξ = 0.6, we consider different values for α − τ ∈ (0, 1], arbitrarily chosen.

18

n

X τ tτ −1 (1 − α) ατ yn+1 (t) = y(0) + n Ξ2 (xn , y n , tn ) + AB(α) AB(α)Γ(α) j=0

Z

tj+1

σ τ −1 (tn+1 − σ)α−1 Ξ2 (x, y, σ)dσ,

(59)

tj

Approximating σ τ −1 Ξi (x, y, σ) ∀ i = 1, 2 in [tj , tj+1 ], we have the following numerical scheme τ tτ −1 (1 − α) τ (∆t)α xn+1 (t) = x(0) + n Ξ1 (xn , y n , tn ) + × AB(α) AB(α)Γ(α + 2) " n   X tτj −1 Ξ1 (xj , y j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − j=0



−1 tτj−1 Ξ1 (xj−1 , y j−1 , tj−1 )

(n − j + 1)

α+1

#  − (n − j) (n − j + 1 + α) , α

τ tτ −1 (1 − α) τ (∆t)α yn+1 (t) = y(0) + n Ξ2 (xn , y n , tn ) + × AB(α) AB(α)Γ(α + 2) " n   X tτj −1 Ξ2 (xj , y j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) −

(60)

j=0

#

  −1 tτj−1 Ξ2 (xj−1 , y j−1 , tj−1 ) (n − j + 1)α+1 − (n − j)α (n − j + 1 + α) . We shown the chaos mechanism in the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model. We consider the following dimensionless parameters set: C = 1, η = 1; 1.5; 2, ς = 1; 2; 3, γ = 1; 0.85; 0.70 and ξ = 0.1; 0.3; 0.6. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1 and y(0) = 1. The numerical results given in Figs. 8a-8i shows the phase portraits for the fractional-order coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model involving the fractal-fractional Atangana-Baleanu-Caputo derivative considering the solution scheme (60) for several values of α − τ ∈ (0, 1], arbitrarily chosen.

2.3

Continuous-time vallis model for El Ni˜ no phenomenon.

In this section, we consider the vallis model for El Ni˜ no phenomenon, this model is a modified form of the former complex model of Saltzman to portray buoyancy-driven convection patterns in the classical rectangular RayleighBernard problem. This model can also be obtained by a simple transformation of the vallis model for El Ni˜ no [35]. Now, we consider this model invoving the Liouville-Caputo and Atangana-Baleanu derivatives, also, we consider fractal-fractional operators with power-law and Mittag-Leffler memories. Liouville-Caputo fractional-order derivative. Considering the Liouville-Caputo fractional-order derivative given by Eq. (11), the fractional-order vallis model for El Ni˜ no phenomenon corresponding to the system (5) is given by C α 0 Dt x(t)

= βy(t) − C(x(t) + p),

C α 0 Dt y(t) = x(t)z(t) − y(t), C α 0 Dt z(t) = −x(t)y(t) − z(t) +

(61) 1.

where α ∈ (0, 1]. The model is subject to initial conditions x(0) (t) = x(0);

y(0) (t) = y(0);

z(0) (t) = z(0).

(62)

The Laplace transform to Liouville-Caputo fractional order derivative (11) is given by [21] L

n

m−1 o X (s) = sα F (s) − sα−k−1 f (k) (0).

C α 0 Dt f (t)

k=0

19

(63)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 8: Numerical simulation for the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model given by Eq. (55). In (a)-(c) C = 1, η = 1, ς = 1, γ = 1, and ξ = 0.1. In (d)-(f) C = 1, η = 1.5, ς = 2, γ = 0.85, and ξ = 0.3. In (g)-(i) C = 1, η = 2, ς = 3, γ = 0.70, and ξ = 0.6, we consider different values for α − τ ∈ (0, 1], arbitrarily chosen.

20

Using the Laplace transform given by Eq. (78) and the inverse Laplace transform on both sides of Eq. (61), have the following solutions n1 h i o x(t) = x(0) + L −1 α L βy(t) − C(x(t) + p) (s) (t), s i o n1 h y(t) = y(0) + L −1 α L x(t)z(t) − y(t) (s) (t), s n1 h i o z(t) = z(0) + L −1 α L − x(t)y(t) − z(t) + 1 (s) (t). (64) s The following iterative formula is then proposed n1 h i o x(n) (t) = x(0) + L −1 α L βy(n−1) (t) − C(x(n−1) (t) + p) (s) (t), s n1 h i o y(n) (t) = y(0) + L −1 α L x(n−1) (t)z(n−1) (t) − y(n−1) (t) (s) (t), (65) s n1 h i o z(n) (t) = z(0) + L −1 α L − x(n−1) (t)y(n−1) (t) − z(n−1) (t) + 1 (s) (t), s where, x(0) (t) = x(0); y(0)(t) = y(0); z(0)(t) = z(0). (66) The approximate solution is assumed to be obtain as a limit when n tend to infinity x(t) = lim x(n) (t);

y(t) = lim y(n) (t);

n→∞

n→∞

z(t) = lim z(n) (t). n→∞

(67)

Theorem 3. The recursive method given by Eq. (65) is stable in the interval [0, T ]. Proof. We assume the following, it is possible to find three positive constants u, v and w such that for all 0 ≤ t ≤ T ≤ ∞,

||x(t)|| < u;

||y(t)|| < v;

||z(t)|| < w.

Now, we consider a subset of C2 ((a, b)(0, T )) defined by Z n o 1 R = ϕ : (a, b)(0, T ) → R, (t − σ)α−1 U (σ)S(σ)dσ < ∞ , Γ(α)

(68)

(69)

we now consider the following operator φ defined as   βy(t) − C(x(t) + p), x(t)z(t) − y(t), φ(x, y, z) =  −x(t)y(t) − z(t) + 1. Then

 < φ(x, y, z) − φ(X, Y, Z), (x − X, y − Y, z − Z) >,    β(y(t) − Y (t)) − C[(x(t) − X(t)) + p], = (x(t) − X(t))(z(t) − Z(t)) − (y(t) − Y (t)),    −(x(t) − X(t))(y(t) − Y (t)) − (z(t) − Z(t)) + 1 >,

where, x(t) 6= X(t);

y(t) 6= Y (t);

z(t) 6= Z(t).

(70)

Thus, applying the norm and the absolute value on both sides, we have

<

        

< φ(x, y, z) − φ(X, Y, Z), (x − X, y − Y, z − Z) >, n h io ||y(t)−Y (t)|| ||x(t)−X(t)|| ||x(t) − X(t)||2 , β ||x(t)−X(t)|| 2 − C ||x(t)−X(t)||2 + p n    o ||x(t)−X(t)|| ||z(t)−Z(t)|| ||y(t)−Y (t)|| − ||y(t) − Y (t)||2 , 2 2 2 (t)|| ||y(t)−Y (t)|| (t)|| n ||y(t)−Y    ||y(t)−Y   o ||x(t)−X(t)|| ||y(t)−Y (t)|| ||z(t)−Z(t)|| 2 − ||z(t)−Z(t)|| 2 ||z(t)−Z(t)||2 − ||z(t)−Z(t)||2 + 1 ||z(t) − Z(t)|| . 21

(71)

where, < φ(x, y, z) − φ(X, Y, Z), (x − X, y − Y, z − Z) >,   M ||x(t) − X(t)||2 , N ||y(t) − Y (t)||2 , <  O||z(t) − Z(t)||2 , with

(72)

h ||x(t) − X(t)|| io n ||y(t) − Y (t)|| −C +p , M= β 2 2 ||x(t) − X(t)|| ||x(t) − X(t)|| n ||x(t) − X(t)||  ||z(t) − Z(t)||   ||y(t) − Y (t)|| o N= − , ||y(t) − Y (t)||2 ||y(t) − Y (t)||2 ||y(t) − Y (t)||2 n  ||x(t) − X(t)||  ||y(t) − Y (t)||   ||z(t) − Z(t)||  o O= − − + 1 . ||z(t) − Z(t)||2 ||z(t) − Z(t)||2 ||z(t) − Z(t)||2

(73)

Also if we consider a given non-null vector (x, y, z), the using the some routine as above case, we obtain < φ(x, y, z) − φ(X, Y, Z), (x − X, y − Y, z − Z) >,   M ||x(t) − X(t)||||x(t)||, N ||y(t) − Y (t)||||y(t)||, <  O||z(t) − Z(t)||||z(t)||,

(74)

from the results obtained in Eqs. (72) and (74), we conclude that the used iterative method is stable. Predictor-corrector Adams-Bashforth-Moulton method. Considering the numerical scheme given by Eq. (26), we can propose a solution for system (61) in the following form Zt n−1 h i k X 1 (k) t x(t) = x(0) + (t − σ)α−1 βy(σ) − C(x(σ) + p) dσ, k! Γ(α) k=0

y(t) =

n−1 X

0 (k) t

y(0)

k=0

z(t) =

n−1 X k=0

k

1 + k! Γ(α)

Zt

h i (t − σ)α−1 x(σ)z(σ) − y(σ) dσ,

(75)

0

z(0)(k)

1 tk + k! Γ(α)

Zt

h i (t − σ)α−1 − x(σ)y(σ) − z(σ) + 1 dσ.

0

We shown the chaos mechanism in the vallis model for El Ni˜ no phenomenon. We consider the following dimensionless parameters set: β = 102, C = 3 and p = 0.5. t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1, y(0) = 1 and z(0) = given in Figs. 9a-9i shows the phase portraits for the fractional-order vallis model for El Liouville-Caputo sense applying the solution scheme (75) for several values of α ∈ (0, 1],

The time simulation was 1. The numerical results Ni˜ no phenomenon in the arbitrarily chosen.

Atangana-Baleanu-Caputo fractional-order derivative. Considering the Atangana-Baleanu-Caputo fractional-order derivative given by Eq. (31), the fractional-order vallis model for El Ni˜ no phenomenon corresponding to the system (5) is given by ABC α Dt x(t) 0

= βy(t) − C(x(t) + p),

ABC α Dt y(t) 0 ABC α Dt z(t) 0

= x(t)z(t) − y(t),

= −x(t)y(t) − z(t) + 1.

22

(76)

20

1.5 =1 =0.95 =0.90 =0.85

15

1 =1 =0.95 =0.90 =0.85

1

0.5

10 0.5 0

z(t)

y(t)

x(t)

5 0

0

-0.5 -0.5

-5

-10

-15

-1.5 0

20

40

60

80

100

120

140

160

180

200

-1.5 0

20

40

60

80

100

t (days)

120

140

160

180

200

0

20

40

60

80

t (days)

(a)

120

140

160

180

200

(c)

1 =1 =0.95 =0.90 =0.85

100

t (days)

(b)

1.5

1

=1 =0.95 =0.90 =0.85

-1

-1

1 =1 =0.95 =0.90 =0.85

0.5

=1 =0.95 =0.90 =0.85

0.5

0.5

z(t)

0

z(t)

y(t)

0 0

-0.5

-0.5

-1

-1

-0.5

-1

-1.5 -15

-10

-5

0

5

10

15

-1.5 -15

20

-10

-5

0

x(t)

(d)

10

15

-1.5 -1.5

20

0

1

1.5

20

10 =1 =0.95 =0.90 =0.85

0

-0.5 -10

x(t)

0

0.5

(f)

10

x(t)

0.5

-0.5

y(t)

20 =1 =0.95 =0.90 =0.85

1

-1

(e)

1.5

y(t)

5

x(t)

=1 =0.95 =0.90 =0.85

0

-10

-1 -1.5 1

-20 2 20

0

10 0

-1

-20 1 1

1 0.5

0

0 -1

-10

z(t)

-2

-20

x(t)

-0.5

y(t)

-2

(g)

-1 -1.5

(h)

z(t)

2

0

1 0

-1 -1

z(t)

-2

-2

(i)

Figure 9: Numerical simulation for the vallis model for El Ni˜ no phenomenon given by Eq. (61).

23

y(t)

where α ∈ (0, 1]. The model is subject to initial conditions x(0) (t) = x(0);

y(0) (t) = y(0);

z(0) (t) = z(0).

(77)

The Laplace transform to Liouville-Caputo fractional order derivative (11) is given by [21] L

m−1 o X α (s) = s F (s) − sα−k−1 f (k) (0).

n

C α 0 Dt f (t)

(78)

k=0

Since this model is nonlinear and it may be hard to obtain the exact solution. The aim of this section is to provide a special solution using an iterative scheme. The technique involves coupling the Sumudu transform and its inverse. The Sumudu transform is derived from the classical Fourier integral. This transform is itself linear, preserves linear functions, and does not change units [42]-[43]. The Sumudu transform (ST ) of the Atangana-Baleanu-Caputo fractional derivative is defined as ST

n

o

ABC α Dt f (t) 0

=

  1 B(α)  uα × [ST (f (t)) − f (0)]. αΓ(α + 1)Eα − 1−α 1−α

(79)

Using the Sumudu transform given by Eq. (79) and the inverse Sumudu transform on both sides of Eq. (76), we have the following solutions ( ! ) h i 1−α −1   · ST βy(t) − C(x(t) + p) (s) (t), x(t) = x(0) + ST 1 B(α)αΓ(α + 1)Eα − 1−α uα ( y(t) = y(0) + ST

 B(α)αΓ(α + 1)Eα −

( z(t) = z(0) + ST

1 α 1−α u

1 α 1−α u

 B(α)αΓ(α + 1)Eα −

(

1 α 1−α u

!

1−α

−1

 B(α)αΓ(α + 1)Eα −

1 α 1−α u







) h i · ST x(t)z(t) − y(t) (s) (t),

) i · ST − x(t)y(t) − z(t) + 1 (s) (t). h

!

) h i · ST βy(n−1) (t) − C(x(n−1) (t) + p) (s) (t),

!

) i · ST x(n−1) (t)z(n−1) (t) − y(n−1) (t) (s) (t),



1−α

−1



!

 B(α)αΓ(α + 1)Eα −

(

z(n) (t) = z(0) + ST

1 α 1−α u

1−α

−1

The following iterative formula is then proposed ( 1−α −1  x(n) (t) = x(0) + ST B(α)αΓ(α + 1)Eα − y(n) (t) = y(0) + ST

!

1−α

−1

(80)

h

) i · ST − x(n−1) (t)y(n−1) (t) − z(n−1) (t) + 1 (s) (t). (81) h

The approximate solution is obtained in the limit when n tends to infinity x(t) = lim x(n) (t); n→∞

y(t) = lim y(n) (t); n→∞

z(t) = lim z(n) (t). n→∞

Let (Ξ, || · ||) be a Banach space and P a self-map of Ξ. Let $n+1 = g(P, $n ) be particular recursive procedure. Suppose that, F (P ) the fixed-point of P has at least one element and that $n converges to a point P ∈ F (P ). Let {υn } ⊆ Ξ and define en = ||υn+1 − g(P, υn )||. If limn→∞ en = 0 implies that limn→∞ υ n = P , then the iteration method $n+1 = g(P, $n ) is said to be P-stable. Without any loss of generality, an assumption is made that, the sequence {υn } contains an upper boundary; else we cannot anticipate the probability of convergence. If every condition is satisfied for υn+1 = P υn which is termed as Picard iteration, therefore the iteration is P-stable. We shall then present the following theorem. 24

Theorem 4. Let (Ξ, || · ||) be a Banach space and P a self-map of Ξ. Let $n+1 = g(P, $n ) ||Pκ − P$ || ≤ C||κ − Pκ || + C||κ − υ||,

(82)

for all κ, υ ∈ Ξ, where 0 ≤ C, 0 ≤ C ≤ 1. Suppose that P is Picard P-stable. We take into consideration the following recursive formula given by Eq. (81) linked with Eq. (76). Theorem 5. Let ∆ be a self-map defined as ( ∆[x(n) (t)] = x(n+1) (t) = x(n) (t) + ST

!

1−α

−1

 B(α)αΓ(α + 1)Eα −

∆[y(n) (t)] = y(n+1) (t) = y(n) (t)+ST

−1

∆[z(n) (t)] = z(n+1) (t) = z(n) (t)+ST

−1

1 α 1−α u

!

1−α

n



B(α)αΓ(α + 1)Eα −

1 α 1−α u

h i o  ·ST x(n−1) (t)z(n−1) (t)−y(n−1) (t) (s) (t), !

1−α

n



) i · ST βy(n−1) (t) − C(x(n−1) (t) + p) , h



B(α)αΓ(α + 1)Eα −

1 α 1−α u

h i o  ·ST −x(n−1) (t)y(n−1) (t)−z(n−1) (t)+1 (s) (t),

is k-stable in L1 (a, b) if ||∆(xn (t)) − ∆(xm (t))|| ≤ ||xn (t) − xm (t)||(1 − λc(π) − µa(π)), ||∆(yn (t)) − ∆(ym (t))|| ≤ ||yn (t) − ym (t)||(1 + λc(π) − (µ + ρ)d(π)), ||∆(zn (t)) − ∆(zm (t))|| ≤ ||zn (t) − zm (t)||(1 + ρf (π) − (µ + γ + d)g(π)). Proof. First considering that ∆ has a fixed point; applying norm on both sides and without loss of generality, we have ||∆(xn (t)) − ∆(xm (t))|| = ( ! ) n h io 1−α −1   ·ST βy(n) (t)−C(x(n) (t)+p)− βy(m) (t)−C(x(m) (t)+p) = xn (t)−xm (t)+ST , 1 α B(α)αΓ(α + 1)E − u α

1−α

( 1−α −1  ≤ ||xn (t)−xm (t)||+ ST B(α)αΓ(α + 1)Eα −

) n  ·ST βy(n) (t)−C(x(n) (t)+p)−βy(m) (t)−C(x(m) (t)+p) , 1 α 1−α u (83) !

or ( 1−α −1  ≤ ||xn (t)−xm (t)||+ ST B(α)αΓ(α + 1)Eα −

) n  ·ST −(βy(n) (t)−C(x(n) (t)+p))+βy(m) (t)−C(x(m) (t)+p) . 1 α 1−α u (84) !

Thus, ||∆(xn (t)) − ∆(xm (t))|| ≤ ||xn (t) − xm (t)||(1 − λc(π) − µa(π)), ( ! ) where c(π) and a(π) are the ST −1

1−α 

1 B(α)αΓ(α+1)Eα − 1−α uα



· ST .

With the same idea we have the following ||∆(yn (t)) − ∆(ym (t))|| ≤ ||yn (t) − ym (t)||(1 + λc(π) − (µ + ρ)d(π)), ||∆(zn (t)) − ∆(zm (t))|| ≤ ||zn (t) − zm (t)||(1 + ρf (π) − (µ + γ + d)g(π)), we conclude that, P is Picard p-stable. 25

(85)

Let (Ξ, || · ||) be a Banach space and P a self-map of Ξ. Let $n+1 = g(P, $n ) be particular recursive procedure. Suppose that, F (P ) the fixed-point of P has at least one element and that $n converges to a point P ∈ F (P ). Theorem 6. We show that the special solution of Eq. (76) using the iteration method (81) is unique. Proof. Consider the following Hilbert space H = L2 ((a, b) × (0, k)) that can be defined as the set of those functions Z Z V : (a, b) × [0, T ] → R, uvdudv < ∞. (86) We now consider the following operator   βy(t) − C(x(t) + p), x(t)z(t) − y(t), η(0, 0), η =  −x(t)y(t) − z(t) + 1. We prove the inner product of (5(x(t) − X(t), y(t) − Y (t), z(t) − Z(t), (ω1 , ω2 , ω3 ))),

(87)

where (x(t) − X(t), y(t) − Y (t), z(t) − Z(t)), are special solutions of the above system. Using the relationship between the norm and the inner function we obtain   − β(y(t) − Y (t)) − C((x(t) − X(t)) + p) ≤ ||β(y(t) − Y (t) − C((x(t) − X(t)) + p)||||ω1 || + µ||x(t) − X(t)||||ω1 ||, 

 − β(y(t) − Y (t)) − C((x(t) − X(t)) + p), ω1 ≤ β||y(t) − Y (t)||||ω1 || − C||(x(t) − X(t))||||ω1 || + ||p||||ω1 ||,   (x(t) − X(t))(z(t) − Z(t)) − (y(t) − Y (t)), ω2 ≤ ||(x(t) − X(t))(z(t) − Z(t))||||||ω2 || − ||(y(t) − Y (t))||||ω2 ||, (88)   −(x(t)−X(t))(y(t)−Y (t))−(z(t)−Z(t))+1, ω3 ≤ −||(x(t)−X(t))(y(t)−Y (t))||||ω3 ||−||(z(t)−Z(t))||||ω3 ||+1||ω3 ||, for large number m, n and k the solutions converge to the exact solution. Using the topology concept, we can find three very small positive parameters (λm , λn and λk ) ||x(t) − X(t)|| <

λm , $

λn , ε λk , ||z(t) − Z(t)|| < π

||y(t) − Y (t)|| <

(89)

where,   $ = 3 − β(y(t) − Y (t)) − C((x(t) − X(t)) + p) ||ω1 ||,   ε = 3 (x(t) − X(t))(z(t) − Z(t)) − (y(t) − Y (t)) ||ω2 ||,   π = 3 − (x(t) − X(t))(y(t) − Y (t)) − (z(t) − Z(t)) + 1 ||ω3 ||.

(90)

Using the topology concept, we conclude that λm < 0, λn < 0 and, where   $ = 3 − ||β(y(t) − Y (t))|| − ||C((x(t) − X(t))|| + p) 6= 0,   ε = 3 ||(x(t) − X(t))(z(t) − Z(t))|| − ||(y(t) − Y (t))|| 6= 0,   π = 3 − ||(x(t) − X(t))(y(t) − Y (t))|| − ||(z(t) − Z(t))|| + 1 6= 0,

(91) (92)

where ||ω1 ||, ||ω2 || and ||ω3 || 6= 0; ||X(t) − x(t)|| = ||Y (t) − y(t)|| = ||Z(t) − z(t)|| = 0; x(t) = X(t), y(t) = Y (t) and z(t) = Z(t). 26

Adams-Moulton method. Using the numerical scheme given by Eq. (45), we have ( ( ! !)) y(n+1) (t) − y(n) (t) x(n+1) (t) − x(n) (t) 1−α n x(n+1) (t) − x(n) (t) = x(0) (t) + β −C +p + B(α) 2 2 ( ! !) ∞ y(ν+1) (t) − y(ν) (t) x(ν+1) (t) − x(ν) (t) α X α ξ β −C +p , + B(α) ν=0 ν 2 2 ( ( ! ! !)) y(n+1) (t) − y(n) (t) x(n+1) (t) − x(n) (t) z(n+1) (t) − z(n) (t) 1−α n y(n+1) (t) − y(n) (t) = y(0) (t)+ − + B(α) 2 2 2 ( ! ! !) ∞ y(ν+1) (t) − y(ν) (t) x(ν+1) (t) − x(ν) (t) z(ν+1) (t) − z(ν) (t) α X α + − , (93) ξ B(α) ν=0 ν 2 2 2 ( ( ! ! ! )) x(n+1) (t) − x(n) (t) y(n+1) (t) − y(n) (t) z(n+1) (t) − z(n) (t) 1−α n z(n+1) (t)−z(n) (t) = z(0) (t)+ − − +1 + B(α) 2 2 2 ( ! ! ! ) ∞ x(ν+1) (t) − x(ν) (t) y(ν+1) (t) − y(ν) (t) z(ν+1) (t) − z(ν) (t) α X α + ξ − − +1 . B(α) ν=0 ν 2 2 2 We shown the chaos mechanism in the vallis model for El Ni˜ no phenomenon. We consider the following dimensionless parameters set: β = 102, C = 3 and p = 0.5. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1, y(0) = 1 and z(0) = 1. The numerical results given in Figs. 10a-10i shows the phase portraits for the fractional-order vallis model for El Ni˜ no phenomenon in the Atangana-Baleanu-Caputo sense applying the solution scheme (93) for several values of α ∈ (0, 1], arbitrarily chosen. Fractal-Fractional power-law derivative. The fractional-fractional vallis model for El Ni˜ no phenomenon involving the power-law derivative corresponding to the system (61) is given by F F −P α,τ 0 Dt x(t) = Ξ1 (x(t), y(t), t), F F −P α,τ 0 Dt y(t)

= Ξ2 (x(t), y(t), t),

F F −P α,τ 0 Dt x(t)

= Ξ3 (x(t), y(t), t),

(94)

where Ξ1 (x(t), y(t), z(t), t) = βy(t) − C(x(t) + p), Ξ2 (x(t), y(t), z(t), t) = x(t)z(t) − y(t), Ξ3 (x(t), y(t), z(t), t) = −x(t)y(t) − z(t) + 1, α, τ ∈ (0, 1] and F F −P0 Dtα,τ is the power-law fractal-fractional derivative given by Eq. (49). The numerical solution of Eq. (94) is obtained considering the numerical scheme given by Eq. (54) n   τ (∆t)α X h τ −1 tj Ξ1 (xj , y j , z j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − Γ(α + 2) j=0  i −1 tτj−1 Ξ1 ((xj−1 , y j−1 , z j−1 , tj−1 ) (n + 1 − j)α+1 − (n − j)α (n − j + 1 + α) ,

xn+1 (t) = x(0) +

n   τ (∆t)α X h τ −1 tj Ξ2 (xj , y j , z j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − Γ(α + 2) j=0  i −1 tτj−1 Ξ2 (xj−1 , y j−1 , z j−1 , tj−1 ) (n + 1 − j)α+1 − (n − j)α (n − j + 1 + α) ,

yn+1 (t) = y(0) +

(95) 27

20

1.5 =1 =0.95 =0.90 =0.85

15

1 =1 =0.95 =0.90 =0.85

1

0.5

10 0.5 0

z(t)

y(t)

x(t)

5 0

0

-0.5 -0.5

-5

-15

-1.5 0

20

40

60

80

100

120

140

160

180

200

=1 =0.95 =0.90 =0.85

-1

-1

-10

-1.5 0

20

40

60

80

100

t (days)

120

140

160

180

200

0

20

40

60

80

t (days)

(a)

100

120

140

160

180

200

t (days)

(b)

(c)

1.5

1 1

=1 =0.95 =0.90 =0.85

1

=1 =0.95 =0.90 =0.85

0.5

=1 =0.95 =0.90 =0.85

0.5

0.5 0

z(t)

z(t)

y(t)

0 0

-0.5

-0.5

-1

-1

-0.5

-1

-1.5 -15

-10

-5

0

5

10

15

-1.5 -15

20

-10

-5

0

x(t)

(d)

10

15

-1.5 -1.5

20

0.5

1

1.5

20 =1 =0.95 =0.90 =0.85

10

x(t)

0

0

(f)

0

=1 =0.95 =0.90 =0.85

10

x(t)

=1 =0.95 =0.90 =0.85

0.5

-0.5

y(t)

20

1

-1

(e)

1.5

y(t)

5

x(t)

0

-0.5 -10

-10

-1 -1.5 1

-20 2 20

0

10 0

-1

-20 1 1

1 0.5

0

0 -1

-10

z(t)

-2

-20

x(t)

-0.5

y(t)

-2

(g)

-1 -1.5

(h)

z(t)

2

0

1 0

-1 -1

z(t)

-2

-2

(i)

Figure 10: Numerical simulation for the vallis model for El Ni˜ no phenomenon given by Eq. (76).

28

y(t)

n   τ (∆t)α X h τ −1 zn+1 (t) = z(0) + tj Ξ3 (xj , y j , z j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − Γ(α + 2) j=0  i −1 tτj−1 Ξ3 (xj−1 , y j−1 , z j−1 , tj−1 ) (n + 1 − j)α+1 − (n − j)α (n − j + 1 + α) .

We shown the chaos mechanism in the vallis model for El Ni˜ no phenomenon. We consider the following dimensionless parameters set: β = 102, C = 3 and p = 0.5. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1, y(0) = 1 and z(0) = 1. The numerical results given in Figs. 11a-11i shows the phase portraits for the fractional-order vallis model for El Ni˜ no phenomenon involving the fractal-fractional Liouville-Caputo derivative considering the solution scheme (95) for several values of α − τ ∈ (0, 1], arbitrarily chosen.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 11: Numerical simulation for the vallis model for El Ni˜ no phenomenon given by Eq. (94).

29

Fractal-Fractional Mittag-Leffler derivative. The fractional-fractional vallis model for El Ni˜ no phenomenon involving the Mittag-Leffler memory corresponding to the system (61) is given by F F −M α,τ 0 Dt x(t) = Ξ1 (x(t), y(t), t), F F −M α,τ 0 Dt y(t)

= Ξ2 (x(t), y(t), t),

F F −M α,τ 0 Dt x(t)

= Ξ3 (x(t), y(t), t),

(96)

where Ξ1 (x(t), y(t), z(t), t) = βy(t) − C(x(t) + p), Ξ2 (x(t), y(t), z(t), t) = x(t)z(t) − y(t), Ξ3 (x(t), y(t), z(t), t) = −x(t)y(t) − z(t) + 1, α, τ ∈ (0, 1] and F F −P0 Dtα,τ is the Mittag-Leffler fractal-fractional derivative given by Eq. (57). The numerical solution of Eq. (96) is obtained considering the numerical scheme given by Eq. (60) τ tτ −1 (1 − α) τ (∆t)α xn+1 (t) = x(0) + n Ξ1 (xn , y n , z n , tn ) + × AB(α) AB(α)Γ(α + 2) " n   X tτj −1 Ξ1 (xj , y j , z j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − j=0

#   −1 tτj−1 Ξ1 (xj−1 , y j−1 , z j−1 , tj−1 ) (n − j + 1)α+1 − (n − j)α (n − j + 1 + α) , τ (∆t)α τ tτ −1 (1 − α) Ξ2 (xn , y n , z n , tn ) + × yn+1 (t) = y(0) + n AB(α) AB(α)Γ(α + 2) " n   X tτj −1 Ξ2 (xj , y j , z n , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) −

(97)

j=0



−1 tτj−1 Ξ2 (xj−1 , y j−1 , z j−1 , tj−1 )

(n − j + 1)

α+1

#  − (n − j) (n − j + 1 + α) , α

τ tτ −1 (1 − α) τ (∆t)α zn+1 (t) = z(0) + n Ξ1 (xn , y n , z n , tn ) + × AB(α) AB(α)Γ(α + 2) " n   X tτj −1 Ξ1 (xj , y j , z j , tj ) (n + 1 − j)α (n − j + 2 + α) − (n − j)α (n − j + 2 + 2α) − j=0 −1 tτj−1 Ξ1 (xj−1 , y j−1 , z j−1 , tj−1 )

#   α+1 α (n − j + 1) − (n − j) (n − j + 1 + α) ,

We shown the chaos mechanism in the vallis model for El Ni˜ no phenomenon. We consider the following dimensionless parameters set: β = 102, C = 3 and p = 0.5. The time simulation was t = 200 [days], step size h = 1 × 10−3 , initial conditions: x(0) = 1, y(0) = 1 and z(0) = 1. The numerical results given in Figs. 12a-12i shows the phase portraits for the fractional-order vallis model for El Ni˜ no phenomenon involving the fractal-fractional Atangana-Baleanu derivative considering the solution scheme (97) for several values of α − τ ∈ (0, 1], arbitrarily chosen.

3

Conclusions

In this paper, we obtain multiple attractors and periodicity using the stochastic approach and differential and integral operators with power-law and Mittag-Leffler law for the coupled dynamical El Ni˜ no/La Ni˜ na-Southern oscillation model and the continuous-time vallis model for el Ni˜ no. Chaos diagram and attractor plots are applied to illustrate the coexisting attractors and complexity of the systems. We prove the existence of a unique set of 30

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 12: Numerical simulation for the vallis model for El Ni˜ no phenomenon given by Eq. (96).

31

the solutions of the new model using the fixed-point theorem. The generalized models with power-law and MittagLeffler kernel were solved numerically via the Adams-Bashforth-Moulton, Adams-Moulton and Atangana-Toufik schemes. In the fractal-fractional cases, we presented numerical simulations for different values of fractional order α and fractal dimension τ . The numerical simulations show that these operators captures more strange behaviors and irregularities that the obtained with local derivative.

Acknowledgments Jos´e Francisco G´ omez Aguilar acknowledges the support provided by CONACyT: C´atedras CONACyT para j´ovenes investigadores 2014 and SNI-CONACyT.

Conflicts of Interest The author declare no conflict of interest.

References [1] D. Kumar, D. Baleanu. Fractional Calculus and its Applications in Physics. Frontiers in Physics, 7, (2019), 1-18. [2] A. Atangana. Fractional discretization: The African’s tortoise walk. Chaos, Solitons & Fractals, 130, (2020), 1-18. [3] A. Atangana, T. Mekkaoui. Trinition the complex number with two imaginary parts: Fractal, chaos and fractional calculus. Chaos, Solitons & Fractals, 128, (2019), 366-381. [4] T. Abdeljawad, Q.M. Al-Mdallal, F. Jarad. Fractional logistic models in the frame of fractional operators generated by conformable derivatives. Chaos, Solitons & Fractals, 119, (2019), 94-101. [5] X.J. Yang, M. Abdel-Aty, C. Cattani. A new general fractional-order derivative with Rabotnov fractionalexponential kernel applied to model the anomalous heat transfer. Thermal Science, 23(3 Part A), (2019), 16771681. [6] E. Bonyah, A. Atangana, M. Chand. Analysis of 3D IS-LM macroeconomic system model within the scope of fractional calculus. Chaos, Solitons & Fractals: X, 1, (2019), 1-18. [7] X.J. Yang, J.A. Tenreiro-Machado. A new fractal nonlinear Burgers’ equation arising in the acoustic signals propagation. Mathematical Methods in the Applied Sciences, 1, (2019), 1-14. ¨ [8] A. Fernandez, M.A. Ozarslan, D. Baleanu. On fractional calculus with general analytic kernels. Applied Mathematics and Computation, 354, (2019), 248-265. [9] M.M. El-Dessoky, M.A. Khan. Application of fractional calculus to combined modified function projective synchronization of different systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(1), (2019), 1-13. [10] X.J. Yang, F. Gao, J.T. Machado, D. Baleanu. A new fractional derivative involving the normalized sinc function without singular kernel. The European Physical Journal Special Topics, 226(16-18), (2017), 3567-3575. [11] X.J. Yang, F. Gao, H.M. Srivastava. New rheological models within local fractional derivative. Rom. Rep. Phys, 69(3), (2017), 1-11. [12] E.F.D. Goufo. On chaotic models with hidden attractors in fractional calculus above power law. Chaos, Solitons & Fractals, 127, (2019), 24-30. [13] X.J. Yang, F. Gao, Y. Ju, H.W. Zhou. Fundamental solutions of the general fractional-order diffusion equations. Mathematical Methods in the Applied Sciences, 41(18), (2018), 9312-9320.

32

[14] X.J. Yang, Y.Y. Feng, C. Cattani, M. Inc. Fundamental solutions of anomalous diffusion equations with the decay exponential kernel. Mathematical Methods in the Applied Sciences, 1, (2019), 1-16. [15] J. Crum, J.A. Levine, A. Gillette. Extending discrete exterior calculus to a fractional derivative. ComputerAided Design, 114, (2019), 64-72. [16] X.J. Yang. New-rheological problems involving general fractional derivatives with nonsingular power-law kernels. Proceedings of the Romanian Academy Series A-Mathematics Physics Technical Sciences Information Science, 6, (2017), 1-8. [17] X.J. Yang, H.M. Srivastava, J.T. Machado. A new fractional derivative without singular kernel. Therm. Sci, 20(2), (2016), 753-756. [18] R. Behinfaraz, S. Ghaemi, S. Khanmohammadi. Adaptive synchronization of new fractional-order chaotic systems with fractional adaption laws based on risk analysis. Mathematical Methods in the Applied Sciences, 42(6), (2019), 1772-1785. [19] A. Bakhet, Y. Jiao, F. He. On the Wright hypergeometric matrix functions and their fractional calculus. Integral Transforms and Special Functions, 30(2), (2019), 138-156. [20] A. Atangana, E.F.D. Goufo. The Caputo-Fabrizio fractional derivative applied to a singular perturbation problem. International Journal of Mathematical Modelling and Numerical Optimisation, 9(3), (2019), 241-253. [21] M. Caputo, F. Mainardi. A new dissipation model based on memory mechanism. Pure Appl. Geophys. 91, (1971), 134-147. [22] M. Caputo, M. Fabricio. A New Definition of Fractional Derivative without Singular Kernel. Progr. Fract. Differ. Appl., 1(2), (2015), 73-85. [23] A. Atangana, D. Baleanu. New Fractional Derivatives with Nonlocal and Non-Singular Kernel. Theory and Application to Heat Transfer Model. Therm Sci. 20(2), (2016), 763-769. [24] A. Atangana. Fractal-fractional differentiation and integration: Connecting fractal calculus and fractional calculus to predict complex system. Chaos, Solitons & Fractals, 102, (2017), 396-406. [25] A. Atangana, S. Qureshi. Modeling attractors of chaotic dynamical systems with fractal-fractional operators. Chaos, Solitons & Fractals, 123, (2019), 320-337. [26] F.F. Jin. An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model. Journal of the atmospheric sciences, 54(7), (1997), 811-829. [27] Y. Zeng. The Laplace-Adomian-Pade technique for the ENSO model. Mathematical Problems in Engineering, 1, (2013), 1-15. [28] J.Q. Mo, W.T. Lin, J. Zhu. The variational iteration solving method for El Ni˜ no/La Ni˜ na-Southern Oscillation model. Adv Math, 35(2), (2006), 232-236 . [29] J.Q. Mo, W.T. Lin. Generalized variation iteration solution of an atmosphere-ocean oscillator model for global climate. J Syst Sci Complex, 24(2), 2011, 271-27 . [30] C.X. Qun, S.J. Qiang, Z.X. Qian, Z.L. Lun, Z.W. Min, Z. Jun. Modified variational iteration method for an El Ni˜ no Southern Oscillation delayed oscillator. Chin Phys B, 21(2), (2012), 1-20. [31] J. Singh, D. Kumar, J.J. Nieto. Analysis of an El Ni˜ no-Southern Oscillation model with a new fractional derivative. Chaos, Solitons & Fractals, 99, (2017), 109-115. [32] G.K. Vallis. El ni˜ no: a chaotic dynamical system?, Science 232, (1986), 243-245. [33] G.K. Vallis. Conceptual models of El Ni˜ no and the southern oscillation. J Geophys Res, 93, (1988), 13979-13991. [34] M. Borghezan, P.C. Rech. Chaos and periodicity in Vallis model for El Ni˜ no. Chaos, Solitons & Fractals, 97, (2017), 15-18. 33

[35] B. Alkahtani, A. Atangana. Chaos on the Vallis model for El Ni˜ no with fractional operators. Entropy, 18(4), (2016), 1-17. [36] A. Atangana, S. Jain. The role of power decay, exponential decay and Mittag-Leffler function’s waiting time distribution: Application of cancer spread. Physica A: Statistical Mechanics and its Applications, 512, (2018), 330-351. [37] L. Changpin, T. Chunxing. On the fractional Adams method. Computers & Mathematics with Applications, 58(8), (2009), 1573-1588. [38] I. Podlubny. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Academic press, New York, 198, (1998). [39] K. Diethelm. The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type. Springer Science & Business Media, (2010). [40] E. Zeidler. Non-linear functional analysis and its application. Springer-Verlag, New York, (1986). [41] B.S.T. Alkahtani. Chua’s circuit model with Atangana-Baleanu derivative with fractional order. Chaos, Solitons & Fractals, 89, (2016), 1-5. [42] J. Singh, Y.S. Shishodia. A New Reliable Approach for Two-Dimensional and Axisymmetric Unsteady Flows Between Parallel Plates. Zeitschrift f¨ ur Naturforschung A, 68(10-11), (2013), 629-634. [43] J. Singh, D. Kumar, M.A. Qurashi, D. Baleanu. A Novel Numerical Approach for a Nonlinear Fractional Dynamical Model of Interpersonal and Romantic Relationships. Entropy, 19(7), (2017), 1-17.

34

Highlights We obtain multiple attractors and periodicity using differential and integral operators with powerlaw and Mittag-Leffler law. Liouville-Caputo and Atangana-Baleanu fractional derivatives are considered. Novel differential and integral operators with fractional order and fractal dimension are considered. Numerical simulations are presented for specific parameters. The results obtained presents more information that were not revealed in the models with local derivative.