First Derivative Approximations and Applications

Abstract

In this paper, we consider constructions of first derivative approximations using the generating function. The weights of the approximations contain the powers of a parameter whose modulus is less than one. The values of the initial weights are determined, and the convergence and order of the approximations are proved. The paper discusses applications of approximations of the first derivative for the numerical solution of ordinary and partial differential equations and proposes an algorithm for fast computation of the numerical solution. Proofs of the convergence and accuracy of the numerical solutions are presented and the performance of the numerical methods considered is compared with the Euler method. The main goal of constructing approximations for integer-order derivatives of this type is their application in deriving high-order approximations for fractional derivatives, whose weights have specific properties. The paper proposes the construction of an approximation for the fractional derivative and its application for numerically solving fractional differential equations. The theoretical results for the accuracy and order of the numerical methods are confirmed by the experimental results presented in the paper.

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}{\mathsf{\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 has 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,\dots ,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 }$. The Grünwald approximation is a second-order shifted approximation of the 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\mathsf{\Gamma }(2-\alpha )\right)$. Their weights satisfy

$${\sigma }_0^{\left(\alpha \right)}>0,{\sigma }_1^{\left(\alpha \right)}<{\sigma }_2^{\left(\alpha \right)}<\dots <{\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\mathsf{\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)}{\mathsf{\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 is 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 the 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^{″}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\mathsf{\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. The main contributions of the paper are the constructions of approximations (3) and (5). The examples in the paper demonstrate that these approximations can be used to construct efficient methods for the numerical solution of differential and fractional differential equations. The method for deriving an approximation (3) can be used to derive approximations of the second derivative and higher-order derivatives [26]. A question for future work is to apply these approximations to the construction of high-order approximations of the fractional derivative that have properties (1) of the weights.

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 $|a|<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^{″}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 obtain

$${\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.$$

(11)

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 the 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 (11). 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).$$

(12)

The numerical solution obtained from (12) 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.$$

(13)

Numerical solution (11) is a special case of (13) when the parameter a is zero, $a=0$. Direct computation of numerical solution (13) 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 us to compute (13) 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 (13) is given below.

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

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

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

(14)

where $x\in [0,1]$. The numerical results for the error and order of numerical solutions (11) and (13) of Equation (14) are given in

Table 1. Error and order of numerical methods (11) and (13) of Equation (14).

h	NS(11)		NS(13), a = −0.1		NS(13), a = 0.9
	Error	Order	Error	Order	Error	Order
$0.0005$	$4.2960\times {10}^{-4}$	$1.00012$	$3.5607\times {10}^{-4}$	$1.00025$	$1.5386\times {10}^{-3}$	$0.99992$
$0.00025$	$2.1479\times {10}^{-4}$	$1.00006$	$1.7802\times {10}^{-4}$	$1.00012$	$7.6933\times {10}^{-4}$	$0.99996$
$0.000125$	$1.0739\times {10}^{-4}$	$1.00003$	$8.9006\times {10}^{-5}$	$1.00006$	$3.8467\times {10}^{-4}$	$0.99998$

. The abbreviation NS for numerical solution is used. The calculation of the numerical solutions was performed using Mathematica 13 software. Although Formula (9) holds asymptotically, the experimental results in Table 1 suggest that numerical method (11) has 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 in 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}}.$$

(15)

By differentiating (15) 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}},$$

(16)

$$\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}}.$$

(17)

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.$$

(18)

$${\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 (15) and (16), we obtain

$$\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.$$

(19)

From (18) and (19), 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)$$

(20)

and holds for every value of n, where $n=2,\dots ,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 obtain

$${\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.$$

(21)

Numerical solution (21) 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 (21) 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 & u_0=y_0;u_1=u_0+h\ast F_1;u_2=u_1+h\ast F_2;\hfill \\ \hfill & A=1-a,H=h/A;S=u_1+a\ast (u_0-u_1)/A;\hfill \\ \hfill & \mathrm{for}n\mathrm{from}2\mathrm{to}N\hfill \\ \hfill & S=A\ast u_{n-1}+a\ast S;\hfill \\ \hfill & u_n=H\ast F_n+S;\hfill \end{array}$$

The algorithm uses ten operations in Step 1 and five operations for calculating each value of the sequence $u_n$ in Step 2 for $n=3,\dots ,N$. The total number of operations for calculating numerical solution (21) 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^{″}h+\mathcal{O}\left(h^2\right),$$

(22)

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].$$

