Detailed Error Analysis for a Fractional Adams Method on Caputo–Hadamard Fractional Differential Equations

Abstract

We consider a predictor–corrector numerical method for solving Caputo–Hadamard fractional differential equation over the uniform mesh $\log t_j=\log a+\left(\log \frac{t_N}{a}\right)\left(\frac{j}{N}\right),j=0,1,2,\dots ,N$ with $a\ge 1$, where $\log a=\log t_0<\log t_1<\dots <\log t_N=\log T$ is a partition of $[\log a,\log T]$. The error estimates under the different smoothness properties of the solution y and the nonlinear function f are studied. Numerical examples are given to verify that the numerical results are consistent with the theoretical results.

1. Introduction

Fractional differential equations have recently become an active research area due to the applications in fields including mechanics, computer science, and biology [1,2,3,4]. There are different fractional derivatives, e.g., Caputo, Riemman–Liouville, Riesz, etc. Such fractional derivatives are well studied in literature. One particular fractional derivative, the Hadamard fractional derivative, is also important and has been used to model physical problems in many fields [5,6,7,8,9,10,11]. The Hadamard fractional derivative was first introduced in early 1892 [12], and the Caputo–Hadamard derivative was suggested by Jarad et al. [8]. In this paper, we will discuss the numerical method for solving a Caputo–Hadamard fractional initial value problem. We will be analyzing the smoothness properties of various aspects of such equations and explain how these properties will affect the convergence order of the numerical method.

Consider the following Caputo–Hadamard fractional differential equation, with $\alpha >0$, [8]

$${}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)=f\left(t,y\left(t\right)\right),{\delta }^{k}y\left(a\right)=y_a^{\left(k\right)},k=0,1,\dots ,\lceil \alpha \rceil -1,$$

(1)

for $1\le a\le t\le T$. Here, $\lceil \alpha \rceil$ denotes the least integer bigger than or equal to $\alpha$. Here, $f(t,y)$ is a nonlinear function with respect to $y\in \mathbb{R}$, and the initial values $y_a^{\left(k\right)}$ are given. We also define $\alpha$ such that $n-1<\alpha <n$, for $n=1,2,3\dots$. Here, the fractional derivative ${}_{CH}{D}_{a,t}^{\alpha }$ denotes the Caputo–Hadamard derivative defined as,

$${}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)=\frac{1}{\Gamma (\lceil \alpha \rceil -\alpha )}{\int }_a^t{\left(\log \frac{t}{s}\right)}^{\lceil \alpha \rceil -\alpha -1}{\delta }^{n}y\left(s\right)\frac{ds}{s},t\ge a\ge 1,$$

(2)

with ${\delta }^{n}y\left(s\right)={\left(s\frac{d}{ds}\right)}^{n}y\left(s\right)$. It is well known that (1) is equivalent to the following Volterra integral equation, with $\alpha >0$ [13,14],

$$y\left(t\right)=\sum _{\nu =0}^{\lceil \alpha \rceil -1}{y}_a^{\left(\nu \right)}\frac{{\left(\log \frac{t}{a}\right)}^{\nu }}{\nu !}+\frac{1}{\Gamma \left(\alpha \right)}{\int }_a^t{\left(\log \frac{t}{s}\right)}^{\alpha -1}f\left(s,y\left(s\right)\right)\frac{\mathrm{d}\mathrm{s}}{s}.$$

(3)

