1 Introduction

Fractional differential equations are a generalization of ordinary differential equations, where integer-order derivatives are replaced by fractional derivatives. Over the last decade, these equations have attracted a lot of attention from researchers from different areas, since fractional 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. (2018, 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 \(\psi \)-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 \(\psi \)-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 \(\psi \)-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), Asanov and Abdujabbarov (2011), Asanov et al. (2016), Delves and Walsh (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 \(\psi \)-Caputo fractional derivative of order \(\alpha >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.

2 Preliminaries on fractional calculus

Let \(\alpha >0\) be a real number and \(x:[a,b]\rightarrow \mathbb R\) a function. Given another function \(\psi \in C^1[a,b]\) such that \(\psi \) is increasing and \(\psi '(t)\not =0\), for all \(t\in [a,b]\), the \(\psi \)-Riemann–Liouville fractional integral of x, of order \(\alpha \), is defined as

$$\begin{aligned} {I_{a+}^{\alpha ,\psi }} x(t)=\frac{1}{\varGamma (\alpha )}\int _a^t \psi '(\tau )(\psi (t)-\psi (\tau ))^{\alpha -1}x(\tau )\,\mathrm{d}\tau , \end{aligned}$$

where \(\varGamma (\cdot )\) is the Gamma function, i.e., \(\varGamma (z)=\int _0^\infty t^{z-1}e^{-t}\mathrm{d}t\), \(z>0\). The \(\psi \)-Riemann–Liouville fractional derivative of x, of order \(\alpha \), is defined as

$$\begin{aligned} {D_{a+}^{\alpha ,\psi }} x(t)= \left( \frac{1}{\psi '(t)}\frac{\mathrm{d}}{\mathrm{d}t}\right) ^n {I_{a+}^{n-\alpha ,\psi }}x(t), \end{aligned}$$

where \(n=[\alpha ]+1\). The other important definition, which will be used in this work, is the \(\psi \)-Caputo fractional derivative of x of order \(\alpha \):

$$\begin{aligned} {^CD_{a+}^{\alpha ,\psi }} x(t)={D_{a+}^{\alpha ,\psi }} \left[ x(t)-\sum _{k=0}^{n-1}\frac{x^{[k]}_\psi (a)}{k!}(\psi (t)-\psi (a))^k\right] , \end{aligned}$$

where

$$\begin{aligned} n=[\alpha ]+1 \quad \text{ for } \, \alpha \notin \mathbb N, \quad n=\alpha \, \text{ otherwise }, \end{aligned}$$

and

$$\begin{aligned} x^{[k]}_\psi (t)=\left( \frac{1}{\psi '(t)}\frac{\mathrm{d}}{\mathrm{d}t}\right) ^kx(t). \end{aligned}$$

We see at once that, if \(\alpha =m\in \mathbb N\), then the \(\psi \)-Caputo fractional derivative coincides with the ordinary derivative

$$\begin{aligned} {^CD_{a+}^{m,\psi }} x(t)=x^{[m]}_\psi (t), \end{aligned}$$

while if \(\alpha \in \mathbb {N}\) and \(x \in C^n[a,b]\), then

$$\begin{aligned} {^CD_{a+}^{\alpha ,\psi }} x(t)=\frac{1}{\varGamma (n-\alpha )}\int _a^t \psi '(\tau )\left( \psi (t)-\psi (\tau )\right) ^{n-\alpha -1}x^{[n]}_\psi (\tau )\,\mathrm{d}\tau . \end{aligned}$$

Important relations between the two fractional operators are as follows (Almeida et al. 2018):

  1. 1.

    If \(x\in C[a,b]\), then

    $$\begin{aligned} {^CD_{a+}^{\alpha ,\psi }} {I_{a+}^{\alpha ,\psi }} x(t)=x(t). \end{aligned}$$
  2. 2.

    If \(x\in C^{n-1}[a,b]\), then

    $$\begin{aligned} {I_{a+}^{\alpha ,\psi }} {^CD_{a+}^{\alpha ,\psi }} x(t)=x(t)-\sum _{k=0}^{n-1}\frac{x^{[k]}_\psi (a)}{k!}\left( \psi (t)-\psi (a)\right) ^k. \end{aligned}$$

3 Fractional differential equations

In this section, we consider the following nonlinear fractional differential equation with the \(\psi \)-Caputo derivative:

$$\begin{aligned} {^CD_{a+}^{\alpha ,\psi }}x(t)=F(t,x(t))+f(t), \quad t\in [a,b], \end{aligned}$$
(1)

subject to the initial conditions

$$\begin{aligned} x(a)=x_{a}, \quad x_{\psi }^{[k]}(a)=x_{a}^{k},\quad k=1,\ldots ,n-1, \end{aligned}$$
(2)

where

  1. 1.

    \(1<\alpha \notin \mathbb N\) and \(n=[\alpha ]+1\),

  2. 2.

    \(x_{a}\) and \(x_{a}^{k}\), for \( k=1,\ldots ,n-1\), are fixed reals,

  3. 3.

    \(F:[a,b]\times \mathbb R\rightarrow \mathbb R\) is continuous and \(F(t,0)=0\), for all \(t\in [a,b]\),

  4. 4.

    \(f:[a,b]\rightarrow \mathbb R\) is continuous.

Regarding the existence and uniqueness of solutions for nonlinear fractional differential equation with the \(\psi \)-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\in 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:

$$\begin{aligned} x(t)=I_{a+}^{\alpha ,\psi }F(t,x(t))+g(t), \quad t\in [a,b], \end{aligned}$$
(3)

where

$$\begin{aligned} I_{a+}^{\alpha ,\psi }F(t,x(t))= & {} \frac{1}{\varGamma (\alpha )}\int _{a}^{t}{(\psi (t)-\psi (\tau ))^{\alpha -1}F(\tau ,x(\tau ))}\,\mathrm{d}\psi (\tau ),\nonumber \\ g(t)= & {} \sum _{k=0}^{n-1}\frac{x_{a}^{k}}{k!}\left( \psi (t)-\psi (a)\right) ^{k}+I_{a+}^{\alpha ,\psi }f(t), \quad t\in [a,b]. \end{aligned}$$
(4)

In what follows we assume that:

  1. (H1)

    Function F is Lipschitz with respect to the second variable, that is, there exists a positive constant L such that

    $$\begin{aligned} \vert F(t,x_{1})-F(t,x_{2})\vert \le L\vert x_{1}- x_{2}\vert , \quad \forall t\in [a,b],\quad \forall x_{1},x_{2}\in \mathbb R. \end{aligned}$$
  2. (H2)

    There exist nonnegative constants \(L_{1}\) and \(\beta \in (0,1]\) such that

    $$\begin{aligned} \vert F(t_{1},x)-F(t_{2},x)\vert \le L_{1}\vert t_{1}- t_{2}\vert ^{\beta } \vert x\vert , \quad \forall t_{1}, t_{2}\in [a,b],\quad \forall x\in \mathbb R. \end{aligned}$$

Lemma 1

Let \(\alpha >1\) and \(s_{1}, s_{2}\in [c,d]\), where \(0\le c<d\). Then, \(\vert s_{1}^{\alpha }-s_{2}^{\alpha }\vert \le c_{0}\vert s_{1}-s_{2}\vert ,\) where \(c_{0}=\alpha d^{\alpha -1}.\)

Proof

The proof follows from the mean value theorem. Let \(\sigma (t)=t^\alpha \), \(t\in [c,d]\). Then,

$$\begin{aligned} \vert s_{1}^{\alpha }-s_{2}^{\alpha }\vert \le \sup _{t\in (c,d)}|\sigma '(t)|\cdot \vert s_{1}-s_{2}\vert \le \alpha d^{\alpha -1}\vert s_{1}-s_{2}\vert . \end{aligned}$$

\(\square \)

For a given function \(\varphi \in C[a,b]\), let \(\Vert \cdot \Vert _C\) denote the usual norm:

$$\begin{aligned} \Vert \varphi \Vert _C=\sup _{t\in [a,b]}\vert \varphi (t)\vert . \end{aligned}$$

Lemma 2

Let \(\alpha >1\), \(f\in C[ a,b]\), and \(x\in C[a,b]\) be the solution to Eq. (3). Then, under assumption (H1) it holds that

$$\begin{aligned} \Vert x\Vert _{C}\le L_{0}\Vert g\Vert _{C}, \end{aligned}$$
(5)

where

$$\begin{aligned} L_{0}=\exp \left\{ \frac{L(b-a)}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert \psi '\Vert _{C}\right\} , \end{aligned}$$

and g is defined by (4).

Proof

Since \(F(t,0)=0\), it follows that

$$\begin{aligned} | F(\tau ,x(\tau ))|= | F(\tau ,x(\tau ))-F(\tau ,0)|\le L | x(\tau )|. \end{aligned}$$

Hence, by (3), we obtain

$$\begin{aligned} \vert x(t)\vert \le \frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert \psi '\Vert _{C}\int _{a}^{t}\vert x(s)\vert \mathrm{d}s+\Vert g\Vert _{C}, \quad t\in [a,b]. \end{aligned}$$

Applying the Grönwall inequality we get (5). \(\square \)

Lemma 3

Let \(\alpha \in (1,2),\) \(s_{1},s_{2}\in [c,d],\) where \(0\le c<d\), \(p>1\) such that \(p(2-\alpha )<1\), and \(\frac{1}{p}+\frac{1}{q}=1\). Then,

$$\begin{aligned} \vert s_{1}^{\alpha -1}-s_{2}^{\alpha -1}\vert \le c_{1}\vert s_{1}-s_{2}\vert ^{\frac{1}{q}}, \end{aligned}$$

where

$$\begin{aligned} c_{1}=(\alpha -1)\left( \frac{1}{1-p(2-\alpha )}\right) ^{\frac{1}{p}}\left( d^{1-p(2-\alpha )}-c^{1-p(2-\alpha )}\right) ^{\frac{1}{p}}. \end{aligned}$$

Proof

Let \(s_{1},s_{2}\in [ c,d],\) with \(s_{2}> s_{1}.\) Then, by the Hölder inequality,

$$\begin{aligned} s_{2}^{\alpha -1}-s_{1}^{\alpha -1}= & {} (\alpha -1)\int _{s_{1}}^{s_{2}}\frac{1}{\tau ^{2-\alpha }}\mathrm{d}\tau \le (\alpha -1)\left( \int _{s_{1}}^{s_{2}}{\tau ^{-p(2-\alpha )}\mathrm{d}\tau }\right) ^{\frac{1}{p}}\left( \int _{s_{1}}^{s_{2}}{\mathrm{d}\tau }\right) ^{\frac{1}{q}}\\= & {} c_{1}\vert s_{2}-s_{1}\vert ^{\frac{1}{q}}. \end{aligned}$$

\(\square \)

Theorem 2

Let \(\alpha >1\) and \(x\in C[a,b]\) be the solution to Eq. (3). Then, under assumption (H1), it holds that

$$\begin{aligned} \vert x(t_{2})-x(t_{1})\vert \le c_{2}\vert t_{2}- t_{1}\vert , \quad \forall t_{1}, t_{2}\in [a,b], \end{aligned}$$

where

$$\begin{aligned} c_{2}= & {} \frac{\Vert \psi '\Vert _{C}}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}[L_{0}L\Vert g\Vert _{C}+\Vert f\Vert _{C}]\\&+\sum _{k=1}^{n-1}{\frac{\vert x_{a}^{k}\vert }{(k-1)!}\left( \psi (b)-\psi (a)\right) ^{k-1}\Vert \psi '\Vert _{C}}, \end{aligned}$$

and \(L_{0}\) as defined in Lemma 2.

Proof

Let \(\alpha >1\) and \(t_{2}>t_{1}\), where \(t_{1},t_{2}\in [a,b]\). Then, by Eqs. (3) and (4) and Lemma 2, we have

$$\begin{aligned} \vert x(t_{2})-x(t_{1})\vert&\le \frac{L_{0}L\Vert g\Vert _{C}}{\varGamma (\alpha )}\int _{a}^{t_{1}}{\left[ (\psi (t_{2})-\psi (\tau ))^{\alpha -1}-(\psi (t_{1})-\psi (\tau ))^{\alpha -1}\right] }\mathrm{d}\psi (\tau )\\&\quad +\frac{LL_{0}\Vert g\Vert _{C}}{\varGamma (\alpha )}\int _{t_{1}}^{t_{2}}{(\psi (t_{2})-\psi (\tau ))^{\alpha -1}}\mathrm{d}\psi (\tau )\\&\quad +\sum _{k=1}^{n-1}{\frac{|x_{a}^{k}|}{k!}\left[ (\psi ( t_{2})-\psi (a))^{k}-(\psi ( t_{1})-\psi (a))^{k}\right] }\\&\quad +\frac{\Vert f\Vert _{C}}{\varGamma (\alpha )}\int _{a}^{t_{1}}{\left[ (\psi ( t_{2})-\psi (\tau ))^{\alpha -1}-(\psi ( t_{1})-\psi (\tau ))^{\alpha -1}\right] \mathrm{d}\psi (\tau )}\\&\quad +\frac{\Vert f\Vert _{C}}{\varGamma (\alpha )}\int _{t_{1}}^{t_{2}}{(\psi (t_{2})-\psi (\tau ))^{\alpha -1}\mathrm{d}\psi (\tau )}. \end{aligned}$$

Computing the integrals and applying the mean value theorem, we get

$$\begin{aligned} \vert x(t_{2})-x(t_{1})\vert&\le \frac{L_{0}L\Vert g\Vert _{C}}{\alpha \varGamma (\alpha )}\left[ (\psi (t_{2})-\psi (a))^{\alpha }-(\psi (t_{1})-\psi (a))^{\alpha }-(\psi ( t_{2})-\psi ( t_{1}))^{\alpha }\right. \\&\quad \left. +(\psi (t_{2})-\psi (t_{1}))^{\alpha }\right] +\frac{\Vert f\Vert _{C}}{\alpha \varGamma (\alpha )}\left[ (\psi (t_{2})-\psi (a))^{\alpha }\right. \\&\quad \left. -(\psi ( t_{1})-\psi (a))^{\alpha }-(\psi ( t_{2})-\psi (t_{1}))^{\alpha }+(\psi ( t_{2})-\psi ( t_{1}))^{\alpha }\right] \\&\quad +\sum _{k=1}^{n-1}{\frac{\vert x_{a}^{k}\vert }{(k-1)!}}(\psi (b)-\psi (a))^{k-1}\Vert \psi '\Vert _{C}(t_{2}-t_{1}). \end{aligned}$$

Finally, by Lemma 1 (with \(d=\psi (b)-\psi (a)\)) and the mean value theorem, we obtain

$$\begin{aligned} \vert x(t_{2})-x(t_{1})\vert&\le \frac{L_{0}L\Vert g\Vert _{C}}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert \psi '\Vert _{C}(t_{2}-t_{1})\\&\quad +\frac{\Vert f\Vert _{C}}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert \psi '\Vert _C(t_{2}-t_{1})\\&\quad +\sum _{k=1}^{n-1}{\frac{\vert x_{a}^{k}\vert }{(k-1)!}}(\psi (b)-\psi (a))^{k-1}\Vert \psi '\Vert _C(t_{2}-t_{1}). \end{aligned}$$

\(\square \)

Theorem 3

Let \(x\in C[a,b] \) be the solution to Eq. (3), \(G=\lbrace (t,\tau ):a\le \tau \le t\le b\rbrace \) and

$$\begin{aligned} k(t,\tau )=\frac{1}{\varGamma (\alpha )}\left( \psi (t)-\psi (\tau )\right) ^{\alpha -1}F(\tau ,x(\tau )). \end{aligned}$$
(6)

Under assumptions (H1) and (H2) the following hold:

  1. 1.

    If \(\alpha >2\) and \(\beta =1\), then

    $$\begin{aligned} \vert k(t,\tau _{1})-k(t,\tau _{2})\vert \le c_{3}\vert \tau _{1}-\tau _{2}\vert , \end{aligned}$$
    (7)

    for all \((t,\tau _{1}),(t,\tau _{2})\in G\), where

    $$\begin{aligned} c_{3}= & {} \frac{1}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -2}\\&\times \left[ (L_{0}L_{1}\Vert g\Vert _{C}+Lc_{2})(\psi (b)-\psi (a))+(\alpha -1)LL_{0}\Vert g\Vert _{C}\Vert \psi '\Vert _{C}\right] . \end{aligned}$$
  2. 2.

    If \(\alpha \in (1,2), p>1\) with \(p(2-\alpha )<1,\frac{1}{p}+\frac{1}{q}=1\), and \(\gamma =\min \lbrace \beta ,\frac{1}{q}\rbrace \), then

    $$\begin{aligned} \vert k(t,\tau _{1})-k(t,\tau _{2})\vert \le c_{4}\vert \tau _{1}-\tau _{2}\vert ^{\gamma }, \end{aligned}$$
    (8)

    for all \(( t,\tau _{1}),(t,\tau _{2})\in G\), where

    • for \(\beta \ge \frac{1}{q}\)

      $$\begin{aligned} c_{4}= & {} \frac{1}{\varGamma (\alpha )}\Big [L_{0}L_{1}(\psi (b)-\psi (a))^{\alpha -1}(b-a)^{\beta -\frac{1}{q}}\Vert g\Vert _{C}\\&+\,Lc_{2}(\psi (b)-\psi (a))^{\alpha -1}(b-a)^{1-\frac{1}{q}}\\&+ \, LL_{0}(\alpha -1)\left( \frac{1}{1-p(2-\alpha )}\right) ^{\frac{1}{p}}(\psi (b)-\psi (a))^{\frac{1}{p}(1-p(2-\alpha ))}\Vert g\Vert _{C}\Vert \psi '(t)\Vert _{c}\Big ]; \end{aligned}$$

    • for \(\beta <\frac{1}{q}\)

      $$\begin{aligned} c_{4}= & {} \frac{1}{\varGamma (\alpha )}\left[ L_{0}L_{1}(\psi (b)-\psi (a))^{\alpha -1}\Vert g\Vert _{C}\right. \\&\left. +Lc_{2}(\psi (b)-\psi (a))^{\alpha -1}(b-a)^{1-\beta }+LL_{0}(\alpha -1)\left( \frac{1}{1-p(2-\alpha )}\right) ^{\frac{1}{p}}\right. \\&\left. \times (\psi (b)-\psi (a))^{\frac{1}{p}(1-p(2-\alpha ))} \Vert g\Vert _{C}\Vert \psi '\Vert _{C}(b-a)^{\frac{1}{q}-\beta }\right] . \end{aligned}$$

Proof

Observe that, for all \((t,\tau _{1}),(t,\tau _{2})\in G\), we have

$$\begin{aligned} k(t,\tau _{2})-k(t,\tau _{1})= & {} \frac{1}{\varGamma (\alpha )}\left[ (\psi (t)-\psi (\tau _{2}))^{\alpha -1} -(\psi (t)-\psi (\tau _{1}))^{\alpha -1}\right] F(\tau _{2},x(\tau _{2}))\nonumber \\&+\frac{1}{\varGamma (\alpha )}(\psi (t)-\psi (\tau _{1}))^{\alpha -1}\left[ F(\tau _{2},x(\tau _{2}))-F(\tau _{1},x( \tau _{2}))\right] \nonumber \\&+\frac{1}{\varGamma (\alpha )}(\psi (t)-\psi (\tau _{1}))^{\alpha -1}\left[ F(\tau _{1},x(\tau _{2}))-F(\tau _{1},x(\tau _{1}))\right] . \end{aligned}$$
(9)

Applying Lemmas 1 and 2, and Theorem 2, we can assert that

$$\begin{aligned}&\left| \frac{1}{\varGamma (\alpha )}(\psi (t)-\psi (\tau _{1}))^{\alpha -1} [F(\tau _{2},x(\tau _{2}))-F(\tau _{1},x(\tau _{2}))] \right| \nonumber \\&\quad \le \frac{L_{0}L_{1}}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert g\Vert _{C}| \tau _{1}-\tau _{2}|^{\beta }, \end{aligned}$$
(10)
$$\begin{aligned}&\left| \frac{1}{\varGamma (\alpha )}(\psi (t)-\psi (\tau _{1}))^{\alpha -1}[F(\tau _{1},x(\tau _{2})) -F(\tau _{1},x(\tau _{1}))] \right| \nonumber \\&\quad \le \frac{L c_{2}}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}|\tau _{1}-\tau _{2}|, \end{aligned}$$
(11)

and

$$\begin{aligned}&\left| \frac{1}{\varGamma (\alpha )}[ (\psi (t)-\psi (\tau _{2}))^{\alpha -1}-(\psi (t)-\psi (\tau _{1}))^{\alpha -1}] F(\tau _{2},x(\tau _{2}))\right| \nonumber \\&\quad \le \frac{LL_{0}(\alpha -1)}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -2}\Vert \psi '\Vert _{C}\Vert g\Vert _{C}| \tau _{1}-\tau _{2}|, \quad \alpha >2. \end{aligned}$$
(12)

If \(\alpha >2\) and \(\beta =1,\) then taking into account (10)–(12), from (9) we obtain (7). On the other hand, if \(\alpha \in (1,2), p>1, p(2-\alpha )<1,\frac{1}{p}+\frac{1}{q}=1,\) then by Lemmas 2 and 3, and Theorem 2 it follows that

$$\begin{aligned}&\left| \frac{1}{\varGamma (\alpha )}\Big [(\psi (t)-\psi (\tau _{1}))^{\alpha -1}-(\psi (t)-\psi (\tau _{2}))^{\alpha -1}\Big ] F(\tau _{2},x(\tau _{2}))\right| \nonumber \\&\quad \le \frac{LL_{0}}{\varGamma (\alpha )}(\alpha -1)\left( \frac{1}{1-p(2-\alpha )}\right) ^{\frac{1}{p}}(\psi (b)-\psi (a))^{\frac{1}{p}(1-p(2-\alpha ))} \Vert g\Vert _{C}\Vert \psi '\Vert _{C}|\tau _{1}-\tau _{2}|^{\frac{1}{q}}.\nonumber \\ \end{aligned}$$
(13)

In this case, taking into account (10)–(11) and (13), from (9) we obtain (8). \(\square \)

4 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).

For \(n\in \mathbb N\), let

$$\begin{aligned} h=\frac{b-a}{2n},\quad t_{k}=a+kh, \quad k=0,1,2,\ldots ,2n. \end{aligned}$$

We substitute \(t=t_{k}\), \(k=0,1,2,\ldots ,2n\) into integral equation (3) and examine the following system of equations:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle x(t_{0})= g(t_{0}), \quad t_{0}=a \\ \displaystyle x(t_{1})=\int _{a}^{ t_{1}}{ k(t_{1}, \tau )\mathrm{d}\psi (\tau )+g (t_{1})} \\ \displaystyle x(t_{2i})=\int _{a}^{ t_{2i}}{ k(t_{2i} , \tau )\mathrm{d}\psi (\tau )+g (t_{2i}), \quad i=1,\ldots ,n}\\ \displaystyle x(t_{2j+1})=\int _{a}^{ t_{2j+1}}{ k(t_{2j+1}, \tau )\mathrm{d}\psi (\tau )+g(t_{2j+1}), } \quad j=1,\ldots ,n-1, \end{array}\right. \end{aligned}$$
(14)

where \(k(t,\tau )\) is defined by (6). Using the generalized midpoint rule for integrals in system (14) we get

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \int _{a}^{ t_{2i}}{ k(t_{2i},\tau )\mathrm{d}\psi (\tau )=\sum _{m=1}^{i}k(t_{2i},t_{2m-1})[\psi (t_{2m})-\psi (t_{2m-2})]}\\ \displaystyle \displaystyle \quad +\sum _{m=1}^{i}{R_{m}^{(2n)}(x,i),}\quad i=1,\ldots ,n \\ \displaystyle \int _{a}^{t_{1}}{k(t_{1},\tau )\mathrm{d}\psi (\tau )=k( t_{1}, t_{0})[\psi ( t_{1})-\psi (t_{0})] +R_{0}^{(0)}(x)} \\ \displaystyle \int _{a}^{ t_{2j+1}}{ k(t_{2j+1}, \tau )\mathrm{d}\psi (\tau )=\int _{a}^{ t_{1}}{k(t_{2j+1}, s)\mathrm{d}\psi (s)}+}\int _{ t_{1}}^{ t_{2j+1}}{ k(t_{2j+1}, \tau )\mathrm{d}\psi (\tau )} \\ \displaystyle = k(t_{2j+1}, t_{0})[ \psi ( t_{1})-\psi ( t_{0})]+R_{j}^{(0)}(x)+\sum _{m=1}^{j}{ k(t_{2j+1}, } t_{2m})\\ \displaystyle \quad \times [ \psi ( t_{2m+1})-\psi ( t_{2m-1})]+\sum _{m=1}^{j}{R_{m}^{(2n-1)}(x,j), } \quad j=1,\ldots ,n-1 ,\\ \end{array}\right. \end{aligned}$$
(15)

where

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle R_{0}^{(0)}(x)=\int _{a}^{ t_{1}}{ [k(t_{1}, \tau )-k( t_{1},t_{0})]\mathrm{d}\psi (\tau )} \\ \displaystyle R_{j}^{(0)}(x)=\int _{a}^{ t_{1}}{ [k(t_{2j+1}, \tau )-k( t_{2j+1},t_{0})]\mathrm{d}\psi (\tau ), \quad j=1,\ldots ,n-1} \\ \displaystyle R_{m}^{(2n)}(x,i)=\int _{ t_{2m-2}}^{ t_{2m}}{ [k(t_{2i}, \tau )-k(t_{2i}, t_{2m-1})]\mathrm{d}\psi (\tau ),\quad m=1,\ldots ,i} \\ \displaystyle R_{m}^{(2n-1)}(x,j)=\int _{ t_{2m-1}}^{ t_{2m+1}}{ [k(t_{2j+1}, \tau )-k(t_{2j+1}, t_{2m})]\mathrm{d}\psi (\tau ), \quad m=1,\ldots ,j.} \\ \end{array}\right. \end{aligned}$$
(16)

If \(\alpha >2\) and \(\beta =1\), then on the strength of Theorem 3 for \(R_{m}^{(2n)}(x,i)\) and \(R_{m}^{(2n-1)}(x,j)\) we obtain the following estimates:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \vert R_{m}^{(2n)}(x,i)\vert \le c_{3}h[ \psi ( t_{2m})-\psi (t_{2m-2})] , \quad m=1,\ldots ,i \\ \displaystyle \vert R_{m}^{(2n-1)}(x,j)\vert \le c_{3}h[ \psi ( t_{2m+1})-\psi (t_{2m-1})] , \quad m=1,\ldots ,j, \\ \end{array}\right. \end{aligned}$$
(17)

where \(i=1,2,\ldots ,n, j=1,2,\ldots ,n-1\) (cf. Corollary 2 in Asanov et al. 2011b). Then taking into account (16) and (17), we have

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \left| \sum _{m=1}^{i}{R_{m}^{(2n)}(x,i) }\right| \le c_{3}h[ \psi (t_{2i})-\psi (a)] , \quad i=1,\ldots ,n\\ \displaystyle \vert R_{0}^{(0)}(x)\vert \le c_{3}h[ \psi ( t_{1})-\psi (a)]\\ \displaystyle \vert R_{j}^{(0)}(x)\vert \le c_{3}h[ \psi ( t_{1})-\psi (a)] \\ \displaystyle \left| \sum _{m=1}^{j}{R_{m}^{(2n-1)}(x,j) }\right| \le c_{3}h[ \psi (t_{2j+1})-\psi ( t_{1})] , \quad j=1,\ldots ,n-1. \\ \end{array}\right. \end{aligned}$$
(18)

If \(\alpha \in (1,2)\), \(p>1\) with \(p(2-\alpha )<1\), \(\frac{1}{p}+\frac{1}{q}=1\), \(\gamma =\min \lbrace \beta \) , \(\frac{1}{q}\rbrace \), then on the strength of Theorem 3, for \(R_{m}^{(2n)}(x,i)\) and \(R_{m}^{(2n-1)}(x,j)\), we obtain the following estimates:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \vert R_{m}^{(2n)}(x,i)\vert \le c_{4}h^{\gamma }[ \psi (t_{2m})-\psi ( t_{2m-2})] , \quad m=1,\ldots ,i \\ \displaystyle \vert R_{m}^{(2n-1)}(x,j)\vert \le c_{4}h^{\gamma }[ \psi (t_{2m+1})-\psi ( t_{2m-1})] , \quad m=1,\ldots ,j ,\\ \end{array}\right. \end{aligned}$$
(19)

for \(i=1,2,\ldots ,n, j=1,2,\ldots ,n-1\) (cf. Corollary 2 in Asanov et al. 2011b). Then, taking into account (16) and (19), we have

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \left| \sum _{m=1}^{i}{R_{m}^{(2n)}(x,i) }\right| \le c_{4}h^{\gamma }[ \psi ( t_{2i})-\psi (a)] ,\quad i=1,\ldots ,n \\ \displaystyle \vert R_{0}^{(0)}(x)\vert \le c_{4}h^{\gamma }[ \psi ( t_{1})-\psi (a)]\\ \displaystyle \vert R_{j}^{(0)}(x)\vert \le c_{4}h^{\gamma }[ \psi ( t_{1})-\psi (a)]\\ \displaystyle \left| \sum _{m=1}^{j}{R_{m}^{(2n-1)}(x,j) }\right| \le c_{4}h^{\gamma }[ \psi ( t_{2j+1})-\psi ( t_{1})] , \quad j=1,\ldots ,n-1. \\ \end{array}\right. \end{aligned}$$
(20)

Combining relations (6) and (15) we can rewrite (14) as

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle x(t_{0})= g(t_{0}), \quad t_{0}=a \\ \displaystyle x(t_{1})=\frac{1}{\varGamma (\alpha )}(\psi ( t_{1})-\psi ( t_{0}))^{\alpha -1}F( t_{0},x( t_{0}))[ \psi ( t_{1})-\psi ( t_{0})] +g(t_{1})+R_{0}^{(0)}(x)\\ \displaystyle x(t_{2i})=\frac{1}{\varGamma (\alpha )}\sum _{m=1}^{i}{(\psi ( t_{2i})-\psi ( t_{2m-1}))^{\alpha -1}}F( t_{2m-1},x( t_{2m-1}))\\ \displaystyle \quad \times [ \psi ( t_{2m})-\psi ( t_{2m-2})] +g(t_{2i})+\sum _{m=1}^{i}{R_{m}^{(2n)}(x,i)}, \quad i=1,\ldots ,n\\ \displaystyle x(t_{2j+1})=\frac{1}{\varGamma (\alpha )}(\psi ( t_{2j+1})-\psi (t_{0}))^{\alpha -1}F( t_{0},x( t_{0}))(\psi ( t_{1})-\psi ( t_{0}))\\ \displaystyle \quad +\frac{1}{\varGamma (\alpha )}\sum _{m=1}^{j}{(\psi ( t_{2j+1})-\psi (t_{2m}))^{\alpha -1}}F( t_{2m},x( t_{2m}))[ \psi ( t_{2m+1})-\psi ( t_{2m-1})] \\ \displaystyle \quad +g (t_{2j+1})+R_{j}^{(0)}(x)+\sum _{m=1}^{j}{R_{m}^{(2n-1)}(x,j),} \quad j=1,\ldots ,n-1. \\ \end{array}\right. \end{aligned}$$
(21)

Omitting

$$\begin{aligned} \sum _{m=1}^{i}{R_{m}^{(2n)}(x,i)},\, R_{0}^{(0)}(x),\, R_{j}^{(0)}(x),\, \sum _{m=1}^{j}{R_{m}^{(2n-1)}(x)} \end{aligned}$$

in equations of system (21) and writing the sought solution x at the nodes \(t_k\), we get the system of equations in terms of \(x_{k}\):

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle x_{0}= g(t_{0}),\quad t_{0}=a \\ \displaystyle x_{1}=\frac{1}{\varGamma (\alpha )}(\psi ( t_{1})-\psi ( t_{0}))^{\alpha -1}F( t_{0},x_{0})[ \psi ( t_{1})-\psi ( t_{0})] +g(t_{1})\\ \displaystyle x_{2i}=\frac{1}{\varGamma (\alpha )}\sum _{m=1}^{i}{(\psi ( t_{2i})-\psi (t_{2m-1}))^{\alpha -1}}F( t_{2m-1}, x_{2m-1})\\ \displaystyle \quad \times [ \psi ( t_{2m})-\psi ( t_{2m-2})] +g (t_{2i}),\quad i=1,\ldots ,n\\ \displaystyle x_{2j+1}=\frac{1}{\varGamma (\alpha )}(\psi ( t_{2j+1})-\psi (t_{0}))^{\alpha -1}F( t_{0}, x_{0})(\psi ( t_{1})-\psi ( t_{0}))\\ \displaystyle \quad +\frac{1}{\varGamma (\alpha )}\sum _{m=1}^{j}{(\psi ( t_{2j+1})-\psi (t_{2m}))^{\alpha -1}}F( t_{2m}, x_{2m})[ \psi ( t_{2m+1})-\psi (t_{2m-1})]\\ \displaystyle \quad +g (t_{2j+1}), \quad j=1,\ldots ,n-1, \\ \end{array}\right. \end{aligned}$$
(22)

where \(x_{0}=x(a)\) and \(x_{i}\approx x(t_{i})\), for \(i=1,2,\ldots ,2n\).

Now, we make the truncation error and convergent order analysis for the proposed numerical scheme. Let x be the solution to integral equation (3). As before, an approximation for \(x(t_k)\) at a node \(t_k\) is denoted by \(x_{k}\).

Theorem 4

Let \(\alpha >2\) and \(\beta =1\). Under assumptions (H1) and (H2) it holds that

$$\begin{aligned} \vert x(t_{k})-x_{k}\vert \le c_{6}h,\quad k=0,1,\ldots ,2n, \end{aligned}$$
(23)

as \(n\rightarrow \infty \), where

$$\begin{aligned} c_{6}= c_{3}[ \psi (b)-\psi (a)] e^{ c_{5}(b-a)},\quad c_{5}=\frac{L}{\varGamma (\alpha )}(\psi (b)- \psi (a))^{\alpha }. \end{aligned}$$

Proof

Let the error be denoted by

$$\begin{aligned} z_{k}=x(t_{k})-x_{k},\quad k=0,1,\ldots ,2n. \end{aligned}$$

Taking into account (21) and (22), we have

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle z_{0}=0 \\ \displaystyle z_{1}=R_{0}^{(0)}(x)\\ \displaystyle z_{2i}=\frac{1}{\varGamma (\alpha )}\sum _{m=1}^{i}{(\psi ( t_{2i})-\psi (t_{2m-1}))^{\alpha -1}}[ F( t_{2m-1},x( t_{2m-1}))\\ \displaystyle \quad -F( t_{2m-1}, x_{2m-1})] (\psi ( t_{2m})-\psi (t_{2m-2}))+\sum _{m=1}^{i}{R_{m}^{(2n)}(x,i),\quad i=1,\ldots ,n} \\ \displaystyle z_{2j+1}=\frac{1}{\varGamma (\alpha )}\sum _{m=1}^{j}{(\psi ( t_{2j+1})-\psi ( t_{2m}))^{\alpha -1}}[ F( t_{2m},x( t_{2m}))\\ \displaystyle \quad -F( t_{2m}, x_{2m})] (\psi ( t_{2m+1})-\psi ( t_{2m-1}))\\ \displaystyle \quad +R_{j}^{(0)}(x)+\sum _{m=1}^{j}{R_{m}^{(2n-1)}(x,j),\quad j=1,\ldots ,n-1.}\\ \end{array}\right. \end{aligned}$$
(24)

From this, by (18), we conclude that

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle \vert z_{1}\vert \le R(h)\\ \displaystyle \vert z_{k}\vert \le R(h)+ c_{5}h\sum \nolimits _{j=1}^{k-1}{\vert z_{j}\vert },\quad k=2,\ldots ,2n, \\ \end{array}\right. \end{aligned}$$
(25)

where

$$\begin{aligned} R(h)= c_{3}[ \psi (b)-\psi (a)] h. \end{aligned}$$
(26)

Now, let us consider the system

$$\begin{aligned} \epsilon _{k}=R(h)+ c_{5}h\sum _{j=1}^{k-1}{\vert \epsilon _{j}\vert },\quad k=2,\ldots ,2n, \end{aligned}$$
(27)

and \(\epsilon _{1}=R(h)\) as an initial condition. It is easily seen that \(\vert z_{k}\vert \le \epsilon _{k}\) for \(k=1,\ldots ,2n.\) Indeed, this can be verified by mathematical induction as follows: for \(k=1\) it is obvious. Let \(\vert z_{j}\vert \le \epsilon _{j}\) for \(j=1,\ldots ,k-1.\) Then, by inequality (25), we get

$$\begin{aligned} \vert z_{k}\vert \le R(h)+ c_{5}h\sum _{j=1}^{k-1}{ \epsilon _{j}=\epsilon _{k}.} \end{aligned}$$

Observe that \(\epsilon _{j}\), \(j=1,\ldots ,2n\), given by

$$\begin{aligned} \epsilon _{j}=R(h)[ 1+ c_{5}h] ^{j-1}, \end{aligned}$$

satisfy system (27). In fact, using the inequality

$$\begin{aligned} (1+\gamma )^{k-1}-1=\gamma \sum _{j=1}^{k-1}{(1+\gamma )^{j-1},\quad k\ge 2}, \end{aligned}$$

and taking \(\gamma = c_{5}h\) we obtain

$$\begin{aligned} R(h)+ c_{5}h\sum _{j=1}^{k-1}\epsilon _{j}= & {} R(h)\left\{ 1+c_{5}h\sum _{j=1}^{k-1}{(1+ c_{5}h)^{j-1}}\right\} \\= & {} R(h)\lbrace 1+[ (1+ c_{5}h)^{k-1}-1] \rbrace = \epsilon _{k}, \quad k\ge 2. \end{aligned}$$

Consequently, we get the following estimates:

$$\begin{aligned} \vert z_{k}\vert \le R(h)(1+ c_{5}h)^{k-1}, \quad k=1,\ldots ,2n. \end{aligned}$$
(28)

Using the fact that \((1+t)^{\frac{1}{t}}\) is increasing and approaches the number e as \(t\rightarrow 0^{+}\) we get

$$\begin{aligned} (1+ c_{5}h)^{k-1}\le (1+ c_{5}h)^{\frac{b-a}{h}}=[ (1+c_{5}h)^{\frac{1}{ c_{5}h}}] ^{ c_{5}(b-a)}\le e^{ c_{5}(b-a)} \end{aligned}$$

for \(k\le \frac{b-a}{h}\). This together with (26) and (28) gives (23). \(\square \)

Theorem 5

Let \(\alpha \in (1,2)\), \(p>1\) with \(p(2-\alpha )<1\),

$$\begin{aligned} \frac{1}{p}+\frac{1}{q}=1, \gamma =\min \left\{ \beta ,\frac{1}{q}\right\} ,\quad \gamma _{0}=\frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha }<1. \end{aligned}$$

Under assumptions (H1) and (H2) it holds that

$$\begin{aligned} \vert x(t_{k})-x_{k}\vert \le c_{7}h^{\gamma },\quad k=0,1,\ldots ,2n, \end{aligned}$$
(29)

as \(n\rightarrow \infty \), where

$$\begin{aligned} c_{7}=\frac{c_{4}}{1-\gamma _{0}}[ \psi (b)-\psi (a)]. \end{aligned}$$

Proof

Using the conditions of Theorem 5 and the estimates (20), from (24) we obtain following inequalities for \(z_{k}:\)

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle |z_{1}|\le c_{4}h^{\gamma }[ \psi (b)-\psi (a)] \\ \displaystyle |z_{2i}|\le \frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\sup _{k=1,\ldots ,2i-1} |z_{k}|\sum _{m=1}^{i}{[ \psi ( t_{2m})-\psi (t_{2m-2})]} \\ \displaystyle \quad +c_{4}h^{\gamma }[ \psi ( t_{2i})-\psi (a)] ,\quad i=1,\ldots ,n \\ \displaystyle |z_{2j+1}|\le \frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\sup _{k=2,\ldots ,2j} |z_{k}|\sum _{m=1}^{j}{[ \psi ( t_{2m+1})-\psi (t_{2m-1})]} \\ \displaystyle \quad +c_{4}h^{\gamma }[ \psi ( t_{2j+1})-\psi ( t_{1})] ,\quad j=1,\ldots ,n-1. \\ \end{array}\right. \end{aligned}$$
(30)

Set \(\Vert z\Vert =\sup _{k=0,\ldots ,2n} |z_{k}|\). Then, by (30), we have

$$\begin{aligned} |z_{2i}|\le \frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert z\Vert [ \psi ( t_{2i})-\psi (a)]+c_{4}h^{\gamma }[ \psi (t_{2i})-\psi (a)] ,\quad i=1,\ldots ,n, \end{aligned}$$

and

$$\begin{aligned} |z_{2j+1}|\le & {} \frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha -1}\Vert z\Vert \left[ \psi ( t_{2j+1})-\psi (a)\right] \\&+c_{4}h^{\gamma }\left[ \psi ( t_{2j+1})-\psi ( t_{1})\right] ,\quad j=1,\ldots ,n-1. \end{aligned}$$

Hence

$$\begin{aligned} \Vert z\Vert \le \frac{L}{\varGamma (\alpha )}(\psi (b)-\psi (a))^{\alpha }\Vert z\Vert +c_{4}h^{\gamma }[ \psi (b)-\psi (a)], \end{aligned}$$

what finishes the proof. \(\square \)

5 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 \(\alpha >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

$$\begin{aligned} E = \max _{i=0,1,\ldots ,2n}|x_{i}- x( t_{i})|, \end{aligned}$$

and the elapsed CPU time in seconds of the proposed numerical scheme (22).

Example 1

Consider the following fractional differential equation:

$$\begin{aligned} {^CD_{0+}^{5/2,\psi }}x(t)=\frac{x^{3}(t)}{1+x^{2}(t)}-\frac{1}{2}, \quad t\in [0,1], \end{aligned}$$
(31)

with the initial conditions and kernel

$$\begin{aligned} x_{0}=1, \, x_{0}^{1}=x_{0}^{2}=0 \quad \text{ and } \quad \psi (t)=t+t\sqrt{t}, \end{aligned}$$
(32)

respectively. Therefore, in integral equation (3) we have \(\alpha =\frac{5}{2}\) and

$$\begin{aligned} F(t,x)=\frac{x^{3}}{1+x^{2}},\quad g(t)=1-\frac{1}{2\varGamma (\frac{7}{2})}\left( t+t\sqrt{t}\right) ^{\frac{5}{2}}. \end{aligned}$$

In this case, conditions (H1) and (H2) are satisfied for \( L=9/8\), \(L_{1}=0\), and \(\beta =1\). Note that, the exact solution to (31) and (32) is \(x(t)=1\), \(t\in [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:

$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)=\frac{x^{5}(t)}{(1+x^{2}(t))^{2}}-\frac{243}{100}, \quad t\in [0,1], \end{aligned}$$
(33)

with the initial conditions and kernel

$$\begin{aligned} x_{0}=3, \quad x_{0}^{1}=0 \quad \text{ and } \quad \psi (t)=2t+t^{\frac{4}{3}}, \end{aligned}$$
(34)

respectively. Therefore, in integral equation (3) we have \(\alpha =\frac{3}{2}\) and

$$\begin{aligned} F(t,x)=\frac{x^{5}}{(1+x^{2})^{2}},\quad g(t)=3-\frac{243}{100\varGamma \left( \frac{5}{2}\right) }\left( 2t+t^{\frac{4}{3}}\right) ^{\frac{3}{2}}. \end{aligned}$$

In this case, conditions (H1) and (H2) are met for \(L=125/108\), \(L_{1}=0\), and \(\beta =1\). Note that the exact solution to (33) and (34) is \(x(t)=3\), \(t\in [0,1]\). The exact and numerical solutions for \(n=20,40,60\), as well as the evolution of the error, are shown in Fig. 2. The maximum of absolute error for different values of n and the elapsed CPU time in seconds are displayed in Table 2.

Fig. 1
figure 1

Comparison of the exact and numerical solutions (left), and errors (right) in Example 1

Full size image
Table 1 Maximum of the absolute error and elapsed CPU time in seconds for the numerical scheme (Example 1)
Full size table
Fig. 2
figure 2

Comparison of the exact and numerical solutions (left), and errors (right) in Example 2

Full size image
Table 2 Maximum of the absolute error and elapsed CPU time in seconds for the numerical scheme (Example 2)
Full size table

Example 3

Consider the following fractional differential equation:

$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)=x(t)+1, \quad t\in [0,5], \end{aligned}$$

with the initial conditions and kernel

$$\begin{aligned} x_{0}=0,\quad x_{0}^{1}=0,\quad \psi (t)=\sqrt{t+1}, \end{aligned}$$

respectively. Therefore, in integral equation (3) we have \(\alpha =\frac{3}{2}\) and

$$\begin{aligned} F(t,x)=x,\quad g(t)=\frac{\left( \sqrt{t+1}-1\right) ^{\frac{3}{2}}}{\varGamma \left( \frac{5}{2}\right) }. \end{aligned}$$

Once \({^CD_{0+}^{3/2,\psi }}E_{3/2}(\psi (t)-\psi (0))^{3/2}=(\psi (t)-\psi (0))^{3/2}\), where \(E_{3/2}(\cdot )\) denotes the Mittag–Leffler function of order \(\alpha =3/2\), we conclude that the solution to this problem is \(x(t)=E_{3/2}(\psi (t)-\psi (0))^{3/2}-1, t\in [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.

Fig. 3
figure 3

Comparison of the exact and numerical solutions (left), and errors (right) in Example 3

Full size image
Table 3 Maximum of the absolute error and elapsed CPU time in seconds for the numerical scheme (Example 3)
Full size table

In the next examples, we analyse the stability of the proposed numerical method.

Example 4

Consider the two fractional differential equations with initial conditions:

$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)=x(t)+1, \quad t\in [0,5], \quad x_{0}=0, \, x_{0}^{1}=0 \end{aligned}$$
(35)

and

$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)=x(t)+1, \quad t\in [0,5], \quad x_{0}=\mu , \quad x_{0}^{1}=0. \end{aligned}$$
(36)

The kernel is given by \(\psi (t)=\sqrt{t+1}\) and \(\mu \) is a real number. Let x and \(\overline{x}\) be solutions of (35) and (36), respectively. In this case, \(F(t,x)=x\) for both equations and function g is given by

$$\begin{aligned} g(t)=\frac{\left( \sqrt{t+1}-1\right) ^{\frac{3}{2}}}{\varGamma \left( \frac{5}{2}\right) } \quad \text{ and } \quad g(t)=\mu +\frac{\left( \sqrt{t+1}-1\right) ^{\frac{3}{2}}}{\varGamma \left( \frac{5}{2}\right) }, \end{aligned}$$

respectively. By the Grönwall inequality (cf. Almeida et al. 2019; Sousa and Oliveira 2019), we conclude that

$$\begin{aligned} |x(t)-\overline{x}(t)|\le |\mu | E_{3/2}\left( \left( \sqrt{t+1}-1\right) ^{3/2}\right) , \quad t\in [0,5]. \end{aligned}$$

Consider the two functions

$$\begin{aligned} X:\, t\mapsto |x(t)-\overline{x}(t)| \quad \text{ and } \quad E: \, t\mapsto |\mu | E_{3/2}\left( (\sqrt{t+1}-1)^{3/2}\right) . \end{aligned}$$

In Fig. 4, we present the plots for two values of \(\mu .\) 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.

Fig. 4
figure 4

Comparison of functions X and E, for \(\mu =0.1\) (left) and \(\mu =0.2\) (right) in Example 4

Full size image

In Table 4, we display approximate values of \(\overline{x}(5)\) for different values of \(\mu \) and n.

Table 4 Approximate values of \(\overline{x}(5)\), where \(\overline{x}\) is the solution to (36), for different values of \(\mu \) and n in Example 4
Full size table

Example 5

Let us now consider the following fractional differential equation:

$$\begin{aligned} {^CD_{0+}^{3/2+\mu _3,\psi }}x(t)=\sin ((\mu _1-1) x(t))+x(t)+\mu _1, \quad t\in [0,5], \end{aligned}$$
(37)

with initial conditions

$$\begin{aligned} x_{0}=\mu _2,\quad x_{0}^{1}=0, \end{aligned}$$

where \(\mu _1,\mu _2,\mu _3\in \mathbb R\) are three parameters. Let \(\psi (t)=\sqrt{t+1}\). In this case, we have

$$\begin{aligned} F(t,x)=\sin ((\mu _1-1) x)+x,\quad g(t)=\mu _2+\frac{\mu _1(\sqrt{t+1}-1)^{\frac{3}{2}+\mu _3}}{\varGamma \left( \frac{5}{2}+\mu _3\right) }. \end{aligned}$$

To analyse the stability of the numerical method, we consider perturbations of parameters \(\mu _1,\mu _2,\mu _3\). First, we fix \(\mu _2=0=\mu _3\) and we analyse the numerical solutions of the two following equations:

Fig. 5
figure 5

Numerical solutions of Eq. (37) for \(n=40\) and with respect to perturbations of parameters: \(\mu _1\) (left), \(\mu _2\) (center) and \(\mu _3\) (right). (Example 5)

Full size image
$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)= & {} \sin ((\mu _1-1) x(t))+x(t)+\mu _1, \quad t\in [0,5],\quad x_{0}=0,\quad x_{0}^{1}=0, \end{aligned}$$
(38)
$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)= & {} x(t)+1, \quad t\in [0,5], \quad x_{0}=0, \quad x_{0}^{1}=0, \end{aligned}$$
(39)

for \(\mu _1\in \{0.8,0.9,1,1.1,1.2\}\). Let x and \(\overline{x}\) be solutions of (38) and (39), respectively. Figure 5 (left) shows x for different values of \(\mu _1\) and \(n=40\). In Table 5, we display \(\max _i\vert x(t_i)-\overline{x}(t_i)\vert \). Next, we assume that \(\mu _1=1\), \(\mu _3=0\), and consider numerical solutions to equations:

$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)= & {} x(t)+1, \quad t\in [0,5],\quad x_{0}=\mu _2,\quad x_{0}^{1}=0, \end{aligned}$$
(40)
$$\begin{aligned} {^CD_{0+}^{3/2,\psi }}x(t)= & {} x(t)+1, \quad t\in [0,5], \quad x_{0}=0, \quad x_{0}^{1}=0, \end{aligned}$$
(41)

for \(\mu _2\in \{0,0.1,0.2,0.3,0.4\}\). Let x and \(\overline{x}\) be solutions of (40) and (41), respectively. Numerical results, for different values of \(\mu _2\), are presented in Fig. 5 (center) and Table 6.

Table 5 Numerical results for equations (38) and (39), with \(n=40\), in Example 5
Full size table
Table 6 Numerical results for equations (40) and (41), with \(n=40\), in Example 5
Full size table

Finally, we fix \(\mu _1=1\), \(\mu _2=0\), and analyse numerical solutions to equations:

$$\begin{aligned}&{^CD_{0+}^{3/2+\mu _3,\psi }}x(t)=x(t)+1, \quad t\in [0,5],\quad x_{0}=0, \quad x_{0}^{1}=0, \end{aligned}$$
(42)
$$\begin{aligned}&{^CD_{0+}^{3/2,\psi }}x(t)=x(t)+1, \quad t\in [0,5], \quad x_{0}=0,\quad x_{0}^{1}=0, \end{aligned}$$
(43)

for \(\mu _3\in \{0,0.1,0.2,0.3,0.4\}\). Let x and \(\overline{x}\) be solutions of (42) and (43), respectively. Numerical results, for different values of \(\mu _3\), are presented in Fig. 5 (right) and Table 7.

Table 7 Numerical results for Eqs. (42) and (43), with \(n=40\), in Example 5
Full size table

6 Conclusions

Recently, the Caputo derivative with respect to a kernel function \(\psi \) 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 \(\psi \)-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 \(\psi \)-Caputo derivative. This important issue will be considered in a forthcoming paper.