Modified Three-Point Fractional Formulas with Richardson Extrapolation

Abstract

In this paper, we introduce new three-point fractional formulas which represent three generalizations for the well-known classical three-point formulas; central, forward and backward formulas. This has enabled us to study the function’s behavior according to different fractional-order values of $\alpha$ numerically. Accordingly, we then introduce a new methodology for Richardson extrapolation depending on the fractional central formula in order to obtain a high accuracy for the gained approximations. We compare the efficiency of the proposed methods by using tables and figures to show their reliability.

1. Introduction

The differential equations with non-integer order derivatives are used to represent the dynamical system, called the fractional-order system. In general, fractional-order dynamics (FoD) is used to describe most practical systems; we can consider a special case for the system called the integer-order representation. We recommend the article in [1,2,3], since it discusses several fractional-order dynamical system problems. In 1675, Issac Newton found the interpolation formula, which was subsequently named after J. L. Lagrange (1736–1813), but it first appeared in 1779, presented by Edward Waring (1736–1798). The interpolation method was mainly studied and written by Lagrange, and his works had a substantial influence on mathematicians in the future. This result was published in 1795. The Lagrange polynomials were studied by Axelsson, O. and V. A. Barker in [4], where they found an explicit representation for it and the error when approximating a function on an interval. These polynomials are frequently used for interpolating tabular data. Note that, in similar cases, the explicit representation of the polynomial is not essential, only the values of the polynomial at specified points are required. Additionally, if the data function is unknown, we cannot use the explicit form. In 1932, E. H. Neville gave in his paper a modification of the Lagrange formula, Neville’s Method. The application of the error term to Lagrange interpolation poses a practical challenge; because of this, it is typically unknown what degree of the polynomial is required for the necessary precision until computations have been made. It is customary to compute the results using different polynomials until the desired agreement is attained. As in many areas, Isaac Newton is prominent in studying differential equations. In 1675, he developed interpolation formulas through means using his notation $\mathsf{\Delta }$ in tables of differences. He produced explicit examples by taking a very general approach to the different formulas, including Lagrange’s formulas, which are called by several other names. Consider that $f\in C^n[a,b]$ and $x_0,x_1,\cdots ,x_n$ are distinct numbers in [a,b]. Then, a number $\xi$ exists in (a,b) with

$$f[x_0,x_1,\cdots ,x_n]=\frac{f^{\left(n\right)}\left(\xi \right)}{n!}.$$

A. Aitken applied the following method in 1926 to accelerate a series rate of convergence in his paper [5], Aitken’s method is based on the assumption that the sequence ${\left\{{\widehat{P}}_n\right\}}_{n=0}^{\infty }$, defined by

$${\widehat{P}}_n=P_n-\frac{{(P_{n+1}-P_n)}^2}{P_{n+2}-2P_{n+1}+P_n},$$

converges more rapidly to P than does the original sequence ${\left\{P_n\right\}}_{n=0}^{\infty }$; the procedure for this method is similar to that used by the Japanese mathematician T. S. Kowa. The forward difference notation $Delta$ developed by Aitken’s $Delt{a}^2$ approach was used to build the Newton forward-difference formula. With this notation,

$$f[x_0,x_1]=\frac{f\left(x_1\right)-f\left(x_0\right)}{x_1-x_0}=\frac{1}{h}(f\left(x_1\right)-f\left(x_0\right))=\frac{1}{h}\mathsf{\Delta }f\left(x_0\right),$$

$$f[x_0,x_1,x_2]=\frac{1}{2h}\left[\frac{\mathsf{\Delta }f\left(x_1\right)-\mathsf{\Delta }f\left(x_0\right)}{h}\right]=\frac{1}{2h^2}{\mathsf{\Delta }}^{2}f\left(x_0\right),$$

and, in general,

$$f[x_0,x_1,\cdots ,x_k]=\frac{1}{k!h^k}{\mathsf{\Delta }}^{k}f\left(x_0\right).$$

Then, the Newton forward-difference formula is given by the following:

$$P_n\left(x\right)=f\left(x_0\right)+\sum _{k=1}^n\left(\genfrac{}{}{0pt}{}{s}{k}\right){\mathsf{\Delta }}^{k}f\left(x_0\right).$$

(1)

The Newton backward-difference formula is given by the following:

$$P_n\left(x\right)=f\left[x_n\right]+\sum _{k=1}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{-s}{k}\right){\nabla }^{k}f\left(x_n\right).$$

(2)

Note that if we take the function $f\left(x\right)$ when x lies near the center of the table, then the Newton backward-difference and forward-difference formulas are inappropriate because neither will allow $x_0$ to be close to x for the highest-order difference. Several divided-difference formulas can be employed in this case, and each one has specific circumstances where it works well. These procedures are referred to as “centered-difference formulas”. We focus on Stirling’s method and consider only one centered difference formula by choosing $x_0$ near the point being approximated and by labeling the nodes directly below $x_0$ as $x_1,x_2,\cdots$ and those directly above as $x_{-1},x_{-2},\cdots$. With this convention, Stirling’s formula is given by the following [6]:

$$\begin{array}{cc}\hfill P_n\left(x\right)=& P_{2m+1}\left(x\right)=f\left[x_0\right]+\frac{sh}{2}(f[x_{-1},x_0]+f[x_0,x_1])\hfill \\ \hfill & +s^2{h}^{2}f[x_{-1},x_0,x_1]{\displaystyle \frac{s(s^2-1)h^3}{2}}(f[x_{-2},x_{-1},x_0,x_1]+f[x_{-1},x_0,x_1,x_2])+\cdots \hfill \\ \hfill & +s^2(s^2-1)(s^2-4)\cdots (s^2-{(m-1)}^2)h^{2m}f[x_{-m},\cdots ,x_m]\hfill \\ \hfill & +{\displaystyle \frac{s(s^2-1)\cdots (s^2-m^2)h^{2m+1}}{2}}(f[x_{-m-1},\cdots ,x_m]+f[x_{-m},\cdots ,x_{m+1}]),\hfill \end{array}$$

(3)

if $n=2m+1$ is odd. If $n=2m$ is even, we use the same formula but delete the last line. In the last quarter of the 17th century, Isaac Newton used and popularized the difference equations, and T. Harriot and H. Briggs developed many of these techniques. The navigation techniques were studied by T. Harriot and saw significant advances. H. Briggs was the most responsible for accepting logarithms as an aid to computation. One of the classical techniques is Richardson extrapolation (RE); see [1,2], which depends on a small parameter to accelerate the convergence of numerical sequences by eliminating the lowest order error term $\left(s\right)$ from the corresponding asymptotic expansion. When the sequence is generated by a numerical method, resolving the initial-value problem

$$y^{\prime }\left(t\right)=f(t,y\left(t\right)),y\left(t_0\right)=y_0,$$