Let us begin by reviewing some different numerical methods for solving (1). Diethelm et al. [15] first introduced the fractional Adams method for the Caputo fractional differential equation as a generalization of the classical numerical method for solving first-order ODEs. They discussed the error estimates for the proposed method, including the convergence orders under the different smoothness assumptions of y and f. This work was then extended by Liu et al. [16] to show that when ${{}_{C}D}_{a,t}^{\alpha }y\left(t\right)\notin C^2[0,T]$, the optimal convergence order is not achieved but can be recovered by implimenting a graded mesh. Gohar et al. [7] studied the existence and uniqueness of the solution of (1), and several numerical methods were introduced including the Euler and predictor–corrector methods. Gohar et al. [13] extended their work by introducing the rectangular, trapezoidal, and predictor–corrector methods for solving (1) with a uniform mesh under the smooth assumptions ${{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)\in C^2[a,T]$ with $\alpha \in (0,1)$. Green et al. [14] extended the results in [16] to the Caputo–Hadamard differential Equation (1). Numerical methods for solving Caputo–Hadamard time fractional partial differential equations can be found in [7,17]. More recent works on numerical methods for solving fractional differential equations can be referred to [17,18,19,20,21].

Caputo–Hadamard type fractional differential equations have the significant interests when they come to real applications due to the logarithmic nature of the integral kernel. This is especially true in mechanics and engineering. An example of this is the use of Hadamard type equations in fracture analysis and the modeling of elasticity [22]. Furthermore, the Caputo–Hadamard fractional derivative has been used in fractional turbulent flow models [23]. In biology, the Hadamard fractional derivative has been used in modeling for cancer treatments by radiotherapy [24]. Many of these models require efficient numerical methods for solving them. In literature, the error estimates of the numeircal methods proposed for solving the Caputo–Hadamard fractional differential Equation (1) are based on the assumptions that the solution y and f are sufficiently smooth [7,13]. In this paper, we will consider the error estimates of the numerical methods under the different smoothness assumptions of y and f.

Diethelm et al. [15] considered a variety of smoothness assumptions of y and f to the fractional Adams method for solving Caputo differential equations. The aim of this work is to extend the ideas in [15] for solving Caputo fractional differential equation to Caputo–Hadamard fractional differential Equation (1).

Let us start by briefly recalling the Adams–Bashforth–Moulton method for the Caputo–Hadamard fractional derivative. To construct such a method, we require an approximation of the integral in (3). We will apply the following product trapezoidal quadrature rule,

$${\int }_a^{t_{k+1}}{\left(\log \left(t_{k+1}\right)-\log \left(z\right)\right)}^{\alpha -1}g\left(z\right)\frac{dz}{z}\approx {\int }_a^{t_{k+1}}{\left(\log \left(t_{k+1}\right)-\log \left(z\right)\right)}^{\alpha -1}{\tilde{g}}_{k+1}\left(z\right)\frac{dz}{z},$$

(4)

where the approximation ${\tilde{g}}_{k+1}\left(z\right)$ is the piecewise linear interpolant for g at $t_j,j=0,1,2,\dots ,k+1$. Using the standard techniques from quadrature theory, we can write the right hand side integral as,

$${\int }_a^{t_{k+1}}{\left(\log \left(t_{k+1}\right)-\log \left(z\right)\right)}^{\alpha -1}{\tilde{g}}_{k+1}\left(z\right)=\sum _{j=0}^{k+1}{a}_{j,k+1}g\left(t_j\right),$$

(5)

where,

$$a_{j,k+1}={\int }_a^{t_{k+1}}{\left(\log \left(t_{k+1}\right)-\log \left(z\right)\right)}^{\alpha -1}{\varphi }_{j,k+1}dz,$$

(6)

and,

$${\varphi }_{j,k+1}\left(z\right)=\left\{\begin{array}{cc}\frac{\log z-\log t_{j-1}}{\log t_j-\log t_{j-1}},\hfill & \mathrm{if}{t}_{j-1}<z<t_j,\hfill \\ \frac{\log t_{j+1}-\log u}{\log t_{j+1}-\log t_j},\hfill & \mathrm{if}{t}_j<z<t_{j+1},\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(7)

Evaluating the integral, we get,

$$a_{j,k+1}=\frac{1}{\alpha (\alpha +1)}\left\{\begin{array}{c}\frac{1}{\log \frac{t_1}{a}}{A}_0,j=0,\hfill \\ \frac{1}{\log \frac{t_{j+1}}{t_j}}{A}_j+\frac{1}{\log \frac{t_{j-1}}{t_j}}{B}_j,j=1,2,\dots ,k,\hfill \\ {\left(\log \frac{t_{k+1}}{t_k}\right)}^{\alpha },j=k+1,\hfill \end{array}\right.$$

(8)

$$\begin{array}{cc}\hfill A_j=& {\left(\log \frac{t_{k+1}}{t_{j+1}}\right)}^{\alpha +1}-{\left(\log \frac{t_{k+1}}{t_j}\right)}^{\alpha +1}+(\alpha +1)\left(\log \frac{t_{j+1}}{t_j}\right){\left(\log \frac{t_{k+1}}{t_j}\right)}^{\alpha },j=0,1,...,k,\hfill \\ \hfill B_j=& {\left(\log \frac{t_{k+1}}{t_j}\right)}^{\alpha +1}-{\left(\log \frac{t_{k+1}}{t_{j-1}}\right)}^{\alpha +1}+(\alpha +1)\left(\log \frac{t_j}{t_{j-1}}\right){\left(\log \frac{t_{k+1}}{t_{j-1}}\right)}^{\alpha },j=1,2,...,k.\hfill \end{array}$$

This gives us the formula for the corrector, known as the fractional one-step Adams–Moulton method,

$$y_{k+1}=y_0+\frac{1}{\Gamma \left(\alpha \right)}\left(\sum _{j=0}^k{a}_{j,k+1}f(t_j,y_j)+a_{k+1,k+1}f(t_{k+1},y_{k+1}^P)\right).$$

(9)

We now must determine the predictor formula required to calculate $y_{k+1}^P$. For this, we will generalize the one-step Adams–Bashforth method for ODEs. To do this, we follow a similar method, but now, we will be replacing the integral on the right-hand side by the product rectangle formula rule,

$${\int }_a^{t_{k+1}}{\left(\log \left(t_{k+1}\right)-\log \left(z\right)\right)}^{\alpha -1}g\left(z\right)\frac{dz}{z}\approx \sum _{j=0}^k{b}_{j,k+1}g\left(t_j\right),$$

(10)

where,

$$b_{j,k+1}=\frac{1}{\Gamma (\alpha +1)}\left[{\left(\log \frac{t_{k+1}}{t_j}\right)}^{\alpha }-{\left(\log \frac{t_{k+1}}{t_{j+1}}\right)}^{\alpha }\right],j=0,1,2,\dots ,k.$$

(11)

Therefore, the predictor can be calculated using the fractional Adams–Bashforth method,

$$y_{k+1}^P=y_0+\sum _{j=0}^k{b}_{j,k+1}f(t_j,y_j).$$

(12)

Finally, we can conclude that our basic method, the fractional Adams–Bashforth–Moulton method, is described with the predictor Equation (12), the corrector Equation (9), and weights (8) and (11). For this method, however, we will be using a uniform mesh defined below.

Let N be a positive integer and let $a=t_0<t_1<\cdots <t_N=T$ be the partition on $[a,T]$. We define the following uniform mesh on $\left[\log \right(a),\log (T\left)\right]$ with $\log a=\log t_0<\log t_1<\cdots <\log t_N=\log T,$

$$\frac{\log t_j-\log a}{\log t_N-\log a}=\frac{j}{N},$$

(13)

which implies that,

$$\log t_j=\log a+jh,\mathrm{with}h=(\log t_N-\log a)/N.$$

(14)

when $j=N$ we have $\log t_N=\log T$ and when $j=0$, we have $\log t_0=\log a$. By applying such a mesh on our weights, we can simplify them to,

$$a_{j,k+1}=\frac{h^{\alpha }}{\Gamma (\alpha +2)}\times \left\{\begin{array}{cc}\left(k^{\alpha +1}-(k-\alpha ){(k+1)}^{\alpha }\right),\hfill & \mathrm{if}j=0,\hfill \\ {(k-j+2)}^{\alpha +1}+{(k-j)}^{\alpha +1}-2{(k-j+1)}^{\alpha },\hfill & \mathrm{if}1\le j\le k,\hfill \\ 1,\hfill & \mathrm{if}j=k+1,\hfill \end{array}\right.$$

(15)

and

$$b_{j,k+1}=\frac{h^{\alpha }}{\Gamma (\alpha +1)}\left({(k+1-j)}^{\alpha }-{(k-j)}^{\alpha }\right).$$

(16)

The fractional Adams–Bashforth–Moulton method has many useful properties. First, we may use this method to solve nonlinear fractional differential equation by transforming the nonlinear equation into a Volterra integral equation with a singular kernel and approximating the corresponding Volterra integral equation with some suitable quadrature formula. Second, the fractional Adams–Bashforth–Moulton method is an explicit method, and as such, it can save computer memory and can be more computationally time-efficient. Finally, this method has been shown to have a convergence order of $O\left(h^{2-\alpha }\right)$ for suitable meshes.

We shall consider the error estimates of the predictor–corrector schemes (9) and (12) under the various smoothness assumptions of f and y in this paper. The paper is organized as follows. In Section 2, we will consider some auxiliary results that will aid us for the remainder of the paper in proving the error estimates. In Section 3, we consider the error estimates of the predictor–corrector method for solving (1) with a uniform mesh under different smoothness conditions. In Section 4, we will provide several numerical examples that support the theoretical conclusions made in Section 2 and Section 3.

2. Auxiliary Results

For the remainder of this work, we will be using the Adams method described by the predictor Equation (12) and the corrector Equation (9) and the uniform mesh described by (13) to solve the fractional initial value problem (1). To begin, we must apply some conditions on the function f, namely that f is continuous and follows the Lipschitz condition with respect to its second argument with the Lipschitz constant, L, on a suitable set G. By forcing these conditions, we may use [25] to show that a unique solution y of the initial value problem exists on the interval $[a,T]$. Our goal is to find a suitable approximation for this unique solution.

As such, we will introduce several auxiliary results on certain smoothness properties to help us in our error analysis of this method. Our first result is taken from [13].

Theorem 1.

(a) Let $\alpha >0$. Assume that $f\in C^2\left(G\right)$ where $G\subset {\mathbb{R}}^2$ is a suitable set. Define $\widehat{v}=\lceil \frac{1}{\alpha }\rceil -1$. Then, there exists a function $\varphi \in C^1[a,T]$ and some constants $c_1,c_2,\dots ,c_{\widehat{v}}\in \mathbb{R}$ such that the solution y of (1) can be expressed in the following form:

$$y\left(t\right)=\varphi \left(t\right)+\sum _{v=1}^{\widehat{v}}{c}_v{\left(\log \frac{t}{a}\right)}^{v\alpha }.$$

(17)

(b) Assume that $f\in C^3\left(G\right)$. Define $\widehat{v}=\lceil \frac{2}{\alpha }\rceil -1$. and $\tilde{v}=\lceil \frac{1}{\alpha }\rceil -1$. Then, there exists a function $\psi \in C^2[a,T]$ and some constants $c_1,c_2,\dots ,c_{\widehat{v}}\in \mathbb{R}$ and $d_1,d_2\dots ,d_{\tilde{v}}\in \mathbb{R}$ such that the solution y of (1) can be expressed in the form,

$$y\left(t\right)=\psi \left(t\right)+\sum _{v=1}^{\widehat{v}}{c}_v{\left(\log \frac{t}{a}\right)}^{v\alpha }+\sum _{v=1}^{\tilde{v}}{d}_v{\left(\log \frac{t}{a}\right)}^{1+v\alpha }.$$

(18)

We can also relate such smoothness properties of the solution with that of the Caputo–Hadamard derivatives. From [13], we have,

Theorem 2.

If $y\in C^m[a,T]$ for some $m\in \mathbb{N}$ and $0<\alpha <m$, then,

$${{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)=\Phi \left(\log \frac{t}{a}\right)+\sum _{l=0}^{m-\lceil \alpha \rceil -1}\frac{{\delta }^{l+\lceil \alpha \rceil }y\left(a\right)}{\Gamma (\lceil a\rceil -\alpha +l+1)}{\left(\log \frac{t}{a}\right)}^{\lceil a\rceil -\alpha +l},$$

(19)

where $\Phi \in C^{m-\lceil \alpha \rceil }[a,T]$ and ${\delta }^{n}y\left(s\right)={\left(s\frac{\mathrm{d}}{\mathrm{ds}}\right)}^{n}y\left(s\right)$ with $n\in \mathbb{N}$. Moreover, the ($m-\lceil \alpha \rceil )$th derivative of g satisfies a Lipschitz condition of order $\lceil \alpha \rceil -\alpha$.

A useful corollary can then be drawn from this theorem and can be used to generalize the classical result for derivatives of integer order.

Corollary 1.

Let $y\in C^m[a,T]$ for some $m\in \mathbb{N}$ and assume that $0<\alpha <m$. Then, ${{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)\in C[a,T]$.

We will now introduce some error estimates for both the product rectangle rule and product trapezoidal rule, which we have implemented as the predictor and corrector of our method. Doing so will aid us in producing error estimates and convergence orders of the fractional Adams method.

Theorem 3.

(a) Let $z\left(t\right)=z\left(\log \frac{t}{a}\right)\in C^1[a,T]$ and $\delta z=\left(t\frac{d}{dt}\right)z$. Let h be defined by (14). Then,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t}{a}\right)|\le \frac{1}{\alpha }{\left|\right|\delta z\left|\right|}_{\infty }{\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha }h.$$

(20)

(b) Let $z\left(\log \frac{t}{a}\right)={\left(\log \frac{t}{a}\right)}^p$ for some $p\in (0,1)$. Then,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t}{a}\right)|\le C_{\alpha ,p}^{Re}{\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha +p-1}h.$$

(21)

Proof. 

Using the construction of the product rectangle formula rule, we can show that in both cases above, the quadrature error can be represented by,

$$\begin{array}{cc}\hfill {\int }_a^{t_{k+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t_j}{a}\right)\hfill \\ \hfill =& \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left(z\left(\log \frac{t}{a}\right)-z\left(\log \frac{t_j}{a}\right)\right)\frac{dt}{t}.\hfill \end{array}$$

(22)

We will begin by proving statement (a). We first apply the Mean Value Theorem on z in the second factor of the integrand and note that $\delta z\left(t\right)=\left(t\frac{d}{dt}\right)z\left(t\right)=z^{\prime }\left(\log \frac{t}{a}\right)$. This gives us, with ${\parallel \delta z\parallel }_{\infty }={max}_{t\in [0,T]}\left|\delta z\left(t\right)\right|$,

$$\begin{array}{cc}\hfill |{\int }_a^{t_{k+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\hfill \\ \hfill \le & {\left|\right|\delta z\left|\right|}_{\infty }\sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}(\log t-\log t_j)\frac{dt}{t}\hfill \\ \hfill =& {\left|\right|\delta z\left|\right|}_{\infty }\frac{h^{1+\alpha }}{\alpha }\sum _{j=0}^k\left(\frac{1}{1+\alpha }\left[{(k+1-j)}^{1+\alpha }-{(k-j)}^{1+\alpha }\right]-{(k-j)}^{\alpha }\right)\hfill \\ \hfill =& {\left|\right|\delta z\left|\right|}_{\infty }\frac{h^{1+\alpha }}{\alpha }\left(\frac{{(k+1)}^{1+\alpha }}{1+\alpha }-\sum _{j=0}^k{j}^{\alpha }\right)={\left|\right|\delta z\left|\right|}_{\infty }\frac{h^{1+\alpha }}{\alpha }\left({\int }_0^{k+1}{t}^{\alpha }-\sum _{j=0}^k{j}^{\alpha }\right).\hfill \end{array}$$

(23)

The final bracket can be interpreted as the remainder of the standard rectangle quadrature formula for the function $t^{\alpha }$ over the interval $[0,k+1]$. We can then apply a standard result from quadrature theory that states that this term is bounded by the total variation of the integrand, ${(k+1)}^{\alpha }$. Therefore, we can conclude that,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\le {\left|\right|\delta z\left|\right|}_{\infty }\frac{h^{1+\alpha }}{\alpha }{(k+1)}^{\alpha }.$$

(24)

Now to prove (b). We use the fact that $z\left(\log \frac{t}{a}\right)={\left(\log \frac{t}{a}\right)}^p$ is monotonic and repeated applications of the Mean Value Theorem. We are required to take cases when $0<\alpha <1$ and $\alpha >1$.

Case 1: $0<\alpha <1$,

$$\begin{array}{cc}\hfill |{\int }_a^{t_{k+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\hfill \\ \hfill \le & \sum _{j=0}^k|z\left(\log \frac{t_{j+1}}{a}\right)-z\left(\log \frac{t_j}{a}\right)|{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\frac{dt}{t}\hfill \\ \hfill =& \frac{h^{\alpha +p}}{\alpha }\sum _{j=0}^k\left({(j+1)}^p-j^p\right)\left({(k+1-j)}^{\alpha }-{(k-j)}^{\alpha }\right)\hfill \\ \hfill \le & \frac{h^{\alpha +p}}{\alpha }\left({(k+1)}^{\alpha }-k^{\alpha }+p\alpha \sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}+{(k+1)}^p-k^p\right)\hfill \\ \hfill \le & \frac{h^{\alpha +p}}{\alpha }\left(\alpha k^{\alpha -1}+p\alpha \sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}+p{k}^{p-1}\right).\hfill \end{array}$$

We shall show that,

$$\sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}\le C{k}^{\alpha +p-1}.$$

(25)

Assuming (25) holds at the moment, we then have,

$$\begin{array}{cc}\hfill I\le & \frac{h^{\alpha +p}}{\alpha }\left(\alpha k^{\alpha -1}+p\alpha C{k}^{\alpha +p-1}+p{k}^{p-1}\right)\le C{h}^{\alpha +p}\left(k^{\alpha +p-1}+C{k}^{\alpha +p-1}+k^{\alpha +p-1}\right)\hfill \\ \hfill \le & C{h}^{\alpha +p}{k}^{\alpha +p-1}=C{h}^{\alpha +p}{\left(\log \frac{t_{k+1}}{\alpha }\right)}^{\alpha +p-1}{h}^{-(\alpha +p-1)}=Ch{\left(\log \frac{t_{k+1}}{\alpha }\right)}^{\alpha +p-1}.\hfill \end{array}$$

It remains to prove (25), which we shall do now.

$$\sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}=1^{p-1}{(k-1)}^{\alpha -1}+2^{p-1}{(k-2)}^{\alpha -1}+\dots +{(k-1)}^{p-1}{\left(1\right)}^{\alpha -1}.$$

Let $F\left(x\right)=x^{p-1}{(k-x)}^{\alpha -1}$ for $0\le x\le k$. Then,

$$F^{\prime }\left(x\right)=x^{p-2}{(k-x)}^{\alpha -2}\left[(2-\alpha -p)x-(1-p)k\right].$$

By letting the parenthesis equal zero, we find that $F\left(x\right)$ has a turning point at $x_p=\frac{(1-p)k}{2-\alpha -p}$, meaning $F\left(x\right)$ is decreasing on $[0,x_p]$ and increasing on $[x_p,k]$. Therefore,

$$\begin{array}{cc}\hfill \sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}\le & {\int }_0^{x_p}{x}^{p-1}{(k-x)}^{\alpha -1}dx+{\int }_{x_p}^k{x}^{p-1}{(k-x)}^{\alpha -1}dx\hfill \\ \hfill \le & {(k-x_p)}^{\alpha -1}{\int }_0^{x_p}{x}^{p-1}dx+x_p^{p-1}{\int }_{x_p}^k{(k-x)}^{\alpha -1}dx\hfill \\ \hfill =& {\left(k-\frac{1-p}{2-\alpha -p}k\right)}^{\alpha -1}\frac{x_p^p}{p}+x_p^{p-1}\frac{{(k-x_p)}^{\alpha }}{\alpha }\hfill \\ \hfill =& C_1{k}^{\alpha -1}{k}^p+C_2{k}^{p-1}{k}^{\alpha }\le C{k}^{p+\alpha -1}.\hfill \end{array}$$

Case 2: $\alpha >1$,

$$\begin{array}{cc}\hfill |{\int }_a^{t_{k+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\hfill \\ \hfill =& \frac{h^{\alpha +p}}{\alpha }\sum _{j=0}^k\left({(j+1)}^p-j^p\right)\left({(k+1-j)}^{\alpha }-{(k-j)}^{\alpha }\right)\hfill \\ \hfill \le & \frac{h^{\alpha +p}}{\alpha }\left({(k+1)}^{\alpha }-k^{\alpha }+p\alpha \sum _{j=1}^{k-1}{j}^{p-1}{(k-j+1)}^{\alpha -1}+{(k+1)}^p-k^p\right)\hfill \\ \hfill \le & \frac{h^{\alpha +p}}{\alpha }\left(\alpha {(k+1)}^{\alpha -1}+p\alpha \sum _{j=1}^{k-1}{j}^{p-1}{(k-j+1)}^{\alpha -1}+{(k+1)}^p-k^p\right).\hfill \end{array}$$

We shall show that,

$$\sum _{j=1}^{k-1}{j}^{p-1}{(k-j+1)}^{\alpha -1}\le C{k}^{\alpha +p-1}.$$

(26)

Assuming (26) holds at the moment, we then have,

$$\begin{array}{cc}\hfill I\le & \frac{h^{\alpha +p}}{\alpha }\left(\alpha {(k+1)}^{\alpha -1}+p\alpha C{k}^{\alpha +p-1}+p{k}^{p-1}\right)\hfill \\ \hfill \le & C{h}^{\alpha +p}\left(C{k}^{\alpha +p-1}+C{k}^{\alpha +p-1}+C{k}^{\alpha +p-1}\right)\hfill \\ \hfill \le & C{h}^{\alpha +p}{k}^{\alpha +p-1}=C{h}^{\alpha +p}{\left(\log \frac{t_{k+1}}{\alpha }\right)}^{\alpha +p-1}{h}^{-(\alpha +p-1)}=Ch{\left(\log \frac{t_{k+1}}{\alpha }\right)}^{\alpha +p-1}.\hfill \end{array}$$

It remains to prove (26), which we shall do now.

$$\sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}=1^{p-1}{\left(k\right)}^{\alpha -1}+2^{p-1}{(k-1)}^{\alpha -1}+\dots +{(k-1)}^{p-1}{\left(2\right)}^{\alpha -1}.$$

Let $F\left(x\right)=x^{p-1}{(k-x)}^{\alpha -1}$ for $0\le x\le k$. Then,

$$F^{\prime }\left(x\right)=x^{p-2}{(k-x)}^{\alpha -2}\left[(2-\alpha -p)x-(1-p)(k+1)\right].$$

By letting the parenthesis equal zero, we find that $F\left(x\right)$ has a turning point at $x_p=\frac{(1-p)(k+1)}{2-\alpha -p}$, meaning F(x) is decreasing on $[0,x_p]$ and increasing on $[x_p,k]$. Therefore,

(i)

If $2-\alpha -p>0$, then $x_p>0$. Therefore,

$$\begin{array}{cc}\hfill \sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}\le & {\int }_0^{x_p}{x}^{p-1}{(k+1-x)}^{\alpha -1}dx+{\int }_{x_k}^p{x}^{p-1}{(k+1-x)}^{\alpha -1}dx\hfill \\ \hfill \le & {(k+1-0)}^{\alpha -1}\frac{x_p^p}{p}+x_p^{p-1}\frac{{(k+1-x_p)}^{\alpha }}{\alpha }\hfill \\ \hfill =& {(k+1)}^{\alpha -1}{C}_1{(k+1)}^p+C_2{(k+1)}^{p-1}{C}_3{(k+1)}^{\alpha }\hfill \\ \hfill \le & C{(k+1)}^{\alpha +p-1}\le C{k}^{\alpha +p-1}.\hfill \end{array}$$

(ii)

If $2-\alpha -p<0$, then $x_p<0$. This means that $F^{\prime }\left(x\right)<0$ and shows that $F\left(x\right)$ is a decreasing function. Therefore,

$$\begin{array}{cc}\hfill \sum _{j=1}^{k-1}{j}^{p-1}{(k-j)}^{\alpha -1}\le & {\int }_0^k{x}^{p-1}{(k+1-x)}^{\alpha -1}dx\hfill \\ \hfill \le & {(k+1-0)}^{\alpha -1}{\int }_0^k{x}^{p-1}dx\le {(k+1)}^{\alpha -1}\frac{k^p}{p}\le C{k}^{\alpha +p-1},\hfill \end{array}$$

which gives us the desired result. □

Next, we come to corresponding results for the trapezoidal formula that has been used for the corrector of our method. The proofs of this theorem are similar to those of the previous theorem.

Theorem 4.

(a) Let $z\left(t\right)=z\left(\log \frac{t}{a}\right)\in C^2[a,T]$ and ${\delta }^{2}z={\left(t\frac{d}{dt}\right)}^{2}z$. Let h be defined by (14). Then, there exists a constant $C_{\alpha }^{Tr}$ depending on α such that,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\le C_{\alpha }^{Tr}\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }{\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha }{h}^2.$$