(23)

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

Table 2. Error and order of numerical method (21) of Equation (23).

h	NS(21), a = 0.5		NS(21), a = −1		NS(21), a = (h − 1)/(h + 1)
	Error	Order	Error	Order	Error	Order
$0.0005$	$3.1026\times {10}^{-4}$	$0.9982$	$3.7475\times {10}^{-7}$	$1.9990$	$4.0785\times {10}^{-7}$	$1.9981$
$0.00025$	$1.5523\times {10}^{-4}$	$0.9991$	$9.3719\times {10}^{-8}$	$1.9995$	$1.0203\times {10}^{-7}$	$1.9991$
$0.000125$	$7.7640\times {10}^{-5}$	$0.9995$	$2.3433\times {10}^{-8}$	$1.9997$	$2.5515\times {10}^{-8}$	$1.9995$

.

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 (13) and (21) of Equation (10). Denote $M={max}_{x\in [0,X]}|f^{″}\left(x\right)|.$

Lemma 1. 

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,$$

(24)

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

(25)

Proof. 

From the Taylor theorem,

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

(26)

for $k=1,\dots ,n$ and ${\theta }_{n-k}\in (x_{n-k},X)$. By substituting $f_{n-k}$ with (26) in Formula (27) for ${\mathcal{A}}_2$, we obtain

$${\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}^{″}\left({\theta }_{n-k}\right)+(1-a)n^2{a}^{n-1}{f}^{″}\left({\theta }_0\right)\right).$$

From (16) and (17), we obtain

$$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)}}.$$

□

Denote $M_1={max}_{x\in [0,X]}|f^{\prime }\left(x\right)|.$ Now, we prove that (9) holds for $n>{log}_{1/a}N$.

Corollary 1. 

Let $f\in C^2[0,X]$ and $n>{log}_{1/a}N$. 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)=f_n^{\prime }+G_{n}h,$$

(27)

where the coefficient $G_n$ satisfies the estimate

$$|G_n|<{\displaystyle \frac{(1+a)M}{2(1-a)}}+{\displaystyle \frac{M_1}{X}}$$

(28)

Proof. 

From (27), it follows that $G_{n}h=E_n^{1}h-a^n{f}_n^{\prime }$. Then,

$$|G_n|\le |E_n^1|+{\displaystyle \frac{a^n|f_n^{\prime }|N}{X}}<|E_n^1|+{\displaystyle \frac{M_1}{X}}<{\displaystyle \frac{(1+a)M}{2(1-a)}}+{\displaystyle \frac{M_1}{X}}$$

because $a^n<a^{{log}_{1/a}N}=1/N$. □

In the next lemma, we prove that approximation (20) holds for all functions in $C^2[0,X]$.

Lemma 2. 

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

(29)

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}^{″}\left({\theta }_{n-k}\right)+(1-2a){(n-1)}^2{a}^{n-2}{f}^{″}\left({\theta }_1\right)+n^2{a}^{n-1}{f}^{″}\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 Lemmas 1 and 2 depends on the step size h, and approximations ${\mathcal{A}}_2$ and ${\mathcal{A}}_3$ can be defined for any interval $[l,X]$. Formula ${\mathcal{A}}_2$ is a first-order approximation of the first derivative for $n>{log}_{1/a}N$, while ${\mathcal{A}}_3$ is an approximation of the first derivative for all $n\ge 2$. The sum of the coefficients of the second derivative in the formula for $E_m^2$ is

$${\displaystyle \frac{1}{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)={\displaystyle \frac{a+1-2a^n}{2(1-a)}}.$$

Therefore, (3) holds for $n>{log}_{1/a}N$. Now, we use the bound for the errors of approximations ${\mathcal{A}}_2$ and ${\mathcal{A}}_3$ to prove the convergence of numerical methods (13) and (21). Let $e_k=y_k-u_k$ be the error of numerical solution (21) at the point $x_k$ for $k=1,\dots ,N$. First, we estimate the errors $e_1$ and $e_2$. From the Taylor Theorem, the solution of Equation (10) satisfies

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

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

Then,

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

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

In Theorems 1 and 2, we derive bounds for the errors of numerical solutions (13) and (21) of Equation (10). The theorems are proved by induction.

Theorem 1. 

The errors of numerical solution (21) satisfy

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

(30)

Proof. 

The estimate (30) for the error of numerical method (21) holds for $n=1$ and $n=2$. Suppose that (30) holds for $n=1,\dots ,m-1$. The solution of Equation (10) and the error of (21) 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,\dots ,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.$$

(31)

By applying (31) 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}|E_k|<{\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 (13).

$$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 2. 

The error of numerical solution (13) 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).$$