the discretization step size $h>0$ can be selected as the parameter in RE. The application of RE to one-step sequences, such as Runge–Kutta methods, are described, for example, in [7,8]. In [9], global (also known as passive) or local (active) versions of RE are implemented with Runge–Kutta sequences. These integrated methods can be used to solve air pollution issues [10] or to improve machine learning [11]. For an overview of real systems (e.g., physical and biological systems) where fractional derivatives are suitable to describe their dynamics, see [12,13,14]; in the article [15], in order to describe the behavior of porous electrodes, the authors presented a method for electrochemical capacitor quality management using fractional order models. For more applications on fractional order differential equations see [16,17,18,19,20,21].

The organization of this paper is arranged as follows: Section 2 aims to recall some basic facts and definitions connected with fractional calculus. Section 3 demonstrates the main results of this work so that it contains the definition of the modified three-point fractional formulas. Section 4 develops the classical Richardson extrapolation to be valid for fractional calculus. Section 5 illustrates numerical results that confirm the theoretical findings of this work, followed by the final section that summarizes the conclusion.

2. Preliminaries

In this section, basic definitions and theorems such as the Riemann–Liouville integral and derivative, the Caputo derivative, etc. are introduced with their references.

Definition 1

([22]). The fractional integral of a function f(Y,t) ∈$C_{\delta }$$(\delta ⩽-1)$ and of order $\mu >0$ initially defined by Riemann–Liouville is presented as follows:

$$\begin{array}{cc}\hfill & J^{\mu }f(Y,t)=\frac{1}{\mathsf{\Gamma }\left(\mu \right)}{\int }_0^{t}f(Y,\tau ){(t-\tau )}^{\mu -1}d\tau ,t>0,\mu >0,\hfill \end{array}$$

(4)

where $Y=(x_1,x_2,x_3,\dots ,x_n)$ is a vector of n variables.

Some of the properties of the Riemann–Liouville integral are as follows:

$$\begin{array}{cc}\hfill \left(1\right)& J^{0}f(Y,t)=f(Y,t).\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \left(2\right)& J^{\mu }{t}^{\gamma }=\frac{\mathsf{\Gamma }(\gamma +1)t^{\alpha +\gamma }}{\mathsf{\Gamma }(\alpha +\gamma +1)}\gamma \ge -1.\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \left(3\right)& J^{\mu }{J}^{\beta }f(Y,t)=J^{\beta }{J}^{\mu }f(Y,t)\alpha ,\beta \ge 0.\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \left(4\right)& J^{\mu }{J}^{\beta }f(Y,t)=J^{\mu +\beta }f(Y,t)\alpha ,\beta \ge 0.\hfill \end{array}$$

(8)

Definition 2

([15]). The Caputo time fractional derivative of order $\alpha >0$ with m being a small number greater than α is defined as follows:

$$\begin{array}{c}\hfill D_t^{\alpha }f(Y,t)=\left\{\begin{array}{cc}\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^t{(t-\tau )}^{-\alpha +n-1}\frac{{\partial }^{m}f(Y,\tau )}{\partial {\tau }^m}d\tau \hfill & ,m-1<\alpha <m\hfill \\ \frac{{\partial }^{m}f(Y,t)}{\partial t^m}\hfill & ,\alpha =m\hfill \end{array}\right.\end{array}$$

(9)

where $Y=(x_1,x_2,x_3,\dots ,x_n)$ is a vector of n variables and for $n-1<\alpha \le n$, $n\in \mathit{N}$, $t>0$, f(Y,t) is a real value function while the parameter α is the order of the derivative.

Some of the characteristics of the Caputo derivative are listed below:

(1) $D_t^{\alpha }c=0$, where c is constant.

(2)

$$D_t^{\alpha }{t}^{\rho }=\left\{\begin{array}{cc}\frac{\mathsf{\Gamma }(\rho +1)}{\mathsf{\Gamma }(\rho -\alpha +1)}{t}^{\rho -\alpha }\hfill & :\rho >\alpha -1\hfill \\ 0\hfill & :Otherwise\hfill \end{array}\right.$$

(3) $D_t^{\alpha }$ is linear, i.e.,

$$\begin{array}{c}\hfill D_t^{\alpha }(\mu f(Y,t)+\omega k(Y,t))=\mu D_t^{\alpha }\left(f(Y,t)\right)+\omega D_t^{\alpha }\left(k\left(t\right)\right),\end{array}$$

where $\mu$ and $\omega$ are constant.

Additionally, we need two of its basic properties. If $n-1<\alpha \le n$, $n\in N$, then

$$D^{\alpha }{J}^{\alpha }f\left(x\right)=f\left(x\right),J^{\alpha }{D}^{\alpha }f\left(x\right)=f\left(x\right)-\sum _{i=1}^n{f}^i\left(0^{+}\right)\frac{x^i}{i!},x>0.$$

(10)

Definition 3

([6]). The Riemann—Liouville fractional derivative. The fractional derivative can be defined using the definition of the fractional integral. To this end, suppose that $v=n-u$, where $0<v<1$ and n is the smallest integral greater than u. Then, the fractional derivative of $f\left(x\right)$ of the order u is

$$D^{u}f\left(x\right)=D^n\left[D^{-v}f\left(x\right)\right].$$

Definition 4

([6]). The linear Lagrange interpolating polynomial through $(x_0,y_0)$ and $(x_1,y_1)$, where

$$L_0\left(x\right)=\frac{x-x_1}{x_0-x_1}and{L}_1\left(x\right)=\frac{x-x_0}{x_1-x_0},$$

is

$$P\left(x\right)=L_0\left(x\right)f\left(x_0\right)+L_1\left(X\right)f\left(x_1\right)=\frac{x-x_1}{x_0-x_1}f\left(x_0\right)+\frac{x-x_0}{x_1-x_0}f\left(x_1\right).$$

Note that

$$L_0\left(x_0\right)=1,L_0\left(x_1\right)=0,L_1\left(x_0\right)=0,L_1\left(x_1\right)=1,$$

which implies that

$$P\left(x_0\right)=1·f\left(x_0\right)+0·f\left(x_1\right)=f\left(x_0\right)=y_0$$

and

$$P\left(x_1\right)=0·f\left(x_0\right)+1·f\left(x_1\right)=f\left(x_1\right)=y_1.$$

Therefore, P is a unique polynomial with a degree of, at most, one that passes through $(x_0,y_0)$ and $(x_1,y_1)$.

Theorem 1

([6]). Consider $x_0$,$x_1$, …, $x_n$ as $n+1$ distinct numbers to define the nth Lagrange interpolating polynomial, where the values of the function f are provided by these numbers, with the one with the highest degree n existing with

$$f\left(x_k\right)=p\left(x_k\right),$$

for each $k=0,1,\dots ,n$. This polynomial is given by the following:

$$p\left(x\right)=f\left(x_0\right)L_{n,0}\left(x\right)+\dots +f\left(x_n\right)L_{n,n}\left(x\right)=\sum _{k=0}^{n}f\left(x_k\right)L_{n,k}\left(x\right),$$

(11)

where, for each $k=0,1,\dots ,n$,

$$Ln,k\left(x\right)=\frac{(x-x_0)(x-x_1)\dots (x-x-k-1)(x-x_{k+1})\dots (x-x_n)}{(x_k-x_0)(x_k-x_1)\dots (x_k-x_{k-1})(x_k-x_{k+1})\dots (x_K-x_n)},$$

$$\prod _{\begin{array}{c}i=0\\ i\ne k\end{array}}^n\frac{(x-x_i)}{(x_k-x_i)}.$$

We write $L_{n,k}\left(x\right)$ simply as $L_k\left(x\right)$ when there is no confusion as to its degree.

Theorem 2

([6]). The interval $[a,b]$ contains the following distinct numbers $x_0,x_1,\cdots ,x_n$ and $f\in C^{n+1}[a,b]$. Then, for each x in $[a,b]$, a number $\xi \left(x\right)$ (generally unknown) between $x_0,x_1,\cdots ,x_n$ and hence in $(a,b)$, exists with

$$f\left(x\right)=P\left(x\right)+\frac{f^{(n+1)}\left(\xi \left(x\right)\right)}{(n+1)!}(x-x_0)(x-x_1)\cdots (x-x_n),$$

where $P\left(x\right)$ is the interpolating polynomial given in (11).

3. Fractional Central Formula

We introduce in this section the fractional central formula to study the function’s behavior according to different numerical or fractional values of $\alpha$.

Theorem 3.

Suppose that $x_0,x_1$, and $x_2$ are distinct points in the interval $[a,b]$ and $f\in C^{2+1}[a,b]$. Then, for each $x\in [a,b]$, a point ξ$\left(x\right)$ (generally unknown) within $x_0,x_1,$ and $x_2$ and, hence, in $(a,b)$ exists with $a=x_0<x_1=x_0+h<x_2=x_0+2h=b$ and $h>0$, so that, by Theorem (2), we have the following:

$$\begin{array}{c}{D}^{\alpha }f\left(x\right)=\frac{x^{2-\alpha }}{h^2\mathsf{\Gamma }(3-\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right))-\frac{x^{1-\alpha }}{2h^2\mathsf{\Gamma }(2-\alpha )}(f\left(x_0\right)(x_1+x_2)-2f\left(x_1\right)(x_0+x_2)\\ +f\left(x_2\right)(x_0+x_1)+{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{6}}\left(\frac{6x^{3-\alpha }}{\mathsf{\Gamma }(4-\alpha }\right)-\frac{2(x_0+x_1+x_2)}{\mathsf{\Gamma }(3-\alpha )}{x}^{2-\alpha }+\frac{(x_0{x}_1+x_0{x}_2+x_1{x}_2)}{\mathsf{\Gamma }(2-\alpha )}{x}^{1-\alpha }).\end{array}$$

(12)

Proof. 

In general, $f\left(x_k\right)=P\left(x_k\right)$ for each $k=0,1,2,\cdots ,n$. This polynomial is given by the following: $P\left(x\right)=f\left(x_0\right)L_{n,0}\left(x\right)+\cdots +f\left(x_n\right)L_{n,n}\left(x\right)={\sum }_{k=0}^{n}f\left(x_k\right)L_{n,k}\left(x\right)$, where for each $k=0,1,2,\cdots ,n$, we have the following:

$$L_{n,k}\left(x\right)=\frac{(x-x_0)(x-x_1)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_n)}{(x_k-x_0)(x_k-x_1)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots (x_k-x_n)}.$$

(13)

Consider $L_{n,k}\left(x\right)$ as $L_k\left(x\right)$ throughout the paper when the degree is clear and unambiguous.

Then, we can write $f\left(x\right)$ in general form as follows:

$$f\left(x\right)=\sum _{k=0}^{m}f\left(x_k\right)L_k\left(x\right)+\frac{(x-x_0)\cdots (x-x_n)}{(n+1)!}f{\left(\xi \left(x\right)\right)}^{n+1}.$$

(14)

For n = 2, we have the following:

$$f\left(x\right)=\sum _{k=0}^{2}f\left(x_k\right)L_k\left(x\right)+\frac{(x-x_0)(x-x_1)(x-x_2)}{3!}{f}^{\left(3\right)}\left(\xi \left(x\right)\right).$$

By simplifying the summation, we have the following:

$$f\left(x\right)=f\left(x_0\right)L_0\left(x\right)+f\left(x_1\right)L_1\left(x\right)+f\left(x_2\right)L_2\left(x\right)+{\displaystyle \frac{1}{6}}{f}^{\left(3\right)}\left(\xi \left(x\right)\right)(x-x_0)(x-x_1)(x-x_2).$$

By using (13), we have the following:

$$\begin{array}{c}f\left(x\right)=f\left(x_0\right)\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+f\left(x_1\right)\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}+f\left(x_2\right)\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}\\ +{\displaystyle \frac{1}{6}}{f}^{\left(3\right)}\left(\xi \left(x\right)\right)(x-x_0)(x-x_1)(x-x_2).\end{array}$$

(15)

Pull out the common factors:

$$\begin{array}{c}f\left(x\right)={\displaystyle \frac{f\left(x_0\right)}{(x_0-x_1)(x_0-x_2)}}(x^2-(x_1+x_2)x+x_1{x}_2)+\frac{f\left(x_1\right)}{(x_1-x_0)(x_1-x_2)}(x^2-(x_0+x_2)x+x_0{x}_2)+\\ {\displaystyle \frac{f\left(x_2\right)}{(x_2-x_0)(x_2-x_1)}}(x^2-(x_0+x_1)x+x_0{x}_1)+\frac{f^{\left(3\right)}\left(\xi \left(x\right)\right)}{6}(x-x_0)(x-x_1)(x-x_2).\end{array}$$

By using $x_1=x_0+h$ and $x_2=x_0+2h$, we can write $f\left(x\right)$ with respect to h:

$$\begin{array}{c}f\left(x\right)=\frac{f\left(x_0\right)}{(-h)(-2h)}(x^2-(x_1+x_2)x+x_1{x}_2)+\frac{f\left(x_1\right)}{\left(h\right)(-h)}(x^2-(x_0+x_2)x+x_0{x}_2)\\ +{\displaystyle \frac{f\left(x_2\right)}{\left(h\right)\left(2h\right)}}(x^2-(x_0+x_1)x+x_0{x}_1)+\frac{f^{\left(3\right)}\left(\xi \left(x\right)\right)}{6}(x^3-(x_0+x_1+x_2)x^2+(x_0{x}_1+x_0{x}_2+x_1{x}_2)x-x_0{x}_1{x}_2).\end{array}$$

The Caputo time fractional derivative can be defined for $f\left(x\right)$ by (2):