(27)

(b) Let $z\left(\log \frac{t}{a}\right)\in C^1[a,T]$ and assume that $\delta z$ fulfils a Lipschitz condition of order μ for some $\mu \in (0,1)$. Then, there exists a positive constant $B_{\alpha ,\mu }^{Tr}$ and $M(z,\mu )$ such that,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\le B_{\alpha ,\mu }^{Tr}M(z,\mu ){\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha }{h}^{1+\mu }.$$

(28)

(c) Let $z\left(\log \frac{t}{a}\right)={\left(\log \frac{t}{a}\right)}^p$ for some $p\in (0,2)$ and $\varsigma :=min(2,p+1)$. Then,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\le C_{\alpha ,p}^{Tr}{\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha +p-\varsigma }{h}^{\varsigma }.$$

(29)

Proof. 

By construction of the product trapezoidal formula rule, we can show that in the first two cases above, the quadrature error can be represented as,

$$\begin{array}{cc}\hfill {\int }_a^{t_{k+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}\left(\log \frac{t_j}{a}\right)\hfill \\ \hfill =& \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left(z\left(\log \frac{t}{a}\right)-P_1\left(t\right)\right)\frac{dt}{t},\hfill \end{array}$$

(30)

where,

$$P_1\left(t\right)=\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}z\left(t_j\right)+\frac{\log \frac{t}{t_j}}{\log \frac{t_j+1}{t_j}}z\left(t_{j+1}\right),t\in [t_j,t_{j+1}].$$

(31)

We will begin by proving statement (a). To find an estimate for (30), we must simplify the second factor of the integrand by applying the Mean Value Theorem. Therefore,

$$\begin{array}{cc}\hfill z(\log & \frac{t}{a})-P_1\left(t\right)\hfill \\ \hfill =& \frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}\left(z\left(\log \frac{t}{a}\right)-z\left(\log \frac{t_j}{a}\right)\right)+\frac{\log \frac{t}{t_j}}{\log \frac{t_j+1}{t_j}}\left(z\left(\log \frac{t}{a}\right)-z\left(\log \frac{t_{j+1}}{a}\right)\right)\hfill \\ \hfill =& \frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}\delta z\left(c_1\right)(\log t-\log t_j)+\frac{\log \frac{t}{t_j}}{\log \frac{t_j+1}{t_j}}\delta z\left(c_2\right)(\log t-\log t_{j+1})\hfill \\ \hfill =& \frac{(\log t-\log t_{j+1})(\log t-\log t_j)}{(\log t_{j+1}-\log t_j)}\left(\delta z\left(c_1\right)-\delta z\left(c_2\right)\right)\hfill \\ \hfill \le & h{\delta }^{2}z\left(c_3\right)(c_1-c_2)\le \left|\right|{\delta }^{2}z{\left|\right|}_{\infty }{h}^2.\hfill \end{array}$$

Applying the above to (30), we get,

