key: cord-0147478-mfq9ngag
authors: Aljhani, Sami; Noorani, Mohd Salmi Md; Douaifia, Redouane; Abdelmalek, Salem
title: A high-order predictor-corrector method for initial value problems with fractional derivative involving Mittag-Leffler kernel: epidemic model case study
date: 2022-05-07
journal: nan
DOI: nan
sha: f9fd5c4f2716ea41d2c0521663dce24be32439e0
doc_id: 147478
cord_uid: mfq9ngag

In this paper, we propose a numerical scheme of the predictor-corrector type for solving nonlinear fractional initial value problems, the chosen fractional derivative is called the Atangana-Baleanu derivative defined in Caputo sense (ABC). This proposed method is based on Lagrangian quadratic polynomials to approximate the nonlinearity implied in the Volterra integral which is obtained by reducing the given fractional differential equation via the properties of the ABC-fractional derivative. Through this technique, we get corrector formula with high accuracy which is implicit as well as predictor formula which is explicit and has the same precision order as the corrective formula. On the other hand, the so-called memory term is computed only once for both prediction and correction phases, which indicates the low cost of the proposed method. Also, the error bound of the proposed numerical scheme is offered. Furthermore, numerical experiments are presented in order to assess the accuracy of the new method on two differential equations. Moreover, a case study is considered where the proposed predictor-corrector scheme is used to obtained approximate solutions of ABC-fractional generalized SI (susceptible-infectious) epidemic model for the purpose of analyzing dynamics of the suggested system as well as demonstrating the effectiveness of the new method to solve systems dealing with real-world problems.

Despite the fact that fractional calculus has a lengthy history in mathematics, it has only recently seen a considerable number of real-world applications. Fractional derivatives are important due to their attractive characteristics (e.g. memory effect), including their wide dynamical range [1, 2, 3] . Many definitions of fractional derivatives and integrals exist, including Riemann-Liouville, Caputo, Grunwald-Letnikov [4] . In 2015, Caputo and Fabrizio (CF) have suggested a new idea of fractional derivatives based on the exponential decay [5] . Following that, Atangana and Baleanu (AB) proposed a novel concept of fractional derivatives with non-singular and non-local kernel based on the Mittag-Leffler function [6] . Due to the difficulty or impossibility of obtaining to obtain exact solution of fractional differential equations, numerical techniques are required to provide approximate solutions. Over the past few decades, numerical methods to solve initial value problems have attracted great interest from the research community, which support various disciplines, including medicine, chemistry, biology, physics, business, and the arts, etc. Where computational algorithms are implemented. Moreover, numerical analysis techniques are the perfect tools to evaluate the behavior of real-life complicated models. Therefore, the implementation of accurate numerical methods is suibtle in the problems representing real world phenomena [6, 7, 8, 9] . Over 2000 years ago, the simplest method in the literature was introduced with the concept of interpolation. In contrast, several other techniques were introduced for nonlinear systems, such as Newton's method, Euler's method, Gaussian elimination, and Lagrange interpolation [10, 11, 12, 13] . Over the past few years, several numerical tools have been developed for solving fractional ordinary differential equations with general nonlinearity [14, 15, 16, 17, 18, 19, 20] . Such as the fractional Adams-Bashforth method, which is the most popular scheme in which the basis of Lagrange interpolation is used [22, 23] and the fractional Adams-Bashforth/Moulton method developed with Newton linear/quadratic interpolations [24, 25] . Baleanu et al. [17] used the definition of the Atangana-Baleanu derivative in the Caputo sense to produce a predictor-corrector approach for fractional differential equation as well as to prove the existence and uniqueness of such problem. They have converted the ABC-fractional ordinary differential equation into a Volterra Integral equation, then they have used the rectangle rule for the predictor and the trapezoid rule for the corrector. On the other hand, based on the techniques and approximations appeared in [18, 19] for finding approximate solutions of fractional differential equations with Caputo and Caputo-Fabrizio fractional derivatives, respectively. We can to propose a new and efficient predictor-corrector method for solving ABC-fractional initial value problems by using quadratic lagrange interpolation, memory term, and second-order Taylor's expansion.

