A universal predictor–corrector algorithm for numerical simulation of generalized fractional differential equations

This study introduces some remarks on generalized fractional integral and differential operators, which generalize some familiar fractional integral and derivative operators, with respect to a given function. We briefly explain how to formulate representations of generalized fractional operators. Then, mainly, we propose a predictor–corrector algorithm for the numerical simulation of initial value problems involving generalized Caputo-type fractional derivatives with respect to another function. Numerical solutions of some generalized Caputo-type fractional derivative models have been introduced to demonstrate the applicability and efficiency of the presented algorithm. The proposed algorithm is expected to be widely used and utilized in the field of simulating fractional-order models.

ory and non-locality characteristics, in many fields [1][2][3][4][5][6][7][8][9][10][11]. Fractional-order derivatives of a given function include the entire history of the function; the next state of a fractional-order model depends not only on its current state but also on all of its historical states. Fractional integral and derivative operators have been treated as the preferred mathematical operators for characterization of memory in physical modeling issue [12][13][14][15][16][17]. In the literature, several fractional operator definitions have been proposed. A family of fractional derivatives, such as Riemann-Liouville, Hadamard, and Caputo derivatives, are introduced using fractional integrals. The Riemann-Liouville fractional integral, which is one of the most studied definitions, of order α > 0 is defined by [1][2][3][4][5][6].
Depending on the above fractional integral definition, the Riemann-Liouville and Caputo fractional derivatives of order α > 0 are defined by respectively, such that m − 1 < α ≤ m and m ∈ IN . An alternative characterization of the Caputo differential operator can be introduced as a modification of Riemann-Liouville differential operator in the form The Caputo differential operator is used extensively to model many physical phenomena in fractional differential calculus applications because it has many features that are similar to ordinary ones. The next type of fractional integral operators, which is related to our study, is the Hadamard fractional integral operator given by [3,7] H I α The associated Hadamard fractional derivative operator of order α > 0 is given by where m = α . Recently, Katugampola [18,19] introduces a new fractional integral operator, which generalizes Riemann-Liouville integral operator, as follows where ρ > 0. Consequently, for m −1 < α ≤ m where m ∈ IN , the Riemann-type and Caputo-type fractional derivative operators corresponding to Katugampola fractional integral operator, given in (7), are introduced in [18,20] and [21] as, and respectively, where m = α and ρ > 0. Furthermore, Osler [22] highlighted a useful extension of fractional integral and derivative operators by initiating fractional-order integrals and derivatives of a given function f related to another function σ as and has been studied in [3], where σ is considered to be any monotone function on [a, b] that has a continuous derivative on (a, b), and in [7], where σ is considered to be any increasing and positive monotone function on [a, b] that has a continuous derivative on (a, b). Then, the fractional derivative concept D α,σ (x) [3,7] under the additional condition that σ (x) = 0 on [a, b]. These conditions, as will be shown in the next section, can be adjusted in a more convenient manner. In [23], Almeida introduces a Caputo-type fractional derivative with respect to the function σ based on the fractional integral given in (10), which has characteristics similar to those of the Caputo derivative given in (3), as On the other hand, the growing interest in applications involving fractional-order differential models has motivated the improvement, development and extension of analytical and numerical methods for solving such models. The predictor-corrector (P-C) approach, presented in [24], is one of the most robust, stable and efficient methods that is widely used for the numerical solution of initial value problems (IVPs) with Caputo fractional derivatives. Some applications and simulations of this approach can been found in [25][26][27][28][29][30][31]. In [32], an optimized algorithm of the P-C approaches has been developed. The suggested algorithm, which implements the first-order Taylor series linearization of the nonlinear equation, provides approximate numerical solutions to Caputo-type IVPs with better accuracy than those obtained using the approach discussed in [24]. In [21], an adaptive algorithm of the P-C approaches to handle IVPs containing the generalized Caputo-type fractional derivative defined in (9) is introduced. The details and the mathematical derivation of the adaptive P-C algorithm are presented in [21].
In the present study, we provide precise conditions on the function σ to define fractional integrals with respect to the function σ , namely generalized fractional integrals. The conditions, that σ is considered to be any strictly increasing function that has a continuous derivative on (a, b), are given in the sense that the generalized fractional integral, for a < x ≤ b, is well defined. Accordingly, we introduce fractional derivative representations with respect to the function σ , namely generalized fractional derivatives, under the additional condition that σ (x) = 0 on [a, b]. We present the formulation of generalized fractional integral and derivative operators with respect to the function σ in a way that they can be considered as generalizations of many well-known fractional operators. Then, basically, we suggest a universal P-C algorithm to simulate nonlinear IVPs containing generalized Caputo-type fractional derivatives with respect to another function numerically. Besides, some test problems are discussed to display the efficiency, performance, features and accuracy of the suggested algorithm. Throughout our work, σ will indicate a strictly increasing and differentiable function on (a, b), where d dx σ (x) = δ(x).