$$\begin{array}{c}{D}^{\alpha }f\left(x\right)=\frac{f\left(x_0\right)}{2h^2}\left[\frac{\mathsf{\Gamma }(2+1)x^{2-\alpha }}{\mathsf{\Gamma }(-\alpha +2+1)}-\frac{(x_1+x_2)\mathsf{\Gamma }(1+1)x^{1-\alpha }}{\mathsf{\Gamma }(-\alpha +1+1)}+0\right]\\ +{\displaystyle \frac{f\left(x_1\right)}{h^2}}\left[\frac{\mathsf{\Gamma }(2+1)x^{2-\alpha }}{\mathsf{\Gamma }(-\alpha +2+1)}-\frac{(x_0+x_2)\mathsf{\Gamma }(1+1)x^{1-\alpha }}{\mathsf{\Gamma }(-\alpha +1+1)}+0\right]\\ +{\displaystyle \frac{f\left(x_2\right)}{2h^2}}\left[\frac{\mathsf{\Gamma }(2+1)x^{2-\alpha }}{\mathsf{\Gamma }(-\alpha +2+1)}-\frac{(x_0+x_1)\mathsf{\Gamma }(1+1)x^{1-\alpha }}{\mathsf{\Gamma }(-\alpha +1+1)}+0\right]\\ +{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{\sigma }}\left[\frac{\mathsf{\Gamma }(3+1)x^{3-\alpha }}{\mathsf{\Gamma }(-\alpha +3+1)}-\frac{(x_0+x_1+x_2)\mathsf{\Gamma }(2+1)x^{2-\alpha }}{\mathsf{\Gamma }(-\alpha +2+1)}\frac{(x_0{x}_1+x_2{x}_1+x_0{x}_2)\mathsf{\Gamma }(1+1)x^{1-\alpha }}{\mathsf{\Gamma }(-\alpha +1+1)}+0\right].\end{array}$$

Simplify the formula.