$$\begin{array}{cc}\hfill \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left(z\left(\log \frac{t}{a}\right)-P_1\left(t\right)\right)\frac{dt}{t}\hfill \\ \hfill \le \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}& {\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left|\right|{\delta }^2{z\left|\right|}_{\infty }{h}^2\frac{dt}{t}=\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }{h}^2\sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\frac{dt}{t}\hfill \\ \hfill =\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }{h}^2& \sum _{j=0}^k{\left[\frac{-1}{\alpha }{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha }\right]}_{t_j}^{t_{j+1}}\hfill \\ \hfill =\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }{h}^2& \frac{-1}{\alpha }\sum _{j=0}^k{\left(\log \frac{t_{k+1}}{t_{j+1}}\right)}^{\alpha }-{\left(\log \frac{t_{k+1}}{t_j}\right)}^{\alpha }\hfill \\ \hfill =\frac{1}{\alpha }\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }& h^2\sum _{j=0}^k{h}^{\alpha }{(k-j+1)}^{\alpha }-h^{\alpha }{(k-j)}^{\alpha }=\frac{1}{\alpha }\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }{h}^2{h}^{\alpha }{(k+1)}^{\alpha }\hfill \\ \hfill =\frac{1}{\alpha }\left|\right|{\delta }^{2}z{\left|\right|}_{\infty }& h^2{\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha }.\hfill \end{array}$$

Now to prove (b). As $z\in C^1[a,T]$, this time, we are unable to apply the mean value theorem for a second time. Instead, we will apply the Lipschitz condition of order $\mu$ on $\delta z$ for some $\mu \in (0,1)$. Therefore,

$$\begin{array}{cc}\hfill z(\log & \frac{t}{a})-P_1\left(t\right)\hfill \\ \hfill =& \frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}\left(z\left(\log \frac{t}{a}\right)-z\left(\log \frac{t_j}{a}\right)\right)+\frac{\log \frac{t}{t_j}}{\log \frac{t_j+1}{t_j}}\left(z\left(\log \frac{t}{a}\right)-z\left(\log \frac{t_{j+1}}{a}\right)\right)\hfill \\ \hfill =& \frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}\delta z\left(c_1\right)(\log t-\log t_j)+\frac{\log \frac{t}{t_j}}{\log \frac{t_j+1}{t_j}}\delta z\left(c_2\right)(\log t-\log t_{j+1})\hfill \\ \hfill =& \frac{(\log t-\log t_{j+1})(\log t-\log t_j)}{(\log t_{j+1}-\log t_j)}\left(\delta z\left(c_1\right)-\delta z\left(c_2\right)\right)\hfill \\ \hfill \le & hM|c_1-c_2{|}^{\mu }\le M{h}^{1+\mu },\hfill \end{array}$$

where M is a constant depending on z and $\mu$. Applying the above to (30), we get,

$$\begin{array}{cc}\hfill & \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left(z\left(\log \frac{t}{a}\right)-P_1\left(t\right)\right)\frac{dt}{t}\hfill \\ \hfill \le & \sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}M{h}^{1+\mu }\frac{dt}{t}=M{h}^{1+\mu }\sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\frac{dt}{t}\hfill \\ \hfill =& M{h}^{1+\mu }\sum _{j=0}^k{\left[\frac{-1}{\alpha }{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha }\right]}_{t_j}^{t_{j+1}}=M{h}^{1+\mu }\frac{-1}{\alpha }\sum _{j=0}^k{\left(\log \frac{t_{k+1}}{t_{j+1}}\right)}^{\alpha }-{\left(\log \frac{t_{k+1}}{t_j}\right)}^{\alpha }\hfill \\ \hfill =& \frac{1}{\alpha }M{h}^{1+\mu }\sum _{j=0}^k{h}^{\alpha }{(k-j+1)}^{\alpha }-h^{\alpha }{(k-j)}^{\alpha }=\frac{1}{\alpha }M{h}^{1+\mu }{h}^{\alpha }{(k+1)}^{\alpha }\hfill \\ \hfill =& \frac{1}{\alpha }M{h}^{1+\mu }{\left(\log \frac{t_{k+1}}{a}\right)}^{\alpha }.\hfill \end{array}$$

Finally, we shall prove (c). We will start by proving the theorem when $0<p<1$.

Case 1: $0<\alpha <1$. Let