Generalized fractional integral operator
In this subsection, we consider the generalization of fractional integral given in formula (10), where σ is assumed to be any strictly increasing function that has a continuous derivative on (a, b), which will be called, in this work, the generalized fractional integral. The considered generalized fractional integral operator can be seen as a generalization of many well-known definitions, such as Riemann-Liouville, Hadamard and Katugampola integral operators. We explain the mathematical derivation approach and then present some properties of the studied generalized integral operator concept. First, let's recall that the Riemann-Liouville integral operator, given in (1), can be seen as a direct extension of the Cauchy integral formula where the integer n is replaced with any real α > 0, using the property that (n + 1) = n!. Now, to derive the examined generalized fractional integral operator consider, for n ∈ IN and strictly increasing differen- In case of n = 2, exchanging the order of integration, using the Dirichlet technique [5], we obtain Repeating the above iteration n − 1 times, we get which suggests the following definition, when the integer n is replaced with any real α > 0.

Definition 1
The generalized fractional integral operator of order α with respect to the function σ , where α > 0, a ≥ −∞, b ≤ ∞, and the function σ is assumed to ensure the convergence of the integral.

Remark 1
The assumption that the function σ must be strictly increasing on [a, b] is sufficient to ensure that σ (x) − σ (ξ) > 0 when a < ξ < x and that σ is injective, and so has an inverse on [a, b].
According to the concept of generalized fractional integral operator, introduced in Definition 1, we can simply draw the following integral operators as special cases 1. In the case of σ (x) = x, the generalized integral operator I α,x a+ reduces fully to the Riemann-Liouville integral operator I α a+ , given in (1). 2. In the case of σ (x) = ln x, the generalized integral operator I α,ln x a+ , where a ≥ 0, reduces to the Hadamard integral operator H I α a+ , given in (5). 3. In the case of σ (x) = x ρ /ρ, the generalized integral operator I α,x ρ /ρ a+ , where ρ > 0 and a ≥ 0, reduces to the Katugampola integral operator I α,ρ a+ f (x), given in (7). 4. In the case of σ (x) = x ρ , the generalized integral operator I α,x ρ a+ , where ρ > 0 and a ≥ 0, reduces to the integral operator Remark 2 The Erdélyi-Kober fractional integral operator, which is defined as [33] (20) with ρ > 0 and η ∈ IR, can be reformulated using the integral operator I α,x ρ a+ , given in (18), as It is known that the Riemann-Liouville fractional integral operator given in (1) satisfies the semigroup property I α a+ I β a+ = I α+β a+ , when α, β > 0. Now, we give a similar semigroup property of the generalized fractional integral operator.
Theorem 1 Let α > 0 and β > 0. Then, for a < x ≤ b, [23] I α,σ (x) Next, we investigate a useful relationship between the generalized fractional integral operator with respect to the function σ , I α,σ (x) a+ , and the Riemann-Liouville integral operator I α a+ . This relationship will be used to formulate classes of generalized fractional derivative operators.
Proof Implementing the change of variables u = σ (ξ) in the generalized integral formula of the function f related to the function σ , given in Definition 1, yields (24) where σ −1 is the inverse function of σ . Using the fact that formula (23) can be, thus, derived. This completes the proof.

Generalized fractional derivative operators
This subsection introduces, based on the representation of the generalized fractional integral operator presented in Definition 1, generalized forms of fractional derivative operators. These forms generalize the Riemann-Liouville, Caputo, Hadamard and Katugampola fractional derivative operators. Some of the generalized derivative operators will be derived based on relation (23). If 0 < α ≤ 1, we may introduce the generalized Riemann-type fractional derivative operator of order α with respect to the function σ , R D α,σ (x) a+ , by replacing the operator I α σ (a) in formula (23) with the composite operator D I 1−α σ (a) . One can easily recognize that, for where δ(x) = 0 on [a, b]. In general, using the con- is characterized as,

Definition 2
The generalized Riemann-type fractional derivative operator of order α with respect to the function σ , R D α,σ (x) a+ , is defined (provided it exists) by Remark 3 We can observe that in the case of σ (x) = x the generalized Riemann-type derivative operator R D α,x a+ reduces to the Riemann-Liouville derivative operator given in (2), in the case of σ (x) = ln x the generalized Riemann-type derivative operator R D α,ln x a+ reduces to the Hadamard derivative operator given in (6) and in the case of σ (x) = x ρ /ρ the generalized Riemann-type derivative operator R D α,x ρ /ρ a+ reduces to the Katugampola derivative operator given in (8). Now, we turn our concern to introduce a Caputo-type of generalized fractional derivative operators. The generalized Caputo-type fractional derivative, C D α,σ (x) a+ , can be obtained by replacing the operator I α σ (a) in formula (23) with the composite operator I m−α σ (a) D m . So, using the conceptual relationship is characterized as,