$$\begin{array}{c}{D}^{\alpha }f\left(x\right)=\frac{f\left(x_0\right)}{2h^2}[\frac{2}{\mathsf{\Gamma }(3-\alpha )}{x}^{2-\alpha }-\frac{(x_1+x_2)}{\mathsf{\Gamma }(2-\alpha )}{x}^{1-\alpha }]-{\displaystyle \frac{f\left(x_1\right)}{h^2}}[\frac{2}{\mathsf{\Gamma }(3-\alpha )}{x}^{2-\alpha }-\frac{(x_0+x_2)}{\mathsf{\Gamma }(2-\alpha )}{x}^{1-\alpha }\\ +{\displaystyle \frac{f\left(x_2\right)}{2h^2}}[\frac{2x^{2-\alpha }}{\mathsf{\Gamma }(3-\alpha )}-\frac{(x_0+x_1)}{\mathsf{\Gamma }(2-\alpha )}{x}^{1-\alpha }]+{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{6}}\left[\frac{6}{\mathsf{\Gamma }(4-\alpha )}{x}^{3-\alpha }-\frac{2(x_0+x_1+x_2)}{\mathsf{\Gamma }(3-\alpha )}{x}^{2-\alpha }+\frac{(x_0{x}_1+x_0{x}_2+x_1{x}_2}{\mathsf{\Gamma }(2-\alpha )}{x}^{1-\alpha }\right].\end{array}$$

Consequently, the formula can be expressed as follows:

$$\begin{array}{c}{D}^{\alpha }f\left(x\right)=\frac{f\left(x_0\right)}{\left(2h^2\right)}\frac{\left(2x^{2-\alpha }\right)}{\mathsf{\Gamma }(3-\alpha )}-\frac{\left(f\left(x_1\right)\right)}{h^2}\frac{\left(2x^{2-\alpha }\right)}{\mathsf{\Gamma }(3-\alpha )}+\frac{f\left(x_2\right)}{2h^2}\frac{\left(2x^{2-\alpha }\right)}{\mathsf{\Gamma }(3-\alpha )}-\frac{\left(f\left(x_0\right)\right)}{\left(2h^2\right)}\frac{\left((x_1+x_2)x^{1-\alpha }\right)}{\mathsf{\Gamma }(2-\alpha )}\\ +{\displaystyle \frac{\left(f\left(x_1\right)\right)}{h^2}}\frac{\left((x_0+x_2)x^{1-\alpha }\right)}{\mathsf{\Gamma }(2-\alpha )}-\frac{\left(f\left(x_2\right)\right)}{\left(2h^2\right)}\frac{\left((x_0+x_1)x^{1-\alpha }\right)}{\mathsf{\Gamma }(2-\alpha )}-\frac{f^{\left(3\right)}\left(\xi \right)}{6}[\frac{6x^{3-\alpha }}{\mathsf{\Gamma }(4-\alpha )}-\frac{(x_0+x_1+x_2)\left(2x^{2-\alpha }\right)}{\mathsf{\Gamma }(3-\alpha )}\\ +{\displaystyle \frac{(x_0{x}_1+x_2{x}_1+x_0{x}_2)x^{1-\alpha }}{\mathsf{\Gamma }(2-\alpha )}}.\end{array}$$

Finally, we obtain the final form:

$$\begin{array}{c}{D}^{\alpha }f\left(x\right)=\frac{x^{2-\alpha }}{h^2\mathsf{\Gamma }(3-\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right))\\ -\frac{x^{1-\alpha }}{2h^2\mathsf{\Gamma }(3-\alpha )}(f\left(x_0\right)(x_1+x_2)-2f\left(x_1\right)(x_0+x_2)+f\left(x_2\right)(x_0+x_1))\\ +{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{6}}\left(\frac{6x^{3-\alpha }}{\mathsf{\Gamma }(4-\alpha )}-\frac{2(x_0+x_1+x_2)}{\mathsf{\Gamma }(3-\alpha )}{x}^{2-\alpha }+\frac{(x_0{x}_1+x_0{x}_2+x_1{x}_2)}{\mathsf{\Gamma }(2-\alpha )}\right).\end{array}$$

□

Lemma 1.

From Theorem (3), we can obtain the classical central formula in the general form given by the following:

$$f^{\prime }\left(x_1\right)=\frac{f\left(x_2\right)-f\left(x_0\right)}{2h}-\frac{h^2}{6}{f}^{\left(3\right)\left(\xi \right)}.$$

Proof. 

To obtain the central formula from Theorem (3), let $x=x_1$; then, we have the following:

$$\begin{array}{c}{D}^{\alpha }f\left(x_1\right)=\frac{x_1^{2-\alpha }}{h^2\mathsf{\Gamma }(3-\alpha )}\left(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right)\right)\\ -\frac{x_1^{1-\alpha }}{2h^2\mathsf{\Gamma }(2-\alpha )}\left(f\left(x_0\right)(x_1+x_2)-2f\left(x_1\right)(x_0+x_2)+f\left(x_2\right)(x_0+x_1)\right)\\ +{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{6}}\left(\frac{6x_1^{3-\alpha }}{\mathsf{\Gamma }(4-\alpha )}-\frac{2(x_0+x_1+x_2)}{\mathsf{\Gamma }(3-\alpha )}{x}_1^{2-\alpha }+\frac{(x_0{x}_1+x_0{x}_2+x_1{x}_2)}{\mathsf{\Gamma }(2-\alpha )}{x}_1^{1-\alpha }\right),\end{array}$$

By applying the power rule for the Caputo derivative, we have the following:

$$\begin{array}{c}{f}^{\prime }\left(x\right)={\displaystyle \frac{x_1}{h^2\mathsf{\Gamma }\left(2\right)}}(f\left(x_0\right)+2f\left(x_1\right)+f\left(x_2\right))-\frac{1}{2h^2}(f\left(x_0\right)(x_1+x_2)\\ +2f\left(x_1\right)(x_0+x_2)+f\left(x_2\right)(x_0+x_1))+\frac{f^{\left(3\right)}\left(\xi \right)}{6}[\frac{6x_1^2}{2}-2(x_0+x_1+x_2)+x_0{x}_1+x_2{x}_1+x_0{x}_2].\end{array}$$

Consequently,

$$\begin{array}{c}{f}^{\prime }\left(x\right)={\displaystyle \frac{1}{h^2}}(x_{1}f\left(x_0\right)+2x_{1}f\left(x_1\right)+x_{1}f\left(x_2\right))-\frac{1}{2}{x}_{1}f\left(x_0\right)-\frac{1}{2}{x}_{2}f\left(x_0\right)+x_{0}f\left(x_1\right)+x_{2}f\left(x_1\right)-\\ \frac{1}{2}{x}_{0}f\left(x_2\right)-\frac{1}{2}{x}_{1}f\left(x_2\right)+\frac{f^{\left(3\right)}\left(\xi \right)}{6}[3x_1^2-2x_0{x}_1-2x_1^2-2x_2{x}_1+x_0{x}_1+x_2{x}_1+x_0{x}_2],\end{array}$$

By substituting $x=x_1$, we have the following:

$$\begin{array}{c}{f}^{\prime }\left(x_1\right)=\frac{1}{h^2}\left[\frac{1}{2}{x}_{1}f\left(x_0\right)-2x_{1}f\left(x_1\right)+\frac{1}{2}{x}_{1}f\left(x_2\right)-\frac{1}{2}{x}_{2}f\left(x_0\right)+x_{0}f\left(x_1\right)+x_{2}f\left(x_1\right)-\frac{1}{2}{x}_{0}f\left(x_2\right)\right]\\ +{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{6}}\left[x_1^2-x_0{x}_1-x_1{x}_2+x_0{x}_2\right].\end{array}$$

Consequently,

$$\begin{array}{c}{f}^{\prime }\left(x_1\right)=\frac{1}{h^2}[\frac{1}{2}{x}_{1}f\left(x_0\right)+2f\left(x_1\right)[-2x_1+x_0+x_2]+\frac{1}{2}{x}_{1}f\left(x_2\right)-\frac{1}{2}{x}_{2}f\left(x_0\right)-\frac{1}{2}{x}_{0}f\left(x_2\right)+\\ \frac{f^{\left(3\right)}\left(\xi \right)}{6}[x_1^2-x_0{x}_1-x_1{x}_2+x_0{x}_2]],\end{array}$$

Simplify:

$$\begin{array}{c}{f}^{\prime }\left(x_1\right)=\frac{1}{h^2}\left[\frac{1}{2}f\left(x_0\right)(x_1-x_2)+\frac{1}{2}f\left(x_2\right)(x_1-x_0)\right]+\frac{f^{\left(3\right)}\left(\xi \right)}{6}\left[(x_1-x_0)(x_1-x_2)\right],\end{array}$$

Finally, we have the following:

$$\begin{array}{c}{f}^{\prime }\left(x_1\right)=\frac{1}{h^2}[\frac{-h}{2}f\left(x_0\right)+\frac{h}{2}f\left(x_2\right)]+\frac{f^{\left(3\right)}\left(\xi \right)}{6}\left(h\right)(-h).\end{array}$$

Then, we have the central series:

$$\begin{array}{c}{f}^{\prime }\left(x_1\right)=\frac{f\left(x_2\right)-f\left(x_0\right)}{2h}+\frac{h^2}{6}{f}^{\left(3\right)}\left(\xi \right).\end{array}$$

□

Remark 1.

Note that if we let $x=x_0$ in the Theorem (3), we obtain the classical form of the forward formula in the same way as in Lemma (1).

Remark 2.

Note that if we let $x=x_2$ in the Theorem (3), we obtain the classical form of the backward formula the same way as in Lemma (1).

4. Richardson Extrapolation

To obtain high-accuracy results for low-order formulas, we use Richardson’s extrapolation. For more about the history and applications of Richardson’s extrapolation method, we recommend the article by Joyce, D. C. [23]. In this section, we introduce our methodology for Richardson extrapolation depending on our fractional central form [21]. Depending on the modified central formula [21], for $x\in [a,b]$, assume that $\alpha =1$ and, without losing generality, we have the following:

$$f^{\prime }\left(x\right)=\frac{f(x+h)-f(x-h)}{2h}-\frac{h^2}{6}{f}^{\left(3\right)}\left(\xi \right).$$

(16)

Note that we denote $f^{\prime }\left(x\right)\mathrm{as}Rv$, $\frac{f(x+h)-f(x-h)}{2h}\mathrm{as}Sv\left(h\right)$ and $c=-\frac{f^{\left(3\right)}\left(\xi \right)}{6}$. Then, we obtain the following equation:

$$Rv=\left(Sv\right)\left(h\right)+c{h}^2.$$

(17)

If we substitute $\frac{h}{2}$ instead of h in (17), we have the following:

$$Rv=\left(Sv\right)\left(\frac{h}{2}\right)+c{\left(\frac{h}{2}\right)}^2.$$

By simplifying the previous equation and substituting with $(-4)$, we have the following:

$$-4(Rv=\left(Av\right)\left(\frac{h}{2}\right)+\frac{1}{4}c{h}^2).$$

(18)

By subtracting (18) from (17), we have the following:

$$-3Rv=\left(Sv\right)\left(h\right)-\left(Sv\right)\left(\frac{h}{2}\right)-3\left(Sv\right)\left(\frac{h}{2}\right),$$

then

$$Rv=\left(Sv\right)\left(\frac{h}{2}\right)+\frac{\left(Sv\right)\left(\frac{h}{2}\right)-\left(Av\right)\left(h\right)}{3}.$$

Suppose that $Rv=M_2\left(h\right)$, $M_1\left(h\right)=\left(Sv\right)\left(h\right)$ and $N_1\left(\frac{h}{2}\right)$. We obtain the following from that:

$$M_2\left(h\right)=M_1\left(\frac{h}{2}\right)+\frac{M_1\left(\frac{h}{2}\right)-M_1\left(h\right)}{3}.$$

If we continue in the same way, we obtain the following

Table 1. The 4th order approximation of Richardson’s extrapolation.

Q(h${}^2$)	Q(h${}^4$)	Q(h${}^6$)	Q(h${}^8$)
M${}_1\left(h\right)$	-	-	-
M${}_1\left(\frac{h}{2}\right)$	M${}_2\left(h\right)$	-	-
M${}_1\left(\frac{h}{4}\right)$	M${}_2\left(\frac{h}{2}\right)$	M${}_3\left(h\right)$	-
M${}_1\left(\frac{h}{8}\right)$	M${}_2\left(\frac{h}{4}\right)$	M${}_3\left(\frac{h}{2}\right)$	M${}_4\left(h\right)$

, which illustrates how we can use Richardson’s extrapolation to build a fourth-order approximation, using four first-order approximations:

Consider for $M_1$

$$M_1\left(h\right)=\frac{f(x+h)-f(x-h)}{2h},$$

$$M_1\left(\frac{h}{2}\right)=\frac{f(x+\frac{h}{2})-f(x-\frac{h}{2})}{h},$$

$$M_1\left(\frac{h}{4}\right)=\frac{f(x+\frac{h}{4})-f(x-\frac{h}{4})}{\frac{h}{2}},$$

$$M_1\left(\frac{h}{8}\right)=\frac{f(x+\frac{h}{8})-f(x-\frac{h}{8})}{\frac{h}{4}}.$$

Consider for $M_2$

$$M_2\left(h\right)=M_1\left(\frac{h}{2}\right)+\frac{M_1\left(\frac{h}{2}\right)-M_1\left(h\right)}{3},$$

$$M_2\left(\frac{h}{2}\right)=M_1\left(\frac{h}{4}\right)+\frac{M_1\left(\frac{h}{4}\right)-M_1\left(\frac{h}{2}\right)}{3},$$

$$M_2\left(\frac{h}{4}\right)=M_1\left(\frac{h}{8}\right)+\frac{M_1\left(\frac{h}{8}\right)-M_1\left(\frac{h}{4}\right)}{3}.$$

Consider for $M_3$

$$M_3\left(h\right)=M_2\left(\frac{h}{2}\right)+\frac{M_2\left(\frac{h}{2}\right)-M_2\left(h\right)}{15},$$

$$M_3\left(\frac{h}{2}\right)=M_2\left(\frac{h}{4}\right)+\frac{M_2\left(\frac{h}{4}\right)-M_2\left(\frac{h}{2}\right)}{15}.$$

Consider that for $M_4$:

$$M_4\left(h\right)=M_3\left(\frac{h}{2}\right)+\frac{M_3\left(\frac{h}{2}\right)-M_3\left(h\right)}{63}.$$

Performing this process further results for each $i=2,3,\cdots$, the $Q\left(h^{2i}\right)$ approximation becomes

$$M_i\left(h\right)=M_{i-1}\left(\frac{h}{2}\right)+\frac{M_{i-1}\left(\frac{h}{2}\right)-M_{i-1}\left(h\right)}{4^{i-1}-1}.$$

(19)

Remark 3.

We can generalize the $Q\left(h^{2i}\right)$ approximation for the fractional case when α is fractional and can obtain the same result as (19).

5. Numerical Examples

In this section, we compare the efficiency between the fractional central formula and the Richardson extrapolation; tables and figures are used to present and compare the outcomes in order to demonstrate the validity of the suggested approaches.

Example 1.

Let the main function be $f\left(x\right)=4x^2-2x-15$, and consider $\alpha =1$; then, the Caputo derivative is given by the following:

$$f^{\prime }\left(x\right)=4\frac{\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }\left(2\right)}x-2\frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left(1\right)}.$$