$$\begin{array}{cc}\hfill & A=|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\hfill \\ \hfill =& \left|\sum _{j=0}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}{\left(\log \frac{t_j}{a}\right)}^p+\frac{\log \frac{t}{t_j}}{\log \frac{t_{j+1}}{t_j}}{\left(\log \frac{t_{j+1}}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill \le & \left|{\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_0}}{\log \frac{t_1}{t_0}}{\left(\log \frac{t_1}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill +& \left|\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}{\left(\log \frac{t_j}{a}\right)}^p+\frac{\log \frac{t}{t_j}}{\log \frac{t_{j+1}}{t_j}}{\left(\log \frac{t_{j+1}}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill =& I+II.\hfill \end{array}$$

For I, we have,

$$\begin{array}{cc}\hfill I=& \left|{\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_0}}{\log \frac{t_1}{t_0}}{\left(\log \frac{t_1}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill =& \left|{\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{a}}{h}{h}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill \le & {\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}{h}^p\frac{dt}{t}=h^p{h}^{\alpha }\left[\frac{{(k+1)}^{\alpha }-k^{\alpha }}{\alpha }\right]\hfill \\ \hfill =& h^{\alpha +p}\frac{\alpha k^{\alpha -1}}{\alpha }=h^{p+\alpha }{k}^{\alpha -1}=C{h}^{p+\alpha }{(k+1)}^{\alpha -1}=C{h}^{p+1}\log {\left(\frac{t_{k+1}}{a}\right)}^{\alpha -1}.\hfill \end{array}$$

For $II$, we have,

$$\begin{array}{cc}\hfill II=& \sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}[{\left(\log \frac{t}{a}\right)}^p-(\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}{\left(\log \frac{t_j}{a}\right)}^p\hfill \\ \hfill & +\frac{\log \frac{t}{t_j}}{\log \frac{t_{j+1}}{t_j}}{\left(\log \frac{t_{j+1}}{a}\right)}^p\left)\right]\frac{dt}{t}\hfill \\ \hfill =& \sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}[\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}\left({\left(\log \frac{t}{a}\right)}^p-{\left(\log \frac{t_j}{a}\right)}^p\right)\hfill \\ \hfill & +\frac{\log \frac{t}{t_j}}{\log \frac{t_{j+1}}{t_j}}\left({\left(\log \frac{t}{a}\right)}^p-{\left(\log \frac{t_{j+1}}{a}\right)}^p\right)]\frac{dt}{t}\hfill \\ \hfill =& \sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}[\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}p{\xi }_2^{p-1}(\log t-\log t_{j+1})\hfill \\ \hfill & +\frac{\log \frac{t}{t_j}}{\log \frac{t_{j+1}}{t_j}}p{\xi }_1^{p-1}(\log t-\log t_j)]\frac{dt}{t}\hfill \\ \hfill =& C\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\frac{\log \frac{t}{t_j}\log \frac{t}{t_{j+1}}}{\log \frac{t_{j+1}}{t_j}}\left[{\xi }_2^{p-1}-{\xi }_1^{p-1}\right]\frac{dt}{t}.\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill \left|II\right|\le & C\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left|\frac{\log \frac{t}{t_j}\log \frac{t}{t_{j+1}}}{\log \frac{t_{j+1}}{t_j}}\right|\left[{\left(\log \frac{t_j}{a}\right)}^{p-1}-{\left(\log \frac{t_{j+1}}{a}\right)}^{p-1}\right]\frac{dt}{t}\hfill \\ \hfill \le & C{h}^p\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[j^{p-1}-{(j+1)}^{p-1}\right]\frac{dt}{t}\hfill \\ \hfill =& C{h}^{p+\alpha }\left[\sum _{j=1}^{k-1}\left[j^{p-1}-{(j+1)}^{p-1}\right]\left[{(k+1-j)}^{\alpha }-{(k-j)}^{\alpha }\right]+\left[k^{p-1}-{(k+1)}^{p-1}\right]\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[\sum _{j=1}^{k-1}\left|(p-1)j^{p-2}\right|\left[\alpha {(k-j)}^{\alpha -1}\right]+k^{p-2}\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[\sum _{j=1}^{k-1}{j}^{p-2}{(k-j)}^{\alpha -1}+k^{p-2}\right].\hfill \end{array}$$

Similar to our previous proof, we now must find the bounds for the summation.

$$\begin{array}{cc}\hfill \sum _{j=1}^{k-1}{j}^{p-2}{(k-j)}^{\alpha -1}\le & {\int }_1^{x_p}{x}^{p-2}{(k-x)}^{\alpha -1}dx+{\int }_{x_p}^k{x}^{p-2}{(k-x)}^{\alpha -1}dx\hfill \\ \hfill \le & {(k-x_p)}^{\alpha -1}{\int }_1^{x_p}{x}^{p-2}dx+x_p^{p-2}{\int }_{x_p}^k{(k-x)}^{\alpha -1}dx\hfill \\ \hfill =& {(k-x_p)}^{\alpha -1}\frac{1}{1-p}+x_p^{p-2}\frac{{(k-x_p)}^{\alpha }}{\alpha }\hfill \\ \hfill =& C_1{k}^{\alpha -1}+C_2{k}^{p-2}·C_3{k}^{\alpha }\le C{k}^{\alpha -1}+C{k}^{\alpha +p-2}.\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill \left|II\right|\le & C{h}^{\alpha +p}\left[\sum _{j=1}^{k-1}{j}^{p-2}{(k-j)}^{\alpha -1}+k^{p-2}\right]\le C{h}^{\alpha +p}\left[C{k}^{\alpha -1}+C{k}^{\alpha +p-2}+k^{p-2}\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[C{k}^{\alpha -1}+C{k}^{\alpha +p-2}+k^{\alpha +p-2}\right]\le C{h}^{\alpha +p}\left[C{k}^{\alpha -1}+C{k}^{\alpha +p-2}\right]\hfill \\ \hfill \le & C{h}^{p+1}\log {\left(\frac{t_k}{a}\right)}^{\alpha -1}+C{h}^2\log {\left(\frac{t_{k+1}}{a}\right)}^{\alpha +p-2}\le C{h}^{p+1}\log {\left(\frac{t_{k+1}}{a}\right)}^{\alpha -1}.\hfill \end{array}$$

Thus, we get the desired result.

Case 2: $1<\alpha <2$. Let

$$\begin{array}{cc}\hfill & A=|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}z\left(\log \frac{t}{a}\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}z\left(\log \frac{t_j}{a}\right)|\hfill \\ \hfill \le & \left|{\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_0}}{\log \frac{t_1}{t_0}}{\left(\log \frac{t_1}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill +& \left|\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_{j+1}}}{\log \frac{t_j}{t_{j+1}}}{\left(\log \frac{t_j}{a}\right)}^p+\frac{\log \frac{t}{t_j}}{\log \frac{t_{j+1}}{t_j}}{\left(\log \frac{t_{j+1}}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill =& I+II.\hfill \end{array}$$

For I, we have,

$$\begin{array}{cc}\hfill I=& \left|{\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{t_0}}{\log \frac{t_1}{t_0}}{\left(\log \frac{t_1}{a}\right)}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill =& \left|{\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{\left(\log \frac{t}{a}\right)}^p-\left(\frac{\log \frac{t}{a}}{h}{h}^p\right)\right]\frac{dt}{t}\right|\hfill \\ \hfill \le & {\int }_{t_0}^{t_1}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}{h}^p\frac{dt}{t}=h^p{h}^{\alpha }\left[\frac{{(k+1)}^{\alpha }-k^{\alpha }}{\alpha }\right]\hfill \\ \hfill =& h^{\alpha +p}\frac{\alpha {(k+1)}^{\alpha -1}}{\alpha }=h^{p+\alpha }{(k+1)}^{\alpha -1}=h^{p+\alpha }{(k+1)}^{\alpha -1}=C{h}^{p+1}{t}_{k+1}^{\alpha -1}.\hfill \end{array}$$

For $II$, we have,

$$\begin{array}{cc}\hfill \left|II\right|\le & C\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left|\frac{\log \frac{t}{t_j}\log \frac{t}{t_{j+1}}}{\log \frac{t_{j+1}}{t_j}}\right|\left[{\left(\log \frac{t_j}{a}\right)}^{p-1}-{\left(\log \frac{t_{j+1}}{a}\right)}^{p-1}\right]\frac{dt}{t}\hfill \\ \hfill \le & C{h}^p\sum _{j=1}^k{\int }_{t_j}^{t_{j+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[j^{p-1}-{(j+1)}^{p-1}\right]\frac{dt}{t}\hfill \\ \hfill =& C{h}^{p+\alpha }\left[\sum _{j=1}^{k-1}\left[j^{p-1}-{(j+1)}^{p-1}\right]\left[{(k+1-j)}^{\alpha }-{(k-j)}^{\alpha }\right]+\left[k^{p-1}-{(k+1)}^{p-1}\right]\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[\sum _{j=1}^{k-1}\left|(p-1)j^{p-2}\right|\left[\alpha {(k-j+1)}^{\alpha -1}\right]+k^{p-2}\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[\sum _{j=1}^{k-1}{j}^{p-2}{(k-j+1)}^{\alpha -1}+k^{p-2}\right].\hfill \end{array}$$

Similar to our previous proof, we now must find the bounds for the summation.

$$\begin{array}{cc}\hfill \sum _{j=1}^{k-1}{j}^{p-2}{(k-j+1)}^{\alpha -1}\le & {\int }_1^{x_p}{x}^{p-2}{(k+1-x)}^{\alpha -1}dx+{\int }_{x_p}^k{x}^{p-2}{(k+1-x)}^{\alpha -1}dx\hfill \\ \hfill =& {(k+1-x_p)}^{\alpha -1}\frac{1}{1-p}+x_p^{p-2}\frac{{(k+1-x_p)}^{\alpha }}{\alpha }\hfill \\ \hfill =& C_1{(k+1)}^{\alpha -1}+C_2{(k+1)}^{p-2}·C_3{(k+1)}^{\alpha }\hfill \\ \hfill \le & C{(k+1)}^{\alpha -1}+C{(k+1)}^{\alpha +p-2}.\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill \left|II\right|\le & C{h}^{\alpha +p}\left[\sum _{j=1}^{k-1}{j}^{p-2}{(k-j)}^{\alpha -1}+k^{p-2}\right]\le C{h}^{\alpha +p}\left[C{(k+1)}^{\alpha -1}+C{(k+1)}^{\alpha +p-2}+k^{p-2}\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[C{(k+1)}^{\alpha -1}+C{(k+1)}^{\alpha +p-2}+{(k+1)}^{\alpha +p-2}\right]\hfill \\ \hfill \le & C{h}^{\alpha +p}\left[C{(k+1)}^{\alpha -1}+C{(k+1)}^{\alpha +p-2}\right]\hfill \\ \hfill \le & C{h}^{p+1}\log {\left(\frac{t_{k+1}}{a}\right)}^{\alpha -1}+C{h}^2\log {\left(\frac{t_{k+1}}{a}\right)}^{\alpha +p-2}\le C{h}^{p+1}\log {\left(\frac{t_{k+1}}{a}\right)}^{\alpha -1}.\hfill \end{array}$$

Thus, we get the desired result. The proof for this theorem when $1<p<2$ is similar to when $0<p<1$. As such, it has been omitted. □

Remark 1.

In the final part of Theorem 4, there is a case when $\alpha ,p<1$. This would result in a $\varsigma =p+1$, and, in turn, this would make the right hand side exponent of $\log (t_{k+1}/a)$ become negative, taking a value of $\alpha -1$. This results in a case in which the error increases as the interval of integration decreases. The explanation for such a situation is that as the interval of integration decreases, so too does the weight function, making the integral more difficult to calculate and resulting in an increase in error.

3. Error Analysis for the Adams Method

In this section, we will be presenting the main error estimates for the Adams method for solving (1). We will be investigating different smoothness conditions for each of y and f and how they affect the error and convergence order.

3.1. A General Result

Using the error estimates established in the previous section, we will present the general convergence order of the Adams–Bashforth–Moulton method relating to the smoothness properties of the given function f and the solution y.

Lemma 1.

Assume that the solution y of (1) is such that,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)\right]\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}\left[{}_{CH}{D}_{a,t}^{\alpha }y\left(t_j\right)\right]|\le C_1{\left(\log \frac{t_{k+1}}{a}\right)}^{{\gamma }_1}{h}^{{\delta }_1},$$

(32)

and,

$$|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)\right]\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}\left[{}_{CH}{D}_{a,t}^{\alpha }y\left(t_j\right)\right]|\le C_1{\left(\log \frac{t_{k+1}}{a}\right)}^{{\gamma }_2}{h}^{{\delta }_2},$$

(33)

with some ${\gamma }_1,{\gamma }_2\ge 0$ and ${\delta }_1,{\delta }_2>0$. Let $y_j$ be the approximate solution of (1). Then, for some suitably chosen $T>1$, we have,

$$\underset{0\le j\le N}{max}|y\left(t_j\right)-y_j|=\mathcal{O}\left(h^q\right),$$

where $q=min\{{\delta }_1+\alpha ,{\delta }_2\}$.

Proof. 

We will show that, for sufficiently small h,

$$|y\left(t_j\right)-y_j|\le C{h}^q,$$

(34)

for all $j\in \{0,1,2,\dots ,N\}$, where C is a suitable constant. This proof will be based on mathematical induction. By using the initial conditions, we can confirm the basis step at $j=0$. We now assume that (34) is true for all $j=0,1,\dots ,k$ for some $k\le N-1$. Finally, we will prove that inequality (34) is true for $j=k+1$. To show this, we must start by finding the error of the predictor $y_{k+1}^P$. By the definition of the predictor, we can show the following:

$$\begin{array}{cc}\hfill \left|y\right(t_{k+1}-y_{k+1}^P|=& \frac{1}{\Gamma \left(\alpha \right)}|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}f(t,y\left(t\right))\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}f(t_j,y_j)y\left(t_j\right)|\hfill \\ \hfill \le & \frac{1}{\Gamma \left(\alpha \right)}|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}{}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)\frac{dt}{t}-\sum _{j=0}^k{b}_{j,k+1}{}_{CH}{D}_{a,t}^{\alpha }y\left(t_j\right)|\hfill \\ \hfill & -\frac{1}{\Gamma \left(\alpha \right)}\sum _{j=0}^k{b}_{j,k+1}|f\left(t_j,y\left(t_j\right)\right)-f(t_j,y_j)|\hfill \\ \hfill \le & \frac{C_1\log {\left(\frac{t_{k+1}}{a}\right)}^{{\gamma }_1}}{\Gamma \left(\alpha \right)}{h}^{{\gamma }_1}+\frac{1}{\Gamma \left(\alpha \right)}\sum _{j=0}^k{b}_{j,k+1}LC{h}^q\hfill \\ \hfill \le & \frac{C_1\log {\left(\frac{t_{k+1}}{a}\right)}^{{\gamma }_1}}{\Gamma \left(\alpha \right)}{h}^{{\gamma }_1}+\frac{CL{\left(\log \frac{T}{a}\right)}^{\alpha }}{\Gamma (\alpha +1)}{h}^q.\hfill \end{array}$$

For this proof, we have used several properties. These include the Lipschitz condition placed on f, the assumption of the error on the rectangle formula, and the understanding of the underlying predictor weights, $b_{j,k+1}>0$ for all j and k and,

$$\sum _{j=0}^k{b}_{j,k+1}\le \frac{1}{\alpha }{\left(\log \frac{T}{a}\right)}^{\alpha }.$$

(35)

Now, we have a bound for the predictor error. We also need to analyze the corrector error. To do so, we will be using a similar argument as with the predictor case as well as using the assumption made for mathematical induction. Note that,

$$\begin{array}{cc}\hfill & |y\left(t_{k+1}\right)-y_{k+1}|\hfill \\ \hfill =& \frac{1}{\Gamma \left(\alpha \right)}|{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}f\left(t,y\left(t\right)\right)\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}f(t_j,y_j)-a_{k+1,k+1}f(t_{k+1},y_{k+1}^P)|\hfill \\ \hfill \le & |{\int }_a^{t_{k+1}}{\left(\log \frac{t_{k+1}}{t}\right)}^{\alpha -1}\left[{}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)\right]\frac{dt}{t}-\sum _{j=0}^k{a}_{j,k+1}\left[{}_{CH}{D}_{a,t}^{\alpha }y\left(t_j\right)\right]|\hfill \\ \hfill +& \frac{1}{\Gamma \left(\alpha \right)}\sum _{j=0}^k{a}_{j,k+1}|f\left(t_j,y\left(t_j\right)\right)-f(t_j,y_j)|\hfill \\ \hfill +& \frac{1}{\Gamma \left(\alpha \right)}{a}_{k+1,k+1}|f\left(t_{k+1},y\left(t_{k+1}\right)\right)-f(t_{k+1},y_{k+1}^P)|\hfill \\ \hfill \le & \frac{C_2\log {\left(\frac{t_{k+1}}{a}\right)}^{{\gamma }_2}}{\Gamma \left(\alpha \right)}{h}^{{\delta }_2}+\frac{CL}{\Gamma \left(\alpha \right)}{h}^q\sum _{j=0}^k{a}_{j,k+1}+a_{k+1,k+1}\frac{L}{\Gamma \left(\alpha \right)}\left(\frac{C_1{\left(\log \frac{T}{a}\right)}^{{\gamma }_1}}{\Gamma \left(\alpha \right)}{h}^{{\gamma }_1}+\frac{CL{\left(\log \frac{T}{a}\right)}^{\alpha }}{\Gamma (\alpha +1)}{h}^q\right)\hfill \\ \hfill \le & \left(\frac{C_2{\left(\log \frac{T}{a}\right)}^{{\gamma }_2}}{\Gamma \left(\alpha \right)}+\frac{CL{\left(\log \frac{T}{a}\right)}^{\alpha }}{\Gamma (\alpha +1)}+\frac{C_{1}L{\left(\log \frac{T}{a}\right)}^{{\gamma }_1}}{\Gamma \left(\alpha \right)\Gamma (\alpha +2)}+\frac{C{L}^2{\left(\log \frac{T}{a}\right)}^{\alpha }}{\Gamma (\alpha +1)\Gamma (\alpha +2)}{h}^{\alpha }\right)h^q.\hfill \end{array}$$

