Fractional differential equations and Volterra–Stieltjes integral equations of the second kind

In this paper, we construct a method to find approximate solutions to fractional differential equations involving fractional derivatives with respect to another function. The method is based on an equivalence relation between the fractional differential equation and the Volterra–Stieltjes integral equation of the second kind. The generalized midpoint rule is applied to solve numerically the integral equation and an estimation for the error is given. Results of numerical experiments demonstrate that satisfactory and reliable results could be obtained by the proposed method.

tional derivatives provide an excellent tool for the description of memory and hereditary properties of various materials and processes (Atanackovic and Stankovic 2009;Carpinteri and Mainardi 2014;Demir et al. 2012;Kulish and Lage 2002;Meerschaert 2011;Podlubny 1999;Rezazadeh et al. 2018;Tariq et al. 2018;Vazquez 2005). Several types of fractional derivatives have been suggested to describe more accurately real-world phenomena, each one with their own advantages and disadvantages (Djida et al. 2017;Kilbas et al. 2006;Osman 2017;Osman et al. 2019;Podlubny 1999;Rezazadeh et al. 2019). A more general unifying perspective to the subject was proposed in Agrawal (2010), Klimek and Lupa (2013), Malinowska et al. (2015), by considering fractional operators depending on general kernels. In this work, we follow the special case of this approach that was developed in Almeida (2017a, b), Almeida et al. (2018Almeida et al. ( , 2019, Garra et al. (2019), Kilbas et al. (2006), Yang and Machado (2017). Namely, we focus on nonlinear fractional differential equations involving a Caputo-type fractional derivative with respect to another function, called ψ-Caputo derivative. These types of equations have been successfully used to model the world population growth (Almeida 2017a, b) and gross domestic product (Almeida 2017b). On the other hand, it is well known that fractional differential equations often have to be solved numerically (Arqub and Maayah 2018;Arqub and Al-Smadi 2018a, b;Diethelm 2010;Ford and Connolly 2006;Lubich 1985;Morgado et al. 2013;Sousa and Oliveira 2019). Therefore, knowing the usefulness of the ψ-Caputo fractional derivative, an important issue is to discuss the numerical methods for differential equations with this derivative.
Motivated by the above discussion, in this paper, we provide a numerical scheme to solve fractional differential equations with the ψ-Caputo derivative. The main idea is to rewrite the considered equation as the Volterra-Stieltjes integral equation and then apply the generalized midpoint rule that was developed by Asanov et al. (2011b). Regarding integral equations, the best general reference is the handbook by Polyanin and Manzhirov (2008). Results on nonclassical Volterra integral equations of the first kind can be found in Apartsyn (2003). In Asanov (1998), problems of regularization, uniqueness and existence of solutions for Volterra integral and operator equations of the first kind are studied. Some properties of Voltera and Volterra-Stieltjes integral operators are given in Bukhgeim (1999), and Banas and Regan (2005). In Banas et al. (2000) and Federson and Bianconi (2001), quadratic integral equations of Urysohn-Stieltjes type and their applications are investigated. Various numerical solution methods for integral equations are presented in Asanov et al. (2011a, b), Abdujabbarov (2011), Asanov et al. (2016), Delves andWalsh (1974), Federson et al. (2002). In particular, the generalized trapezoid rule and the generalized midpoint rule to evaluate the Stieltjes integral approximately by employing the notion of derivative of a function by means of a strictly increasing function (Asanov et al. 2011a, b;Asanov 2001), the generalized trapezoid rule for linear Volterra-Stieltjes integral equations of the second kind (Asanov et al. 2016), and the generalized midpoint rule for linear Fredholm-Stieltjes integral equations of the second kind (Asanov and Abdujabbarov 2011).
The paper is organized as follows. First, in Sect. 2, we review some necessary concepts and results on fractional calculus. In Sect. 3, we state an initial value problem with the ψ-Caputo fractional derivative of order α > 1 and prove several results concerning this problem which will be needed in the forthcoming sections. Then, in Sect. 4, using the generalized midpoint rule, we exhibit a numerical procedure to solve the Volterra-Stieltjes integral equation that corresponds to the given fractional differential equation. An upper bound formula for the error in our approximation is derived. Results of numerical experiments are presented in Sect. 5. We verify the accuracy and analyse the stability of the numerical scheme. Section 6 provides conclusions and suggestions for future work and thus completes this work.

Preliminaries on fractional calculus
Let α > 0 be a real number and x : [a, b] → R a function. Given another function ψ ∈ C 1 [a, b] such that ψ is increasing and ψ (t) = 0, for all t ∈ [a, b], the ψ-Riemann-Liouville fractional integral of x, of order α, is defined as where n = [α] + 1. The other important definition, which will be used in this work, is the ψ-Caputo fractional derivative of x of order α: We see at once that, if α = m ∈ N, then the ψ-Caputo fractional derivative coincides with the ordinary derivative Important relations between the two fractional operators are as follows (Almeida et al. 2018): 2. If x ∈ C n−1 [a, b], then

Fractional differential equations
In this section, we consider the following nonlinear fractional differential equation with the ψ-Caputo derivative: subject to the initial conditions where 1. 1 < α / ∈ N and n = [α] + 1, 2. x a and x k a , for k = 1, . . . , n − 1, are fixed reals, Regarding the existence and uniqueness of solutions for nonlinear fractional differential equation with the ψ-Caputo derivative of type (1), we refer the reader to Almeida et al. (2018). For stability results, we suggest the work of Almeida et al. (2019).
We shall prove some preparatory results providing a basis for the later development of a numerical method for solving the initial value problem (1) and (2). First, let us recall that problem (1) and (2) can be rewritten as the Volterra-Stieltjes integral equation, i.e., a Volterra-type integral equation involving the Riemann-Stieltjes integral.

Theorem 1 Almeida et al. (2018) A function x ∈ C n−1 [a, b] is a solution to problem (1) and
(2) if and only if x satisfies the following Volterra-Stieltjes integral equation: where I α,ψ In what follows we assume that: (H1) Function F is Lipschitz with respect to the second variable, that is, there exists a positive constant L such that (H2) There exist nonnegative constants L 1 and β ∈ (0, 1] such that Lemma 1 Let α > 1 and s 1 , Proof The proof follows from the mean value theorem. Let σ (t) = t α , t ∈ [c, d]. Then, For a given function ϕ ∈ C [a, b], let · C denote the usual norm: a, b], and x ∈ C [a, b] be the solution to Eq. (3). Then, under assumption (H1) it holds that where and g is defined by (4).

Theorem 2 Let α > 1 and x ∈ C[a, b] be the solution to Eq. (3). Then, under assumption (H1), it holds that
and L 0 as defined in Lemma 2.

Numerical method and error analysis
The purpose of this section is to construct a method to find approximate solutions to fractional initial value problems of type (1) and (2). We use an equivalence relation between problem (1), (2) and integral equation (3). Then, the approximation routine is based on the generalized midpoint rule that was first developed by Asanov et al. (2011b). They proposed an approximation of the Stieltjes integral with the use of the notion of the derivative of a function with respect to the strictly increasing function (see Asanov 2001). The generalized midpoint rule summarizes the midpoint rule (Kalitkin 1978), and here we use this method to solve numerically the integral equation (3).
Theorem 4 Let α > 2 and β = 1. Under assumptions (H1) and (H2) it holds that as n → ∞, where Proof Let the error be denoted by Taking into account (21) and (22), we have From this, by (18), we conclude that where Now, let us consider the system and 1 = R(h) as an initial condition. It is easily seen that |z k | ≤ k for k = 1, . . . , 2n. Indeed, this can be verified by mathematical induction as follows: for k = 1 it is obvious. Let |z j | ≤ j for j = 1, . . . , k − 1. Then, by inequality (25), we get Observe that j , j = 1, . . . , 2n, given by satisfy system (27). In fact, using the inequality and taking γ = c 5 h we obtain Consequently, we get the following estimates: Using the fact that (1 + t) 1 t is increasing and approaches the number e as t → 0 + we get This together with (26) and (28) gives (23).

Under assumptions (H1) and (H2) it holds that
as n → ∞, where Proof Using the conditions of Theorem 5 and the estimates (20), from (24) we obtain following inequalities for z k : Set z = sup k=0,...,2n |z k |. Then, by (30), we have what finishes the proof.

Numerical examples
In the following examples, to show the efficiency of the proposed numerical method, we approximate the solutions for some fractional differential equations of order α > 1. In all examples, corresponding systems of algebraic equations of type (22) are solved using the command fsolve in Maple. In tables, we present the absolute error and the elapsed CPU time in seconds of the proposed numerical scheme (22).
Example 1 Consider the following fractional differential equation: with the initial conditions and kernel respectively. Therefore, in integral equation (3) we have α = 5 2 and In this case, conditions (H1) and (H2) are satisfied for L = 9/8, L 1 = 0, and β = 1. Note that, the exact solution to (31) and (32) is x(t) = 1, t ∈ [0, 1]. The exact and numerical solutions for n = 20, 40, 60, as well as the evolution of the error, are shown in Fig. 1. The maximum of absolute error for different values of n and the elapsed CPU time in seconds are displayed in Table 1. As expected, as n increases, the error decreases.
Example 2 Consider the following fractional differential equation: with the initial conditions and kernel respectively. Therefore, in integral equation (3) we have α = 3 2 and     Table 2.
Example 3 Consider the following fractional differential equation:  with the initial conditions and kernel respectively. Therefore, in integral equation (3) we have α = 3 2 and (0)) 3/2 , where E 3/2 (·) denotes the Mittag-Leffler function of order α = 3/2, we conclude that the solution to this problem is x(t) = E 3/2 (ψ(t) − ψ(0)) 3/2 − 1, t ∈ [0, 5]. The exact and numerical solutions for n = 20, 40, 60, as well as the evolution of the error, are shown in Fig. 3. The maximum of absolute error for different values of n and the elapsed CPU time in seconds are displayed in Table 3. In the next examples, we analyse the stability of the proposed numerical method.
Example 4 Consider the two fractional differential equations with initial conditions: The kernel is given by ψ(t) = √ t + 1 and μ is a real number. Let x and x be solutions of (35) and (36), respectively. In this case, F(t, x) = x for both equations and function g is given by respectively. By the Grönwall inequality (cf. Almeida et al. 2019; Sousa and Oliveira 2019), we conclude that Consider the two functions In Fig. 4, we present the plots for two values of μ. It can be observed that the proposed numerical scheme preserves the underlying structural stability of the initial value problem, with respect to small perturbation of the initial conditions. In Table 4, we display approximate values of x(5) for different values of μ and n. Example 5 Let us now consider the following fractional differential equation: with initial conditions where μ 1 , μ 2 , μ 3 ∈ R are three parameters. Let ψ(t) = √ t + 1. In this case, we have To analyse the stability of the numerical method, we consider perturbations of parameters μ 1 , μ 2 , μ 3 . First, we fix μ 2 = 0 = μ 3 and we analyse the numerical solutions of the two following equations:  Table 6 Numerical results for equations (40) and (41), with n = 40, in Example 5  Table 7 Numerical results for Eqs. (42) and (43), with n = 40, in Example 5 for μ 1 ∈ {0.8, 0.9, 1, 1.1, 1.2}. Let x and x be solutions of (38) and (39), respectively. Figure 5 (left) shows x for different values of μ 1 and n = 40. In Table 5, we display max i |x(t i ) − x(t i )|. Next, we assume that μ 1 = 1, μ 3 = 0, and consider numerical solutions to equations: for μ 2 ∈ {0, 0.1, 0.2, 0.3, 0.4}. Let x and x be solutions of (40) and (41), respectively. Numerical results, for different values of μ 2 , are presented in Fig. 5 (center) and Table 6. Finally, we fix μ 1 = 1, μ 2 = 0, and analyse numerical solutions to equations: for μ 3 ∈ {0, 0.1, 0.2, 0.3, 0.4}. Let x and x be solutions of (42) and (43), respectively. Numerical results, for different values of μ 3 , are presented in Fig. 5 (right) and Table 7.

Conclusions
Recently, the Caputo derivative with respect to a kernel function ψ was proposed and applied to some real-world processes (Almeida 2017a, b; Voyiadjis and Sumelka 2019). As usual, an important issue is to develop numerical methods for fractional differential equations with this new type of derivative. In this paper, such a numerical scheme is presented and its error bound is investigated. The idea is based on an equivalence relation between the fractional differential equation with ψ-Caputo derivative and the Volterra-Stieltjes integral equation.
To the latter equation, the generalized midpoint rule is applied. The numerical simulations show that the proposed numerical scheme preserves the underlying structural stability of the initial value problem, with respect to small perturbation of the initial data, and satisfactory and reliable results could be obtained. It means that the approximation routine presented in (Asanov et al. 2011b) for Stieltjes integral can also be successfully applied to solve fractional differential equations, where the fractional derivative operator depends on an increasing function. Nevertheless, as mentioned in Introduction, there exist various numerical methods for integral equations. Therefore, we expect that some of those methods could be adopted for fractional differential equations with ψ-Caputo derivative. This important issue will be considered in a forthcoming paper.
Yang XJ, Machado JA (2017) A new fractional operator of variable order: application in the description of anomalous diffusion. Phys A 481:276-283 Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.