Take the point of evaluation $x_0=2$; $f^{\prime }\left(2\right)=14$. Then, the exact values of the derivatives (the classical central formula) equal the Caputo derivative. By taking the step initial value $h=0.01$ and the order of extrapolation $N=8$, we obtain the results shown in

Table 2. The numerical result when $\alpha =1$ using our central formula with Richardson extrapolation.

h	Q(${\mathit{h}}^2$)	Q(${\mathit{h}}^4$)	Q(${\mathit{h}}^6$)	Q(${\mathit{h}}^8$)	Q(${\mathit{h}}^{10}$)	Q(${\mathit{h}}^{12}$)	Q(${\mathit{h}}^{14}$)	Q(${\mathit{h}}^{16}$)
0.0100	14.0800	0	0	0	0	0	0	0
0.0050	14.0400	14.0267	0	0	0	0	0	0
0.0025	14.0200	14.0133	14.0124	0	0	0	0	0
0.0013	14.0100	14.0067	14.0062	14.0061	0	0	0	0
0.0006	14.0050	14.0033	14.0031	14.0031	14.0030	0	0	0
0.0003	14.0025	14.0017	14.0016	14.0015	14.0015	14.0015	0	0
0.0002	14.0012	14.0008	14.0008	14.0008	14.0008	14.0008	14.0008	0
0.0001	14.00064	14.0004	14.0004	14.0004	14.0004	14.0004	14.0004	14.0004

.

Consider the fractional case when $\alpha =0.75$ and for $x_0=2$:

$$D^{0.75}f\left(2\right)=4\frac{\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }\left(2.25\right)}{x}^{1.25}-2\frac{\mathsf{\Gamma }\left(2\right)}{\mathsf{\Gamma }\left(1.25\right)x^{0.25}}=14.1697,$$

and the Richardson extrapolation in view of the fractional central formula formulated by using Caputo derivative as shown in

Table 3. The numerical result when $\alpha =0.75$ using our central formula with Richardson extrapolation.

h	Q(${\mathit{h}}^2$)	Q(${\mathit{h}}^4$)	Q(${\mathit{h}}^6$)	Q(${\mathit{h}}^8$)	Q(${\mathit{h}}^{10}$)	Q(${\mathit{h}}^{12}$)	Q(${\mathit{h}}^{14}$)	Q(${\mathit{h}}^{16}$)
0.0100	14.2714	0	0	0	0	0	0	0
0.0050	14.2205	14.2036	0	0	0	0	0	0
0.0025	14.1951	14.1866	14.1855	0	0	0	0	0
0.0013	14.1824	14.1782	14.1776	14.1775	0	0	0	0
0.0006	14.1760	14.1739	14.1736	14.1736	14.1736	0	0	0
0.0003	14.1729	14.1718	14.1717	14.1716	14.1716	14.1716	0	0
0.0002	14.1713	14.1707	14.1707	14.1707	14.1707	14.1707	14.1707	0
0.0001	14.1705	14.1702	14.1702	14.1702	14.17024	14.1702	14.1702	14.1702