The following is a breakdown of the paper's structure. Basic definitions and notations, including the Atangana-Baleanu fractional derivatives and integral, are introduced in section 2. In section 3, a new numerical method is proposed for solving fractional initial value problems involving the Atangana-Baleanu fractional derivative in Caputo sense. The schemes' error estimates are obtained in section 4. Finally, in section 5, we present numerical examples to demonstrate the efficacy of the proposed scheme, and we also compare the solutions with other methods [8, 17] . Also, the stability analysis of the ABC-fractional-order SI (Susceptible-Infection) epidemiological model and its comparison with the numerical results obtained via the proposed scheme, which investigates the significance of disease prevalence in the community, is discussed in the last section.

The Atangana-Baleanu fractional derivative in the Caputo sense (ABC) is defined as (cf. [6] ):

where α ∈ (0, 1), AB(α) > 0 is a normalization function obeying AB(0) = AB(1) = 1 e.g.

, E α (.) denotes the Mittag-Leffler function of order α defined by

, z, α ∈ C, and Re(α) > 0,

and Γ(.) denotes Gamma function, defined as

The fractional integral for the ABC, which is newly defined with a nonlocal kernel and does not have singularities at t = s, is defined as follows (cf. [27] ):

We consider the initial-value problem with Atangana-Baleanu derivative

where f is a smooth nonlinear function that guarantees the existence of a unique solution for (5) , with the fractional order α ∈ (0, 1), and y 0 ∈ R. A continuous function y(t) is the solution of (5) if and only if it is the solution of following Voltera-integral equation:

At point t n+1 = h(n + 1), n = 0, 1, . . . , N ∈ N with h = T N , we get

where, y lag n+1 and y inc n+1 are called lag term and increment term respectively, which take the following forms:

and

To start presenting the proposed numerical method to solve the initial-value problem (5), we need the following lemma:

is the space of all polynomials of degree less than or equal to two. Let ψ k , k = 0, ...., N be the restricted value of ψ(t) on t k ( 0 ≤ k ≤ N ). Then there

Now, we use quadratic Lagrange polynomial of f (s, y(s)) over the intervals [t i−1 , t i+1 ] , 1 ≤ i ≤ N − 1:

where

But, on [t 0 , t 1 ] we can interpolate f (s, y(s)) with the points t 0 , t 1 2 , t 1 , then we get

where

.

The approximation of y(t n+1 ) denoting by y n+1 , according to (11)- (14) , y n+1 takes the following form:

where

and

where,

The predictor term can be approximated as follows:

where

Throughout this section, we need the following lemmas: be non-negative sequences with second one is monotonic increasing and satisfy that

where, M > 0 is independent of h > 0, and 0 < θ ≤ 1. Then,

Lemma 5. ( [18] ) There exist K 1 , K 2 > 0 such that for α ∈ (0, 1), and j = 0, 1, 2, we have

Let T P n+1 be the truncation error of prediction at point t n+1 , defined by

Theorem 1. Assume that f (., y(.)) ∈ C 3 ([0, T ]). Then, there exists C > 0 (independent of all grid parameters) such that:

Proof 1. We set the notation,

, then

where,

Thanks to lemma 2, there exists C 1 > 0, such that

The Taylor expansion of f around t n gives:

where,

According to lemma 2 and (33), we have

it follows that

According to (30) and (35), we have

and

We thus obtain the estimate

Consequently, there exists C 2 > 0 such that

According to (30) , (34) and (42), proof of theorem is achieved.

Error) Assume that f (., y(.)) ∈ C 3 ([0, T ]) and is Lipschitz continuous in its second argument, i.e.

Then, there exist K 1 , K 2 > 0 such that, the global predictor error satisfies

.

it follows that

Hence, according to the lemma 3 there exist C 1 , K 1 > 0, such that

On the other hand, for α ∈ (0, 1) we have

Consequently, there exists K 2 > 0 such that the estimate (48), becomes

. Then there exists positive constants C 1 and C 2 (independent of grid parameters) such that:

where,

with, AB g = α AB(α)Γ(α) .

Proof 3. It follows immediately that,