Definition 3
The generalized Caputo-type fractional derivative operator of order α with respect to the function σ , C D α,σ (x) a+ , is defined (provided it exists) by where m − 1 < α ≤ m and δ(x) = 0 on [a, b].
Remark 4 In particular, if 0 < α < 1, then, for a < x ≤ b, the generalized Caputo-type fractional derivative reduces to In view of the following results, which have been proved in [23], we can observe that the properties of the generalized Caputo-type fractional derivative appear closer to those of ordinary derivatives than those of the generalized Riemann-type fractional derivative. and

The universal predictor-corrector algorithm
In this section, we develop an extended P-C algorithm, namely the universal P-C algorithm, to conveniently, effectively and accurately offer numerical solutions for IVPs involving generalized Caputo-type fractional derivatives. We will derive the main steps of the suggested algorithm to deal with the IVP where C D α,σ (x) a+ is the generalized Caputo-type fractional derivative operator of order α with respect to the function σ given in Definition 3, a ≥ 0 and δ(x) = 0 on [a, b]. At first, for m − 1 < α ≤ m, m ∈ IN and y ∈ C m [a, b], IVP (31) is equivalent, referring to Theorem 3, to the integral equation where Let the functions f and σ be assumed so that a unique solution to IVP (31)

exits in the interval [a, b].
To describe our algorithm, let us define a grid in the interval [a, b] with N + 1 non-equispaced nodes x k , k = 0, 1, . . . N , given by where h = σ (b)−σ (a) N and σ −1 is the inverse function of σ . Now, we are going to generate the approximation y k ≈ y(x k ) over the non-uniform grid {x k , k = 1, 2, . . . N } to simulate numerically the solution of IVP (31). To do this, assuming that the numerical values of y(x) at x 0 , x 1 , . . . , x k have been evaluated which are denoted as y 0 , y 1 , . . . , y k where y 0 = y(a), we need to calculate the approximation y k+1 ≈ y(x k+1 ) by means of the integral Performing the substitution yields or equivalently Next, by applying the trapezoidal rule with respect to the weight function (σ (x k+1 ) − .) α−1 to the integrals given on the right-hand side of Eq. (38) (i.e., replacing the function f σ −1 (u), y(σ −1 (u)) by its piecewise linear interpolating polynomial with nodes chosen at the σ (x n ), n = 0, 1, . . . k + 1), we get y(x n+1 )) .

(39)
Hence, substituting the approximations given in Eq. (39) in to Eq. (38) and collecting like terms, we obtain the corrector formula for y(x k+1 ), k = 0, 1, . . . , N − 1, as where a n,k+1 = The principal step of our algorithm is to exchange the quantity y(x k+1 ) appeared on the right-hand side of Eq. (40) with the predictor value y P (x k+1 ). The predictor value y P (x k+1 ) can be gained by implementing the Adams-Bashforth method to the integral Eq. (37). In this case, if we replace the function f σ −1 (u), y(σ −1 (u)) by the quantity f (x n , y(x n )) at each integral in Eq. (38), we obtain That is Consequently, the universal P-C algorithm, for providing the numerical approximation y k+1 ≈ y(x k+1 ), is finally outlined by the rule k n=0 a n,k+1 f (x n , y n ) where y n ≈ y(x n ), n = 0, 1, . . . k, and the predicted value y P k+1 ≈ y P (x k+1 ) can be evaluated as expressed in Eq. (43) with the weights a n,k+1 being determined according to Eq. (41). m −1 < α ≤ m, m ∈ IN , and δ(x) = 0 on [a, b]. In the case of 0 < α ≤ 1 and δ(x) = 0 on (a, b), the same P-C algorithm for IVP (31) can be derived in a similar manner, using (28), as given in Eq. (44) where υ(x k+1 ) = y(a).

Remark 6
One can simply observe that in the case of σ (x) = x, our P-C algorithm reduces to the P-C approach presented in [24] and in the case of σ (x) = x ρ /ρ, where (0 < ρ ≤ 1 when m − 1 < α ≤ m and m ∈ IN ) or (ρ > 0 when 0 < α ≤ 1), our P-C algorithm reduces to the adaptive P-C algorithm presented in [21].