Due to both ${\gamma }_1$ and ${\gamma }_2$ being non-negative and the relations ${\delta }_2\le q$ and ${\delta }_1+\alpha \le q$, we may choose a sufficiently small T such that the second summand in the parentheses is bounded above by $C/2$. Then, by fixing the value of T, we may bound the rest of the terms by $C/2$ as well, given a small enough value of h and by choosing a large C value. Finally, we can state that the error estimate for the corrector is now bounded above by $C{h}^q$. □

3.2. Error Estimates with Smoothness Assumptions on the Solution

In this subsection, we will be introducing some error estimates given certain smoothness assumptions being placed on the solution of our inital value problem. To do this, we will be using the general error estimate introduced above as well as the auxiliary results demonstrated in Section 2. For our first case, we will assume that y is sufficiently differentiable. This means the outcome is dependent on $\alpha$—more specifically, when $\alpha <1$ and when $\alpha \ge 1$.

Theorem 5.

Let $\alpha >0$ and assume ${}_{CH}{D}_{a,t}^{\alpha }y\in C^2[a,T]$ for some suitable T. Then,

$$\underset{0\le j\le N}{max}|y\left(t_j\right)-y_j|=\left\{\begin{array}{cc}\mathcal{O}\left(h^2\right)\hfill & \alpha \ge 1,\hfill \\ \mathcal{O}\left(h^{1+\alpha }\right)\hfill & \alpha <1.\hfill \end{array}\right.$$

(36)

We can note from this theorem that the order of convergence depends on $\alpha$. More specifically, a larger value of $\alpha$ gives a better order of convergence. The reason for this is due to the fact that we have chosen to discretize an integral operator that will behave more smoothly as the value of $\alpha$ increases, thus creating a higher order. We could have chosen to discretize the differential operator directly from the inital value problem. It can be shown that the opposite result occurs compared to the integral operator. As $\alpha$ increases, the smoothness of the operator deteriorates, and the convergence order diminishes. More specifically, it has been shown that when $\alpha \ge 2$, no convergence is achieved. This means that such a method is effective when $\alpha$ is small but will be ineffective for larger $\alpha$ values. This is one of the main advantages of this method for the Caputo–Hadamard derivative, as convergence is not only achieved but is most optimal at larger values of $\alpha$.

Proof. 

By applying Theorems 3 and 4 and Lemma 1 with ${\gamma }_1={\gamma }_2=\alpha >0$, ${\delta }_1=1$, and ${\delta }_2=2$, we can show that the an $\mathcal{O}\left(h^q\right)$ error bound where,

$$q=min\{1+\alpha ,2\}=\left\{\begin{array}{cc}2\hfill & \mathrm{if}\alpha \ge 1,\hfill \\ 1+\alpha \hfill & \mathrm{if}\alpha <1.\hfill \end{array}\right.$$

(37)

□

We can note from this theorem that this is for an optimal situation. The function we must approximate when applying the Adams method is the function $f(·,y(·))={}_{CH}{D}_{a,t}^{\alpha }y$. Therefore, the error bound is heavily dependent on obtaining a good approximations of f. To obtain such good error bounds, we can look to quadrature theory, which gives a well-known condition for obtaining this, namely the function $f\in C^2$ over the interval of the integral. As we can see from the above theorem, that is the case that we are looking at here. Therefore, this can be considered as an optimal situation. However, this can be also considered as an unusual situation. Rarely are we given enough information such that we can determine the smoothness of y or, in this case, ${}_{CH}{D}_{a,t}^{\alpha }y$, so we cannot rely on such a theorem. In general, we must formulate our theorems using the data given on the function f, which will be considered in the next subsection.

We next move on to giving some theorems that now give the smoothness assumptions based on the unknown solution y instead of ${}_{CH}{D}_{a,t}^{\alpha }y$. Theorem 2 implies that the smoothness of y often implies the nonsmoothness of ${}_{CH}{D}_{a,t}^{\alpha }y$; we are expecting some difficulties in finding the error estimates. However, Theorem 2 also gives us the form of ${}_{CH}{D}_{a,t}^{\alpha }y$ under such smoothness conditions and gives us information about the singularities and smoothness of ${}_{CH}{D}_{a,t}^{\alpha }y$. As such, we can use this information to give the following result.

Theorem 6.

Let $\alpha >1$ and assume that $y\in C^{1+\lceil \alpha \rceil }[a,T]$ for some suitable T. Then,

$$\underset{0\le j\le N}{max}|y\left(t_j\right)-y_j|=\mathcal{O}\left(h^{1+\lceil \alpha \rceil -\alpha }\right).$$

(38)

Proof. 

By applying Theorem 2, we have that ${}_{CH}{D}_{a,t}^{\alpha }y\left(t\right)=C{\left(\log \frac{t}{a}\right)}^{\lceil \alpha \rceil }+g\left(\log \frac{t}{a}\right)$ where $g\in C^1[a,T]$ and $\delta g$ follows the Lipschitz condition of order $\lceil \alpha \rceil -\alpha$. Therefore, by applying Theorems 3 and 4 and Lemma 1, we get that ${\gamma }_1=0,{\gamma }_2=\alpha -1$, ${\delta }_1=1$ and ${\delta }_2=1+\lceil \alpha \rceil -\alpha$. With $\alpha >1$, we can say that ${\delta }_1+\alpha =1+\alpha >2>{\delta }_2$. Therefore, $min\{{\delta }_1+\alpha ,{\delta }_2\}={\delta }_2$, so an error bound of $\mathcal{O}\left(h_2^{\delta }\right)$ is finally given. □

From the above theorem, we can draw some conclusions, namely that the fractional part of $\alpha$ plays a huge role in the overall convergence order. More specifically, as $\lceil \alpha \rceil -\alpha$ increases or the fractional part of $\alpha$ decreases, the convergence order improves. This means that the convergence order no longer forms a monotonic function of $\alpha$ but rather oscillates between 1 and 2 depending on $\alpha$. However, this theorem does show that under such smoothness properties, this method converges for all $\alpha >0$.

Theorem 7.

Let $0<\alpha <1$ and assume that $y\in C^2[a,T]$ for some suitable T. Then, for $1\le j\le N$,

$$|y\left(t_j\right)-y_j|\le C{\left(\log \frac{t_j}{a}\right)}^{\alpha -1}\times \left\{\begin{array}{cc}{h}^{1+\alpha }\hfill & \mathrm{if}0<\alpha <\frac{1}{2},\hfill \\ h^{2-\alpha }\hfill & \mathrm{if}\frac{1}{2}\le \alpha <1,\hfill \end{array}\right.$$

(39)

where C is a constant independent of j.

From this, we can take the following corollary.

Corollary 2.

Under the assumption of Theorem 7, we have,

$$\underset{0\le j\le N}{max}|y\left(t_j\right)-y_j|=\left\{\begin{array}{cc}\mathcal{O}\left(h^{2\alpha }\right)\hfill & \mathrm{if}0<\alpha <\frac{1}{2},\hfill \\ \mathcal{O}\left(h\right)\hfill & \mathrm{if}\frac{1}{2}\le \alpha <1.\hfill \end{array}\right.$$

(40)

Moreover, for every $ϵ\in (a,T)$, we have,

$$\underset{t_j\in [ϵ,T]}{max}|y\left(t_j\right)-y_j|=\left\{\begin{array}{cc}\mathcal{O}\left(h^{1+\alpha }\right)\hfill & \mathrm{if}0<\alpha <\frac{1}{2},\hfill \\ \mathcal{O}\left(h^{2-\alpha }\right)\hfill & \mathrm{if}\frac{1}{2}\le \alpha <1.\hfill \end{array}\right.$$

(41)

Proof of Theorem 7.

For this proof, we would be following that of Theorem 6. However, as $0<\alpha <1$, we have that ${\gamma }_2=\alpha -1<0$, meaning that we can no longer apply Lemma 1. We will need to adapt the proof of this lemma in order to apply it to this case. To do so, we will retain the main structure of the proof and application of mathematical induction but will now change the induction hypothesis to that of (39). With such a change in hypothesis, we must now find estimates for the terms ${\sum }_{j=1}^{k-1}{b}_{j,k+1}{\left(\log \frac{t_j}{a}\right)}^{{\delta }_2}$ and ${\sum }_{j=1}^{k-1}{a}_{j,k+1}{\left(\log \frac{t_j}{a}\right)}^{{\delta }_2}$. By applying the Mean Value Theorem, it is known that $0\le b_{j,k+1}\le h^{\alpha }{(k-j)}^{\alpha -1}$ and $0\le a_{j,k+1}\le c{h}^{\alpha }{(k-j)}^{\alpha -1}$ for $1\le j\le k-1$ and c is independent of j and k. Applying these bounds for the weights, we reduce the problem to finding bounds for the sum $S_k:={\sum }_{j=1}^{k-1}{j}^{{\gamma }_2}{(k-j)}^{\alpha -1}$. Looking back at Theorems 3 and 4, we have shown similar results, and it is easily shown that $S_k=C{k}^{{\gamma }_2+\alpha }$ when both the indices ${\gamma }_2$ and $\alpha -1$ are in the interval $(0,1)$. By applying this, we can complete this proof using structure of Lemma 1. □

3.3. Error Estimates with Smoothness Assumptions on the Given Data

In this final subsection, we will be looking at how changing the smoothness assumptions of the given function f can change the error and convergence order for this method. We will be looking at both when $\alpha <1$ and when $\alpha >1$.

Theorem 8.

Let $\alpha >1$. Then, with $f\in C^3\left(G\right)$,

$$\underset{0\le j\le N}{max}|y\left(t_j\right)-y_j|=\mathcal{O}\left(h^2\right).$$

(42)

Proof. 

We will split this proof into when $\alpha \ge 2$ an when $1<\alpha <2$. When $\alpha \ge 2$, we can adapt a result from Miller and Feldstein [26] to apply here which shows that $y\in C^2[a,T]$. Then, given the smoothness assumpions on f, and applying the chain rule, ${}_{CH}{D}_{a,t}^{\alpha }y\in C^2[a,T]$. This then fulfills the conditions of Theorem 5, which gives the desired result.

Now, for when $1<\alpha <2$, we wish to apply Lemma 1. To do this, we must find the constants ${\gamma }_1,{\gamma }_2,{\delta }_1$ and ${\delta }_2$ in the lemma. As in the case of our previous theorems, we must determine the behavior and smoothness of y. We find this information by applying Theorem 1b, which tells us that y is in the form,

$$y\left(t\right)=\psi \left(t\right)+\sum _{v=1}^{\widehat{v}}{c}_v{\left(\log \frac{t}{a}\right)}^{v\alpha }+\sum _{v=1}^{\tilde{v}}{d}_v{\left(\log \frac{t}{a}\right)}^{1+v\alpha }.$$

(43)

As $\psi \in C^2[a,T]$, this implies that $y\in C^1[a,T]$. Similar to the case of $\alpha \ge 2$, we can deduce ${}_{CH}{D}_{a,t}^{\alpha }y\in C^1[a,T]$. This then fulfills the conditions of Theorem 3a, giving us that ${\gamma }_1=\alpha$ and ${\delta }_1=1$. Furthermore, we may apply Theorem 4a,c to find the remaining values such that ${\gamma }_2=min\{\alpha ,2\alpha -2\}=2\alpha -2\ge 0$ and ${\delta }_2=2$. By applying Lemma 1, the required result is achieved. □

4. Numerical Examples

In this section, we will be considering some numerical examples to confirm the theoretical results presented in the previous sections. We will be presenting examples below with $0<\alpha <2$ as the applications for $\alpha \ge 2$ are not currently of interest. We will be solving the following examples for $a=1$ and $T=2$. The following examples and graphs were created in blueMATLAB using the following algorithm.

• Give an initial value $y_0$.
• Find $y_1^P$ using (11) and (12) such that,

  $$y_1^P=y_0+b_{0,1}f(t_0,y_0),$$

  (44)

  where,

  $$b_{0,1}=\frac{1}{\Gamma (\alpha +1)}\left[{\left(\log \frac{t_1}{t_0}\right)}^{\alpha }-{\left(\log \frac{t_1}{t_1}\right)}^{\alpha }\right]=\frac{1}{\Gamma (\alpha +1)}{\left(\log \frac{t_1}{a}\right)}^{\alpha }.$$

  (45)
• Find $y_1$ using (8) and (9) such that,

  $$y_1=y_0+\frac{1}{\Gamma \left(\alpha \right)}\left(a_{0,1}f(t_0,y_0)+a_{1,1}f(t_1,y_1^P)\right).$$

  (46)
• Repeat steps 2–3 to find $y_2,y_3,\dots ,y_N$.

Example 1

(${}_{CH}{D}_{a,t}^{\alpha }y$ is smooth). Consider the following nonlinear fractional differential equation, with $0<\alpha <2$,

$$\left\{\begin{array}{c}{{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)=f(t,y),1\le a<t\le T,\hfill \\ y\left(a\right)=0,\delta y\left(a\right)=0,\hfill \end{array}\right.$$

(47)

where,

$$\begin{array}{cc}\hfill f(t,y)& =\frac{\Gamma \left(9\right)}{\Gamma (9-\alpha )}{\left(\log \frac{t}{a}\right)}^{8-\alpha }-3\frac{\Gamma (5+\alpha /2)}{\Gamma (5-\alpha /2)}{\left(\log \frac{t}{a}\right)}^{4-\alpha /2}+\frac{9}{4}\Gamma (\alpha +1)\hfill \\ \hfill & +{\left(\frac{3}{2}{\left(\log \frac{t}{a}\right)}^{\alpha /2}-{\left(\log \frac{t}{a}\right)}^4\right)}^3-y^{3/2}.\hfill \end{array}$$

(48)

This example has a nonlinear and non smooth f. However, due to the form of this equation, it is well known that the solution y is given as,

$$y\left(t\right)={\left(\log \frac{t}{a}\right)}^8-3{\left(\log \frac{t}{a}\right)}^{4+\alpha /2}+\frac{9}{4}{\left(\log \frac{t}{a}\right)}^{\alpha }.$$

(49)

As such we can say,

$${{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)=\frac{\Gamma \left(9\right)}{\Gamma (9-\alpha )}{\left(\log \frac{t}{a}\right)}^{8-\alpha }-\frac{\Gamma (5+\alpha /2)}{\Gamma (5-\alpha /2)}{\left(\log \frac{t}{a}\right)}^{4-\alpha /2}+\frac{9}{4}\Gamma (\alpha +1).$$

(50)

For $\alpha \le 4$, we have that ${{}_{CH}D}_{a,t}^{\alpha }y\in C^2[a,T]$; this fulfills the requirements for Theorem 5. As such, we can show the theorem holds for such an example. Let N be a positive integer and let $\log a=\log t_0<\log t_1<\dots <\log t_N=\log T$ be the uniform mesh on the interval $[\log a,\log T]$. such that $\log t_j=\log a+jh$ for $j=0,1,2,\dots ,N$ and $h=\left(\log \frac{T}{a}\right)/N$. Therefore, we have by Theorem 5,

$$\left|\right|e_N{\left|\right|}_{\infty }=\underset{0\le j\le N}{max}|y\left(t_j\right)-y_j|=\left\{\begin{array}{cc}\mathcal{O}\left(h^2\right)\hfill & \alpha \ge 1,\hfill \\ \mathcal{O}\left(h^{1+\alpha }\right)\hfill & \alpha <1.\hfill \end{array}\right.$$

(51)

In

Table 1. Table showing the maximum absolute error and EOC for solving (47) using the predictor–corrector method for $0<\alpha <1$.

N	$\mathit{\alpha }=0.4$	EOC	$\mathit{\alpha }=0.6$	EOC	$\mathit{\alpha }=0.8$	EOC
10	6.839 × 10${}^{-2}$		2.607 × 10${}^{-2}$		1.350 × 10${}^{-2}$
20	2.064 × 10${}^{-2}$	1.728	7.433 × 10${}^{-3}$	1.810	3.542 × 10${}^{-3}$	1.930
40	6.576 × 10${}^{-3}$	1.650	2.209 × 10${}^{-3}$	1.751	9.534 × 10${}^{-4}$	1.894
80	2.193 × 10${}^{-3}$	1.584	6.768 × 10${}^{-4}$	1.706	2.610 × 10${}^{-4}$	1.869
160	7.561 × 10${}^{-4}$	1.536	2.120 × 10${}^{-4}$	1.675	7.225 × 10${}^{-5}$	1.853
320	2.671 × 10${}^{-4}$	1.502	6.744 × 10${}^{-5}$	1.653	2.015 × 10${}^{-5}$	1.842
640	9.600 × 10${}^{-5}$	1.476	2.169 × 10${}^{-5}$	1.637	5.652 × 10${}^{-6}$	1.834
1280	3.497 × 10${}^{-5}$	1.457	7.026 × 10${}^{-6}$	1.626	1.591 × 10${}^{-6}$	1.829

and

Table 2. Table showing the maximum absolute error and EOC for solving (47) using the predictor–corrector method for $\alpha >1$.

N	$\mathit{\alpha }=1.25$	EOC	$\mathit{\alpha }=1.50$	EOC	$\mathit{\alpha }=1.75$	EOC
10	6.166 × 10${}^{-3}$		5.475 × 10${}^{-3}$		5.369 × 10${}^{-3}$
20	1.463 × 10${}^{-3}$	2.076	1.319 × 10${}^{-3}$	2.053	1.319 × 10${}^{-3}$	2.026
40	3.497 × 10${}^{-4}$	2.065	3.207 × 10${}^{-4}$	2.040	3.259 × 10${}^{-4}$	2.017
80	8.408 × 10${}^{-5}$	2.056	7.853 × 10${}^{-5}$	2.030	8.091 × 10${}^{-5}$	2.010
160	2.033 × 10${}^{-5}$	2.049	1.934 × 10${}^{-5}$	2.022	2.014 × 10${}^{-5}$	2.006
320	4.935 × 10${}^{-6}$	2.042	4.783 × 10${}^{-6}$	2.016	5.022 × 10${}^{-6}$	2.004
640	1.203 × 10${}^{-6}$	2.036	1.187 × 10${}^{-6}$	2.011	1.253 × 10${}^{-6}$	2.002
1280	2.943 × 10${}^{-7}$	2.031	2.951 × 10${}^{-7}$	2.008	3.131 × 10${}^{-7}$	2.001

, we can see the maximum absolute error and experimental order of convergence (EOC) for the predictor–corrector method at varying $\alpha$ and N values. For our different $0<\alpha <2$, we have chosen N values as $N=10\times 2^l,l=0,1,2,\dots ,7$.The maximum absolute errors $\left|\right|e_N{\left|\right|}_{\infty }$ were obtained as shown above with respect to N, and we calculate the experimental order of convergence or EOC as $\log \left(\frac{\left|\right|e_N{\left|\right|}_{\infty }}{\left|\right|e_{2N}{\left|\right|}_{\infty }}\right).$ As we can see, the EOCs for the above example follow that of Theorem 5.

In Figure 1, we have plotted the order of convergence for Example 1. From Equation (51), we have for $\alpha <1$,

$$\left({\log }_2\left|\right|e_N\left|\right|\right)\le \left({\log }_{2}C\right)+\left({\log }_2{h}^{(1+\alpha )}\right)\le \left({\log }_{2}C\right)+\left(1+\alpha \right)\left({\log }_{2}h\right).$$

(52)

We can now plot this graph such that $y=\left({\log }_2\left|\right|e_N\left|\right|\right)$ and let $x=\left({\log }_{2}h\right)$ and $h=\frac{1}{5\times 2^l},l=0,1,\dots ,7$. Doing this, we get that the gradient of the graph would equal the EOC. To compare this to the theoretical order of convergence, we have also plotted the straight line $y=(1+\alpha )x$. For Figure 1, we use $\alpha =0.8$. We can observe that the two lines drawn are parallel. Therefore, we can conclude that the order of convergence of this predictor–corrector method is $O\left(h^{1+\alpha }\right)$ for when $\alpha <1$. A similar result can be obtained for when $\alpha >1$. Figure 2 shows the same graph but for $\alpha =1.75$. However, now we can see that the line is parallel to the straight line $y=2x$, which is what we expected as the convergence order should tend to 2 for $\alpha >1$.

Example 2

(y is smooth). Consider the following nonlinear fractional differential equation, with $0<\alpha <2$,

$$\left\{\begin{array}{c}{{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)=f(t,y),1\le a<t\le T,\hfill \\ y\left(a\right)=0,\delta y\left(a\right)=1,\hfill \end{array}\right.$$

(53)

where,

$$\begin{array}{cc}\hfill f(t,y)& =\left\{\begin{array}{cc}\frac{2}{\Gamma (3-\alpha )}{\left(\log \frac{t}{a}\right)}^{2-\alpha }-y+{\left(\log \frac{t}{a}\right)}^2-\log \frac{t}{a},\hfill & \mathrm{for}\alpha >1,\hfill \\ \frac{2}{\Gamma (3-\alpha )}{\left(\log \frac{t}{a}\right)}^{2-\alpha }-\frac{1}{\Gamma (2-\alpha )}{\left(\log \frac{t}{a}\right)}^{1-\alpha }-y+{\left(\log \frac{t}{a}\right)}^2-\log \frac{t}{a},\hfill & \mathrm{for}\alpha \le 1.\hfill \end{array}\right.\hfill \end{array}$$

(54)

The exact solution of this equation is,

$$y\left(t\right)={\left(\log \frac{t}{a}\right)}^2-\log \frac{t}{a}.$$

(55)

In

Table 3. Table showing the maximum absolute error and EOC for solving (53) using the predictor–corrector method for $0<\alpha \le 0.5$.

N	$\mathit{\alpha }=0.1$	EOC	$\mathit{\alpha }=0.3$	EOC	$\mathit{\alpha }=0.5$	EOC
10	2.225 × 10${}^{-2}$		9.375 × 10${}^{-3}$		5.123 × 10${}^{-3}$
20	1.261 × 10${}^{-2}$	0.819	3.345 × 10${}^{-3}$	1.487	1.622 × 10${}^{-3}$	1.660
40	5.625 × 10${}^{-3}$	1.164	1.196 × 10${}^{-3}$	1.484	5.255 × 10${}^{-4}$	1.626
80	2.387 × 10${}^{-3}$	1.237	4.343 × 10${}^{-4}$	1.461	1.739 × 10${}^{-4}$	1.595
160	1.004 × 10${}^{-3}$	1.249	1.606 × 10${}^{-4}$	1.436	5.856 × 10${}^{-5}$	1.571
320	4.242 × 10${}^{-4}$	1.243	6.033 × 10${}^{-5}$	1.412	1.998 × 10${}^{-5}$	1.552
640	1.805 × 10${}^{-4}$	1.233	2.299 × 10${}^{-5}$	1.392	6.884 × 10${}^{-6}$	1.537
1280	7.742 × 10${}^{-5}$	1.221	8.864 × 10${}^{-6}$	1.375	2.389 × 10${}^{-6}$	1.527

,

Table 4. Table showing the maximum absolute error and EOC for solving (53) using the predictor–corrector method for $0.5<\alpha <1$.

N	$\mathit{\alpha }=0.7$	EOC	$\mathit{\alpha }=0.9$	EOC
10	5.507 × 10${}^{-3}$		1.162 × 10${}^{-2}$
20	1.931 × 10${}^{-3}$	1.512	5.100 × 10${}^{-3}$	1.188
40	7.031 × 10${}^{-4}$	1.457	2.300 × 10${}^{-3}$	1.151
80	2.635 × 10${}^{-4}$	1.416	1.051 × 10${}^{-3}$	1.129
160	1.008 × 10${}^{-4}$	1.386	4.847 × 10${}^{-4}$	1.116
320	3.920 × 10${}^{-5}$	1.363	2.247 × 10${}^{-4}$	1.109
640	1.541 × 10${}^{-5}$	1.347	1.045 × 10${}^{-4}$	1.105
1280	6.110 × 10${}^{-6}$	1.335	4.863 × 10${}^{-5}$	1.103

and

Table 5. Table showing the maximum absolute error and EOC for solving (53) using the predictor–corrector method for $\alpha >1$.

N	$\mathit{\alpha }=1.25$	EOC	$\mathit{\alpha }=1.50$	EOC	$\mathit{\alpha }=1.75$	EOC
10	6.166 × 10${}^{-3}$		5.475 × 10${}^{-3}$		5.369 × 10${}^{-3}$
20	1.463 × 10${}^{-3}$	2.076	1.319 × 10${}^{-3}$	2.053	1.319 × 10${}^{-3}$	2.026
40	3.497 × 10${}^{-4}$	2.065	3.207 × 10${}^{-4}$	2.040	3.259 × 10${}^{-4}$	2.017
80	8.408 × 10${}^{-5}$	2.056	7.853 × 10${}^{-5}$	2.030	8.091 × 10${}^{-5}$	2.010
160	2.033 × 10${}^{-5}$	2.049	1.934 × 10${}^{-5}$	2.022	2.014 × 10${}^{-5}$	2.006
320	4.935 × 10${}^{-6}$	2.042	4.783 × 10${}^{-6}$	2.016	5.022 × 10${}^{-6}$	2.004
640	1.203 × 10${}^{-6}$	2.036	1.187 × 10${}^{-6}$	2.011	1.253 × 10${}^{-6}$	2.002
1280	2.943 × 10${}^{-7}$	2.031	2.951 × 10${}^{-7}$	2.008	3.131 × 10${}^{-7}$	2.001

, we can see the maximum absolute error and experimental order of convergence (EOC) for the predictor–corrector method at varying $\alpha$ and N values. For our different $0<\alpha <2$, we have chosen N values as $N=10\times 2^l,l=0,1,2,\dots ,7$. As $y\in C^2[a,T]$ we can apply Theorems 6 and 7 and more specifically, Corollary 2. As we can see, the EOCs for the above example verifies these theorems. When $\alpha \le 0.5$, we find that the EOC is close to $1+\alpha$ and when $0.5<\alpha <1$, we find that the EOC is close to $2-\alpha$. Finally, when $\alpha >1$, we have that the EOC is close to $1+\lceil \alpha \rceil -\alpha$.

Example 3

(f is smooth). Consider the following nonlinear fractional differential equation, with $0<\alpha <2$,

$$\left\{\begin{array}{c}{{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)=-y\left(t\right),1\le a<t\le T,\hfill \\ y\left(a\right)=1,\delta y\left(a\right)=0.\hfill \end{array}\right.$$

(56)

The exact solution of this equation is $y\left(t\right)=E_{\alpha ,1}\left(-{\left(\log t\right)}^{\alpha }\right)$. Therefore, ${{}_{CH}D}_{a,t}^{\alpha }y\left(t\right)=-E_{\alpha ,1}\left(-{\left(\log t\right)}^{\alpha }\right)$, where $E_{\alpha ,\gamma }\left(z\right)$ is defined as the Mittag–Leffler function,

$$E_{\alpha ,\gamma }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\gamma )},\alpha ,\gamma >0.$$

(57)

Therefore, f is smooth. As such, this equation fulfils the conditions of Theorem 8. In

Table 6. Table showing the maximum absolute error and EOC for solving Example 3 using the predictor–corrector method for $\alpha <1$.

N	$\mathit{\alpha }=0.3$	EOC	$\mathit{\alpha }=0.6$	EOC	$\mathit{\alpha }=0.9$	EOC
10	1.353 × 10${}^{-3}$		6.520 × 10${}^{-4}$		3.414 × 10${}^{-4}$
20	4.324 × 10${}^{-4}$	1.646	1.937 × 10${}^{-4}$	1.751	8.704 × 10${}^{-5}$	1.972
40	1.466 × 10${}^{-4}$	1.560	6.019 × 10${}^{-5}$	1.686	2.270 × 10${}^{-5}$	1.939
80	5.159 × 10${}^{-5}$	1.507	1.919 × 10${}^{-5}$	1.649	5.990 × 10${}^{-6}$	1.922
160	1.863 × 10${}^{-5}$	1.470	6.211 × 10${}^{-6}$	1.628	1.591 × 10${}^{-6}$	1.913
320	6.863 × 10${}^{-6}$	1.441	2.027 × 10${}^{-6}$	1.615	4.239 × 10${}^{-7}$	1.908
640	2.570 × 10${}^{-6}$	1.417	6.652 × 10${}^{-7}$	1.608	1.132 × 10${}^{-7}$	1.905
1280	9.752 × 10${}^{-7}$	1.398	2.189 × 10${}^{-7}$	1.604	3.026 × 10${}^{-8}$	1.904

and

Table 7. Table showing the maximum absolute error and EOC for solving Example 3 using the predictor–corrector method for $\alpha >1$.

N	$\mathit{\alpha }=1.25$	EOC	$\mathit{\alpha }=1.50$	EOC	$\mathit{\alpha }=1.75$	EOC
10	2.457 × 10${}^{-4}$		2.108 × 10${}^{-4}$		1.507 × 10${}^{-4}$
20	5.847 × 10${}^{-5}$	2.071	5.141 × 10${}^{-5}$	2.036	3.723 × 10${}^{-5}$	2.017
40	1.400 × 10${}^{-5}$	2.062	1.261 × 10${}^{-5}$	2.027	9.236 × 10${}^{-6}$	2.011
80	3.374 × 10${}^{-6}$	2.053	3.109 × 10${}^{-6}$	2.020	2.298 × 10${}^{-6}$	2.007
160	8.213 × 10${}^{-7}$	2.039	7.693 × 10${}^{-7}$	2.015	5.728 × 10${}^{-7}$	2.004
320	2.016 × 10${}^{-7}$	2.027	1.909 × 10${}^{-7}$	2.011	1.430 × 10${}^{-7}$	2.003
640	4.975 × 10${}^{-8}$	2.019	4.748 × 10${}^{-8}$	2.008	3.570 × 10${}^{-8}$	2.002
1280	1.233 × 10${}^{-8}$	2.013	1.183 × 10${}^{-8}$	2.005	8.919 × 10${}^{-9}$	2.001

, we can see the maximum absolute error and experimental order of convergence (EOC) for the predictor–corrector method at varying $\alpha$ and N values. For our different $0<\alpha <2$, we have chosen N values as $N=10\times 2^l,l=0,1,2,\dots ,7$. As f is sufficiently smooth, we can apply Theorem 8. As we can see, the EOCs for the above example verifies these theorems. When $0<\alpha <1$, we have that the EOC is close to $1+\alpha$. For $\alpha >1$, we have the EOC is close to 2.

5. Conclusions

In this paper, we proposed a predictor–corrector method for solving Caputo–Hadamard fractional differential equations. Both the initial data f and the unknown solution y were investigated to see how the different smoothness conditions affected the convergence order. It was found under optimal conditions and with $\alpha \ge 1$ that we can obtain a theoretical convergence order of 2. However, under certain smoothness conditions, a suboptimal convergence order is achieved, often depending on the fractional part of $\alpha$. Several numerical simulations are given to support the theoretical results obtained above in terms of the error estimates.