Thus, according to the lemma 2, there exist C, C > 0 such that

The desired estimate (52) can be derived by direct consideration of the last inequality (54). Then, we have

where C > 0 (independent of grid parameters), given E 1 ,

Proof 4. By taking into account the previous results, there exist C j > 0 (with j = 1, . . . , 10), such that

Consequently, by the discrete Gronwall's inequality (i.e. lemma 4) the desired result holds.

Since the global error indicated in the Theorem 4 is dependent on that of the start-up (E 1 2 , E 1 , and E 2 ). We suggest employing the start-up scheme described in Appendix A to generate approximate solutions for the first stages (i.e. y 1 2 , y 1 , and y 2 ).

In this section, we give some numerical experiments through the proposed predictor-corrector scheme (PPC) (15)-(21), the predictor-corrector method introduced by Baleanu-Jajarmi-Hajipour (BJH-PC) [17] , and the explicit numerical scheme introduced by Toufik-Atangana (TAE) [8] , to show the efficiency and accuracy of our new method. The experimental order of convergence (EOC) is computed by

where AE is the absolute error, which takes the following form:

Example 1. Consider the following fractional Initial value problem (IVP): where α ∈ (0, 1), and n ∈ N. On account of [17] , the exact solution of (59) is given by

In addition, the problem (59) is numerically solved. (21) are roughly 3. Also, we note that the approximate solutions obtained by our proposed scheme get closer to the exact solutions by the increase in α. From Figure 1 and Figure 2 , we note that the approximate solution obtained by PPC (15)-(21) almost matches with the exact solution with small step size compared to its counterparts. These figures and tables, indicate the efficacy of the current predictor-corrector numerical method (PPC) (15)-(21).

Example 2. Consider the following Atangana-Beleanu-fractional differential equation: where α ∈ (0, 1). According to [17] , the exact solution of (61) is given by

with E α,β (.) denotes the Mittag-Leffler function of two parameters α and β, which defined by

, z, α, β ∈ C, and Re(α), Re(β) > 0,

and the exact solution (62) is calculated using the algorithm mlf.m (see [28] ) evaluated with accuracy 10 −12 . Furthermore, the problem (61) is numerically solved. Table 2 and table 2 show the comparison of absolute error of different numerical methods ((PPC) (15)- (21) , BJH-PC [17] , and TAE [8] ) for (61) with AB(α) = 1, AB(α) = 1 − α + α Γ(α) , t ∈ [0, 1] and various values of fractional order α ∈ {0.5, 0.55, 0.7, 0.9, 0.95}, where we notice that our method is superior to them in terms of accuracy (i.e. PPC (15)-(21) achieves a lower error than TAE and BJH-PC). From Figure 3 and Figure 4 , we note that the approximate solution obtained by PPC (15)-(21) almost matches with the exact solution with small step size compared to its counterparts. These figures and tables, again confirm the efficacy of the current predictor-corrector numerical method (PPC) (15)- (21) .

In the previous examples, we have considered some differential equations with known exact solutions. Let us now analyze a realistic ABC-fractional epidemic model using qualitative analysis (Lyapunov theory) of solutions and validate the theoretical results numerically by means of the proposed predictor-corrector method. We consider the epidemic phenomenon proposed in [29] , which is an extended version of the SI epidemic model with a nonlinear incidence ϕ, and we include the fractional derivative which involving Mittag-Leffer kernel. Thus, the novel system is described as follows:

where u(t), v(t) represent the numbers of susceptible individuals, and infective individuals at time t, respectively. Λ is the birth rate of the population, γ is the disease transmission rate, µ is the natural death rate, and σ := σ + µ, with σ is the death rate due to disease, which all are positive real numbers, and α ∈ (0, 1). The incidence function ϕ(v) introduces a nonlinear relation between the susceptible individuals and infective individuals. We suppose that ϕ ∈ C 1 (R + ; R + ), and satisfies:

According to the above condition and Lemma 1 in [30] , it is classical task to prove the existence of unique solution for system (64) (see e.g. [31, 32] ). As well as, the system (64) has a disease-free equilibrium

