First Derivative Approximations and Applications

1. Introduction

The fractional derivatives are a generalization of the integer-order derivatives and provide powerful tools for analyzing functions and systems. Many of the methods for studying integer-order derivatives also apply to fractional derivatives. The Caputo fractional derivative of order $\alpha$, where $0<\alpha <1$ is defined as

$$f^{\left(\alpha \right)}\left(x\right)=D^{\alpha }f\left(x\right)={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_l^x{\displaystyle \frac{f^{\prime }\left(\xi \right)}{{(x-\xi )}^{\alpha }}}d\xi .$$

Fractional calculus is an active research area and have been used to describe phenomena that cannot be adequately explained by integer-order derivatives. Fractional differential equations have found applications in various fields, including physics, engineering, chemistry and biology [1,2,3,4,5,6]. Numerical methods are used to analyze models that involve fractional differential equations. The finite difference schemes for solving numerically fractional differential equations use approximations of the fractional derivative. Consider the interval $[l,x]$ and a uniform grid with a step size of $h=(x-l)/N$, where N is a positive integer. Denote $x_n=l+nh$ and $f_n=f\left(x_n\right)$. The approximations of fractional derivative $f_n^{\left(\alpha \right)}$ have the form $h^{-\alpha }{\sum }_{k=0}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}$ where ${\sigma }_k^{\left(\alpha \right)}$ are the weights of the approximation and $n=1,\cdots ,N$. The generating function of an approximation of the fractional derivative is defined as $G\left(x\right)={\sum }_{k=0}^{\infty }{\sigma }_k^{\left(\alpha \right)}{x}^k$. Two important approximations of the fractional derivative are the Grünwald difference approximation and the L1 approximation. The Grünwald difference approximation of the fractional derivative of order $\alpha$ has a first-order accuracy and a generating function ${(1-x)}^{\alpha }$. Grünwald approximation is a second-order shifted approximation of he fractional derivative with a shift parameter $\alpha /2$. L1 approximation has an order $2-\alpha$ and a generating function ${(x-1)}^2{\mathrm{Li}}_{\alpha -1}\left(x\right)/\left(x\mathrm{\Gamma }(2-\alpha )\right)$. Their weights satisfy

$${\sigma }_0^{\left(\alpha \right)}>0,{\sigma }_1^{\left(\alpha \right)}<{\sigma }_2^{\left(\alpha \right)}<\cdots <{\sigma }_{m-1}^{\left(\alpha \right)}<0,\sum _{k=0}^m{\sigma }_k^{\left(\alpha \right)}=0.$$

(1)

The properties of the weights of an approximation of the fractional derivative (1) allow for an effective analysis of the convergence of numerical schemes for fractional differential equations [7,8,9,10]. The construction of approximations for the fractional derivative and schemes for numerically solving fractional differential equations is an area of ongoing research. Second-order approximations of the fractional derivative are constructed by Arshad et al. [11], Nasir and Nafa [12]. Approximations of order $3-\alpha$ are constructed by Alikhanov [13], Gao et al. [14], Xing and Yan [15]. High-order approximations and numerical schemes for fractional differential equations are studied in [16,17,18,19,20]. Implicit ADI schemes for fractional differential equations are constructed by Nasrollahzadeh and Hosseini [21], Wang et al. [22]. Denote $s_n={\sum }_{k=1}^{n-1}{k}^{-\alpha }-n^{1-\alpha }/(1-\alpha )-\zeta \left(\alpha \right)$. In [23] we obtain an approximation of the fractional derivative and its asymptotic formula

$${\displaystyle \frac{1}{2\mathrm{\Gamma }(1-\alpha )h^{\alpha }}}\sum _{k=0}^m{\omega }_k^{\left(\alpha \right)}{f}_{n-k}=f_m^{\left(\alpha \right)}+{\displaystyle \frac{\zeta \left(\alpha \right)}{\mathrm{\Gamma }(1-\alpha )}}{f}_n^{\prime }{h}^{1-\alpha }+\mathcal{O}\left(h^{2-\alpha }\right),$$

(2)

where

$${\omega }_0^{\left(\alpha \right)}=1,{\omega }_1^{\left(\alpha \right)}=2^{-\alpha },{\omega }_k^{\left(\alpha \right)}={(k+1)}^{-\alpha }-{(k-1)}^{-\alpha },$$

$${\sigma }_{m-1}^{\left(\alpha \right)}=-{(n-1)}^{-\alpha }+2s_n,{\sigma }_n^{\left(\alpha \right)}=-{(n-2)}^{-\alpha }-2s_n.$$

Approximation (2) has an order $1-\alpha$ and a generating function $(1-x^2){\mathrm{Li}}_{-\alpha }\left(x\right)/\left(2x\right)$. The main class of approximations of integer-order derivatives are the finite-difference approximations which have first-order accuracy and second-order at the midpoint. Finite-difference approximations are local numerical differentiation methods and are a special case of the Grünwald difference approximation, where the order is an integer. Methods for global numerical differentiation are studied in [24,25]. The methods for constructing approximations of fractional derivatives using a generating function also apply to the construction of integer-order derivative approximations. In [26,27] we construct approximations of the first and second derivatives whose generating functions are transformations of the exponential and logarithmic functions. In [26,28] we construct second-order approximations of the fractional derivative using the asymptotic formula of the L1 approximation and second derivative approximations. In this paper we apply the method from [26,27] for constructing an approximation of first derivative and its second-order asymptotic formula

$${\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{f}_{n-k}-{\displaystyle \frac{a^{n-2}(1-2a)}{1-a}}{f}_1-{\displaystyle \frac{a^{n-1}}{1-a}}{f}_0\right)=f_n^{\prime }-{\displaystyle \frac{C_a}{2}}{f}_n^{\prime \prime }h+\mathcal{O}\left(h^2\right),$$

(3)

where the coefficient $C_a$ has a value $C_a=(1+a)/(1-a)$. Approximation (3) has a generating function $G\left(x\right)=(1-a)(1-x)/(1-ax)$ and is a global approximation of the first derivative. The structure of the paper is as follows. In section 2 we construct approximation (3) and the following first-order approximation

$${\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{f}_{n-k}-a^{n-1}{f}_0\right)=(1-a^n)f_n^{\prime }+\mathcal{O}\left(h\right).$$

(4)

The weights of approximations (3) and (4) satisfy conditions (1). We consider the application of these approximations for numerically solving an ordinary differential equation and propose an algorithm that computes the numerical solutions using $\mathcal{O}\left(N\right)$ arithmetic operations. The computational time and accuracy of the numerical methods are compared with Euler’s method. In section 3 we derive estimates for the errors of approximations (3) and (4) and the corresponding numerical solutions. In Section 4, we construct numerical solutions for first-order ODEs which use approximation (3) of the first derivative. We determine the values of the parameter a such that the numerical methods achieve an arbitrary order of accuracy in the interval $(0,2]$. In section 5 we construct a finite difference scheme for the heat diffusion equation which uses approximation (3) for the first partial derivative. The stability and convergence of the scheme are analyzed. In section 6, we consider an application of approximation (3) for constructing an approximation of the Caputo fractional derivative. By applying approximation (3) with the parameter value $\alpha /2$ to the first derivative $f_n^{\prime }$ in the asymptotic formula (2), we obtain the following approximation of a fractional derivative of order $2-\alpha$.

$${\displaystyle \frac{1}{2\mathrm{\Gamma }(1-\alpha )h^{\alpha }}}\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+\mathcal{O}\left(h^{2-\alpha }\right),$$

(5)

where

$${\sigma }_0^{\left(\alpha \right)}=1-(2-\alpha )\zeta \left(\alpha \right),{\sigma }_1^{\left(\alpha \right)}=2^{-\alpha }+2{(1-\alpha /2)}^2\zeta \left(\alpha \right),$$

$${\sigma }_k^{\left(\alpha \right)}={(k+1)}^{-\alpha }-{(k-1)}^{-\alpha }+2{(1-\alpha /2)}^2{(\alpha /2)}^{k-1}\zeta \left(\alpha \right),(2\le k\le n-2),$$

$${\sigma }_{n-1}^{\left(\alpha \right)}=-{(n-1)}^{-\alpha }+2s_n+2(1-\alpha ){(\alpha /2)}^{n-2}\zeta \left(\alpha \right),$$

$${\sigma }_n^{\left(\alpha \right)}=-{(n-2)}^{-\alpha }-2s_n+2{(\alpha /2)}^{n-1}\zeta \left(\alpha \right).$$

We prove that the weights of approximation (5) satisfy properties (1). The numerical experiments presented in the paper validate the theoretical results on the accuracy and error of the numerical methods.

2. Asymptotic Formula

In this section we use the method from [26,27] to construct an approximation for the first derivative and its asymptotic formula by specifying the generating function of the approximation. Let

$$G\left(x\right)={\displaystyle \frac{(1-a)(1-x)}{1-ax}},$$

where $\left|a\right|<1$.The generating function $G\left(x\right)$ has properties $G\left(1\right)=0,G^{\prime }\left(1\right)=-1$. Denote

$$H\left(x\right)=G\left(e^{-x}\right)={\displaystyle \frac{(1-a)(1-e^{-x})}{1-a{e}^{-x}}}.$$

The functions $G\left(x\right)$ and $H\left(x\right)$ have Maclaurin series

$$G\left(x\right)=(1-a)-{(1-a)}^2\sum _{k=1}^{\infty }{a}^{k-1}{x}^k,$$

(6)

$$H\left(x\right)=x-{\displaystyle \frac{(1+a)x^2}{2(1-a)}}+{\displaystyle \frac{(1+4a+a^2)x^3}{6{(1-a)}^2}}+\mathcal{O}\left(x^4\right).$$

(7)

Consider a uniform grid on the interval $[0,X]$ with step size of $h=X/N$ where N is a positive integer and let $f_n=f\left(x_n\right)=f\left(nh\right)$. From (6) and (7) we obtain

$${\mathcal{A}}_1\left[f_n\right]={\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{f}_{n-k}\right)=f_n^{\prime }-{\displaystyle \frac{1+a}{2(1-a)}}{f}_n^{\prime \prime }h+\mathcal{O}\left(h^2\right),$$

(8)

where the function $f\in C^2[0,X]$ and satisfies the condition $f\left(0\right)=0$. Now we extend asymptotic formula (8) to all functions in the class $C^2[0,X]$. Let

$${\mathcal{A}}_2\left[f_n\right]={\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{f}_{n-k}+c_0{f}_0\right)=f_n^{\prime }+\mathcal{O}\left(h\right),$$

where ${\mathcal{A}}_2$ satisfies the condition ${\mathcal{A}}_2\left[1\right]=0$.

$$1-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}+c_0=0,1-(1-a^{n-1})+c_0=0.$$

The coefficient $c_0$ has a value $c_0=-a^{n-1}$ and

$${\mathcal{A}}_2\left[f_n\right]={\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{f}_{n-k}-a^{n-1}{f}_0\right)=f_n^{\prime }+\mathcal{O}\left(h\right).$$

(9)