(32)

Proof. 

The estimate (32) for the error of numerical method (13) holds for $n=1$ and $n=2$. Suppose that (32) holds for $n=1,\dots ,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,\dots ,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}$$

(33)

By applying (33) 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}|E_k|$ and $F={max}_{x\in [0,X]}|F\left(x\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 2. 

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.$$

(34)

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

(35)

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 (34) 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 (34) 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).$$

(36)

The numerical solution (36) is calculated with $\mathcal{O}\left(N\right)$ arithmetic operations using an algorithm that is similar to the algorithms for calculating (13) and (21). From Formula (36), 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 (34) 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).$$

(37)

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.$$

(38)

Experimental results for the error and order of numerical solutions (35) and (37) of Equation (38) and positive values of L are given in

Table 3. Error and order of numerical methods (35) and (37) of Equation (38).

h	NS(35), L = 2		NS(37), L = 2, a = −0.9		NS(37), L = 20, a = 0.5
	Error	Order	Error	Order	Error	Order
$0.0005$	$2.1521\times {10}^{-4}$	$0.9998$	$1.1346\times {10}^{-5}$	$1.0023$	$9.7015\times {10}^{-5}$	$0.9990$
$0.00025$	$1.0761\times {10}^{-4}$	$0.9999$	$5.6688\times {10}^{-6}$	$1.0011$	$4.8524\times {10}^{-5}$	$0.9995$
$0.000125$	$5.3809\times {10}^{-5}$	$0.9999$	$2.8333\times {10}^{-6}$	$1.0006$	$2.4266\times {10}^{-5}$	$0.9997$

. Now, we prove the convergence of numerical method (37) when the value of L is positive. Denote by $e_n=y_n-u_n$ the error of (37) at the point $x_n$. First, we estimate the errors $e_1$ and $e_2$. The solution of Equation (34) satisfies

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

(39)

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

(40)

From the 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}^{″}\left({\theta }_1\right)\hfill \\ \hfill & ={\displaystyle \frac{h}{1-a}}(y_0^{\prime }-y_1^{\prime })+{\displaystyle \frac{h^2}{2}}{y}^{″}\left({\theta }_1\right)={\displaystyle \frac{h^2}{2}}{y}^{″}\left({\theta }_1\right)-{\displaystyle \frac{h^2}{1-a}}{y}^{″}\left({\theta }_2\right).\hfill \end{array}$$

Hence,