, and a unique endemic equilibrium E * := (u * , v * ) if R 0 := Λγ µσ ϕ (0) > 1 (cf. [29] ). In the remainder, the following notation R + := [0, +∞) will be useful. The feasible region of the suggested model (64) is given by

Now, we investigate the global asymptotic stability of the equilibrium points E 0 and E * , which is determined by the reproduction number R 0 , as a result, the cases R 0 1 and R 0 > 1 are treated independently. Firstly, we provide the global stability result of the disease-free equilibrium E 0 .

Theorem 5. Let α ∈ (0, 1) and R 0 1. Assume that (65) holds. Then, the disease-free equilibrium E 0 is globally asympotatically stable in the feasible region Υ.

Proof 5. We begin by rewriting the first equation of the system (64), with taking into account this relation-

We consider the following candidate Lyapunov function:

such that

Thanks to the property of ABC-fractional derivative given in Lemma 3.1 from [33] , we get

Then, by employing (64)-(66), and Lemma 1 [30] , we obtain

Therefore, R 0 1, guarantees that ABC 0 D α t V ≤ 0 for α ∈ (0, 1). Consequently, by the Lyapunov stability theorem (see Theorem 2.7 from [33] and [34] ), we conclude that the disease-free equilibrium E 0 is globally asymptotically stable in Υ.

Secondly, we provide the global stability result of the endemic equilibrium E * . Before do that, we state a lemma which will be useful (see Lemma 2 from [29] ): Lemma 6. Assume that the function ϕ satisfied (65) and

We thus get the following inequality

where v * is the second component of the equilibrium point E * .

Theorem 6. Let α ∈ (0, 1) and R 0 > 1. Suppose that (65) holds. Then, the endemic equilibrium E * is globally asymptotically stable in interior of the feasible region Υ.

Proof 6. The first two equation of system (64) at E * , clearly give

It follows that

We consider the following candidate Lyapunov function:

such that

We now apply the property of ABC-fractional derivative given in Lemma 3.2 from [33] , to obtain

From (75) and by immediately computation, we get

Therefore, the nonnegativity of function F (defined in (72)) on (0, +∞), and the lemma 6 ensure that ABC 0 D α t L ≤ 0 for α ∈ (0, 1). Consequently, by the Lyapunov stability theorem (see Theorem 2.7 and Theorem 2.8 [33] as well as [34] ), we get the desired result.

Next, we present some numerical simulation of the general model (64) through several sets of parameters (see table 3 ) to examine the effectiveness of the proposed predictor-corrector scheme in the current work (i.e. PPC (15)-(21)), as well as the feasibility of the theoretical findings in this example. The following is a description of the results:

• According to the first set of parameters in table 3, with ϕ(v) = v resp. ϕ(v) = v 1+0.01v , and (u 0 , v 0 ) = 0.52, Λ µ − 0.52 , the Theorem 5 guarantees (for α ∈ (0, 1) and any positive initial data belongs to Υ) the global asymptotic stability of the disease-free equilibrium E 0 = (u 0 , 0) = (0.6818, 0), which is confirmed by the numerical results depicted in Figure 5 (resp. Figure 7) with various values of fractional order α ∈ {0.8, 0.85, 0.9, 0.99}.

• According to the second set of parameters in table 3, with ϕ(v) = v resp. ϕ(v) = v 1+0.01v , where (u 0 , v 0 ) = 0.6, Λ µ − 0.6 , the Theorem 6 guarantees (for α ∈ (0, 1) and any positive initial data belongs to Υ) the global asymptotic stability of the endemic equilibrium E * = (u * , v * ) = (0.125, 0.5087) (resp. E * = (0.1256, 0.5083)), which is confirmed by the numerical results depicted in Figure 6 (resp. Figure  8) with various values of fractional order α ∈ {0.8, 0.85, 0.9, 0.99}.

• According to the third set of parameters in table 3, with ϕ(v) = v resp. ϕ(v) = v 1+0.01v , and α = 0.99, the Theorem 5 guarantees (for any positive initial data belongs to Υ) the global asymptotic stability of the disease-free equilibrium E 0 = (u 0 , 0) = (0.9375, 0), which is confirmed by the numerical results depicted in Figure 9 (resp. Figure 11) with various values of initial data.