.

We introduce the absolute error for Richardson extrapolation in

Table 4. The absolute error for Richardson extrapolation between $\alpha =1$ and $\alpha =0.75$.

h	Q$\left({\mathit{h}}^{2\mathit{i}}\right)$	Absolute Error for $\mathit{\alpha }=1$	Absolute Error for $\mathit{\alpha }=0.75$
0.0100	Q($h^2$)	$8\times {10}^{-2}$	$1.017\times {10}^{-1}$
0.0050	Q($h^4$)	$2.67\times {10}^{-2}$	$3.39\times {10}^{-2}$
0.0025	Q($h^6$)	$1.24\times {10}^{-2}$	$1.58\times {10}^{-2}$
0.0013	Q($h^8$)	$6.1\times {10}^{-3}$	$7.8\times {10}^{-3}$
0.0006	Q($h^{10}$)	$3\times {10}^{-3}$	$3.9\times {10}^{-3}$
0.0003	Q($h^{12}$)	$1.5\times {10}^{-3}$	$1.9\times {10}^{-3}$
0.0002	Q($h^{14}$)	$8\times {10}^{-4}$	$1\times {10}^{-3}$
0.0001	Q($h^{16}$)	$4\times {10}^{-4}$	$5\times {10}^{-4}$

by tacking different values of α to see the efficiency of the proposed methods.

In what follows, we introduce (Figure 1) some graphical comparisons between the results generated using the analytical and fractional central formulas.

Example 2.

Let the main function $f\left(x\right)=\frac{1}{4}\sqrt{x^3}-8\sqrt{x}$, and consider $\alpha =1$; then, the Caputo derivative is given by the following:

$$f^{\prime }\left(x\right)=\frac{1}{4}\frac{\mathsf{\Gamma }\left(\frac{5}{2}\right)}{\mathsf{\Gamma }(\frac{5}{2}-1)}\sqrt{x}-8\frac{\mathsf{\Gamma }\left(\frac{3}{2}\right)}{\mathsf{\Gamma }(\frac{3}{2}-1)}\frac{1}{\sqrt{x}}.$$

Take the point of evaluation $x_0=\frac{1}{2}$: $f^{\prime }\left(\frac{1}{2}\right)=-5.3917$; then, the exact values of the derivatives (the classical central formula) equal the Caputo derivative. By taking the step initial value $h=0.01$ and the order of extrapolation $N=8$, we get the results shown in

Table 5. The numerical result when $\alpha =1$ using modified central formula with Richardson extrapolation.

h	Q(${\mathit{h}}^2$)	Q(${\mathit{h}}^4$)	Q(${\mathit{h}}^6$)	Q(${\mathit{h}}^8$)	Q(${\mathit{h}}^{10}$)	Q(${\mathit{h}}^{12}$)	Q(${\mathit{h}}^{14}$)	Q(${\mathit{h}}^{16}$)
0.0100	−5.3336	0	0	0	0	0	0	0
0.0050	−5.3624	−5.3720	0	0	0	0	0	0
0.0025	−5.3770	−5.3818	−5.3825	0	0	0	0	0
0.0013	−5.3843	−5.3868	−5.3871	−5.3872	0	0	0	0
0.0006	−5.3880	−5.3892	−5.3894	−5.3894	−5.3894	0	0	0
0.0003	−5.3898	−5.3905	−5.3905	−3.3906	−5.3906	−5.3906	0	0
0.0002	−5.3908	−5.3911	−5.3911	−5.3911	−5.3911	−5.3911	−5.3911	0
0.0001	−5.3912	−5.3914	−5.3914	−5.3914	−5.3914	−5.3914	−5.3914	−5.3914

.

Consider the fractional case when $\alpha =0.98$ and for $x_0=\frac{1}{2}$:

$$D^{0.98}f\left(\frac{1}{2}\right)=\frac{1}{4}\frac{\mathsf{\Gamma }\left(\frac{5}{2}\right)}{\mathsf{\Gamma }(\frac{5}{2}-0.98)}{x}^{\frac{3}{2}-0.98}-8\frac{\mathsf{\Gamma }\left(\frac{3}{2}\right)}{\mathsf{\Gamma }(\frac{3}{2}-0.98)}{x}^{\frac{1}{2}-0.98}=-5.5355,$$

and the Richardson extrapolation in view of the fractional central formula formulated by using Caputo derivative as shown in

Table 6. The numerical result when $\alpha =0.98$ using our central formula with Richardson extrapolation.

h	Q(${\mathit{h}}^2$)	Q(${\mathit{h}}^4$)	Q(${\mathit{h}}^6$)	Q(${\mathit{h}}^8$)	Q(${\mathit{h}}^{10}$)	Q(${\mathit{h}}^{12}$)	Q(${\mathit{h}}^{14}$)	Q(${\mathit{h}}^{16}$)
0.0100	−5.3790	0	0	0	0	0	0	0
0.0050	−5.4069	−5.4162	0	0	0	0	0	0
0.0025	−5.4211	−5.4258	−5.4264	0	0	0	0	0
0.0013	−5.4282	−5.4306	−5.4309	−5.4310	0	0	0	0
0.0006	−5.4318	−5.4330	−5.4331	−5.4332	−5.4332	0	0	0
0.0003	−5.4336	−5.4342	−5.4342	−5.4343	−5.4343	−5.4343	0	0
0.0002	−5.4345	−5.4348	−5.4348	−5.4348	−5.4348	−5.4348	−5.4348	0
0.0001	−5.4349	−5.4351	−5.4351	−5.4351	−5.4351	−5.4351	−5.4351	−5.4351

.

We introduce the absolute error for Richardson extrapolation in

Table 7. The absolute error for Richardson extrapolation between $\alpha =1$ and $\alpha =0.98$.

h	Q$\left({\mathit{h}}^{2\mathit{i}}\right)$	Absolute Error for $\mathit{\alpha }=1$	Absolute Error for $\mathit{\alpha }=0.98$
0.0100	$Q\left(h^2\right)$	$5.81\times {10}^{-2}$	$1.565\times {10}^{-1}$
0.0050	$Q\left(h^4\right)$	$1.97\times {10}^{-2}$	$1.193\times {10}^{-1}$
0.0025	$Q\left(h^6\right)$	$9.2\times {10}^{-3}$	$1.091\times {10}^{-1}$
0.0013	$Q\left(h^8\right)$	$4.5\times {10}^{-3}$	$1.045\times {10}^{-1}$
0.0006	$Q\left(h^{10}\right)$	$2.3\times {10}^{-3}$	$1.023\times {10}^{-1}$
0.0003	$Q\left(h^{12}\right)$	$1.1\times {10}^{-3}$	$1.012\times {10}^{-1}$
0.0002	$Q\left(h^{14}\right)$	$6\times {10}^{-4}$	$1.007\times {10}^{-1}$
0.0001	$Q\left(h^{16}\right)$	$3\times {10}^{-4}$	$1.004\times {10}^{-1}$