$$|{\epsilon }_1|\le {\displaystyle \frac{h^2}{2}}|y^{″}\left({\theta }_1\right)|+{\displaystyle \frac{h^2}{1-a}}|y^{″}\left({\theta }_2\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,|e_1|\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 (37) of Equation (38) 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 3. 

Let $L>0$. Then,

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

(41)

Proof. 

Inequality (41) is satisfied for $n=1$ and $n=2$. Assume that (41) 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}}|e_{m-1}|+{\displaystyle \frac{(|E_m^2|+a|E_{m-1}^2|h^2}{(1+LH)(1-a)}}.$$

From Lemma 2

$$|e_m|<\left(1-{\displaystyle \frac{(1-a)LH}{1+LH}}\right)|e_{m-1}|+{\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 (41) holds for all $n=1,\dots ,N$. □

We will compare the performance of numerical methods (35) and (37) and negative values of the parameter L of Equation (38). Experimental results for the error and order of numerical solutions (35) and (37) and $L=-20,-30,-50$ are given in

Table 4. Error and order of numerical method (35) of Equation (38).

h	NS(35), L = −20		NS(35), L = −30		NS(35), L = −50
	Error	Order	Error	Order	Error	Order
$1.5625\times {10}^{-5}$	$200.1160$	$1.0045$	$2.8992\times {10}^6$	$1.0101$	$8.4295\times {10}^{14}$	$1.0282$
$7.8125\times {10}^{-6}$	$99.9028$	$1.0022$	$1.4445\times {10}^6$	$1.0051$	$4.1738\times {10}^{14}$	$1.0141$
$3.9062\times {10}^{-6}$	$49.9120$	$1.0011$	$7.2099\times {10}^5$	$1.0025$	$2.0767\times {10}^{14}$	$1.0071$

,

Table 5. Error and order of numerical method (37) of Equation (38).

h	NS(37), L = −20		NS(37), L = −30		NS(37), L = −50
	a = 0.99995		a = 0.99999		a = 0.99997
	Error	Order	Error	Order	Error	Order
$1.5625\times {10}^{-5}$	$0.0333$	$0.5963$	$0.0453$	$0.1657$	$0.0179$	$0.4295$
$7.8125\times {10}^{-6}$	$0.0192$	$0.7976$	$0.0366$	$0.3072$	$0.0113$	$0.6594$
$3.9062\times {10}^{-6}$	$0.0103$	$0.8943$	$0.0255$	$0.5211$	$0.0064$	$0.8317$

and

Table 6. Error and order of numerical method (37) of Equation (38) and $a=1-h^{0.8}$.

h	NS(37), L = −20		NS(37), L = −30		NS(37), L = −50
	Error	Order	Error	Order	Error	Order
$1.5625\times {10}^{-5}$	$0.0140$	$0.1779$	$0.0092$	$0.1783$	$0.0054$	$0.1785$
$7.8125\times {10}^{-6}$	$0.0124$	$0.1805$	$0.0081$	$0.1808$	$0.0048$	$0.1811$
$3.9062\times {10}^{-6}$	$0.0109$	$0.1828$	$0.0071$	$0.1831$	$0.0042$	$0.1833$

. 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 (35) 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 (37) has an order $0.2$. Although the order of (37) in Table 5 and Table 6 is smaller than one, the errors of the method are significantly smaller than the corresponding errors of numerical solution (35) in Table 4.

Example 3. 

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.$$

(42)

The numerical solutions (43) and (45) for Equation (42), 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,$$

(43)

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

(44)

In a similar way as in Example 2, a recursive formula for calculating the numerical solution of an Equation (42) is derived from (44)

$$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.$$

(45)

The first-order liner ODE

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

(46)

has a solution $y\left(x\right)=sinx-cosx$. The experimental results of numerical solutions (43) and (45) of Equation (46) are given in

Table 7. Error and order of numerical methods (43) and (45) of Equation (50).

h	NS(43)		NS(45), a = −0.3		NS(45), a = −0.9
	Error	Order	Error	Order	Error	Order
$0.0005$	$8.4224\times {10}^{-5}$	$0.9995$	$4.5348\times {10}^{-5}$	$0.9994$	$4.4203\times {10}^{-6}$	$0.9954$
$0.00025$	$4.2119\times {10}^{-5}$	$0.9997$	$2.2679\times {10}^{-5}$	$0.9997$	$2.2136\times {10}^{-6}$	$0.9977$
$0.000125$	$2.1061\times {10}^{-5}$	$0.9999$	$1.1341\times {10}^{-5}$	$0.9999$	$1.1077\times {10}^{-6}$	$0.9989$

.

Numerical solution (45) 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 (45) has an order $1+p$. When $a=1-s{h}^p$, the coefficient $C_a$ satisfies

$$C_{a}h={\displaystyle \frac{1}{s}}(2-s{h}^p)h^{1-p}={\displaystyle \frac{2}{s}}{h}^{1-p}-h.$$

Numerical solution (45) has an order $1-p$. A summary of these results is given below.

• • Numerical method (45) 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}$$

(47)

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. Error and order of numerical method (45) of Equation (50) and $a=(s{h}^p-1)/(s{h}^p+1)$.

h	NS(45), s = 1, p = −0.3		NS(45), s = 2, p = 0.7		NS(45), s = 3, p = 1
	Error	Order	Error	Order	Error	Order
$0.0005$	$8.2264\times {10}^{-4}$	$0.6986$	$8.1004\times {10}^{-7}$	$1.6943$	$1.1272\times {10}^{-7}$	$2.0000$
$0.00025$	$5.0670\times {10}^{-4}$	$0.6991$	$2.5012\times {10}^{-7}$	$1.6954$	$2.8180\times {10}^{-8}$	$2.0000$
$0.000125$	$3.1203\times {10}^{-4}$	$0.6995$	$7.7183\times {10}^{-8}$	$1.6962$	$7.0449\times {10}^{-9}$	$2.0000$

.

• • (45) has an order $1-p$, when $s>0$ and $0<p<1$,

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

(48)

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

Table 9. Error and order of numerical method (45) of Equation (50) and $a=1-s{h}^p$.

h	NS(45), s = 1, p = 0.2		NS(45), s = 2, p = 0.5		NS(45), s = 3, p = 0.7
	Error	Order	Error	Order	Error	Order
$0.0005$	$6.8543\times {10}^{-4}$	$0.7721$	$3.6584\times {10}^{-3}$	$0.4823$	$1.1166\times {10}^{-2}$	$0.2866$
$0.00025$	$4.0015\times {10}^{-4}$	$0.7764$	$2.6093\times {10}^{-3}$	$0.4875$	$9.1303\times {10}^{-3}$	$0.2903$
$0.000125$	$2.3303\times {10}^{-4}$	$0.7800$	$1.8564\times {10}^{-3}$	$0.4912$	$7.4526\times {10}^{-3}$	$0.2929$

.

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

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

(49)

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

Table 10. Error and order of numerical method (45) of Equation (50) and $a=-1+s{h}^p$.

h	NS(45), s = 1, p = 0.25		NS(45), s = 2, p = 0.75		NS(45), s = 3, p = 1.5
	Error	Order	Error	Order	Error	Order
$0.0005$	$6.7943\times {10}^{-6}$	$1.2702$	$2.6885\times {10}^{-7}$	$1.7395$	$1.7315\times {10}^{-7}$	$1.9696$
$0.00025$	$2.8227\times {10}^{-6}$	$1.2672$	$8.0463\times {10}^{-8}$	$1.7404$	$4.3926\times {10}^{-9}$	$1.9789$
$0.000125$	$1.1748\times {10}^{-6}$	$1.2646$	$2.4064\times {10}^{-8}$	$1.7414$	$1.1094\times {10}^{-9}$	$1.9853$

.

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

Example 4. 

Consider the logistic equation

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

(50)

Equation (50) has a solution

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

The explicit schemes of Equation (50) 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,$$

(51)

$$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.$$

(52)

Numerical method (52) 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})).$$