• According to the fourth set of parameters in table 3, with ϕ(v) = v resp. ϕ(v) = v 1+0.01v , and α = 0.99, the Theorem 6 guarantees (for any positive initial data belongs to Υ) the global asymptotic stability of the endemic equilibrium E * = (u * , v * ) = (0.125, 0.52) (resp. E * = (0.1256, 0.0.5196)), which is confirmed by the numerical results depicted in Figure 10 (resp. Figure 12) with various values of initial data. 

A novel scheme for solving Caputo time-fractional nonlinear equations: theory and application

Asymptotic stability conditions for autonomous time-fractional reaction-diffusion systems

Numerical solution of fractional-order HIV model using homotopy method. Discrete Dynamics in Nature and Society

Fractional calculus in bioengineering

A new definition of fractional derivative without singular kernel

New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model

Numerical solutions of certain new models of the time-fractional Gray-Scott

New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models

Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus with optimal control analysis with a case study

Polynomial interpolation: Lagrange versus newton. Mathematics of computation

A note on convergence of Newton interpolating polynomials

Visualizing and understanding the components of Lagrange and Newton interpolation

Efficient algorithms for polynomial interpolation and numerical differentiation

The FracPECE subroutine for the numerical solution of differential equations of fractional order. Forschung und wissenschaftliches Rechnen

Detailed error analysis for a fractional Adams method

Numerical computation of a fractional derivative with non-local and nonsingular kernel

On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag-Leffler kernel

A high-order predictor-corrector method for solving nonlinear differential equations of fractional order. Fractional Calculus and applied analysis

A fast and high-order numerical method for nonlinear fractional-order differential equations with non-singular kernel

Numerical approaches to fractional calculus and fractional ordinary differential equation

Weakly singular discrete Gronwall inequalities

Numerical analysis for the fractional diffusion and fractional Buckmaster equation by the two-step Laplace Adam-Bashforth method

The decoupled Crank-Nicolson/Adams-Bashforth scheme for the Boussinesq equations with nonsmooth initial data

New numerical scheme with newton polynomial: theory, methods, and applications

A Newton interpolation based predictor-corrector numerical method for fractional differential equations with an activator-inhibitor case study

A predictor-corrector approach for the numerical solution of fractional differential equations

Chaos in a simple nonlinear system with Atangana-Baleanu derivatives with fractional order

Mittag-Leffler function

Global and local asymptotic stability of an epidemic reaction-diffusion model with a nonlinear incidence

On the global asymptotic stability of solutions to a generalised Lengyel-Epstein system

Existence and uniqueness results for a smoking model with determination and education in the frame of non-singular derivatives

Mathematical analysis of an extended SEIR model of COVID-19 using the ABCfractional operator

Stability and Lyapunov functions for systems with Atangana-Baleanu Caputo derivative: an HIV/AIDS epidemic model

The proof of Lyapunov asymptotic stability theorems for Caputo fractional order systems

The authors are grateful to Pr. Bongsoo Jang and Mr. Junseo Lee for sharing the algorithms file related to the methods developed in [18] . The work of S. Aljhani is supported by Taibah University -Saudi Arabia.

To find a desired accuracy for y 1 and y 2 , we find the approximate solutions at points t 1 4 and t 1 2 using the constant, linear and quadratic interpolation. Let I 0 (f a ) be the constant interpolation of f at point t = a that is I 0 (f a ) = f (a). Let I 1 (f a , f b ) and I 2 (f a , f b , f c ) be linear and quadratic interpolation of f with grids (a, b) and (a, b, c) respectively. In the algorithm below we describe how to find the approximate solution of y(t) at points t 1 4 , t 1 2 , t 1 , and t 2 (i.e. y 1 4 , y 1 2 , y 1 , and y 2 ) using a predictor-corrector scheme: • Approximate solution of y(t) at point t 1 4 :• Approximate solution of y(t) at point t 1 2 :• Approximate solution of y(t) at point t 1 :• Approximate solution of y(t) at point t 2 :where AB f = 1 − α AB(α) and AB g = α AB(α)Γ(α).