by tacking different values of α to see the efficiency of the proposed methods.

In what follows, we introduce (Figure 2) some graphical comparisons between the results generated using the analytical and fractional central formulas.

Example 3.

Let the main function be $f\left(x\right)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\frac{x^8}{8!}-\frac{x^{1}0}{10!}$, and consider $\alpha =1$; then, the Caputo derivative of the function is given by the following:

$$f^{\prime }\left(x\right)=-\frac{\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }\left(2\right)}\frac{1}{2!}x+\frac{\mathsf{\Gamma }\left(5\right)}{\mathsf{\Gamma }\left(4\right)}\frac{1}{4!}{x}^3-\frac{\mathsf{\Gamma }\left(7\right)}{\mathsf{\Gamma }\left(6\right)}\frac{1}{6!}{x}^5+\frac{\mathsf{\Gamma }\left(9\right)}{\mathsf{\Gamma }\left(8\right)}\frac{1}{8!}{x}^7-\frac{\mathsf{\Gamma }\left(11\right)}{\mathsf{\Gamma }\left(10\right)}\frac{1}{10!}{x}^9.$$

Take the point of evaluation $x_0=\frac{\pi }{6}$: $f^{\prime }\left(\frac{\pi }{6}\right)=-0.5000$; then, the exact values of the derivatives (the classical central formula) equal the Caputo derivative. By taking the step initial value $h=0.01$ and the order of extrapolation $N=8$, we obtain the results shown in

Table 8. The numerical result for $\alpha =1$ using our central formula with Richardson extrapolation.

h	Q(${\mathit{h}}^2$)	Q(${\mathit{h}}^4$)	Q(${\mathit{h}}^6$)	Q(${\mathit{h}}^8$)	Q(${\mathit{h}}^{10}$)	Q(${\mathit{h}}^{12}$)	Q(${\mathit{h}}^{14}$)	Q(${\mathit{h}}^{16}$)
0.0100	−0.5086	0	0	0	0	0	0	0
0.0050	−0.5043	−0.5029	0	0	0	0	0	0
0.0025	−0.5022	−0.5014	−0.5013	0	0	0	0	0
0.0013	−0.5011	−0.5007	−0.5007	−0.5007	0	0	0	0
0.0006	−0.5005	−0.5004	−0.5003	−0.5003	−0.5003	0	0	0
0.0003	−0.5003	−0.5002	−0.5002	−0.5002	−0.5002	−−0.5002	0	0
0.0002	−0.5001	−0.5001	−0.5001	−0.5001	−0.5001	−0.5001	−0.5001	0
0.0001	−0.5001	−0.5000	−0.5000	−0.5000	−0.5000	−0.5000	−0.5000	−0.5000

.

Consider the fractional case when $\alpha =0.35$ and for $x_0=\frac{\pi }{6}$:

$$D^{0.35}f\left(\frac{\pi }{6}\right)=-\frac{\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }(3-0.35)}\frac{1}{2!}{x}^{2-0.35}+\frac{\mathsf{\Gamma }\left(5\right)}{\mathsf{\Gamma }(5-0.35)}\frac{1}{4!}{x}^{4-0.35}-\frac{\mathsf{\Gamma }\left(7\right)}{\mathsf{\Gamma }(7-0.35)}\frac{1}{6!}{x}^{6-0.35}$$

$$+{\displaystyle \frac{\mathsf{\Gamma }\left(9\right)}{\mathsf{\Gamma }(9-0.35)}}\frac{1}{8!}{x}^{8-0.35}-\frac{\mathsf{\Gamma }\left(11\right)}{\mathsf{\Gamma }(11-0.35)}\frac{1}{10!}{x}^{10-0.35}=-0.2250,$$

and the Richardson extrapolation in view of the fractional central formula formulated by using Caputo derivative as shown in

Table 9. The numerical result when $\alpha =0.35$ using our central formula with Richardson extrapolation.

h	Q(${\mathit{h}}^2$)	Q(${\mathit{h}}^4$)	Q(${\mathit{h}}^6$)	Q(${\mathit{h}}^8$)	Q(${\mathit{h}}^{10}$)	Q(${\mathit{h}}^{12}$)	Q(${\mathit{h}}^{14}$)	Q(${\mathit{h}}^{16}$)
0.0100	−0.2420	0	0	0	0	0	0	0
0.0050	−0.2382	−0.2369	0	0	0	0	0	0
0.0025	−0.2363	−0.2357	−0.2356	0	0	0	0	0
0.0013	−0.2354	−0.2351	−0.2350	−0.2350	0	0	0	0
0.0006	−0.2349	−0.2348	−0.2347	−0.2347	−0.2347	0	0	0
0.0003	−0.2347	−0.2346	−0.2346	−0.2346	−0.2346	−0.2346	0	0
0.0002	−0.2346	−0.2345	−0.2345	−0.2345	−0.2345	−0.2345	−0.2345	0
0.0001	−0.2345	−0.2345	−0.2345	−0.2345	−0.2345	−0.2345	−0.2345	−0.2345

.

We introduce the absolute error for Richardson extrapolation in

Table 10. The absolute error for Richardson extrapolation between $\alpha =1$ and $\alpha =0.35$.

h	Q$\left({\mathit{h}}^{2\mathit{i}}\right)$	Absolute Error for $\mathit{\alpha }=1$	Absolute Error for $\mathit{\alpha }=0.35$
0.0100	Q(h${}^2$)	$8.6\times {10}^{-2}$	$1.7\times {10}^{-2}$
0.0050	Q(h${}^4$)	$2.9\times {10}^{-2}$	$1.81\times {10}^{-2}$
0.0025	Q(h${}^6$)	$1.3\times {10}^{-2}$	$1.94\times {10}^{-2}$
0.0013	Q(h${}^8$)	$7\times {10}^{-3}$	$2\times {10}^{-2}$
0.0006	Q(h${}^{10}$)	$3\times {10}^{-3}$	$2.03\times {10}^{-2}$
0.0003	Q(h${}^{12}$)	$2\times {10}^{-3}$	$2.04\times {10}^{-2}$
0.0002	Q(h${}^{14}$)	$1\times {10}^{-3}$	$2.05\times {10}^{-2}$
0.0001	Q(h${}^{16}$)	0	$2.05\times {10}^{-2}$

by tacking different values of α to see the efficiency of the proposed methods.

In what follows, we introduce (Figure 3) some graphical comparisons between the results generated using the analytical and fractional central formulas.

6. Conclusions

In this work, we generalized the classical central formula to introduce the modified three-point fractional formula; the function’s behavior for our classical central formula is the same as the main function, and we show this in Section 5 from the figures and tables. Additionally, we made a new methodology for the Richardson extrapolation. In our next work, we will generate our proposed scheme to establish a new generalization of the second-order central formula in the fractional central case.