(53)

Experimental results for the error and order of numerical solutions (51) and (53) of Equation (50) are given in

Table 11. Error and order of numerical methods (51) and (53) of Equation (50) and $L=1,y_0=2$.

h	NS(51)		NS(53), a = 0.3		NS(53), a = 0.5
	Error	Order	Error	Order	Error	Order
$0.0005$	$1.8404\times {10}^{-4}$	$1.0008$	$2.6296\times {10}^{-5}$	$1.0011$	$1.8392\times {10}^{-4}$	$0.9998$
$0.00025$	$9.1994\times {10}^{-5}$	$1.0004$	$1.3143\times {10}^{-5}$	$1.0005$	$9.1966\times {10}^{-5}$	$0.9999$
$0.000125$	$4.5991\times {10}^{-5}$	$1.0002$	$6.5704\times {10}^{-6}$	$1.0002$	$4.5984\times {10}^{-5}$	$0.9999$

.

5. Numerical Solution 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 are proved, and experimental results for the error and order of the method are presented.

Example 5. 

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.$$

(54)

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 $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}$. Now, we construct the finite difference scheme of Equation (54), 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 (54) 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}$$

(55)

The coefficient $K_n^m$ satisfies $|K_n^m|<M_4/12$. From Lemma 2 the coefficient $S_n^m$ satisfies the estimate $|S_n^m|<{\displaystyle \frac{(1+a)M_2}{2(1-a)}}$, 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|.$$

Formula (55) 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}$$

(56)

From Formulas (55) and (56), 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 & \hspace{1em}-{\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 (54) 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}$$

(57)

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 (54) on row m of the grid $\mathcal{J}$, which corresponds to (57), 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}$$

(58)

The system of linear Equation (58) is written in matrix form as

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

(59)