Asymptotic formula (9) holds for all functions $f\in C^2[0,X]$. When the parameter a is zero, ${\mathcal{A}}_2$ is a first-order backward difference approximation of the first derivative. We will consider an application of formula (9) for the numerical solution of an ordinary differential equation and compare the performance of the method with the Euler method.

Example 1. 

Consider the following first order ordinary differential equation

$$y^{\prime }\left(x\right)=F\left(x\right),y\left(0\right)=y_0,x\in [0,X].$$

(10)

By approximating the first derivative at the point $x_n$ using a first-order backward difference approximation, we get

$${\displaystyle \frac{y_n-y_{n-1}}{h}}=F_n+\mathcal{O}\left(h\right).$$

The numerical solution of equation (10) satisfies $u_0=y_0$ and

$$u_n=u_{n-1}+h{F}_n.$$

(NS1)

The value of $u_n$ is computed with one addition and one multiplication, assuming that the values $F_n$ of the function F at the nodes $x_n$ of the grid are known. The total number of arithmetic operations of Euler method is $2N$. Note that in many cases, the computational time required to evaluate the values of $F_n$ exceeds the time needed to compute the numerical solution NS1. Now we construct the numerical solution of equation (10) which uses asymptotic formula (9). By substituting $y_n^{\prime }$ with ${\mathcal{A}}_2$ we find

$${\displaystyle \frac{1-a}{h}}\left(y_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{y}_{n-k}-a^{n-1}{y}_0\right)=F_n+\mathcal{O}\left(h\right).$$

(11)

The numerical solution obtained from (11) by truncating the error term satisfies

$${\displaystyle \frac{1-a}{h}}\left(u_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{u}_{n-k}-a^{n-1}{u}_0\right)=F_n.$$

The sequence ${\left\{u_n\right\}}_{n=0}^N$ has initial conditions $u_0=y_0,u_1=u_0+h{F}_1$ and satisfies

$$u_n={\displaystyle \frac{h{F}_n}{1-a}}+(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{u}_{n-k}+a^{n-1}{u}_0.$$

(NS2)

Numerical solution NS1 is a special case of NS2 when the parameter a is zero, $a=0$. Direct computation of numerical solution NS2 using the above formula requires $\mathcal{O}\left(N^2\right)$ operations. The weights of ${\mathcal{A}}_2$ consist of powers of the parameter a, which allows to compute NS2 with $\mathcal{O}\left(N\right)$ operations. The sequence $u_n$ satisfies

$$u_n={\displaystyle \frac{h{F}_n}{1-a}}+S_n,$$

where

$$S_n=(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{u}_{n-k}+a^{n-1}{u}_0.$$

The sequence $S_n$ satisfies

$$S_n=(1-a)u_{n-1}+a\left((1-a)\sum _{k=2}^{n-1}{a}^{k-2}{u}_{n-k}+a^{n-2}{u}_0\right),$$

$$S_n=(1-a)u_{n-1}+a{S}_{n-1}.$$

Both sequences $S_n$ and $u_n$ are computed with $\mathcal{O}\left(N\right)$ operations. The pseudocode for calculating the numerical solution NS2 is given below.

$$\begin{array}{cc}\hfill \mathit{Step}\mathit{1}\text{:}& A=1-a;H=h/A;S=u_0=y_0;u_1=u_0+h*F_1;\hfill \\ \hfill \mathit{Step}\mathit{2}\text{:}& \mathrm{for}n\mathrm{from}2\mathrm{to}N\hfill \\ \hfill & S=A*u_{n-1}+a*S;\hfill \\ \hfill & u_n=H*F_n+S;\hfill \end{array}$$

The algorithm uses four operations on Step 1 and five operations for calculating each value of the sequence $u_n$ for $n=2,\cdots ,N$ in Step 2. Total number of operations for computation of numerical solution NS2 is $5N-1$. Consider the ordinary differential equation

$$y^{\prime }\left(x\right)=e^x,y\left(0\right)=1,$$

(12)

where $x\in [0,1]$. The numerical results for the error and order of numerical solutions 11 and NS1 of equation (NS2) are given in Table 1. The calculation of the numerical solutions was performed with Mathematica software. Although formula (9) is valid only asymptotically, the experimental results in Table 1 suggest that numerical method NS1 has a first-order accuracy. This fact is proven in the next section. The time for computation of the values of $F_n=e^{nh}$ is greater than the computational time of the numerical solutions. We observed a small difference of the computational time of the numerical solutions in Table 1, which is a fraction of a second. The finite geometric series satisfies

$$\sum _{k=0}^{n-1}{a}^k={\displaystyle \frac{1-a^n}{1-a}}.$$

(13)

By differentiating (13) twice we find

$$\sum _{k=1}^{n-1}k{a}^{k-1}={\displaystyle \frac{1+(n-1)a^n-n{a}^{n-1}}{{(1-a)}^2}},$$

(14)

$$\sum _{k=1}^{n-1}{k}^2{a}^{k-1}={\displaystyle \frac{1+a-n^2{a}^{n-1}+(2n^2-2n-1)a^n-{(n-1)}^2{a}^{n+1}}{{(1-a)}^3}}.$$

(15)

Now we construct an approximation of the first derivative in the form

$${\mathcal{A}}_3\left[f_n\right]={\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{f}_{n-k}-a^{n-2}{c}_1{f}_1-a^{n-1}{c}_0{f}_0\right)=f_n^{\prime }+\mathcal{O}\left(h\right)$$

which satisfies the conditions ${\mathcal{A}}_3\left[1\right]=0$ and ${\mathcal{A}}_3\left[x\right]=1$. The values of the coefficients $c_0$ and $c_1$ are determined from the two conditions.

$${\mathcal{A}}_3\left[1\right]={\displaystyle \frac{1-a}{h}}\left(1-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}-a^{n-2}{c}_1-a^{n-1}{c}_0\right)=0,$$

$$1-(1-a^{n-2})-a^{n-2}{c}_1-a^{n-1}{c}_0=0.$$

The numbers $c_0$ and $c_1$ satisfy

$$c_1+a{c}_0=1.$$

(16)

$${\mathcal{A}}_3\left[x\right]=(1-a)\left(n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}(n-k)-a^{n-2}{c}_1\right)=1.$$

Using formulas (13) and (14) we get

$$\sum _{k=1}^{n-2}(n-k)a^{k-1}=-{\displaystyle \frac{1+a^{n-2}-2a^{n-1}-n+na}{{(1-a)}^2}}.$$

Then

$$(1-a)n+1+a^{n-2}-2a^{n-1}-n+na-(1-a)a^{n-2}{c}_1=1,$$

$$a^{n-2}-2a^{n-1}-(1-a)a^{n-2}{c}_1=1,(1-a)c_1=1-2a.$$

(17)

From (16) and (17) we find

$$c_1={\displaystyle \frac{1-2a}{1-a}},c_0={\displaystyle \frac{1}{1-a}}.$$

Approximation ${\mathcal{A}}_3$ has the form

$${\mathcal{A}}_3\left[f_n\right]={\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{f}_{n-k}-{\displaystyle \frac{a^{n-2}(1-2a)}{1-a}}{f}_1-{\displaystyle \frac{a^{n-1}}{1-a}}{f}_0\right)=f_n^{\prime }+\mathcal{O}\left(h\right)$$

and holds for every value of n, where $n=2,\cdots ,N$. Now we construct the numerical solution of equation (10) which uses approximation ${\mathcal{A}}_3$ for the first derivative. By replacing the first derivative $y_n^{\prime }$ with ${\mathcal{A}}_3$ we get

$${\displaystyle \frac{1-a}{h}}\left(u_n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{u}_{n-k}-{\displaystyle \frac{a^{n-2}(1-2a)}{1-a}}{u}_1-{\displaystyle \frac{a^{n-1}}{1-a}}{u}_0\right)=F_n,$$