Remark 7
Firstly, it is clear, from the P-C rule given in Eq. (44), that the approximation y k+1 is based on the complete data values (y 0 , y 1 , · · · , y k ). This means that the generalized Caputo-type fractional derivative, introduced in Definition 3, holds memory effects. Secondly, in view of the results presented in [34], we can notice that the universal P-C algorithm features are similar to those of the classic Adams-Bashforth-Moulton method and the P-C algorithm of [24]. Hence, our P-C algorithm behaves satisfactorily according to its numerical stability.

Applications
This section deals with generalized Caputo-type fractional differential equations that can be viewed as fractional extensions of some known equations. Further, it investigates the efficiency and performance of the universal P-C algorithm, introduced in the previous section, using some cases of the function σ . Numerical solutions to the considered problems are provided by implementing the proposed P-C algorithm. The numerical simulation results are performed using the symbolic software Mathematica.

Example 1 Our first problem assumes the IVP
is the generalized Caputo-type fractional derivative operator, introduced in Definition 3, A ∈ IR, 0 < α ≤ 1 and σ (x) = x ρ such that 0 < ρ ≤ 1. To simulate numerically the IVP (45) over the interval [a, b] = [0, T ], on the basis of our P-C algorithm, we have 0, 1, . . . , N − 1) and h = T ρ N , for some N ∈ IN . Hence, the approximation y k+1 ≈ y(x k+1 ) can be described by the formula k n=0 a n,k+1 with the weights a n,k+1 being determined according to Eq. (41), and Here, our P-C algorithm aims at searching numerical solutions to IVP (45) at the nodes (kh) 1/ρ , k = 1, 2, . . . , N . Tables 1 and 2 show approximate solutions to IVP (45), when A = 0, using our universal P-C algorithm for different values of α, ρ and N at t = 0.5 and t = 2, respectively. In case of α = 1 and ρ = 1, the exact values of the solution to IVP (45), when A = 0, at t = 0.5 and t = 2, are y(0.5) = 0.75601439 and y(2) = 2.35777165, respectively. In this case, from the numerical values displayed in Tables 1 and 2, we can easily see that the numerical solutions provided using our algorithm are definitely highly compatible with the exact solution. For other cases, from the convergence of numerical values presented in Tables 1 and 2, we can realize the numerical stability characteristic of the proposed P-C algorithm. Certainly, the accuracy of the approximate solutions presented using our P-C algorithm can be enhanced when N becomes large.
In order to describe the dynamic behavior of the model given in Eq. (45) according to the parameters α and ρ versus the variable x, we display in Fig. 1 numerical solutions to IVP (45) with T = 2 and N = 1000 computed using our P-C algorithm, when A = 0, for some specific values of α and ρ. y(x) + βy(x) + γ y 3 (x) = 0, is the generalized Caputo-type fractional derivative operator, introduced in Definition 3, β, γ ∈ IR, 1 < α ≤ 2 and σ (x) = exp(λx) such that λ > 0. To simulate numerically IVP (48) over the interval [a, b] = [0, T ], according to our P-C algorithm, we have x 0 = 0, x k+1 = 1 λ ln exp(λx k ) + h (k = 0, 1, . . . , N − 1) and h = exp(λT )−1 N , for some N ∈ IN . Hence, the approximation y k+1 ≈ y(x k+1 ) can be described by the formula with the weights a n,k+1 being determined according to Eq. (41), and Now, our P-C algorithm aims at searching numerical solutions to IVP (48) at the nodes 1 λ ln kh + 1 , k = 1, 2, . . . , N . Tables 3 and 4 show approximate solutions to IVP (48), when β = 2 and γ = 1, using our universal P-C algorithm for different values of α, λ and N at t = 0.5 and t = 2, respectively. As concluded in the previous example, from the convergence of numerical values presented in Tables 3 and 4, we can observe that our P-C algorithm behaves acceptably with respect to its numerical stability.
The dynamic behavior of the model given in Eq. (48) according to the parameters α and λ versus the variable x is described in Fig. 2. This figure displays numerical solutions to IVP (48) with T = 2.5 and N = 1000 computed using our P-C algorithm, when β = 2 and

Concluding remarks
This paper introduces a survey of formations of generalized fractional integrals and derivatives with respect to another function. A novel P-C algorithm to pro-duce numerical solutions to IVPs with generalized Caputo-type fractional differential equations is developed. The proposed algorithm is universal in the sense that it works with any possible choice of the function σ under the prescribed conditions. The obtained attractive numerical results recommend that the proposed P-C algorithm works successfully in handling IVPs with generalized Caputo-type fractional derivatives conveniently, directly and accurately. Because we have many possible options for the function σ within the proposed algorithm, we hope to draw the attention of researchers to consider and employ the proposed P-C algorithm in providing approximate solutions to many nonlinear models including generalized fractional derivatives.