where ${\mathcal{U}}_m$ and ${\mathcal{F}}_m$ are vectors of dimension $N-1$ which have entries $U_n^m$ and ${\overline{F}}_n^m/(1-a)$ for $n=1,\dots ,N-1$ and ${\mathcal{I}}_m={\left(u_1\left(m\tau \right),0,\dots ,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 (59) is computed with $\mathcal{O}\left(MN\right)$ operations because $I+\eta \Delta$ is a tridiagonal M-matrix [37,38].

A Crank–Nicolson scheme for the heat equation has second-order accuracy [39,40]. The convergence of finite difference schemes for the heat equation are studied in [41,42]. Now, we derive a bound for the error of finite difference scheme (59). 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,\dots ,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 3, we estimate the errors of difference scheme (59) on the first row of the grid $\mathcal{J}$.

Lemma 3. 

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}$$

(60)

Proof. 

The solution of Equation (54) 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|C_4|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}$$

(61)

for $n=1,\dots ,N-1$ and $e_0^1=e_N^1=0$. The errors of numerical method (59) 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,$$

(62)

where ${\mathcal{W}}_1$ is the vector with elements ${\epsilon }_n^1$. From (61), 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,$$

(63)

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,\dots ,N-1$$

and eigenvectors ${\mathcal{V}}_k$ for $k=1,\dots ,N-1$, which have components $v_k^i=sin(k\pi i/N)$, where $i=1,\dots ,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.$$

(64)

The matrix P has eigenvalues $(I+a\eta {\lambda }_k)/(I+\eta {\lambda }_k)$ for $k=1,\dots ,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 4, we prove that finite difference scheme (59) is unconditionally stable and has an accuracy $O(\tau +h^2)$.

Theorem 4. 

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,\dots M$.

Proof. 

By applying (64) 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 3, 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).$$

□

Consider the following heat equation

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

(65)

where $(x,t)\in [0,1]\times [0,1]$. Equation (65) has a solution $u(x,t)=e^{-x-t}$. The graphs of numerical solution (59) and its error for $h=\tau =0.01$ are shown in Figure 1. Numerical method (59) has an accuracy $\mathcal{O}(\tau +h^2)$ for $a=-0.5$. When the parameter $a=(\tau -1)/(\tau +1)$, numerical method (59) has a second-order accuracy $\mathcal{O}({\tau }^2+h^2)$. Experimental results for the error and order of the numerical solution (59) of Equation (65) with $a=-0.5$ and $a=(\tau -1)/(\tau +1)$ are given in

Table 12. Error and order of numerical solution (59) and $N=1000$.

$\mathit{\tau }$	$\mathit{\alpha }=-0.5$		$\mathit{a}=(\mathit{\tau }-1)/(\mathit{\tau }+1)$
	Error	Order	Error	Order
$1/8$	$5.0573\times {10}^{-4}$	$-$	18,947 $\times {10}^{-4}$	$-$
$1/16$	$2.3493\times {10}^{-4}$	$1.1061$	39,067 $\times {10}^{-5}$	$2.2779$
$1/32$	$1.1768\times {10}^{-4}$	$0.9974$	$9.4822\times {10}^{-6}$	$2.0427$
$1/64$	$5.9357\times {10}^{-5}$	$0.9873$	$2.3558\times {10}^{-6}$	$2.0089$
$1/128$	$2.9850\times {10}^{-5}$	$0.9917$	$5.9294\times {10}^{-7}$	$1.9902$

and

Table 13. Error and order of numerical solution (59) and $M=1000$.

h	$\mathit{\alpha }=-0.5$		$\mathit{a}=(\mathit{\tau }-1)/(\mathit{\tau }+1)$
	Error	Order	Error	Order
$1/5$	$2.3007\times {10}^{-4}$	$-$	$2.2645\times {10}^{-4}$	$-$
$1/10$	$5.8200\times {10}^{-5}$	$1.9829$	$5.7235\times {10}^{-5}$	$1.9842$
$1/20$	$1.4924\times {10}^{-5}$	$1.9633$	$1.4425\times {10}^{-5}$	$1.9883$
$1/40$	$3.8741\times {10}^{-6}$	$1.9457$	$3.6144\times {10}^{-6}$	$1.9967$
$1/80$	$1.0286\times {10}^{-6}$	$1.9131$	$9.1106\times {10}^{-7}$	$1.9881$

.