$$u_n={\displaystyle \frac{h{F}_n}{1-a}}+(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{u}_{n-k}+{\displaystyle \frac{a^{n-2}(1-2a)}{1-a}}{u}_1+{\displaystyle \frac{a^{n-1}}{1-a}}{u}_0.$$

(NS3)

Numerical solution NS3 has initial conditions $u_0=y_0,u_1=u_0+h{F}_1,u_2=u_1+h{F}_2$. The algorithm and pseudocode for computing the numerical solution NS3 are given below:

$$u_n={\displaystyle \frac{h{F}_n}{1-a}}+S_n,$$

where

$$S_n=(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{u}_{n-k}+a^{n-1}{u}_0.$$

The sequence $S_n$ satisfies

$$S_n=(1-a)u_{n-1}+a\left((1-a)\sum _{k=2}^{n-1}{a}^{k-2}{u}_{n-k}+a^{n-2}{u}_0\right),$$

$$S_n=(1-a)u_{n-1}+a{S}_{n-1}.$$

The sequences $S_n$ and $u_n$ are calculated in $\mathcal{O}\left(N\right)$ time with the following pseudocode.

$$\begin{array}{cc}\hfill \mathrm{Step}1\text{:}& u_0=y_0;u_1=u_0+h*F_1;u_2=u_1+h*F_2;\hfill \\ \hfill & A=1-a,H=h/A;S=u_1+a*(u_0-u_1)/A;\hfill \\ \hfill \mathrm{Step}2\text{:}& \mathrm{for}n\mathrm{from}2\mathrm{to}N\hfill \\ \hfill & S=A*u_{n-1}+a*S;\hfill \\ \hfill & u_n=H*F_n+S;\hfill \end{array}$$

The algorithm uses ten operations on Step 1 and five operations for calculating each value of the sequence $u_n$ in Step 2 for $n=3,\cdots ,N$. The total number of operations for calculating numerical solution NS3 is $5N$.

From (8) approximation ${\mathcal{A}}_3$ has a second-order asymptotic expansion formula

$${\displaystyle \frac{1-a}{h}}\left(y_n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{y}_{n-k}-{\displaystyle \frac{a^{n-2}(1-2a)}{1-a}}{y}_1-{\displaystyle \frac{a^{n-1}}{1-a}}{y}_0\right)=f_n^{\prime }-{\displaystyle \frac{C_a}{2}}{f}_n^{\prime \prime }h+\mathcal{O}\left(h^2\right),$$

(18)

where $C_a=(1+a)/(1-a)$. Approximation ${\mathcal{A}}_3$ has second-order accuracy, when the coefficient $C_a$ is equal to zero, $C_a=0$ or $C_a=h$.

$${\displaystyle \frac{1+a}{1-a}}=h,1+a=h-ah,a={\displaystyle \frac{h-1}{h+1}}.$$

Approximation ${\mathcal{A}}_3$ is second-order accurate when $a=-1$ and $a=(h-1)/(h+1)$. Consider the following ordinary differential equation

$$y^{\prime }\left(x\right)=sinx+cosx,y\left(0\right)=-1,x\in [0,1].$$

(19)

Equation (19) has a solution $y\left(x\right)=sinx-cosx$. The experimental results of numerical solution NS3 for equation (19) and values of the parameter $a=0.5,a=-1$ and $a=(h-1)/(h+1)$ are given in Table 2.

3. Error Estimates and Convergence

In this section we derive estimates for the errors of approximations ${\mathcal{A}}_2$ and ${\mathcal{A}}_3$ of the first derivative and the errors of numerical solutions NS2 and NS3 of equation (10). Denote $M={max}_{x\in [0,X]}\left|f^{\prime \prime }\left(x\right)\right|.$

Lemma 2. 

let $f\in C^2[0,X]$. Then

$${\mathcal{A}}_2\left[f_n\right]={\displaystyle \frac{1-a}{h}}\left(f_n-(1-a)\sum _{k=1}^{n-1}{a}^{k-1}{f}_{n-k}-a^{n-1}{f}_0\right)=(1-a^n)f_n^{\prime }+E_n^{1}h,$$

(20)

where the coefficient $E_n^1$ satisfies the estimate

$$|E_n^1|<{\displaystyle \frac{(1+a)M}{2(1-a)}},(2\le n\le N).$$

(21)

Proof. 

From Taylor theorem

$$\begin{array}{c}\hfill f_{n-k}=f_n-kh{f}_n^{\prime }+k^2{\displaystyle \frac{h^2}{2}}{f}^{\prime \prime }\left({\theta }_{n-k}\right)\end{array}$$

(22)

for $k=1,\cdots ,n$ and ${\theta }_{n-k}\in (x_{n-k},X)$. By substituting $f_{n-k}$ with (22) in formula (20) for ${\mathcal{A}}_2$ we get

$${\mathcal{A}}_2\left[f\right]=K_0{f}_n+K_1{f}_n^{\prime }+E_n^{1}h,$$

where

$$K_0=1-a-{(1-a)}^2\sum _{k=1}^{n-1}{a}^{k-1}-(1-a)a^{n-1}=0,$$

$$K_1={(1-a)}^2\sum _{k=1}^{n-1}k{a}^{k-1}+(1-a)n{a}^{n-1},$$

$$E_n=-{\displaystyle \frac{1}{2}}\left({(1-a)}^2\sum _{k=1}^{n-1}{k}^2{a}^{k-1}{f}^{\prime \prime }\left({\theta }_{n-k}\right)+(1-a)n^2{a}^{n-1}{f}^{\prime \prime }\left({\theta }_0\right)\right).$$

From (14) and (15) we get

$$K_1=1+(n-1)a^n-n{a}^{n-1}+n{a}^{n-1}-n{a}^n=1-a^n,$$

$$|E_n|\le {\displaystyle \frac{M}{2}}\left({(1-a)}^2\sum _{k=1}^{n-1}{k}^2{a}^{k-1}+(1-a)n^2{a}^{n-1}\right),$$

$$|E_n|\le {\displaystyle \frac{M}{2}}\left({\displaystyle \frac{1+a-n^2{a}^{n-1}+(2n^2-2n-1)a^n-{(n-1)}^2{a}^{n+1}}{1-a}}+(1-a)n^2{a}^{n-1}\right),$$

$$|E_n|\le {\displaystyle \frac{M}{2(1-a)}}\left(1+a-(2n+1)a^n+(2n-1)a^{n+1}\right)<{\displaystyle \frac{M(1+a)}{2(1-a)}}.$$

□

Lemma 3. 

let $f\in C^2[0,X]$. Then

$${\mathcal{A}}_3\left[f_n\right]={\displaystyle \frac{1}{h}}\left((1-a)f_n-{(1-a)}^2\sum _{k=1}^{n-2}{a}^{k-1}{f}_{n-k}-a^{n-2}(1-2a)f_1-a^{n-1}{f}_0\right)=f_n^{\prime }+E_n^{2}h,$$

where the coefficient $E_n^2$ satisfies the estimate

$$|E_n^2|<{\displaystyle \frac{(1+a)M}{2(1-a)}},(2\le n\le N).$$

(23)

Proof. 

Approximation ${\mathcal{A}}_3$ is written in the form

$${\mathcal{A}}_3\left[f_n\right]=K_0{f}_n+K_1{f}_n^{\prime }+E_n^{2}h,$$

where the coefficients $K_0,K_1$ and $E_n^2$ satisfy

$$\begin{array}{cc}\hfill K_0& =1-a-{(1-a)}^2\sum _{k=1}^{n-2}{a}^{k-1}-(1-2a)a^{n-2}-a^{n-1}\hfill \\ \hfill & =1-a-(1-a)(1-a^{n-2})-(1-2a)a^{n-2}-a^{n-1}=0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill K_1& ={(1-a)}^2\sum _{k=1}^{n-2}k{a}^{k-1}+(1-2a)(n-1)a^{n-2}+n{a}^{n-1}\hfill \\ \hfill & =1+(n-2)a^{n-1}-(n-1)a^{n-2}+(1-2a)(n-1)a^{n-2}+n{a}^{n-1}=1,\hfill \end{array}$$

$$E_n^2=-{\displaystyle \frac{1}{2}}\left({(1-a)}^2\sum _{k=1}^{n-2}{k}^2{a}^{k-1}{f}^{\prime \prime }\left({\theta }_{n-k}\right)+(1-2a){(n-1)}^2{a}^{n-2}{f}^{\prime \prime }\left({\theta }_1\right)+n^2{a}^{n-1}{f}^{\prime \prime }\left({\theta }_0\right)\right),$$

$$\begin{array}{cc}\hfill |E_n^2|& \le {\displaystyle \frac{M}{2}}\left({(1-a)}^2\sum _{k=1}^{n-2}{k}^2{a}^{k-1}+(1-2a){(n-1)}^2{a}^{n-2}+n^2{a}^{n-1}\right)\hfill \\ \hfill & \le {\displaystyle \frac{M}{2}}{\displaystyle \frac{1+a-2a^n}{1-a}}<{\displaystyle \frac{M(1+a)}{2(1-a)}}.\hfill \end{array}$$

□

The estimate for the coefficients $E_n^1$ and $E_n^2$ derived in Lemma 2 and Lemma 3 depends on the step size h and approximations ${\mathcal{A}}_2$ and ${\mathcal{A}}_3$ can be defined for any interval $[l,X]$. Asymptotic formula ${\mathcal{A}}_2$ is a first-order approximation of the first derivative for $n>lnN$ while ${\mathcal{A}}_3$ is an approximation of the first derivative for all $n\ge 2$. Now we use the bound for the errors of approximations ${\mathcal{A}}_2$ and ${\mathcal{A}}_3$ to prove the convergence of numerical methods NS2 and NS3. Let $e_k=y_k-u_k$ be the error of numerical solution NS3 at the point $x_k$ for $k=1,\cdots ,N$. First we estimate the errors $e_1$ and $e_2$. From Taylor Theorem the solution of equation (10) satisfies

$$y_1=y_0+h{F}_1-{\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_1\right),$$

$$y_2=y_1+h{F}_2-{\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_2\right).$$

Then

$$e_1=-{\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_1\right),|e_1|\le {\displaystyle \frac{h^2}{2}}\left|y^{\prime \prime }\left({\theta }_1\right)\right|\le {\displaystyle \frac{M{h}^2}{2}},$$

$$e_2=e_1-{\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_2\right),|e_2|\le |e_1|+{\displaystyle \frac{h^2}{2}}\left|y^{\prime \prime }\left({\theta }_2\right)\right|\le M{h}^2.$$

In Theorem 4 and Theorem 5 we derive bounds for the errors of numerical solutions NS2 and NS3 of equation (10). The theorems are proved by induction.

Theorem 4. 

The errors of numerical solution NS3 satisfies

$$|e_n|<{\displaystyle \frac{MXh}{{(1-a)}^2}},(1\le n\le N).$$

(24)

Proof. 

The estimate (24) for the error of numerical method NS3 holds for $n=1$ and $n=2$. Suppose that (24) holds for $n=1,\cdots ,m-1$. The solution of equation (10) and the error of NS3 satisfy

$$y_m={\displaystyle \frac{h{F}_m}{1-a}}+(1-a)\sum _{k=1}^{m-2}{a}^{k-1}{y}_{m-k}+{\displaystyle \frac{a^{m-2}(1-2a)}{1-a}}{y}_1+{\displaystyle \frac{a^{m-1}}{1-a}}{y}_0+{\displaystyle \frac{E_m^2{h}^2}{1-a}},$$

$$e_m=(1-a)\sum _{k=1}^{m-2}{a}^{k-1}{e}_{m-k}+{\displaystyle \frac{a^{m-2}(1-2a)}{1-a}}{e}_1+E_m{h}^2,$$

where $E_k=E_k^2/(1-a)$ for $k=1,\cdots ,N$. The error $e_m$ is expressed with $e_{m-1}$ as

$$e_m=(1-a)e_{m-1}+a(1-a)\sum _{k=2}^{m-2}{a}^{k-2}{e}_{m-k}+{\displaystyle \frac{a^{m-2}(1-2a)}{1-a}}{e}_1+E_m{h}^2,$$

$$e_m=(1-a)e_{m-1}+a\left(e_{m-1}-{\displaystyle \frac{a^{m-3}(1-2a)}{1-a}}{e}_1-E_{m-1}{h}^2\right)+{\displaystyle \frac{a^{m-2}(1-2a)}{1-a}}{e}_1+E_m{h}^2,$$

$$e_m=e_{m-1}+(E_m-a{E}_{m-1})h^2.$$

(25)

By applying (25) successively $m-3$ times, we obtain

$$e_m=e_2+(E_m-a{E}_2)h^2+(1-a)h^2\sum _{k=3}^{m-1}{E}_k.$$

Denote

$$E=\underset{1\le k\le N}{max}\left|E_k\right|<{\displaystyle \frac{M(1+a)}{2{(1-a)}^2}}.$$

The error $e_m$ satisfies the estimate.

$$\begin{array}{cc}\hfill |e_m|& \le |e_2|+(m-2+a)E{h}^2<{\displaystyle \frac{M{h}^2}{{(1-a)}^2}}+{\displaystyle \frac{(1+a)M{h}^2}{2{(1-a)}^2}}(m-2+a)\hfill \\ \hfill & <{\displaystyle \frac{M{h}^2}{{(1-a)}^2}}+{\displaystyle \frac{M{h}^2}{{(1-a)}^2}}(m-2+a)<{\displaystyle \frac{Mm{h}^2}{{(1-a)}^2}}\le {\displaystyle \frac{MXh}{{(1-a)}^2}}.\hfill \end{array}$$

□

Now we estimate the error $e_2$ of numerical solution NS2.

$$y_2={\displaystyle \frac{f_{2}h}{1-a}}+(1-a)y_1+a{y}_0+{\displaystyle \frac{E_2^1{h}^2}{1-a}},$$

$$u_2={\displaystyle \frac{f_{2}h}{1-a}}+(1-a)u_1+a{u}_0.$$

The error $e_2$ satisfies

$$e_2=(1-a)e_1+{\displaystyle \frac{E_2^1{h}^2}{1-a}},$$

$$|e_2|\le (1-a)|e_1|+{\displaystyle \frac{|E_2^1|h^2}{1-a}}\le {\displaystyle \frac{M(1-a)h^2}{2}}+{\displaystyle \frac{M(1+a)h^2}{2{(1-a)}^2}},$$

$$|e_2|\le {\displaystyle \frac{M{h}^2}{2{(1-a)}^2}}\left(2-2a+3a^2-a^3\right)<{\displaystyle \frac{M{h}^2}{{(1-a)}^2}}.$$

Theorem 5. 

The error of numerical solution NS2 satisfies

$$|e_n|<\left({\displaystyle \frac{2a^{2}F}{1-a}}+{\displaystyle \frac{M(1+a)X}{2{(1-a)}^2}}\right)h,(1\le n\le N).$$

(26)

Proof. 

The estimate (26) for the error of numerical method NS2 holds for $n=1$ and $n=2$. Suppose that (26) holds for $n=1,\cdots ,m-1$. The solution of equation (10) satisfies

$$y_m={\displaystyle \frac{(1-a^m)F_{m}h}{1-a}}+(1-a)\sum _{k=1}^{m-1}{a}^{k-1}{y}_{m-k}+a^{m-1}{y}_0+{\displaystyle \frac{E_m^1{h}^2}{1-a}}.$$

Then

$$e_m=-{\displaystyle \frac{a^m{F}_{m}h}{1-a}}+(1-a)\sum _{k=1}^{m-1}{a}^{k-1}{e}_{m-k}+E_m{h}^2,$$

where $E_k=E_k^1/(1-a)$ for $k=1,\cdots ,N$. The error $e_m$ is expressed with $e_{m-1}$ as

$$e_m=-{\displaystyle \frac{a^m{F}_{m}h}{1-a}}+(1-a)e_{m-1}+a(1-a)\sum _{k=2}^{m-1}{a}^{k-2}{e}_{m-k}+E_m{h}^2,$$

$$e_m=-{\displaystyle \frac{a^n{F}_{m}h}{1-a}}+(1-a)e_{m-1}+a\left(e_{m-1}+{\displaystyle \frac{a^{m-1}h{F}_{m-1}}{1-a}}-E_{m-2}{h}^2\right)+E_m{h}^2,$$

$$\begin{array}{c}\hfill e_m=e_{m-1}+{\displaystyle \frac{a^n(F_{m-1}-F_m)h}{1-a}}+(E_m-a{E}_{m-1})h^2.\end{array}$$

(27)

By applying (27) successively $m-2$ times, we obtain

$$e_m=e_1+\left(a^2{F}_1-a^m{F}_n+(a-1)\sum _{k=2}^{m-1}{a}^k{F}_k\right){\displaystyle \frac{h}{1-a}}+(E_m-a{E}_2)h^2+(1-a)h^2\sum _{k=3}^{m-1}{E}_k.$$

Denote $E={max}_{1\le k\le N}\left|E_k\right|$ and $F={max}_{x\in [0,X]}\left|F\left(x\right)\right|$. Then

$$\begin{array}{cc}\hfill |e_m|& \le |e_1|+\left(a^2+a^m+(a-1)\sum _{k=2}^{m-1}{a}^k\right){\displaystyle \frac{Fh}{1-a}}+E{h}^2(1+a)+(1-a)(m-3))\hfill \\ \hfill & <{\displaystyle \frac{M{h}^2}{2}}+\left(a^2+a^m+a^2(1-a^{m-2}\right){\displaystyle \frac{Fh}{1-a}}+{\displaystyle \frac{(1+a)M{h}^2}{2{(1-a)}^2}}(m-2+a),\hfill \end{array}$$

$$|e_m|<\left({\displaystyle \frac{2a^{2}F}{1-a}}+{\displaystyle \frac{M(1+a)X}{2{(1-a)}^2}}\right)h.$$

□

4. Numerical Solutions of First Order ODEs

In this section we construct the numerical methods for first order ODEs which use approximation ${\mathcal{A}}_3$ of the first derivative and we obtain recursive formulas for computation of the numerical solution. The performance of the numerical methods is compared with the performance of the corresponding Euler methods. We show that the numerical methods can achieve an arbitrary order in the interval $(0,2]$ by using appropriate choices of the parameter of approximation ${\mathcal{A}}_3$.

Example 6. 

Consider the first-order linear ODE with constant coefficients

$$y^{\prime }\left(x\right)+Ly\left(x\right)=F\left(x\right),y\left(0\right)=y_0.$$

(28)

By approximating $y_n^{\prime }$ with first order backward difference approximation we obtain

$${\displaystyle \frac{y_n-y_{n-1}}{h}}+L{y}_n=F_n+\mathcal{O}\left(h\right).$$

The numerical solution is computed as $u_0=y_0$ and

$$u_n={\displaystyle \frac{1}{1+Lh}}\left(h{F}_n+u_{n-1}\right).$$

(NS4)

By approximating $y_n^{\prime }$ with approximation ${\mathcal{A}}_3$ we obtain

$${\displaystyle \frac{1-a}{h}}\left(y_n-(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{y}_{n-k}-{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{y}_1-{\displaystyle \frac{a^{n-1}}{1-a}}{y}_0\right)+L{y}_n=y_n^{\prime }+E_n^{2}h,$$

$$\left(1+{\displaystyle \frac{Lh}{1-a}}\right)y_n={\displaystyle \frac{h{F}_n}{1-a}}+(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{y}_{n-k}+{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{y}_1+{\displaystyle \frac{a^{n-1}}{1-a}}{y}_0+{\displaystyle \frac{E_n^2{h}^2}{1-a}}.$$

Denote $H=h/(1-a)$. The solution of equation (28) satisfies

$$y_n={\displaystyle \frac{1}{1+LH}}\left(F_{n}H+(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{y}_{n-k}+{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{y}_1+{\displaystyle \frac{a^{n-1}}{1-a}}{y}_0+{\displaystyle \frac{E_n^2{h}^2}{1-a}}\right).$$

The numerical solution of equation (28) is computed as

$$u_n={\displaystyle \frac{1}{1+LH}}\left(F_{n}H+(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{u}_{n-k}+{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{u}_1+{\displaystyle \frac{a^{n-1}}{1-a}}{u}_0\right).$$

(29)

The numerical solution (29) is calculated with $\mathcal{O}\left(N\right)$ arithmetic operations using an algorithm that is similar to the algorithms for calculating NS2 and NS3. From formula (29), the value of $u_n$ is expressed in terms of $u_{n-1}$ as follows.

$$(1+LH)u_n=F_{n}H+(1-a)u_{n-1}+a\left((1-a)\sum _{k=2}^{n-2}{a}^{k-1}{u}_{n-k}+{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{u}_1+{\displaystyle \frac{a^{n-1}}{1-a}}{u}_0\right),$$

$$(1+LH)u_n=F_{n}H+(1-a)u_{n-1}+a\left(u_{n-1}-F_{n-1}H\right).$$

The numerical solution of equation (28) using the approximation ${\mathcal{A}}_3$ for the first derivative is computed recursively as $u_0=y_0$ and

$$u_n={\displaystyle \frac{1}{1+LH}}\left((1+aLH)u_{n-1}+H(F_n-a{F}_{n-1})\right).$$

(NS5)

Consider the following first-order linear ODE

$$y^{\prime }\left(x\right)+Ly\left(x\right)=(1+L)e^x,y\left(0\right)=y_0.$$

(30)

Experimental results for the error and order of numerical solutions NS4 and NS5 of equation (30) and positive values of L are are given in Table 3. Now we prove the convergence of numerical method NS5 when the value of L is positive. Denote by $e_n=y_n-u_n$ the error of NS5 at the point $x_n$. First, we estimate the errors $e_1$ and $e_2$. The solution of equation (28) satisfies

$$y_1={\displaystyle \frac{1}{1+LH}}\left((1+aLH)y_0+H(F_1-a{F}_0)\right)+{\epsilon }_1,$$

(31)

$$y_2={\displaystyle \frac{1}{1+LH}}\left((1+aLH)y_1+H(F_2-a{F}_1)\right)+{\epsilon }_2.$$

(32)

From Taylor theorem the truncation error ${\epsilon }_1$ satisfies

$$\begin{array}{cc}\hfill (1+LH){\epsilon }_1& =(1+LH)y_1-H{F}_1-(1+aLH)y_0-aH{F}_0\hfill \\ \hfill & =y_1-H{y}_1^{\prime }-y_0+aH{y}_0^{\prime }=h{y}_0^{\prime }-H{y}_1^{\prime }+aH{y}_0^{\prime }+{\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_1\right)\hfill \\ \hfill & ={\displaystyle \frac{h}{1-a}}(y_0^{\prime }-y_1^{\prime })+{\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_1\right)={\displaystyle \frac{h^2}{2}}{y}^{\prime \prime }\left({\theta }_1\right)-{\displaystyle \frac{h^2}{1-a}}{y}^{\prime \prime }\left({\theta }_2\right).\hfill \end{array}$$

Hence

$$|{\epsilon }_1|\le {\displaystyle \frac{h^2}{2}}|y^{\prime \prime }\left({\theta }_1\right)|+{\displaystyle \frac{h^2}{1-a}}\left|y^{\prime \prime }\left({\theta }_2\right)\right|\le {\displaystyle \frac{(3-a)M{h}^2}{2(1-a)}}.$$

In a similar way, we show that $|{\epsilon }_2|\le {\displaystyle \frac{(3-a)M{h}^2}{2(1-a)}}$. The errors $e_1$ and $e_2$ satisfy the estimates

$$e_1={\epsilon }_1,\left|e_1\right|\le {\displaystyle \frac{(3-a)M{h}^2}{2(1-a)}},$$

$$e_2={\displaystyle \frac{1+aLH}{1+LH}}{e}_1+{\epsilon }_2,|e_2|<|e_1|+|{\epsilon }_2|<{\displaystyle \frac{(3-a)M{h}^2}{1-a}}.$$

We will prove the convergence of the numerical method NS5 of equation (30) for positive values of the parameter L and derive an estimate for the error of the method. Denote

$$C=max\left\{{\displaystyle \frac{M{(1+a)}^2}{2L{(1-a)}^2}},{\displaystyle \frac{(3-a)Mh}{2(1-a)}}\right\}.$$

Theorem 7. 

Let $L>0$. Then

$$|e_n|<Ch,(1\le n\le N).$$

(33)

Proof. 

Inequality (33) is satisfied for $n=1$ and $n=2$. Assume that (33) holds for $n\le m-1$. The error $e_m$ satisfies

$$e_m={\displaystyle \frac{1+aLH}{1+LH}}{e}_{m-1}+{\displaystyle \frac{(E_m^2-a{E}_{m-1}^2)h^2}{(1+LH)(1-a)}},$$

$$|e_m|<{\displaystyle \frac{1+aLH}{1+LH}}\left|e_{m-1}\right|+{\displaystyle \frac{\left(\right|E_m^2|+a|E_{m-1}^2|h^2}{(1+LH)(1-a)}}.$$

From Lemma 3:

$$|e_m|<\left(1-{\displaystyle \frac{(1-a)LH}{1+LH}}\right)\left|e_{m-1}\right|+{\displaystyle \frac{{(1+a)}^{2}M{h}^2}{2{(1-a)}^2(1+LH)}}.$$

Using the assumption of induction

$$|e_m|<Ch-{\displaystyle \frac{Lh}{1+LH}}Ch+{\displaystyle \frac{{(1+a)}^{2}M{h}^2}{2{(1-a)}^2(1+LH)}}\le Ch.$$

By the principle of induction, the bound (33) holds for all $n=1,\cdots ,N$. □

We will compare the performance of numerical methods NS4 and NS5 and negative values of the parameter L of equation (30). Experimental results for the error and order of numerical solutions NS4 and NS5 and $L=-20,-30,-50$ are given in Tables 4-6. The results of the numerical experiments presented in this section refer to ODSs that are defined in the interval $[0,1]$. The numerical results in Table 4 show that method NS4 has a first-order accuracy and its error can be very large for small values of the step size h of the method. In Table 5 the values of the parameter a are fixed numbers close to one. In Table 6 the parameter $a=1-h^{0.8}$ and numerical solution NS5 has an order $0.2$. Although the order of NS5 in Tables 5-6 is smaller than one, the errors of the method are significantly smaller than the corresponding errors of numerical solution NS4 in Table 4.

Example 8. 

Consider the first-order linear ODE

$$y^{\prime }\left(x\right)+G\left(x\right)y\left(x\right)=F\left(x\right),y\left(0\right)=y_0.$$

(34)

The numerical solutions NS6 and NS7 for equation (34) which use backward difference approximation and approximation ${\mathcal{A}}_3$ for first derivative are computed as

$$u_n={\displaystyle \frac{1}{1+h{G}_n}}(h{F}_n+u_{n-1}),u_0=y_0,$$

(NS6)

$$u_n={\displaystyle \frac{1}{1-a+h{G}_n}}\left(h{F}_n+{(1-a)}^2\sum _{k=1}^{n-2}{a}^{k-1}{u}_{n-k}+(1-2a)a^{n-2}{u}_1+a^{n-1}{u}_0\right).$$

(35)

In a similar way as in Example 7, a recursive formula for calculating the numerical solution of an equation (34) is derived from (35)

$$u_n={\displaystyle \frac{1}{1-a+h{G}_n}}\left((1-a+ah{G}_{n-1})u_{n-1}+h(F_n-a{F}_{n-1})\right),u_0=y_0.$$

(NS7)

The first-order liner ODE

$$y^{\prime }\left(x\right)+tanxy\left(x\right)=secx,y\left(0\right)=-1$$

(36)

has a solution $y\left(x\right)=sinx-cosx$. The experimental results of numerical solutions NS6 and NS7 of equation (36) are given in Table 7.

Numerical solution NS7 depends on a parameter a. The method can have any order in the interval $(0,2]$ with a suitable choice of the parameter. The coefficient $C_a$ in an asymptotic formula (3) has a value $C_a=(1+a)/(1-a)$. The order of NS7 is $1+p$ when the coefficient $C_a$ is equal to $s{h}^p$.

$${\displaystyle \frac{1+a}{1-a}}=s{h}^p,1+a=s{h}^p-as{h}^p,a={\displaystyle \frac{s{h}^p-1}{s{h}^p+1}}.$$

When $a=-1+s{h}^p$ the coefficient $C_a=s{h}^p/(2-s{h}^p)$ and NS7 has an order $1+p$. When $a=1-s{h}^p$ the coefficient $C_a=(2-s{h}^p)h^{-p}/s$ and NS7 has an order $1-p$. A summary of these results is given below.

• Numerical method NS7 has an order $1+p$ when $s>0$ and $-1<p\le 1$,

$$\begin{array}{c}\hfill a={\displaystyle \frac{s{h}^p-1}{s{h}^p+1}}.\end{array}$$

(37)

Values of the parameter p greater than one, $p>1$ lead to a second-order method. The experimental results for $p=-0.3,p=0.7$ and $p=1$ are given in Table 8.

•NS7 has an order $1-p$, when $s>0$ and $0<p<1$,

$$\begin{array}{c}\hfill a=1-s{h}^p.\end{array}$$

(38)

The experimental results for $p=0.2,p=0.5$ and $p=0.7$ are given in Table 9.

•NS7 has an order $1+p$, when $s>0$ and $0<p\le 1$,

$$\begin{array}{c}\hfill a=-1+s{h}^p.\end{array}$$

(39)

When $p>1$, numerical solution NS7 has a second-order accuracy. The experimental results for $p=0.25,p=0.75$ and $p=1.5$ are given in Table 10.

The logistic equation is a first-order nonlinear ordinary differential equation used for modeling population dynamics [29,30,31].

Example 9. 

Consider the logistic equation

$$y^{\prime }\left(x\right)=Ly\left(x\right)(1-y\left(x\right)),y\left(0\right)=y_0.$$

(40)

Equation (40) has a solution

$$y={\displaystyle \frac{y_0}{y_0+(1-y_0)e^{-Lx}}}.$$

The explicit schemes of equation (40) which use first-order backward difference and approximation ${\mathcal{A}}_3$ of the first derivative satisfy

$$u_n=u_{n-1}+Lh{u}_{n-1}(1-u_{n-1}),u_0=y_0,$$

(NS8)

$$u_n={\displaystyle \frac{Lh}{1-a}}{u}_{n-1}(1-u_{n-1})+(1-a)\sum _{k=1}^{n-2}{a}^{k-1}{u}_{n-k}+{\displaystyle \frac{(1-2a)a^{n-2}}{1-a}}{u}_1+{\displaystyle \frac{a^{n-1}}{1-a}}{u}_0.$$

(41)

Numerical method (41) has initial conditions $u_0=y_0,u_1=u_0+Lh{u}_0(1-u_0)$ and is computed recursively as

$$u_n=u_{n-1}+{\displaystyle \frac{Lh}{1-a}}(u_{n-1}(1-u_{n-1})-a{u}_{n-2}(1-u_{n-2})).$$

(NS9)

Experimental results for the error and order of numerical solutions NS8 and NS9 of equation (40) are given in Table 11.

5. Numerical Solutions of Heat Equation

The heat equation is a partial differential equation that models the distribution of temperature in a one-dimensional object. Heat conduction a classical problem in physics and engineering that can be generalized to higher dimensions and has a wide range of applications. The analytical and numerical solutions of the heat equation are extensively studied in [32,33,34,35,36]. In this section we construct a finite difference scheme for the heat equation which uses approximation ${\mathcal{A}}_3$ for the partial derivative in time and second order difference approximation for the partial derivative in space. The stability and convergence of the method is proved and its performance is compared with the standard method using first order difference approximation for the partial derivative in time.

Example 10. 

Consider the following heat equation

$$\left\{\begin{array}{c}{u}_t(x,z)=D{u}_{xx}(x,t)+F(x,t),(x,t)\in \mathcal{R},\hfill \\ u(x,0)=u_0\left(x\right),u(0,t)=u_1\left(t\right),u(1,t)=u_2\left(t\right),\hfill \end{array}\right.$$

(42)

where D is the thermal diffusivity of the material and $\mathcal{R}=[0,X]\times [0,T]$. Denote $h=X/N,\tau =T/M$, where M and N are positive integers and by $u_n^m=u(nh,m\tau )$ the values of the solution at the nodes $(x_n,y_m)=(nh,m\tau )$ of a rectangular grid $\mathcal{J}$ on $\mathcal{R}$. By approximating the partial derivative in time by first-order backward difference approximation and the partial derivative in space by second-order central difference approximation we obtain

$${\displaystyle \frac{u_n^m-u_n^{m-1}}{\tau }}=D{\displaystyle \frac{u_{n-1}^m-2u_n^m+u_{n+1}^m}{h^2}}+F_n^m+C_n^m\tau +K_n^m{h}^2,$$

$$-{\eta }_0{u}_{n-1}^m+(1+2{\eta }_0)u_n^m-{\eta }_0{u}_{n+1}^m=u_n^{m-1}+\tau F_n^m+\tau (C_n^m\tau +K_n^m{h}^2).$$

(43)

where ${\eta }_0=D\tau /h^2$ and $C_n^m\tau +K_n^m{h}^2$ is the truncation error of the approximations of the partial derivatives at the point $(x_n,y_m)$ of the grid $\mathcal{J}$. The coefficients $C_n^m$ and $K_n^m$ satisfy $|C_n^m|<M_2/2,\left|K_n^m\right|<M_4/12$, where

$$M_2=\underset{(x,t)\in \mathcal{R}}{max}\left|{\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial x^2}}\right|,M_4=\underset{(x,t)\in \mathcal{R}}{max}\left|{\displaystyle \frac{{\partial }^{4}u(x,t)}{\partial x^4}}\right|.$$

The numerical solution $\left\{U_n^m\right\}$ of heat equation (42) on layer m of the grid $\mathcal{J}$ is a solution of the system of linear equations

$$-{\eta }_0{U}_{n-1}^m+(1+2{\eta }_0)U_n^m-{\eta }_0{U}_{n+1}^m=U_n^{m-1}+\tau F_n^m$$

(44)

and has initial conditions $U_n^0=u_0\left(nh\right)$ for $n=0,\cdots ,N$ and boundary conditions

$$U_0^m=u_1\left(m\tau \right),U_N^m=u_2\left(m\tau \right),(1\le m\le M).$$

The system of linear equations (44) is written in matrix form as

$$(I+\eta \Delta ){\mathcal{U}}_m={\mathcal{U}}_{m-1}+\tau {\mathcal{F}}_m+{\eta }_0{\mathcal{I}}_m,$$

(NS10)

where ${\mathcal{U}}_m$ and ${\mathcal{F}}_m$ are vectors of dimension $N-1$ which have entries $U_n^m$ and $F_n^m$ for $n=1,\cdots ,N-1$ and ${\mathcal{I}}_m={\left(u_1\left(m\tau \right),0,\cdots ,0,u_2\left(m\tau \right)\right)}^T$. The matrix $\Delta$ is the tridiagonal matrix of dimension $N-1$ with main diagonal entries equal to 2 and entries above and below the main diagonal equal to $-1$. Numerical solution NS10 has an accuracy $O(\tau +h^2)$ and it is computed with $\mathcal{O}\left(MN\right)$ operations because $I+\eta \Delta$ is a tridiagonal M-matrix [37,38]. Now we construct the finite difference scheme of equation (42) which uses approximation ${\mathcal{A}}_3$ for the partial derivative in time and central difference approximation for the partial derivative in space. The solution of heat equation (42) satisfies

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1-a}{\tau }}\left(u_n^m-(1-a)\sum _{k=2}^{m-2}{a}^{k-1}{u}_n^{m-k}-{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{u}_n^1-{\displaystyle \frac{a^{n-1}}{1-a}}{u}_n^0\right)\hfill \\ \hfill & =D{\displaystyle \frac{u_{n-1}^m-2u_n^m+u_{n+1}^m}{h^2}}+S_n^m\tau +K_n^m{h}^2+F_n^m.\hfill \end{array}$$

Denote $\eta =D\tau /\left((1-a)h^2\right)$.

$$\begin{array}{cc}\hfill & -\eta u_{n-1}^m+(1+2\eta )u_n^m-\eta u_{n+1}^m=(1-a)\sum _{k=1}^{m-2}{a}^{k-1}{u}_n^{m-k}+{\displaystyle \frac{1-2a}{1-a}}{a}^{n-2}{u}_n^1+{\displaystyle \frac{a^{n-1}}{1-a}}{u}_n^0\hfill \\ \hfill & +{\displaystyle \frac{\tau F_n^m}{1-a}}+{\displaystyle \frac{\tau }{1-a}}(S_n^m\tau +K_n^m{h}^2).\hfill \end{array}$$

(45)

From Lemma 3 the coefficient $S_n^m$ satisfies the estimate $|S_n^m|<{\displaystyle \frac{(1+a)M_2}{2(1-a)}}$. The formula (45) is written in the form

$$\begin{array}{cc}\hfill & -\eta u_{n-1}^m+(1+2\eta )u_n^m-\eta u_{n+1}^m=a\left[(1-a)\sum _{k=2}^{m-2}{a}^{k-2}{u}_n^{m-k}+{\displaystyle \frac{1-2a}{1-a}}{a}^{n-3}{u}_n^1+{\displaystyle \frac{a^{n-2}}{1-a}}{u}_n^0\right]\hfill \\ \hfill & +(1-a)u_n^{m-1}+{\displaystyle \frac{\tau F_n^m}{1-a}}+{\displaystyle \frac{\tau }{1-a}}(S_n^m\tau +K_n^m{h}^2).\hfill \end{array}$$

(46)

From formulas (45) and (46) we obtain

$$\begin{array}{cc}\hfill & -\eta u_{n-1}^m+(1+2\eta )u_n^m-\eta u_{n+1}^m=a\left[-\eta u_{n-1}^{m-1}+(1+2\eta )u_n^{m-1}-\eta u_{n+1}^{m-1}\right]+(1-a)u_n^{m-1}\hfill \\ \hfill & -{\displaystyle \frac{a\tau F_n^{m-1}}{1-a}}-{\displaystyle \frac{a\tau }{1-a}}(S_n^{m-1}\tau +K_n^{m-1}{h}^2)+{\displaystyle \frac{\tau F_n^m}{1-a}}+{\displaystyle \frac{\tau }{1-a}}(S_n^m\tau +K_n^m{h}^2).\hfill \end{array}$$

The solution of heat equation (51) satisfies

$$\begin{array}{cc}\hfill & -\eta u_{n-1}^m+(1+2\eta )u_n^m-\eta u_{n+1}^m=-a\eta u_{n-1}^{m-1}+(1+2a\eta )u_n^{m-1}-a\eta u_{n+1}^{m-1}\hfill \\ \hfill & +{\displaystyle \frac{\tau }{1-a}}\left({\overline{F}}_n^m+{\overline{S}}_n^m\tau +{\overline{K}}_n^m{h}^2\right),\hfill \end{array}$$

(47)

where ${\overline{F}}_n^m=F_n^m-a{F}_n^{m-1},{\overline{S}}_n^m=S_n^m-a{S}_n^{m-1},{\overline{K}}_n^m=K_n^m-a{K}_n^{m-1}$. The numerical solution of heat equation (42) on row m of the grid $\mathcal{J}$ which corresponds to (47) satisfies

$$\begin{array}{c}\hfill -\eta U_{n-1}^m+(1+2\eta )U_n^m-\eta U_{n+1}^m=-a\eta U_{n-1}^{m-1}+(1+2a\eta )U_n^{m-1}-a\eta U_{n+1}^{m-1}+{\displaystyle \frac{\tau {\overline{F}}_n^m}{1-a}}.\end{array}$$

(48)

The system of liner equations (58) is written in matrix form as

$$(I+\eta \Delta ){\mathcal{U}}_m=(I+a\eta \Delta ){\mathcal{U}}_{m-1}+\tau {\overline{\mathcal{F}}}_m+\eta {\mathcal{I}}_m,$$

(NS11)

where ${\overline{\mathcal{F}}}_m$ is the vector of dimension $N-1$ with elements ${\overline{F}}_n^m/(1-a)$. Consider the following heat equation

$$\left\{\begin{array}{c}{u}_t(x,z)=3u_{xx}(x,t)-2e^{x+t},\hfill \\ u(x,0)=e^x,u(0,t)=e^t,u(1,t)=e^{1+t}.\hfill \end{array}\right.$$

(49)

where $(x,t)\in [0,1]\times [0,1]$. Experimental results for error and order of numerical methods 44 and NS11 of equation (49) and $h=\tau$ are given in Table 12. The orders of the methods with respect to $\tau$ and h are computed with formulas ${log}_2|e_{\tau }/e_{2\tau }|$ and ${log}_2|(e_h-2e_{2h})/(e_{2h}-2e_{4h})|$, where $e_{\tau }=e_h$ is he error of the method for the given values of $\tau =h$ in Table 12. The two methods, 44 and NS11, have similar performance with respect to computational time and accuracy. Numerical solution NS11 has an order $1+p$ when the parameter $a=({\tau }^p-1)/({\tau }^p+1)$. The numerical results for $p=-0.5,0.5,1$ and $h=2\tau$ are given in Table 13.

A Crank-Nicolson scheme for heat equation has second order accuracy [39,40]. The convergence of finite difference schemes for heat equation are studied in [41,42]. Now we derive a bound for the error of finite difference scheme NS11. Let $e_n^m=U_n^m-u_n^m$ be the error of 59 at the point $(nh,m\tau )$ and ${\mathcal{E}}_m$ be the vector with the errors $e_n^m$ on row m of the grid $\mathcal{J}$, where $n=1,\cdots ,N-1$. The infinity norm of a vector is the maximum absolute value of its components and the infinity norm of a matrix is the maximum of the absolute row sums. In Lemma 11 we estimate the errors of difference scheme NS11 on the first row of the grid $\mathcal{J}$.

Lemma 11. 

Let $u\in C^4[0,X]\times C^2[0,T]$. Then

$$\begin{array}{c}\hfill \parallel {\mathcal{E}}_1\parallel <{\displaystyle \frac{(1+a)\tau }{12(1-a)}}\left(6M_2\tau +M_4{h}^2\right).\end{array}$$

(50)

Proof. 

The solution of equation (42) on the first row of the grid $\mathcal{J}$ satisfies

$$\begin{array}{cc}\hfill -\eta u_{n-1}^1+(1+2\eta )u_n^1-\eta u_{n+1}^1=-a\eta u_{n-1}^0+(1+2a\eta )& u_n^0-a\eta u_{n+1}^0\hfill \\ \hfill & +{\displaystyle \frac{\tau }{1-a}}\left(F_n^1-a{F}_n^0\right)+{\epsilon }_n^1,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\displaystyle \frac{(1-a)(u_n^1-u_n^0)}{\tau }}={\displaystyle \frac{D}{h^2}}\left(u_{n-1}^1-2u_n^1+u_{n+1}^1\right)-& {\displaystyle \frac{aD}{h^2}}\left(u_{n-1}^0-2u_n^0+u_{n+1}^0\right)\hfill \\ \hfill & +F_n^1-a{F}_n^0+{\displaystyle \frac{1-a}{\tau }}{\epsilon }_n^1.\hfill \end{array}$$

The first order difference approximation satisfies

$${\displaystyle \frac{u_n^1-u_n^0}{\tau }}=u_t(x_n,y_1)+C_1\tau =u_t(x_n,y_0)+C_2\tau ,$$

where $|C_i|<M_2/2$ for $i=1,2$. Then

$$\begin{array}{cc}\hfill u_t(x_n,y_1)+C_1\tau -a{u}_t(x_n,y_0)-a{C}_2\tau =D{u}_{xx}& (x_n,y_1)+C_3{h}^2-aD{u}_{xx}(x_n,y_0)-a{C}_4{h}^2\hfill \\ \hfill & +F(x_n,y_1)-aF(x_n,y_0)+{\displaystyle \frac{1-a}{\tau }}{\epsilon }_n^1,\hfill \end{array}$$

where $|C_i|<M_4/12$ for $i=3,4$. The error ${\epsilon }_n^1$ satisfies

$$\begin{array}{c}\hfill {\displaystyle \frac{1-a}{\tau }}{\epsilon }_n^1=-C_1\tau +a{C}_2\tau +C_3{h}^2-a{C}_4{h}^2,\end{array}$$

$$\begin{array}{c}\hfill |{\epsilon }_n^1|\le {\displaystyle \frac{\tau }{1-a}}\left(|C_1|\tau +a|C_2|\tau +|C_3|h^2+a\left|C_4\right|h^2\right).\end{array}$$

Hence

$$\begin{array}{c}\hfill |{\epsilon }_n^1|<{\displaystyle \frac{(1+a)\tau }{12(1-a)}}\left(6M_2\tau +M_4{h}^2\right)\end{array}$$

(51)

for $n=1,\cdots ,N-1$ and $e_0^1=e_N^1=0$. The errors of numerical method NS11 on the first row of the grid $\mathcal{J}$ satisfy

$$\begin{array}{c}\hfill -\eta e_{n-1}^1+(1+2\eta )e_n^1-\eta e_{n+1}^1={\epsilon }_n^1.\end{array}$$

The elements of ${\mathcal{E}}_0$ are equal to zero, $e_n^0=0$ and the vector ${\mathcal{E}}_1$ with the errors of 59 on the first row of the grid $\mathcal{J}$ satisfies

$$(I+\eta \Delta ){\mathcal{E}}_1={\mathcal{W}}_1,$$

(52)

where ${\mathcal{W}}_1$ is the vector with elements ${\epsilon }_n^1$. From (62) the infinity norm of ${\mathcal{W}}_1$ satisfies

$$\parallel {\mathcal{W}}_1\parallel <{\displaystyle \frac{(1+a)\tau }{12(1-a)}}\left(6M_2\tau +M_4{h}^2\right).$$

The matrix $I+\eta \Delta$ is a diagonally dominant M-matrix and its infinity norm satisfies [43,44]

$$\parallel {(I+\eta \Delta )}^{-1}\parallel <1.$$

Therefore

$$\parallel {\mathcal{E}}_1\parallel \le \parallel {(I+\eta \Delta )}^{-1}\parallel ·\parallel {\mathcal{W}}_1\parallel <{\displaystyle \frac{(1+a)\tau }{12(1-a)}}\left(6M_2\tau +M_4{h}^2\right).$$

□

The vector ${\mathcal{E}}_m$ with the errors of difference scheme 59 on row m of the grid $\mathcal{J}$ satisfies

$$(I+\eta \Delta ){\mathcal{E}}_m=(I+a\eta \Delta ){\mathcal{E}}_{m-1}+\tau {\mathcal{W}}_m,$$

(53)

where ${\mathcal{W}}_m$ is the vector with entries ${\overline{S}}_n^m\tau +{\overline{K}}_n^m{h}^2$. The infinity norm of ${\mathcal{W}}_m$ satisfies

$$\parallel {\mathcal{W}}_m\parallel \le M_0(\tau +h^2)$$

where

$$M_0=max\left\{{\displaystyle \frac{{(1+a)}^2{M}_2}{2(1-a)}},{\displaystyle \frac{(1+a)M_4}{12}}\right\}.$$

The matrix $\Delta$ has eigenvalues [45]

$${\lambda }_k=4{sin}^2\left({\displaystyle \frac{k\pi }{2N}}\right),k=1,\cdots ,N-1$$

and eigenvectors ${\mathcal{V}}_k$ for $k=1,\cdots ,N-1$ which have components $v_k^i=sin(k\pi i/N)$, where $i=1,\cdots ,N-1$. Denote by P and Q the following matrices

$$P={(I+\eta \Delta )}^{-1}(I+a\eta \Delta ),Q={(I+\eta \Delta )}^{-1}.$$

The sequence of vectors ${\mathcal{E}}_m$ satisfies

$${\mathcal{E}}_m=P{\mathcal{E}}_{m-1}+\tau Q{\mathcal{W}}_m.$$

(54)

The matrix P has eigenvalues $(I+a\eta {\lambda }_k)/(I+\eta {\lambda }_k)$ for $k=1.\cdots ,N-1$. The spectral radius of the matrix P is smaller than one, $\rho \left(P\right)<1$ and its powers $P^n$ converge to zero. Therefore ${lim}_{n\to \infty }{\parallel P\parallel }^n=0$. Denote

$$p=\underset{n\ge 1}{max}{\parallel P\parallel }^n,q=\parallel Q\parallel .$$

In Theorem 12 we prove that finite difference scheme NS11 is unconditionally stable and has an accuracy $O(\tau +h^2)$.

Theorem 12. 

Let $u\in C^4[0,X]\times C^2[0,T]$. Then

$$\begin{array}{c}\hfill \parallel {\mathcal{E}}_m\parallel <C\left(\tau +h^2\right),\end{array}$$

where $C=p(1+Tq)M_0$ and all $m=2,\cdots M$.

Proof. 

By applying (54) successively $m-1$ times we obtain

$${\mathcal{E}}_m=P^{m-2}{\mathcal{E}}_1+\tau \sum _{k=0}^{m-2}{P}^{k}Q{\mathcal{W}}_{m-k}.$$

Then

$$\parallel {\mathcal{E}}_m\parallel \le \parallel P^{m-2}\parallel ·\parallel {\mathcal{E}}_1\parallel +\tau \sum _{k=0}^{m-2}\parallel P^k\parallel ·\parallel Q\parallel ·\parallel {\mathcal{W}}_{m-k}\parallel .$$

From Lemma 11 the norm of the vector ${\mathcal{E}}_1$ satisfies $\parallel {\mathcal{E}}_1\parallel <M_0(\tau +h^2)$. Hence

$$\parallel {\mathcal{E}}_m\parallel <p{M}_0(\tau +h^2)+pq{M}_0\tau (\tau +h^2)(m-1)<p(1+Tq)M_0(\tau +h^2).$$

□

6. Approximations of the Fractional Derivative

In this section, we discuss the construction of an approximation (5) for the fractional derivative, using the asymptotic formula (2) and the approximation ${\mathcal{A}}_3$ for the first derivative. Two approximations of the fractional derivative of order $2-\alpha$ whose weights satisfy (1) are given below. Let $h=(x-l)/N$. The L1 approximation of Caputo fractional derivative has the form

$${\displaystyle \frac{1}{\mathrm{\Gamma }(2-\alpha )h^{\alpha }}}\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+\mathcal{O}\left(h^{2-\alpha }\right),$$

(55)

where ${\sigma }_0^{\left(\alpha \right)}=1,{\sigma }_n^{\left(\alpha \right)}={(n-1)}^{1-\alpha }-n^{1-\alpha }$,

$${\sigma }_k^{\left(\alpha \right)}={(k+1)}^{1-\alpha }-2k^{1-\alpha }+{(k-1)}^{1-\alpha },(1\le k\le n-1).$$

The L1 approximation and approximation (5) have an order $2-\alpha$ and and their weights satisfy (1). In [7,46] we use the properties of the weights (1) of approximations of the fractional derivative for analyzing the convergence and order of the numerical solutions of two-term ordinary fractional differential equation and the time-fractional Black-Scholes equation. In [26,28] we construct second order approximations of Caputo derivative using the asymptotic formula of L1 approximation. In [23] we obtain an approximation of the Caputo fractional derivative of order $2-\alpha$ by substituting the first derivative in (2) with a first-order backward difference approximation.

$${\displaystyle \frac{1}{2\mathrm{\Gamma }(1-\alpha )h^{\alpha }}}\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+\mathcal{O}\left(h^{2-\alpha }\right),$$

(56)

where ${\sigma }_0^{\left(\alpha \right)}=1-2\zeta \left(\alpha \right),{\sigma }_1^{\left(\alpha \right)}={\displaystyle \frac{1}{2^{\alpha }}}+2\zeta \left(\alpha \right)$,

$${\sigma }_k^{\left(\alpha \right)}={\displaystyle \frac{1}{{(k+1)}^{\alpha }}}-{\displaystyle \frac{1}{{(k-1)}^{\alpha }}},(2\le k\le n-2),$$

$${\sigma }_{n-1}^{\left(\alpha \right)}=-{\displaystyle \frac{1}{{(n-2)}^{\alpha }}}-2s_n,{\sigma }_n^{\left(\alpha \right)}=-{\displaystyle \frac{1}{{(n-1)}^{\alpha }}}+2s_n.$$

and $s_n={\sum }_{k=1}^{n-1}{k}^{-\alpha }-n^{1-\alpha }/(1-\alpha )-\zeta \left(\alpha \right)$. Now we apply approximation ${\mathcal{A}}_3$ for constructing an approximation of the Caputo derivative. By substituting the first derivative $f_n^{\prime }$ in formula (2) with ${\mathcal{A}}_3$ we obtain an approximation of order $2-\alpha$.

$${\displaystyle \frac{1}{2\mathrm{\Gamma }(1-\alpha )h^{\alpha }}}\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+\mathcal{O}\left(h^{2-\alpha }\right),$$

(57)

where ${\sigma }_0^{\left(\alpha \right)}=1-2(1-b)\zeta \left(\alpha \right),{\sigma }_1^{\left(\alpha \right)}={\displaystyle \frac{1}{2^{\alpha }}}+2{(1-b)}^2\zeta \left(\alpha \right)$,

$${\sigma }_k^{\left(\alpha \right)}={\displaystyle \frac{1}{{(k+1)}^{\alpha }}}-{\displaystyle \frac{1}{{(k-1)}^{\alpha }}}+2{(1-b)}^2{b}^{k-1}\zeta \left(\alpha \right),(2\le k\le n-2),$$

$${\sigma }_{n-1}^{\left(\alpha \right)}=-{\displaystyle \frac{1}{{(n-2)}^{\alpha }}}-2s_n+2(1-2b)b^{n-2}\zeta \left(\alpha \right),$$

$${\sigma }_n^{\left(\alpha \right)}=-{\displaystyle \frac{1}{{(n-1)}^{\alpha }}}+2s_n+2b^{n-1}\zeta \left(\alpha \right).$$

The value of zeta function $\zeta \left(\alpha \right)$ is negative when $0<\alpha <1$ and the parameter b of approximation ${\mathcal{A}}_3$ has a modulus smaller than one, $\left|b\right|<1$. Therefore ${\sigma }_0^{\left(\alpha \right)}>0$ and ${\sigma }_k^{\left(\alpha \right)}<0$ for $k>1$. In Lemma 13 we show that the coefficient ${\sigma }_1^{\left(\alpha \right)}$ of approximation (58) is negative when the parameter b has a value $b=\alpha /2$.

Lemma 13. 

Let $b=a/2$. Then

$$\begin{array}{c}\hfill {\sigma }_1^{\left(\alpha \right)}<0.\end{array}$$

(58)

Proof. 

$${\sigma }_1^{\left(\alpha \right)}={\displaystyle \frac{1}{2^{\alpha }}}+2{\left(1-{\displaystyle \frac{\alpha }{2}}\right)}^2\zeta \left(\alpha \right)={\displaystyle \frac{1}{2^{\alpha }}}-{\displaystyle \frac{1}{2}}{\left(2-a\right)}^2\left|\zeta \left(\alpha \right)\right|.$$

Inequality (58) holds when

$${\displaystyle \frac{1}{2^{\alpha }}}-{\displaystyle \frac{1}{2}}{\left(2-\alpha \right)}^2\left|\zeta \left(\alpha \right)\right|<0,$$

$$\left|\zeta \left(\alpha \right)\right|>{\displaystyle \frac{2}{2^{\alpha }{(2-\alpha )}^2}}.$$

Taking the logarithm of both sides, we get

$$ln\left|\zeta \right(\alpha \left)\right|>ln2-\alpha ln2-2ln(2-\alpha ).$$

Let $f\left(\alpha \right)=ln\left|\zeta \right(\alpha \left)\right|-ln2+\alpha ln2+2ln(2-\alpha )$. The function f satisfies $f\left(0\right)=0$ and has first derivative

$$f^{\prime }\left(\alpha \right)={\displaystyle \frac{{\zeta }^{\prime }\left(\alpha \right)}{\zeta \left(\alpha \right)}}+ln2-{\displaystyle \frac{2}{2-\alpha }}.$$

It is sufficient to prove that $f^{\prime }$ is positive. The zeta function and its derivative have Laurent series expansions

$$\begin{array}{cc}\hfill & \zeta \left(\alpha \right)=-{\displaystyle \frac{1}{1-a}}+\sum _{k=0}^{\infty }{\displaystyle \frac{{\gamma }_k}{k!}}{(1-a)}^k,\hfill \\ \hfill & {\zeta }^{\prime }\left(\alpha \right)=-{\displaystyle \frac{1}{{(1-\alpha )}^2}}-\sum _{k=0}^{\infty }{\displaystyle \frac{{\gamma }_{k+1}}{k!}}{(1-a)}^k,\hfill \end{array}$$

where ${\gamma }_k$ are Stieltjes constants

$${\gamma }_0=0,5772,{\gamma }_1=-0,0728,{\gamma }_2=-0,0097,{\gamma }_3=0,0021,{\gamma }_4=0,0023,\cdots .$$

Hence

$$\left|\zeta \left(\alpha \right)\right|<{\displaystyle \frac{1}{1-\alpha }},\left|{\zeta }^{\prime }\left(\alpha \right)\right|>{\displaystyle \frac{1}{{(1-\alpha )}^2}}+{\gamma }_1+{\gamma }_2(1-\alpha ),$$

$$|{\zeta }^{\prime }\left(\alpha \right)|>{\displaystyle \frac{1}{{(1-\alpha )}^2}}-0,08-0,01(1-\alpha ).$$

The derivative of the function f satisfies

$$f^{\prime }\left(\alpha \right)>{\displaystyle \frac{1}{1-\alpha }}-0,08(1-\alpha )-0,01{(1-\alpha )}^2+ln2-{\displaystyle \frac{2}{2-\alpha }}.$$

Denote

$$g\left(\alpha \right)={\displaystyle \frac{1}{1-\alpha }}-0,08(1-\alpha )-0,01{(1-\alpha )}^2+ln2-{\displaystyle \frac{2}{2-\alpha }}.$$

The function g has first, second and third order derivatives

$$g^{\prime }\left(\alpha \right)={\displaystyle \frac{1}{{(1-\alpha )}^2}}+0,08+0,02(1-\alpha )-{\displaystyle \frac{2}{{(2-\alpha )}^2}},$$

$$g^{\prime \prime }\left(\alpha \right)={\displaystyle \frac{2}{{(1-\alpha )}^3}}-0,02-{\displaystyle \frac{4}{{(2-\alpha )}^3}},$$

$$g^{\prime \prime \prime }\left(\alpha \right)={\displaystyle \frac{6}{{(1-\alpha )}^4}}-{\displaystyle \frac{12}{{(2-\alpha )}^4}}.$$

The values of the function g and its derivatives $g^{\prime },g^{\prime \prime },g^{\prime \prime \prime }$ at zero are positive

$$g\left(0\right)=ln2-0,09>0,g^{\prime }\left(0\right)=0,6,g^{\prime \prime }\left(0\right)=1,48,g^{\prime \prime \prime }\left(0\right)=5,25.$$

The third derivative $g^{\prime \prime \prime }$ is positive on the interval $[0,1]$ because

$${\left({\displaystyle \frac{2-\alpha }{1-\alpha }}\right)}^4>2,1+{\displaystyle \frac{1}{1-\alpha }}>\sqrt[4]{2},{\displaystyle \frac{1}{1-\alpha }}>1>\sqrt[4]{2}-1.$$

Therefore the functions $g,g^{\prime },g^{\prime \prime }$ are positive and increasing and $f^{\prime }\left(\alpha \right)>g\left(\alpha \right)>0$. The function f is increasing and positive because $f\left(0\right)=0$. □

Approximation (5) of the fractional derivative is obtained from (58) and value of the parameter $b=\alpha /2$. The two-term and multi-term ordinary fractional differential equations are an important class of equations in fractional calculus which generalize the linear ordinary differential equations with constant coefficients. Their analytical and numerical solutions are studied in [47,48,49,50,51]. In the next examples, we will consider an application of approximation (57) for construction of numerical schemes for the two-term ordinary fractional differential equation and the fractional subdiffusion equation.

Example 14. 

Consider the following two-term equation

$$y^{\left(\alpha \right)}\left(x\right)+Ly\left(x\right)=F\left(x\right),y\left(0\right)=y_0.$$

(59)

The numerical solution ${\left\{u_n\right\}}_{n=1}^N$ of equation (59) which uses L1 approximation (56) of the fractional derivative is computed as

$$u_m={\displaystyle \frac{1}{{\sigma }_0^{\left(\alpha \right)}+L\mathrm{\Gamma }(2-\alpha )h^{\alpha }}}\left(L\mathrm{\Gamma }(2-\alpha )h^{\alpha }{F}_m-\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{u}_{m-k}\right).$$

(NS12)

The numerical solution of equation (59) which uses approximations (56) and (57) of the fractional derivative $y_n^{\left(\alpha \right)}$ satisfies

$$u_m={\displaystyle \frac{1}{{\sigma }_0^{\left(\alpha \right)}+2L\mathrm{\Gamma }(1-\alpha )h^{\alpha }}}\left(2L\mathrm{\Gamma }(1-\alpha )h^{\alpha }{F}_m-\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{u}_{m-k}\right)$$

(NS13)

and has initial conditions

$$u_0=y_0,u_1={\displaystyle \frac{h^{\alpha }\mathrm{\Gamma }(2-\alpha )F_1+u_0}{L{h}^{\alpha }\mathrm{\Gamma }(2-\alpha )+1}},u_2={\displaystyle \frac{{\left(2h\right)}^{\alpha }\mathrm{\Gamma }(2-\alpha )F_2+u_0}{L{\left(2h\right)}^{\alpha }\mathrm{\Gamma }(2-\alpha )+1}}.$$

(60)

The initial conditions (Section 6) are obtained with the approximation of the fractional derivative

$${\displaystyle \frac{f(\ell )-f\left(0\right)}{\mathrm{\Gamma }(2-\alpha ){\ell }^{\alpha }}}=f^{\left(\alpha \right)}(\ell )+O\left({\ell }^{2-\alpha }\right).$$

(61)

Consider the two-term equation

$$y^{\left(\alpha \right)}\left(x\right)+Ly\left(x\right)=x^{1-\alpha }{E}_{1,2-\alpha }\left(x\right)+L{e}^x,y\left(0\right)=1.$$

(62)

Equation (62) has a solution $y\left(x\right)=e^x$, where $E_{1,2-\alpha }\left(x\right)$ is the Mittag-Leffler function

$$E_{1,2-\alpha }\left(x\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{x^k}{\mathrm{\Gamma }(k+2-\alpha )}}.$$

The experimental results of numerical solutions NS12 and NS13 of equation (62) on the interval $[0,1]$ are given in Table 14 and Table 15. The second column of Table 14 pertains to NS13 using approximation (57), while the results in the third column relate to NS13 using approximation (58) with a parameter $b=-0.5$. The numerical results in Table 15 were obtained using approximation (5) for the fractional derivative. The sequence $s_n$ and the coefficients ${\sigma }_{n-1}^{\left(\alpha \right)}$ and ${\sigma }_n^{\left(\alpha \right)}$ of approximations (56) and (57) for $n=3,\cdots ,N$ are computed recursively in $\mathcal{O}\left(N\right)$ arithmetic operations. The NS12 and NS13 methods have similar computational time and accuracy. The approximation (56) is obtained from the formula (2) by replacing the first derivative $f_n^{\prime }$ with a first-order difference approximation. Replacing the second derivative and derivatives of higher order with finite difference approximations in the asymptotic formulas leads to approximations whose weights do not satisfy (1). The method of constructing approximation (5) using the approximation ${\mathcal{A}}_3$ of the first derivative has the advantage that it can be generalized to constructions of high-order fractional derivative approximations which have properties (1).

The time fractional subdiffusion equation is a generalization of heat diffusion equation using a time-fractional derivative of order $\alpha$ with $a\in (0,1)$. The analytical and numerical solutions of fractional subdiffusion equations are studied in [52,53,54,55,56].

Example 15. 

Consider the following fractional subdiffusion equation

$$\left\{\begin{array}{c}{D}_t^{\alpha }u(x,t)=D{u}_{xx}(x,t)+F(x,t),\hfill \\ u(x,0)=u_0\left(x\right),u(0,t)=u_1\left(t\right),u(1,t)=u_2\left(t\right).\hfill \end{array}\right.$$

(63)

where $(x,t)\in [0,1]\times [0,1]$ and $D_t^{\alpha }u(x,t)$ denotes the Caputo derivative of the solution $u(x,t)$ in time. Let $\tau =1/M$ and $h=1/N$, where M and N are positive integers and $\mathcal{J}$ be a rectangular grid with nodes $(x_n,y_n)=(nh,mh)$. The numerical solution $U_n^m$ on the first three rows of $\mathcal{J}$ is computed with the approximation (61) of the partial fractional derivative and second-order central difference approximation of the partial derivative in space

$${\displaystyle \frac{U_n^m-U_n^0}{\mathrm{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }}}={\displaystyle \frac{D}{h^2}}(U_{n-1}^m-2U_n^m+U_{n+1}^m)+F_n^m,$$

$$U_n^m-U_n^0={\displaystyle \frac{D\mathrm{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }}{h^2}}(U_{n-1}^m-2U_n^m+U_{n+1}^m)+\mathrm{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }{F}_n^m.$$

Denote ${\eta }_m=D\mathrm{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }/h^2$. The numerical solution satisfies the system of linear equations

$$-{\eta }_m{U}_{n-1}^m+(1+2{\eta }_m)U_n^m-{\eta }_m{U}_{n+1}^m=U_n^0+\mathrm{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }{F}_n^m$$

and has initial conditions $U_n^0=u_0\left(nh\right)$ for $n=0,1,\cdots ,N$ and boundary conditions $U_0^m=u_1\left(m\tau \right),U_N^m=u_2\left(m\tau \right)$ for $m=1,2,3$. The numerical solution on row m of the grid $\mathcal{J}$, where $m>3$, which uses approximation (5) of the fractional derivative satisfies

$${\displaystyle \frac{1}{2\mathrm{\Gamma }(1-\alpha ){\tau }^{\alpha }}}\sum _{k=0}^m{\sigma }_k^{\left(\alpha \right)}{U}_n^{m-k}={\displaystyle \frac{D}{h^2}}(U_{n-1}^m-2U_n^m+U_{n+1}^m)+F_n^m,$$

$${\sigma }_0^{\left(\alpha \right)}{U}_n^{m-k}+\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{U}_n^{m-k}={\displaystyle \frac{2D\mathrm{\Gamma }(1-a){\tau }^{\alpha }}{h^2}}(U_{n-1}^m-2U_n^m+U_{n+1}^m)+2\mathrm{\Gamma }(1-a){\tau }^{\alpha }{F}_n^m.$$

Denote $\eta =2D\mathrm{\Gamma }(1-a){\tau }^{\alpha }/h^2$. The numerical solution on row m satisfies the tridiagonal system of linear equations

$$-\eta U_{n-1}^m+({\sigma }_0^{\left(\alpha \right)}+2\eta )U_n^m-\eta U_{n+1}^m=\mathrm{\Gamma }(2-\alpha ){\tau }^{\alpha }{F}_n^m-\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{U}_n^{m-k}.$$

(NS14)

Consider the fractional subdiffusion equation

$$\left\{\begin{array}{c}{D}_t^{\alpha }u(x,t)=D{u}_{xx}(x,t)+e^x(t^{1-\alpha }{E}_{1,2-\alpha }\left(t\right)-D{e}^t),\hfill \\ u(x,0)=e^x,u(0,t)=e^t,u(1,t)=e^{1+t}.\hfill \end{array}\right.$$

(64)

Equation (64) has a solution $u(x,t)=e^{x+t}$. Experimental results for the error and order of numerical solution NS14 of equation (64) are given in Table 16.

7. Conclusions

In this paper, we construct approximations of the first derivative whose weights contain powers of a parameter a, where $\left|a\right|<1$. The approximations are applied for the numerical solution of ordinary differential equations and the heat equation. The properties of the weights allow the computation of the numerical solution of ordinary differential equations with $\mathcal{O}\left(N\right)$ operations. The numerical solutions of partial differential equations are computed with $\mathcal{O}\left(MN\right)$ operations, where M is the number of nodes in the time and N is the number of nodes in space. The numerical methods are compared with the Euler method using a first-order backward approximation of the first derivative. In the examples given in the paper, we observe similar computational time and accuracy for the numerical methods, with the difference in computational time being less than a second. The numerical methods for ODEs and PDEs using approximation (3) of the first derivative can achieve an arbitrary order in $(0,2]$ using values of the parameter (37), (38), and (39). The numerical methods constructed in the paper are easy to implement and have a better performance compared to the Euler method for equation (30) with negative values of L and values of the parameter a close to one, as shown in Tables 4-6. The main goal of constructing the approximations of the first derivative (3) and (4) is to apply them to the construction of approximations of the fractional derivative, whose weights have property (1). An example of a construction of an approximation of the fractional derivative of order $2-\alpha$ which uses asymptotic formula (2) and approximation (3) of the first derivative is given in section 6. In future work, we will extend the method for constructing approximations of the second and higher-order derivatives. We will study constructions of high-order approximations of the fractional derivative that satisfy (1), and we will investigate the properties and convergence of the finite difference schemes for fractional differential equations.