The order of the numerical method (59) in time is obtained by fixing $N=1000$ and varying the value of M using the formula ${log}_2|E_{\tau }/E_{2\tau }|$. Similarly, the order in space is computed as ${log}_2|E_h/E_{2h}|$ with the parameter $M=1000$ kept fixed.

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 the Caputo fractional derivative has the form

$${\displaystyle \frac{1}{\mathsf{\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),$$

(66)

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 the Caputo derivative using the asymptotic formula of the 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\mathsf{\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),$$

(67)

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\mathsf{\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),$$

(68)

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, $|b|<1$. Therefore, ${\sigma }_0^{\left(\alpha \right)}>0$ and ${\sigma }_k^{\left(\alpha \right)}<0$ for $k>1$. In Lemma 4, we show that the coefficient ${\sigma }_1^{\left(\alpha \right)}$ of approximation (68) is negative when the parameter b has a value $b=\alpha /2$.

Lemma 4. 

Let $b=a/2$. Then,

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

(69)

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|\zeta \left(\alpha \right)|.$$

Inequality (69) holds when

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

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

Taking the logarithm of both sides, we obtain

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

Let $f\left(\alpha \right)=ln|\zeta (\alpha )|-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,\dots .$$

Hence,

$$|\zeta \left(\alpha \right)|<{\displaystyle \frac{1}{1-\alpha }},|{\zeta }^{\prime }\left(\alpha \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^{″}\left(\alpha \right)={\displaystyle \frac{2}{{(1-\alpha )}^3}}-0.02-{\displaystyle \frac{4}{{(2-\alpha )}^3}},$$

$$g^{‴}\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^{″},g^{‴}$ at zero are positive:

$$g\left(0\right)=ln2-0.09>0,g^{\prime }\left(0\right)=0.6,g^{″}\left(0\right)=1.48,g^{‴}\left(0\right)=5.25.$$

The third derivative $g^{‴}$ 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^{″}$ 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 (68) 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 following examples, we will consider an application of approximation (68) for the construction of numerical schemes for the two-term ordinary fractional differential equation and the fractional subdiffusion equation.

Example 6. 

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.$$

(70)

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

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

(71)

The numerical solution of Equation (70), which uses approximations (67) and (68) of the fractional derivative $y_n^{\left(\alpha \right)}$, satisfies

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

(72)

and has initial conditions

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

(73)

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

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

(74)

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.$$

(75)

Equation (75) 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}{\mathsf{\Gamma }(k+2-\alpha )}}.$$

The experimental results of numerical solutions (71) and (72) of Equation (75) on the interval $[0,1]$ are given in

Table 14. Error and order of numerical solutions (71) and (72) of Equation (75).

h	NS(71), $\mathit{\alpha }=0.25,\mathit{L}=1$		NS(72), $\mathit{\alpha }=0.5,\mathit{L}=10$		$\mathit{\alpha }=0.75,\mathit{L}=5,\mathit{b}=-0.5$
	Error	Order	Error	Order	Error	Order
$0.0005$	$2.9218\times {10}^{-7}$	$1.7225$	$7.9883\times {10}^{-7}$	$1.4919$	$2.3232\times {10}^{-6}$	$1.2443$
$0.00025$	$8.8244\times {10}^{-7}$	$1.7273$	$2.8355\times {10}^{-7}$	$1.4943$	$9.7923\times {10}^{-7}$	$1.2463$
$0.000125$	$2.6579\times {10}^{-8}$	$1.7312$	$1.0053\times {10}^{-7}$	$1.4959$	$4.1236\times {10}^{-7}$	$1.2477$

and

Table 15. Error and order of numerical solution (72) of Equation (75) and $b=\alpha /2$.

h	NS(72),$\mathit{\alpha }=0.25,\mathit{L}=1$		NS(72),$\mathit{\alpha }=0.5,\mathit{L}=10$		NS(72),$\mathit{\alpha }=0.75,\mathit{L}=5$
	Error	Order	Error	Order	Error	Order
$0.0005$	$6.4339\times {10}^{-7}$	$1.7248$	$1.5529\times {10}^{-6}$	$1.4953$	$3.1981\times {10}^{-5}$	$1.2487$
$0.00025$	$1.9405\times {10}^{-7}$	$1.7292$	$5.5026\times {10}^{-7}$	$1.4968$	$1.3453\times {10}^{-5}$	$1.2493$
$0.000125$	$5.8383\times {10}^{-8}$	$1.7328$	$1.9484\times {10}^{-7}$	$1.4978$	$5.6578\times {10}^{-6}$	$1.2496$

. The second column of Table 14 pertains to (72) using approximation (67), while the results in the third column relate to (72) using approximation (68) 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 (67) and (68) for $n=3,\dots ,N$ are computed recursively in $\mathcal{O}\left(N\right)$ arithmetic operations. The (71) and (72) methods have similar computational time and accuracy. The approximation (67) is obtained from 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 the 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 7. 

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.$$

(76)

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 let $\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 (74) 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}{\mathsf{\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\mathsf{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }}{h^2}}(U_{n-1}^m-2U_n^m+U_{n+1}^m)+\mathsf{\Gamma }(2-\alpha ){\left(m\tau \right)}^{\alpha }{F}_n^m.$$

Denote ${\eta }_m=D\mathsf{\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+\mathsf{\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,\dots ,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\mathsf{\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\mathsf{\Gamma }(1-a){\tau }^{\alpha }}{h^2}}(U_{n-1}^m-2U_n^m+U_{n+1}^m)+2\mathsf{\Gamma }(1-a){\tau }^{\alpha }{F}_n^m.$$

Denote $\eta =2D\mathsf{\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=\mathsf{\Gamma }(2-\alpha ){\tau }^{\alpha }{F}_n^m-\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{U}_n^{m-k}.$$

(77)

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.$$

(78)

Equation (78) has a solution $u(x,t)=e^{x+t}$. The graphs of numerical solution (77) of Equation (78) and its error for $h=\tau =0.01$ are shown in Figure 2. Experimental results for the error and order of numerical solution (77) of fractional subdiffusion Equation (78) and $\alpha =0.5, 0.75$ are given in

Table 16. Error and order of numerical solution (77) and $N=1000$.

$\mathit{\tau }$	$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.75$
	Error	Order	Error	Order
$1/8$	$5.6355\times {10}^{-3}$	$-$	$1.6687\times {10}^{-2}$	$-$
$1/16$	$2.1509\times {10}^{-3}$	$1.3896$	$7.5333\times {10}^{-3}$	$1.1473$
$1/32$	$7.9545\times {10}^{-4}$	$1.4351$	$3.2885\times {10}^{-3}$	$1.1958$
$1/64$	$2.8898\times {10}^{-4}$	$1.4607$	$1.4105\times {10}^{-3}$	$1.2213$
$1/128$	$1.0393\times {10}^{-4}$	$1.4753$	$5.9934\times {10}^{-4}$	$1.2347$

and

Table 17. Error and order of numerical solution (77) and $N=1000$.

h	$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.75$
	Error	Order	Error	Order
$1/5$	$1.6057\times {10}^{-6}$	$-$	$1.4085\times {10}^{-6}$	$-$
$1/10$	$3.5567\times {10}^{-4}$	$1.9845$	$3.5532\times {10}^{-4}$	$1.9555$
$1/20$	$8.9727\times {10}^{-5}$	$1.9869$	$8.9352\times {10}^{-5}$	$1.9915$
$1/40$	$2.2870\times {10}^{-5}$	$1.9721$	$2.2487\times {10}^{-5}$	$1.9903$
$1/80$	$6.0154\times {10}^{-6}$	$1.9267$	$5.6314\times {10}^{-6}$	$1.9975$

.

7. Conclusions

In this paper, we construct approximations of the first derivative whose weights contain powers of a parameter a, where $|a|<1$. The approximations are applied to 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 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 time difference 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 (47), (48), and (49). The numerical methods constructed in the paper are easy to implement and perform better compared to the Euler method for Equation (38) with negative values of L and values of the parameter a close to one, as shown in Table 4, Table 5 and Table 6. The main goal of constructing the first derivative approximations (3) and (4) is to use them in constructing fractional derivative approximations, 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. The method for constructing the approximation (5) can be generalized to develop high-order approximations of the fractional derivative. In future work, we plan to extend the results presented in this paper. We will present proofs of the convergence of the numerical methods discussed in the paper and apply the method to construct approximations for the second-order and higher-order derivatives. We aim to study the construction of high-order approximations of the fractional derivative with weights that satisfy (1), as well as to investigate the properties and convergence of finite difference schemes for fractional differential equations.