Open AccessArticle

First Derivative Approximations and Applications

Yuri Dimitrov

^1,2,*

Slavi Georgiev

^2,3

and

Venelin Todorov

Department of Mathematics, Physics and Informatics, University of Forestry, 1756 Sofia, Bulgaria

Department of Informational Modeling, Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria

Department of Applied Mathematics and Statistics, University of Ruse, 7017 Ruse, Bulgaria

Author to whom correspondence should be addressed.

Fractal Fract. 2024, 8(10), 608; https://doi.org/10.3390/fractalfract8100608

Submission received: 6 September 2024 / Revised: 8 October 2024 / Accepted: 16 October 2024 / Published: 18 October 2024

(This article belongs to the Special Issue Advances in Fractional Modeling and Computation)

Download

Browse Figures

Versions Notes

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.

Keywords:

fractional derivative; approximation; numerical solution; convergence

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

α

, where

0 < α < 1

, is defined as

f^{(α)} (x) = D^{α} f (x) = \frac{1}{Γ (1 - α)} \int_{l}^{x} \frac{f^{'} (ξ)}{{(x - ξ)}^{α}} d ξ .

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 + n h

and

f_{n} = f (x_{n})

. The approximations of fractional derivative

f_{n}^{(α)}

have the form

h^{- α} \sum_{k = 0}^{n} σ_{k}^{(α)} f_{n - k}

, where

σ_{k}^{(α)}

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 (x) = \sum_{k = 0}^{\infty} σ_{k}^{(α)} 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

α

has a first-order accuracy and a generating function

{(1 - x)}^{α}

. The Grünwald approximation is a second-order shifted approximation of the fractional derivative with a shift parameter

α / 2

. L1 approximation has an order

2 - α

and a generating function

{(x - 1)}^{2} {Li}_{α - 1} (x) / (x Γ (2 - α))

. Their weights satisfy

σ_{0}^{(α)} > 0, σ_{1}^{(α)} < σ_{2}^{(α)} < \dots < σ_{m - 1}^{(α)} < 0, \sum_{k = 0}^{m} σ_{k}^{(α)} = 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 - α

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^{- α} - n^{1 - α} / (1 - α) - ζ (α)

. In [23], we obtain an approximation of the fractional derivative and its asymptotic formula

\frac{1}{2 Γ (1 - α) h^{α}} \sum_{k = 0}^{m} ω_{k}^{(α)} f_{n - k} = f_{m}^{(α)} + \frac{ζ (α)}{Γ (1 - α)} f_{n}^{'} h^{1 - α} + O (h^{2 - α}),

(2)

where

ω_{0}^{(α)} = 1, ω_{1}^{(α)} = 2^{- α}, ω_{k}^{(α)} = {(k + 1)}^{- α} - {(k - 1)}^{- α},

σ_{m - 1}^{(α)} = - {(n - 1)}^{- α} + 2 s_{n}, σ_{n}^{(α)} = - {(n - 2)}^{- α} - 2 s_{n} .

Approximation (2) has an order

1 - α

and a generating function

(1 - x^{2}) {Li}_{- α} (x) / (2 x)

. 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

\frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} f_{n - k} - \frac{a^{n - 2} (1 - 2 a)}{1 - a} f_{1} - \frac{a^{n - 1}}{1 - a} f_{0}) = f_{n}^{'} - \frac{C_{a}}{2} f_{n}^{″} h + O (h^{2}),

(3)

where the coefficient

C_{a}

has a value

C_{a} = (1 + a) / (1 - a)

. Approximation (3) has a generating function

G (x) = (1 - a) (1 - x) / (1 - a x)

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:

\frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} f_{n - k} - a^{n - 1} f_{0}) = (1 - a^{n}) f_{n}^{'} + O (h) .

(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

O (N)

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

α / 2

to the first derivative

f_{n}^{'}

in the asymptotic Formula (2), we obtain the following approximation of a fractional derivative of order

2 - α

\frac{1}{2 Γ (1 - α) h^{α}} \sum_{k = 0}^{n} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + O (h^{2 - α}),

(5)

where

σ_{0}^{(α)} = 1 - (2 - α) ζ (α), σ_{1}^{(α)} = 2^{- α} + 2 {(1 - α / 2)}^{2} ζ (α),

σ_{k}^{(α)} = {(k + 1)}^{- α} - {(k - 1)}^{- α} + 2 {(1 - α / 2)}^{2} {(α / 2)}^{k - 1} ζ (α), (2 \leq k \leq n - 2),

σ_{n - 1}^{(α)} = - {(n - 1)}^{- α} + 2 s_{n} + 2 (1 - α) {(α / 2)}^{n - 2} ζ (α),

σ_{n}^{(α)} = - {(n - 2)}^{- α} - 2 s_{n} + 2 {(α / 2)}^{n - 1} ζ (α) .

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 (x) = \frac{(1 - a) (1 - x)}{1 - a x},

where

| a | < 1

. The generating function

G (x)

has properties

G (1) = 0, G^{'} (1) = - 1

. Denote

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

The functions

G (x)

and

H (x)

have Maclaurin series

G (x) = (1 - a) - {(1 - a)}^{2} \sum_{k = 1}^{\infty} a^{k - 1} x^{k},

(6)

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

(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 (x_{n}) = f (n h)

. From (6) and (7), we obtain

A_{1} [f_{n}] = \frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} f_{n - k}) = f_{n}^{'} - \frac{1 + a}{2 (1 - a)} f_{n}^{″} h + O (h^{2}),

(8)

where the function

f \in C^{2} [0, X]

and satisfies the condition

f (0) = 0

. Now, we extend asymptotic Formula (8) to all functions in the class

C^{2} [0, X]

. Let

A_{2} [f_{n}] = \frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} f_{n - k} + c_{0} f_{0}) = f_{n}^{'} + O (h),

where

A_{2}

satisfies the condition

A_{2} [1] = 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

A_{2} [f_{n}] = \frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} f_{n - k} - a^{n - 1} f_{0}) = f_{n}^{'} + O (h) .

(9)

Asymptotic Formula (9) holds for all functions

f \in C^{2} [0, X]

. When the parameter a is zero,

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^{'} (x) = F (x), y (0) = 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

\frac{y_{n} - y_{n - 1}}{h} = F_{n} + O (h) .

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

2 N

. 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}^{'}

with

A_{2}

, we find

\frac{1 - a}{h} (y_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} y_{n - k} - a^{n - 1} y_{0}) = F_{n} + O (h) .

(12)

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

\frac{1 - a}{h} (u_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} u_{n - k} - a^{n - 1} u_{0}) = F_{n} .

The sequence

{u_{n}}_{n = 0}^{N}

has initial conditions

u_{0} = y_{0}, u_{1} = u_{0} + h F_{1}

and satisfies

u_{n} = \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

O (N^{2})

operations. The weights of

A_{2}

consist of powers of the parameter a, which allows us to compute (13) with

O (N)

operations. The sequence

u_{n}

satisfies

u_{n} = \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 ((1 - a) \sum_{k = 2}^{n - 1} a^{k - 2} u_{n - k} + a^{n - 2} u_{0}),

S_{n} = (1 - a) u_{n - 1} + a S_{n - 1} .

Both sequences

S_{n}

and

u_{n}

are computed with

O (N)

operations. The pseudocode for calculating the numerical solution (13) is given below.

\begin{matrix} A = 1 - a; H = h / A; S = u_{0} = y_{0}; u_{1} = u_{0} + h * F_{1}; \\ for n from 2 to N \\ S = A * u_{n - 1} + a * S; \\ u_{n} = H * F_{n} + S; \end{matrix}

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

5 N - 1

. Consider the ordinary differential equation

y^{'} (x) = e^{x}, y (0) = 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. 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^{n h}

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} = \frac{1 - a^{n}}{1 - a} .

(15)

By differentiating (15) twice, we find

\sum_{k = 1}^{n - 1} k a^{k - 1} = \frac{1 + (n - 1) a^{n} - n a^{n - 1}}{{(1 - a)}^{2}},

(16)

\sum_{k = 1}^{n - 1} k^{2} a^{k - 1} = \frac{1 + a - n^{2} a^{n - 1} + (2 n^{2} - 2 n - 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

A_{3} [f_{n}] = \frac{1 - a}{h} (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}) = f_{n}^{'} + O (h)

which satisfies the conditions

A_{3} [1] = 0

and

A_{3} [x] = 1

. The values of the coefficients

c_{0}

and

c_{1}

are determined from the two conditions.

A_{3} [1] = \frac{1 - a}{h} (1 - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} - a^{n - 2} c_{1} - a^{n - 1} c_{0}) = 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)

A_{3} [x] = (1 - a) (n - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} (n - k) - a^{n - 2} c_{1}) = 1 .

Using Formulas (15) and (16), we obtain

\sum_{k = 1}^{n - 2} (n - k) a^{k - 1} = - \frac{1 + a^{n - 2} - 2 a^{n - 1} - n + n a}{{(1 - a)}^{2}} .

Then,

(1 - a) n + 1 + a^{n - 2} - 2 a^{n - 1} - n + n a - (1 - a) a^{n - 2} c_{1} = 1,

a^{n - 2} - 2 a^{n - 1} - (1 - a) a^{n - 2} c_{1} = 1, (1 - a) c_{1} = 1 - 2 a .

(19)

From (18) and (19), we find

c_{1} = \frac{1 - 2 a}{1 - a}, c_{0} = \frac{1}{1 - a} .

Approximation

A_{3}

has the form

A_{3} [f_{n}] = \frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} f_{n - k} - \frac{a^{n - 2} (1 - 2 a)}{1 - a} f_{1} - \frac{a^{n - 1}}{1 - a} f_{0}) = f_{n}^{'} + O (h)

(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

A_{3}

for the first derivative. By replacing the first derivative

y_{n}^{'}

with

A_{3}

, we obtain

\frac{1 - a}{h} (u_{n} - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} u_{n - k} - \frac{a^{n - 2} (1 - 2 a)}{1 - a} u_{1} - \frac{a^{n - 1}}{1 - a} u_{0}) = F_{n},

u_{n} = \frac{h F_{n}}{1 - a} + (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} u_{n - k} + \frac{a^{n - 2} (1 - 2 a)}{1 - a} u_{1} + \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} = \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 ((1 - a) \sum_{k = 2}^{n - 1} a^{k - 2} u_{n - k} + a^{n - 2} u_{0}),

S_{n} = (1 - a) u_{n - 1} + a S_{n - 1} .

The sequences

S_{n}

and

u_{n}

are calculated in

O (N)

time with the following pseudocode.

\begin{matrix} u_{0} = y_{0}; u_{1} = u_{0} + h * F_{1}; u_{2} = u_{1} + h * F_{2}; \\ A = 1 - a, H = h / A; S = u_{1} + a * (u_{0} - u_{1}) / A; \\ for n from 2 to N \\ S = A * u_{n - 1} + a * S; \\ u_{n} = H * F_{n} + S; \end{matrix}

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

5 N

From (8), approximation

A_{3}

has a second-order asymptotic expansion formula

\frac{1 - a}{h} (y_{n} - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} y_{n - k} - \frac{a^{n - 2} (1 - 2 a)}{1 - a} y_{1} - \frac{a^{n - 1}}{1 - a} y_{0}) = f_{n}^{'} - \frac{C_{a}}{2} f_{n}^{″} h + O (h^{2}),

(22)

where

C_{a} = (1 + a) / (1 - a)

. Approximation

A_{3}

has second-order accuracy, when the coefficient

C_{a}

is equal to zero,

C_{a} = 0

C_{a} = h

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

Approximation

A_{3}

is second-order accurate when

a = - 1

and

a = (h - 1) / (h + 1)

. Consider the following ordinary differential equation:

y^{'} (x) = sin x + cos x, y (0) = - 1, x \in [0, 1] .

(23)

Equation (23) has a solution

y (x) = sin x - cos x

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

3. Error Estimates and Convergence

In this section, we derive estimates for the errors of approximations

A_{2}

and

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^{″} (x) | .

Lemma 1.

Let

f \in C^{2} [0, X]

. Then,

A_{2} [f_{n}] = \frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} f_{n - k} - a^{n - 1} f_{0}) = (1 - a^{n}) f_{n}^{'} + E_{n}^{1} h,

(24)

where the coefficient

E_{n}^{1}

satisfies the estimate

| E_{n}^{1} | < \frac{(1 + a) M}{2 (1 - a)}, (2 \leq n \leq N) .

(25)

Proof.

From the Taylor theorem,

\begin{matrix} f_{n - k} = f_{n} - k h f_{n}^{'} + k^{2} \frac{h^{2}}{2} f^{″} (θ_{n - k}) \end{matrix}

(26)

for

k = 1, \dots, n

and

θ_{n - k} \in (x_{n - k}, X)

. By substituting

f_{n - k}

with (26) in Formula (27) for

A_{2}

, we obtain

A_{2} [f] = K_{0} f_{n} + K_{1} f_{n}^{'} + 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} = - \frac{1}{2} ({(1 - a)}^{2} \sum_{k = 1}^{n - 1} k^{2} a^{k - 1} f^{″} (θ_{n - k}) + (1 - a) n^{2} a^{n - 1} f^{″} (θ_{0})) .

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} | \leq \frac{M}{2} ({(1 - a)}^{2} \sum_{k = 1}^{n - 1} k^{2} a^{k - 1} + (1 - a) n^{2} a^{n - 1}),

| E_{n} | \leq \frac{M}{2} (\frac{1 + a - n^{2} a^{n - 1} + (2 n^{2} - 2 n - 1) a^{n} - {(n - 1)}^{2} a^{n + 1}}{1 - a} + (1 - a) n^{2} a^{n - 1}),

| E_{n} | \leq \frac{M}{2 (1 - a)} (1 + a - (2 n + 1) a^{n} + (2 n - 1) a^{n + 1}) < \frac{M (1 + a)}{2 (1 - a)} .

□

Denote

M_{1} = {max}_{x \in [0, X]} | f^{'} (x) | .

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,

A_{2} [f_{n}] = \frac{1 - a}{h} (f_{n} - (1 - a) \sum_{k = 1}^{n - 1} a^{k - 1} f_{n - k} - a^{n - 1} f_{0}) = f_{n}^{'} + G_{n} h,

(27)

where the coefficient

G_{n}

satisfies the estimate

| G_{n} | < \frac{(1 + a) M}{2 (1 - a)} + \frac{M_{1}}{X}

(28)

Proof.

From (27), it follows that

G_{n} h = E_{n}^{1} h - a^{n} f_{n}^{'}

. Then,

| G_{n} | \leq | E_{n}^{1} | + \frac{a^{n} | f_{n}^{'} | N}{X} < | E_{n}^{1} | + \frac{M_{1}}{X} < \frac{(1 + a) M}{2 (1 - a)} + \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,

A_{3} [f_{n}] = \frac{1}{h} ((1 - a) f_{n} - {(1 - a)}^{2} \sum_{k = 1}^{n - 2} a^{k - 1} f_{n - k} - a^{n - 2} (1 - 2 a) f_{1} - a^{n - 1} f_{0}) = f_{n}^{'} + E_{n}^{2} h,

where the coefficient

E_{n}^{2}

satisfies the estimate

| E_{n}^{2} | < \frac{(1 + a) M}{2 (1 - a)}, (2 \leq n \leq N) .

(29)

Proof.

Approximation

A_{3}

is written in the form

A_{3} [f_{n}] = K_{0} f_{n} + K_{1} f_{n}^{'} + E_{n}^{2} h,

where the coefficients

K_{0}, K_{1}

and

E_{n}^{2}

satisfy

\begin{matrix} K_{0} & = 1 - a - {(1 - a)}^{2} \sum_{k = 1}^{n - 2} a^{k - 1} - (1 - 2 a) a^{n - 2} - a^{n - 1} \\ = 1 - a - (1 - a) (1 - a^{n - 2}) - (1 - 2 a) a^{n - 2} - a^{n - 1} = 0, \end{matrix}

\begin{matrix} K_{1} & = {(1 - a)}^{2} \sum_{k = 1}^{n - 2} k a^{k - 1} + (1 - 2 a) (n - 1) a^{n - 2} + n a^{n - 1} \\ = 1 + (n - 2) a^{n - 1} - (n - 1) a^{n - 2} + (1 - 2 a) (n - 1) a^{n - 2} + n a^{n - 1} = 1, \end{matrix}

E_{n}^{2} = - \frac{1}{2} ({(1 - a)}^{2} \sum_{k = 1}^{n - 2} k^{2} a^{k - 1} f^{″} (θ_{n - k}) + (1 - 2 a) {(n - 1)}^{2} a^{n - 2} f^{″} (θ_{1}) + n^{2} a^{n - 1} f^{″} (θ_{0})),

\begin{matrix} | E_{n}^{2} | & \leq \frac{M}{2} ({(1 - a)}^{2} \sum_{k = 1}^{n - 2} k^{2} a^{k - 1} + (1 - 2 a) {(n - 1)}^{2} a^{n - 2} + n^{2} a^{n - 1}) \\ \leq \frac{M}{2} \frac{1 + a - 2 a^{n}}{1 - a} < \frac{M (1 + a)}{2 (1 - a)} . \end{matrix}

□

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

A_{2}

and

A_{3}

can be defined for any interval

[l, X]

. Formula

A_{2}

is a first-order approximation of the first derivative for

n > {log}_{1 / a} N

, while

A_{3}

is an approximation of the first derivative for all

n \geq 2

. The sum of the coefficients of the second derivative in the formula for

E_{m}^{2}

\frac{1}{2} ({(1 - a)}^{2} \sum_{k = 1}^{n - 2} k^{2} a^{k - 1} + (1 - 2 a) {(n - 1)}^{2} a^{n - 2} + n^{2} a^{n - 1}) = \frac{a + 1 - 2 a^{n}}{2 (1 - a)} .

Therefore, (3) holds for

n > {log}_{1 / a} N

. Now, we use the bound for the errors of approximations

A_{2}

and

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} - \frac{h^{2}}{2} y^{″} (θ_{1}),

y_{2} = y_{1} + h F_{2} - \frac{h^{2}}{2} y^{″} (θ_{2}) .

Then,

e_{1} = - \frac{h^{2}}{2} y^{″} (θ_{1}), | e_{1} | \leq \frac{h^{2}}{2} | y^{″} (θ_{1}) | \leq \frac{M h^{2}}{2},

e_{2} = e_{1} - \frac{h^{2}}{2} y^{″} (θ_{2}), | e_{2} | \leq | e_{1} | + \frac{h^{2}}{2} | y^{″} (θ_{2}) | \leq 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} | < \frac{M X h}{{(1 - a)}^{2}}, (1 \leq n \leq 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} = \frac{h F_{m}}{1 - a} + (1 - a) \sum_{k = 1}^{m - 2} a^{k - 1} y_{m - k} + \frac{a^{m - 2} (1 - 2 a)}{1 - a} y_{1} + \frac{a^{m - 1}}{1 - a} y_{0} + \frac{E_{m}^{2} h^{2}}{1 - a},

e_{m} = (1 - a) \sum_{k = 1}^{m - 2} a^{k - 1} e_{m - k} + \frac{a^{m - 2} (1 - 2 a)}{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}

e_{m} = (1 - a) e_{m - 1} + a (1 - a) \sum_{k = 2}^{m - 2} a^{k - 2} e_{m - k} + \frac{a^{m - 2} (1 - 2 a)}{1 - a} e_{1} + E_{m} h^{2},

e_{m} = (1 - a) e_{m - 1} + a (e_{m - 1} - \frac{a^{m - 3} (1 - 2 a)}{1 - a} e_{1} - E_{m - 1} h^{2}) + \frac{a^{m - 2} (1 - 2 a)}{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 = max_{1 \leq k \leq N} | E_{k} | < \frac{M (1 + a)}{2 {(1 - a)}^{2}} .

The error

e_{m}

satisfies the estimate.

\begin{matrix} | e_{m} | & \leq | e_{2} | + (m - 2 + a) E h^{2} < \frac{M h^{2}}{{(1 - a)}^{2}} + \frac{(1 + a) M h^{2}}{2 {(1 - a)}^{2}} (m - 2 + a) \\ < \frac{M h^{2}}{{(1 - a)}^{2}} + \frac{M h^{2}}{{(1 - a)}^{2}} (m - 2 + a) < \frac{M m h^{2}}{{(1 - a)}^{2}} \leq \frac{M X h}{{(1 - a)}^{2}} . \end{matrix}

□

Now, we estimate the error

e_{2}

of numerical solution (13).

y_{2} = \frac{f_{2} h}{1 - a} + (1 - a) y_{1} + a y_{0} + \frac{E_{2}^{1} h^{2}}{1 - a},

u_{2} = \frac{f_{2} h}{1 - a} + (1 - a) u_{1} + a u_{0} .

The error

e_{2}

satisfies

e_{2} = (1 - a) e_{1} + \frac{E_{2}^{1} h^{2}}{1 - a},

| e_{2} | \leq (1 - a) | e_{1} | + \frac{| E_{2}^{1} | h^{2}}{1 - a} \leq \frac{M (1 - a) h^{2}}{2} + \frac{M (1 + a) h^{2}}{2 {(1 - a)}^{2}},

| e_{2} | \leq \frac{M h^{2}}{2 {(1 - a)}^{2}} (2 - 2 a + 3 a^{2} - a^{3}) < \frac{M h^{2}}{{(1 - a)}^{2}} .

Theorem 2.

The error of numerical solution (13) satisfies

| e_{n} | < (\frac{2 a^{2} F}{1 - a} + \frac{M (1 + a) X}{2 {(1 - a)}^{2}}) h, (1 \leq n \leq 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} = \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} + \frac{E_{m}^{1} h^{2}}{1 - a} .

Then,

e_{m} = - \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}

e_{m} = - \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} = - \frac{a^{n} F_{m} h}{1 - a} + (1 - a) e_{m - 1} + a (e_{m - 1} + \frac{a^{m - 1} h F_{m - 1}}{1 - a} - E_{m - 2} h^{2}) + E_{m} h^{2},

\begin{matrix} e_{m} = e_{m - 1} + \frac{a^{n} (F_{m - 1} - F_{m}) h}{1 - a} + (E_{m} - a E_{m - 1}) h^{2} . \end{matrix}

(33)

By applying (33) successively

m - 2

times, we obtain

e_{m} = e_{1} + (a^{2} F_{1} - a^{m} F_{n} + (a - 1) \sum_{k = 2}^{m - 1} a^{k} F_{k}) \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 \leq k \leq N} | E_{k} |

and

F = {max}_{x \in [0, X]} | F (x) |

. Then,

\begin{matrix} | e_{m} | & \leq | e_{1} | + (a^{2} + a^{m} + (a - 1) \sum_{k = 2}^{m - 1} a^{k}) \frac{F h}{1 - a} + E h^{2} (1 + a) + (1 - a) (m - 3) \\ < \frac{M h^{2}}{2} + (a^{2} + a^{m} + a^{2} (1 - a^{m - 2}) \frac{F h}{1 - a} + \frac{(1 + a) M h^{2}}{2 {(1 - a)}^{2}} (m - 2 + a), \end{matrix}

| e_{m} | < (\frac{2 a^{2} F}{1 - a} + \frac{M (1 + a) X}{2 {(1 - a)}^{2}}) h .

□

4. Numerical Solutions of First-Order ODEs

In this section, we construct the numerical methods for first-order ODEs which use approximation

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

A_{3}

Example 2.

Consider the first-order linear ODE with constant coefficients

y^{'} (x) + L y (x) = F (x), y (0) = y_{0} .

(34)

By approximating

y_{n}^{'}

with first-order backward difference approximation, we obtain

\frac{y_{n} - y_{n - 1}}{h} + L y_{n} = F_{n} + O (h) .

The numerical solution is computed as

u_{0} = y_{0}

and

u_{n} = \frac{1}{1 + L h} (h F_{n} + u_{n - 1}) .

(35)

By approximating

y_{n}^{'}

with approximation

A_{3}

, we obtain

\frac{1 - a}{h} (y_{n} - (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} y_{n - k} - \frac{1 - 2 a}{1 - a} a^{n - 2} y_{1} - \frac{a^{n - 1}}{1 - a} y_{0}) + L y_{n} = y_{n}^{'} + E_{n}^{2} h,

(1 + \frac{L h}{1 - a}) y_{n} = \frac{h F_{n}}{1 - a} + (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} y_{n - k} + \frac{1 - 2 a}{1 - a} a^{n - 2} y_{1} + \frac{a^{n - 1}}{1 - a} y_{0} + \frac{E_{n}^{2} h^{2}}{1 - a} .

Denote

H = h / (1 - a)

. The solution of Equation (34) satisfies

y_{n} = \frac{1}{1 + L H} (F_{n} H + (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} y_{n - k} + \frac{1 - 2 a}{1 - a} a^{n - 2} y_{1} + \frac{a^{n - 1}}{1 - a} y_{0} + \frac{E_{n}^{2} h^{2}}{1 - a}) .

The numerical solution of Equation (34) is computed as

u_{n} = \frac{1}{1 + L H} (F_{n} H + (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} u_{n - k} + \frac{1 - 2 a}{1 - a} a^{n - 2} u_{1} + \frac{a^{n - 1}}{1 - a} u_{0}) .

(36)

The numerical solution (36) is calculated with

O (N)

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 + L H) u_{n} = F_{n} H + (1 - a) u_{n - 1} + a ((1 - a) \sum_{k = 2}^{n - 2} a^{k - 1} u_{n - k} + \frac{1 - 2 a}{1 - a} a^{n - 2} u_{1} + \frac{a^{n - 1}}{1 - a} u_{0}),

(1 + L H) u_{n} = F_{n} H + (1 - a) u_{n - 1} + a (u_{n - 1} - F_{n - 1} H) .

The numerical solution of Equation (34) using the approximation

A_{3}

for the first derivative is computed recursively as

u_{0} = y_{0}

and

u_{n} = \frac{1}{1 + L H} ((1 + a L H) u_{n - 1} + H (F_{n} - a F_{n - 1})) .

(37)

Consider the following first-order linear ODE:

y^{'} (x) + L y (x) = (1 + L) e^{x}, y (0) = 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. 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} = \frac{1}{1 + L H} ((1 + a L H) y_{0} + H (F_{1} - a F_{0})) + ε_{1},

(39)

y_{2} = \frac{1}{1 + L H} ((1 + a L H) y_{1} + H (F_{2} - a F_{1})) + ε_{2} .

(40)

From the Taylor theorem, the truncation error

ε_{1}

satisfies

\begin{matrix} (1 + L H) ε_{1} & = (1 + L H) y_{1} - H F_{1} - (1 + a L H) y_{0} - a H F_{0} \\ = y_{1} - H y_{1}^{'} - y_{0} + a H y_{0}^{'} = h y_{0}^{'} - H y_{1}^{'} + a H y_{0}^{'} + \frac{h^{2}}{2} y^{″} (θ_{1}) \\ = \frac{h}{1 - a} (y_{0}^{'} - y_{1}^{'}) + \frac{h^{2}}{2} y^{″} (θ_{1}) = \frac{h^{2}}{2} y^{″} (θ_{1}) - \frac{h^{2}}{1 - a} y^{″} (θ_{2}) . \end{matrix}

Hence,

| ε_{1} | \leq \frac{h^{2}}{2} | y^{″} (θ_{1}) | + \frac{h^{2}}{1 - a} | y^{″} (θ_{2}) | \leq \frac{(3 - a) M h^{2}}{2 (1 - a)} .

In a similar way, we show that

| ε_{2} | \leq \frac{(3 - a) M h^{2}}{2 (1 - a)}

. The errors

e_{1}

and

e_{2}

satisfy the estimates

e_{1} = ε_{1}, | e_{1} | \leq \frac{(3 - a) M h^{2}}{2 (1 - a)},

e_{2} = \frac{1 + a L H}{1 + L H} e_{1} + ε_{2}, | e_{2} | < | e_{1} | + | ε_{2} | < \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 \{\frac{M {(1 + a)}^{2}}{2 L {(1 - a)}^{2}}, \frac{(3 - a) M h}{2 (1 - a)}\} .

Theorem 3.

Let

L > 0

. Then,

| e_{n} | < C h, (1 \leq n \leq N) .

(41)

Proof.

Inequality (41) is satisfied for

n = 1

and

n = 2

. Assume that (41) holds for

n \leq m - 1

. The error

e_{m}

satisfies

e_{m} = \frac{1 + a L H}{1 + L H} e_{m - 1} + \frac{(E_{m}^{2} - a E_{m - 1}^{2}) h^{2}}{(1 + L H) (1 - a)},

| e_{m} | < \frac{1 + a L H}{1 + L H} | e_{m - 1} | + \frac{(| E_{m}^{2} | + a | E_{m - 1}^{2} | h^{2}}{(1 + L H) (1 - a)} .

From Lemma 2

| e_{m} | < (1 - \frac{(1 - a) L H}{1 + L H}) | e_{m - 1} | + \frac{{(1 + a)}^{2} M h^{2}}{2 {(1 - a)}^{2} (1 + L H)} .

Using the assumption of induction,

| e_{m} | < C h - \frac{L h}{1 + L H} C h + \frac{{(1 + a)}^{2} M h^{2}}{2 {(1 - a)}^{2} (1 + L H)} \leq C h .

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, Table 5 and Table 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 (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^{'} (x) + G (x) y (x) = F (x), y (0) = y_{0} .

(42)

The numerical solutions (43) and (45) for Equation (42), which use backward difference approximation and approximation

A_{3}

for first derivative, are computed as

u_{n} = \frac{1}{1 + h G_{n}} (h F_{n} + u_{n - 1}), u_{0} = y_{0},

(43)

u_{n} = \frac{1}{1 - a + h G_{n}} (h F_{n} + {(1 - a)}^{2} \sum_{k = 1}^{n - 2} a^{k - 1} u_{n - k} + (1 - 2 a) a^{n - 2} u_{1} + a^{n - 1} u_{0}) .

(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} = \frac{1}{1 - a + h G_{n}} ((1 - a + a h G_{n - 1}) u_{n - 1} + h (F_{n} - a F_{n - 1})), u_{0} = y_{0} .

(45)

The first-order liner ODE

y^{'} (x) + tan x y (x) = sec x, y (0) = - 1

(46)

has a solution

y (x) = sin x - cos x

. The experimental results of numerical solutions (43) and (45) of Equation (46) are given in Table 7.

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}

\frac{1 + a}{1 - a} = s h^{p}, 1 + a = s h^{p} - a s h^{p}, a = \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 = \frac{1}{s} (2 - s h^{p}) h^{1 - p} = \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 \leq 1$ ,

\begin{matrix} a = \frac{s h^{p} - 1}{s h^{p} + 1} . \end{matrix}

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

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

\begin{matrix} a = 1 - s h^{p} . \end{matrix}

(48)

The experimental results for

p = 0.2, p = 0.5

and

p = 0.7

are given in Table 9.

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

\begin{matrix} a = - 1 + s h^{p} . \end{matrix}

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

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^{'} (x) = L y (x) (1 - y (x)), y (0) = y_{0} .

(50)

Equation (50) has a solution

y = \frac{y_{0}}{y_{0} + (1 - y_{0}) e^{- L x}} .

The explicit schemes of Equation (50) which use first-order backward difference and approximation

A_{3}

of the first derivative satisfy

u_{n} = u_{n - 1} + L h u_{n - 1} (1 - u_{n - 1}), u_{0} = y_{0},

(51)

u_{n} = \frac{L h}{1 - a} u_{n - 1} (1 - u_{n - 1}) + (1 - a) \sum_{k = 1}^{n - 2} a^{k - 1} u_{n - k} + \frac{(1 - 2 a) a^{n - 2}}{1 - a} u_{1} + \frac{a^{n - 1}}{1 - a} u_{0} .

(52)

Numerical method (52) has initial conditions

u_{0} = y_{0}, u_{1} = u_{0} + L h u_{0} (1 - u_{0})

and is computed recursively as

u_{n} = u_{n - 1} + \frac{L h}{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.

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

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:

\{\begin{matrix} u_{t} (x, z) = D u_{x x} (x, t) + F (x, t), (x, t) \in R, \\ u (x, 0) = u_{0} (x), u (0, t) = u_{1} (t), u (1, t) = u_{2} (t), \end{matrix}

(54)

where D is the thermal diffusivity of the material and

R = [0, X] \times [0, T]

. Denote

h = X / N, τ = T / M

, where M and N are positive integers, and

u_{n}^{m} = u (n h, m τ)

, the values of the solution at the nodes

(x_{n}, y_{m}) = (n h, m τ)

of a rectangular grid

J

R

. Now, we construct the finite difference scheme of Equation (54), which uses approximation

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{matrix} \frac{1 - a}{τ} (u_{n}^{m} - (1 - a) \sum_{k = 2}^{m - 2} a^{k - 1} u_{n}^{m - k} - \frac{1 - 2 a}{1 - a} a^{n - 2} u_{n}^{1} - \frac{a^{n - 1}}{1 - a} u_{n}^{0}) \\ = D \frac{u_{n - 1}^{m} - 2 u_{n}^{m} + u_{n + 1}^{m}}{h^{2}} + S_{n}^{m} τ + K_{n}^{m} h^{2} + F_{n}^{m} . \end{matrix}

Denote

η = D τ / ((1 - a) h^{2})

\begin{matrix} - η u_{n - 1}^{m} + (1 + 2 η) u_{n}^{m} - η u_{n + 1}^{m} = (1 - a) \sum_{k = 1}^{m - 2} a^{k - 1} u_{n}^{m - k} + \frac{1 - 2 a}{1 - a} a^{n - 2} u_{n}^{1} + \frac{a^{n - 1}}{1 - a} u_{n}^{0} \\ + \frac{τ F_{n}^{m}}{1 - a} + \frac{τ}{1 - a} (S_{n}^{m} τ + K_{n}^{m} h^{2}) . \end{matrix}

(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} | < \frac{(1 + a) M_{2}}{2 (1 - a)}

, where

M_{2} = max_{(x, t) \in R} |\frac{\partial^{2} u (x, t)}{\partial x^{2}}|, M_{4} = max_{(x, t) \in R} |\frac{\partial^{4} u (x, t)}{\partial x^{4}}| .

Formula (55) is written in the form

\begin{matrix} - η u_{n - 1}^{m} + (1 + 2 η) u_{n}^{m} - η u_{n + 1}^{m} = a [(1 - a) \sum_{k = 2}^{m - 2} a^{k - 2} u_{n}^{m - k} + \frac{1 - 2 a}{1 - a} a^{n - 3} u_{n}^{1} + \frac{a^{n - 2}}{1 - a} u_{n}^{0}] \\ + (1 - a) u_{n}^{m - 1} + \frac{τ F_{n}^{m}}{1 - a} + \frac{τ}{1 - a} (S_{n}^{m} τ + K_{n}^{m} h^{2}) . \end{matrix}

(56)

From Formulas (55) and (56), we obtain

\begin{matrix} - η u_{n - 1}^{m} + (1 + 2 η) u_{n}^{m} - η u_{n + 1}^{m} = a [- η u_{n - 1}^{m - 1} + (1 + 2 η) u_{n}^{m - 1} - η u_{n + 1}^{m - 1}] + (1 - a) u_{n}^{m - 1} \\ - \frac{a τ F_{n}^{m - 1}}{1 - a} - \frac{a τ}{1 - a} (S_{n}^{m - 1} τ + K_{n}^{m - 1} h^{2}) + \frac{τ F_{n}^{m}}{1 - a} + \frac{τ}{1 - a} (S_{n}^{m} τ + K_{n}^{m} h^{2}) . \end{matrix}

The solution of heat Equation (54) satisfies

\begin{matrix} - η u_{n - 1}^{m} + (1 + 2 η) u_{n}^{m} - η u_{n + 1}^{m} = - a η u_{n - 1}^{m - 1} + (1 + 2 a η) u_{n}^{m - 1} - a η u_{n + 1}^{m - 1} \\ + \frac{τ}{1 - a} ({\bar{F}}_{n}^{m} + {\bar{S}}_{n}^{m} τ + {\bar{K}}_{n}^{m} h^{2}), \end{matrix}

(57)

where

{\bar{F}}_{n}^{m} = F_{n}^{m} - a F_{n}^{m - 1}, {\bar{S}}_{n}^{m} = S_{n}^{m} - a S_{n}^{m - 1}, {\bar{K}}_{n}^{m} = K_{n}^{m} - a K_{n}^{m - 1}

. The numerical solution of heat Equation (54) on row m of the grid

J

, which corresponds to (57), satisfies

\begin{matrix} - η U_{n - 1}^{m} + (1 + 2 η) U_{n}^{m} - η U_{n + 1}^{m} = - a η U_{n - 1}^{m - 1} + (1 + 2 a η) U_{n}^{m - 1} - a η U_{n + 1}^{m - 1} + \frac{τ {\bar{F}}_{n}^{m}}{1 - a} . \end{matrix}

(58)

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

(I + η Δ) U_{m} = (I + a η Δ) U_{m - 1} + τ F_{m} + η I_{m},

(59)

where

U_{m}

and

F_{m}

are vectors of dimension

N - 1

which have entries

U_{n}^{m}

and

{\bar{F}}_{n}^{m} / (1 - a)

for

n = 1, \dots, N - 1

and

I_{m} = {(u_{1} (m τ), 0, \dots, 0, u_{2} (m τ))}^{T}

. The matrix

Δ

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

O (M N)

operations because

I + η Δ

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

(n h, m τ)

and

E_{m}

be the vector with the errors

e_{n}^{m}

on row m of the grid

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

J

Lemma 3.

Let

u \in C^{4} [0, X] \times C^{2} [0, T]

. Then,

\begin{matrix} ∥ E_{1} ∥ < \frac{(1 + a) τ}{12 (1 - a)} (6 M_{2} τ + M_{4} h^{2}) . \end{matrix}

(60)

Proof.

The solution of Equation (54) on the first row of the grid

J

satisfies

\begin{matrix} - η u_{n - 1}^{1} + (1 + 2 η) u_{n}^{1} - η u_{n + 1}^{1} = - a η u_{n - 1}^{0} + (1 + 2 a η) & u_{n}^{0} - a η u_{n + 1}^{0} \\ + \frac{τ}{1 - a} (F_{n}^{1} - a F_{n}^{0}) + ε_{n}^{1}, \end{matrix}

\begin{matrix} \frac{(1 - a) (u_{n}^{1} - u_{n}^{0})}{τ} = \frac{D}{h^{2}} (u_{n - 1}^{1} - 2 u_{n}^{1} + u_{n + 1}^{1}) - & \frac{a D}{h^{2}} (u_{n - 1}^{0} - 2 u_{n}^{0} + u_{n + 1}^{0}) \\ + F_{n}^{1} - a F_{n}^{0} + \frac{1 - a}{τ} ε_{n}^{1} . \end{matrix}

The first-order difference approximation satisfies

\frac{u_{n}^{1} - u_{n}^{0}}{τ} = u_{t} (x_{n}, y_{1}) + C_{1} τ = u_{t} (x_{n}, y_{0}) + C_{2} τ,

where

| C_{i} | < M_{2} / 2

for

i = 1, 2

. Then,

\begin{matrix} u_{t} (x_{n}, y_{1}) + C_{1} τ - a u_{t} (x_{n}, y_{0}) - a C_{2} τ = D u_{x x} & (x_{n}, y_{1}) + C_{3} h^{2} - a D u_{x x} (x_{n}, y_{0}) - a C_{4} h^{2} \\ + F (x_{n}, y_{1}) - a F (x_{n}, y_{0}) + \frac{1 - a}{τ} ε_{n}^{1}, \end{matrix}

where

| C_{i} | < M_{4} / 12

for

i = 3, 4

. The error

ε_{n}^{1}

satisfies

\begin{matrix} \frac{1 - a}{τ} ε_{n}^{1} = - C_{1} τ + a C_{2} τ + C_{3} h^{2} - a C_{4} h^{2}, \end{matrix}

\begin{matrix} | ε_{n}^{1} | \leq \frac{τ}{1 - a} (| C_{1} | τ + a | C_{2} | τ + | C_{3} | h^{2} + a | C_{4} | h^{2}) . \end{matrix}

Hence,

\begin{matrix} | ε_{n}^{1} | < \frac{(1 + a) τ}{12 (1 - a)} (6 M_{2} τ + M_{4} h^{2}) \end{matrix}

(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

J

satisfy

\begin{matrix} - η e_{n - 1}^{1} + (1 + 2 η) e_{n}^{1} - η e_{n + 1}^{1} = ε_{n}^{1} . \end{matrix}

The elements of

E_{0}

are equal to zero,

e_{n}^{0} = 0

, and the vector

E_{1}

with the errors of (59) on the first row of the grid

J

satisfies

(I + η Δ) E_{1} = W_{1},

(62)

where

W_{1}

is the vector with elements

ε_{n}^{1}

. From (61), the infinity norm of

W_{1}

satisfies

∥ W_{1} ∥ < \frac{(1 + a) τ}{12 (1 - a)} (6 M_{2} τ + M_{4} h^{2}) .

The matrix

I + η Δ

is a diagonally dominant M-matrix and its infinity norm satisfies [43,44]

∥ {(I + η Δ)}^{- 1} ∥ < 1 .

Therefore,

∥ E_{1} ∥ \leq ∥ {(I + η Δ)}^{- 1} ∥ \cdot ∥ W_{1} ∥ < \frac{(1 + a) τ}{12 (1 - a)} (6 M_{2} τ + M_{4} h^{2}) .

□

The vector

E_{m}

with the errors of difference scheme (59) on row m of the grid

J

satisfies

(I + η Δ) E_{m} = (I + a η Δ) E_{m - 1} + τ W_{m},

(63)

where

W_{m}

is the vector with entries

{\bar{S}}_{n}^{m} τ + {\bar{K}}_{n}^{m} h^{2}

. The infinity norm of

W_{m}

satisfies

∥ W_{m} ∥ \leq M_{0} (τ + h^{2})

where

M_{0} = max \{\frac{{(1 + a)}^{2} M_{2}}{2 (1 - a)}, \frac{(1 + a) M_{4}}{12}\} .

The matrix

Δ

has eigenvalues [45]

λ_{k} = 4 {sin}^{2} (\frac{k π}{2 N}), k = 1, \dots, N - 1

and eigenvectors

V_{k}

for

k = 1, \dots, N - 1

, which have components

v_{k}^{i} = sin (k π i / N)

, where

i = 1, \dots, N - 1

. Denote by P and Q the following matrices:

P = {(I + η Δ)}^{- 1} (I + a η Δ), Q = {(I + η Δ)}^{- 1} .

The sequence of vectors

E_{m}

satisfies

E_{m} = P E_{m - 1} + τ Q W_{m} .

(64)

The matrix P has eigenvalues

(I + a η λ_{k}) / (I + η λ_{k})

for

k = 1, \dots, N - 1

. The spectral radius of the matrix P is smaller than one,

ρ (P) < 1

, and its powers

P^{n}

converge to zero. Therefore,

{lim}_{n \to \infty} {∥ P ∥}^{n} = 0

. Denote

p = max_{n \geq 1} {∥ P ∥}^{n}, q = ∥ Q ∥ .

In Theorem 4, we prove that finite difference scheme (59) is unconditionally stable and has an accuracy

O (τ + h^{2})

Theorem 4.

Let

u \in C^{4} [0, X] \times C^{2} [0, T]

. Then,

\begin{matrix} ∥ E_{m} ∥ < C (τ + h^{2}), \end{matrix}

where

C = p (1 + T q) M_{0}

and all

m = 2, \dots M

Proof.

By applying (64) successively

m - 1

times, we obtain

E_{m} = P^{m - 2} E_{1} + τ \sum_{k = 0}^{m - 2} P^{k} Q W_{m - k} .

Then,

∥ E_{m} ∥ \leq ∥ P^{m - 2} ∥ \cdot ∥ E_{1} ∥ + τ \sum_{k = 0}^{m - 2} ∥ P^{k} ∥ \cdot ∥ Q ∥ \cdot ∥ W_{m - k} ∥ .

From Lemma 3, the norm of the vector

E_{1}

satisfies

∥ E_{1} ∥ < M_{0} (τ + h^{2})

. Hence,

∥ E_{m} ∥ < p M_{0} (τ + h^{2}) + p q M_{0} τ (τ + h^{2}) (m - 1) < p (1 + T q) M_{0} (τ + h^{2}) .

□

Consider the following heat equation

\{\begin{matrix} u_{t} (x, z) = 3 u_{x x} (x, t) - 4 e^{- x - t}, \\ u (x, 0) = e^{- x}, u (0, t) = e^{- t}, u (1, t) = e^{- 1 - t} . \end{matrix}

(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 = τ = 0.01

are shown in Figure 1. Numerical method (59) has an accuracy

O (τ + h^{2})

for

a = - 0.5

. When the parameter

a = (τ - 1) / (τ + 1)

, numerical method (59) has a second-order accuracy

O (τ^{2} + h^{2})

. Experimental results for the error and order of the numerical solution (59) of Equation (65) with

a = - 0.5

and

a = (τ - 1) / (τ + 1)

are given in Table 12 and Table 13.

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_{τ} / E_{2 τ} |

. Similarly, the order in space is computed as

{log}_{2} | E_{h} / E_{2 h} |

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

A_{3}

for the first derivative. Two approximations of the fractional derivative of order

2 - α

whose weights satisfy (1) are given below. Let

h = (x - l) / N

. The L1 approximation of the Caputo fractional derivative has the form

\frac{1}{Γ (2 - α) h^{α}} \sum_{k = 0}^{n} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + O (h^{2 - α}),

(66)

where

σ_{0}^{(α)} = 1, σ_{n}^{(α)} = {(n - 1)}^{1 - α} - n^{1 - α}

σ_{k}^{(α)} = {(k + 1)}^{1 - α} - 2 k^{1 - α} + {(k - 1)}^{1 - α}, (1 \leq k \leq n - 1) .

The L1 approximation and approximation (5) have an order

2 - α

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 - α

by substituting the first derivative in (2) with a first-order backward difference approximation.

\frac{1}{2 Γ (1 - α) h^{α}} \sum_{k = 0}^{n} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + O (h^{2 - α}),

(67)

where

σ_{0}^{(α)} = 1 - 2 ζ (α), σ_{1}^{(α)} = \frac{1}{2^{α}} + 2 ζ (α)

σ_{k}^{(α)} = \frac{1}{{(k + 1)}^{α}} - \frac{1}{{(k - 1)}^{α}}, (2 \leq k \leq n - 2),

σ_{n - 1}^{(α)} = - \frac{1}{{(n - 2)}^{α}} - 2 s_{n}, σ_{n}^{(α)} = - \frac{1}{{(n - 1)}^{α}} + 2 s_{n} .

and

s_{n} = \sum_{k = 1}^{n - 1} k^{- α} - n^{1 - α} / (1 - α) - ζ (α)

. Now, we apply approximation

A_{3}

for constructing an approximation of the Caputo derivative. By substituting the first derivative

f_{n}^{'}

in Formula (2) with

A_{3}

, we obtain an approximation of order

2 - α

\frac{1}{2 Γ (1 - α) h^{α}} \sum_{k = 0}^{n} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + O (h^{2 - α}),

(68)

where

σ_{0}^{(α)} = 1 - 2 (1 - b) ζ (α), σ_{1}^{(α)} = \frac{1}{2^{α}} + 2 {(1 - b)}^{2} ζ (α)

σ_{k}^{(α)} = \frac{1}{{(k + 1)}^{α}} - \frac{1}{{(k - 1)}^{α}} + 2 {(1 - b)}^{2} b^{k - 1} ζ (α), (2 \leq k \leq n - 2),

σ_{n - 1}^{(α)} = - \frac{1}{{(n - 2)}^{α}} - 2 s_{n} + 2 (1 - 2 b) b^{n - 2} ζ (α),

σ_{n}^{(α)} = - \frac{1}{{(n - 1)}^{α}} + 2 s_{n} + 2 b^{n - 1} ζ (α) .

The value of zeta function

ζ (α)

is negative when

0 < α < 1

and the parameter b of approximation

A_{3}

has a modulus smaller than one,

| b | < 1

. Therefore,

σ_{0}^{(α)} > 0

and

σ_{k}^{(α)} < 0

for

k > 1

. In Lemma 4, we show that the coefficient

σ_{1}^{(α)}

of approximation (68) is negative when the parameter b has a value

b = α / 2

Lemma 4.

Let

b = a / 2

. Then,

\begin{matrix} σ_{1}^{(α)} < 0 . \end{matrix}

(69)

Proof.

σ_{1}^{(α)} = \frac{1}{2^{α}} + 2 {(1 - \frac{α}{2})}^{2} ζ (α) = \frac{1}{2^{α}} - \frac{1}{2} {(2 - a)}^{2} | ζ (α) | .

Inequality (69) holds when

\frac{1}{2^{α}} - \frac{1}{2} {(2 - α)}^{2} | ζ (α) | < 0,

| ζ (α) | > \frac{2}{2^{α} {(2 - α)}^{2}} .

Taking the logarithm of both sides, we obtain

ln | ζ (α) | > ln 2 - α ln 2 - 2 ln (2 - α) .

Let

f (α) = ln | ζ (α) | - ln 2 + α ln 2 + 2 ln (2 - α)

. The function f satisfies

f (0) = 0

and has first derivative

f^{'} (α) = \frac{ζ^{'} (α)}{ζ (α)} + ln 2 - \frac{2}{2 - α} .

It is sufficient to prove that

f^{'}

is positive. The zeta function and its derivative have Laurent series expansions

\begin{matrix} ζ (α) = - \frac{1}{1 - a} + \sum_{k = 0}^{\infty} \frac{γ_{k}}{k!} {(1 - a)}^{k}, \\ ζ^{'} (α) = - \frac{1}{{(1 - α)}^{2}} - \sum_{k = 0}^{\infty} \frac{γ_{k + 1}}{k!} {(1 - a)}^{k}, \end{matrix}

where

γ_{k}

are Stieltjes constants

γ_{0} = 0.5772, γ_{1} = - 0.0728, γ_{2} = - 0.0097, γ_{3} = 0.0021, γ_{4} = 0.0023, \dots .

Hence,

| ζ (α) | < \frac{1}{1 - α}, | ζ^{'} (α) | > \frac{1}{{(1 - α)}^{2}} + γ_{1} + γ_{2} (1 - α),

| ζ^{'} (α) | > \frac{1}{{(1 - α)}^{2}} - 0.08 - 0.01 (1 - α) .

The derivative of the function f satisfies

f^{'} (α) > \frac{1}{1 - α} - 0.08 (1 - α) - 0.01 {(1 - α)}^{2} + ln 2 - \frac{2}{2 - α} .

Denote

g (α) = \frac{1}{1 - α} - 0.08 (1 - α) - 0.01 {(1 - α)}^{2} + ln 2 - \frac{2}{2 - α} .

The function g has first-, second- and third-order derivatives

g^{'} (α) = \frac{1}{{(1 - α)}^{2}} + 0.08 + 0.02 (1 - α) - \frac{2}{{(2 - α)}^{2}},

g^{″} (α) = \frac{2}{{(1 - α)}^{3}} - 0.02 - \frac{4}{{(2 - α)}^{3}},

g^{‴} (α) = \frac{6}{{(1 - α)}^{4}} - \frac{12}{{(2 - α)}^{4}} .

The values of the function g and its derivatives

g^{'}, g^{″}, g^{‴}

at zero are positive:

g (0) = ln 2 - 0.09 > 0, g^{'} (0) = 0.6, g^{″} (0) = 1.48, g^{‴} (0) = 5.25 .

The third derivative

g^{‴}

is positive on the interval

[0, 1]

because

{(\frac{2 - α}{1 - α})}^{4} > 2, 1 + \frac{1}{1 - α} > \sqrt[4]{2}, \frac{1}{1 - α} > 1 > \sqrt[4]{2} - 1 .

Therefore, the functions

g, g^{'}, g^{″}

are positive and increasing and

f^{'} (α) > g (α) > 0

. The function f is increasing and positive because

f (0) = 0

. □

Approximation (5) of the fractional derivative is obtained from (68) and value of the parameter

b = α / 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^{(α)} (x) + L y (x) = F (x), y (0) = y_{0} .

(70)

The numerical solution

{u_{n}}_{n = 1}^{N}

of Equation (70), which uses the L1 approximation (66) of the fractional derivative, is computed as

u_{m} = \frac{1}{σ_{0}^{(α)} + L Γ (2 - α) h^{α}} (L Γ (2 - α) h^{α} F_{m} - \sum_{k = 1}^{m} σ_{k}^{(α)} u_{m - k}) .

(71)

The numerical solution of Equation (70), which uses approximations (67) and (68) of the fractional derivative

y_{n}^{(α)}

, satisfies

u_{m} = \frac{1}{σ_{0}^{(α)} + 2 L Γ (1 - α) h^{α}} (2 L Γ (1 - α) h^{α} F_{m} - \sum_{k = 1}^{m} σ_{k}^{(α)} u_{m - k})

(72)

and has initial conditions

u_{0} = y_{0}, u_{1} = \frac{h^{α} Γ (2 - α) F_{1} + u_{0}}{L h^{α} Γ (2 - α) + 1}, u_{2} = \frac{{(2 h)}^{α} Γ (2 - α) F_{2} + u_{0}}{L {(2 h)}^{α} Γ (2 - α) + 1} .

(73)

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

\frac{f (ℓ) - f (0)}{Γ (2 - α) ℓ^{α}} = f^{(α)} (ℓ) + O (ℓ^{2 - α}) .

(74)

Consider the two-term equation

y^{(α)} (x) + L y (x) = x^{1 - α} E_{1, 2 - α} (x) + L e^{x}, y (0) = 1 .

(75)

Equation (75) has a solution

y (x) = e^{x}

, where

E_{1, 2 - α} (x)

is the Mittag–Leffler function

E_{1, 2 - α} (x) = \sum_{k = 0}^{\infty} \frac{x^{k}}{Γ (k + 2 - α)} .

The experimental results of numerical solutions (71) and (72) of Equation (75) on the interval

[0, 1]

are given in Table 14 and Table 15. 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

σ_{n - 1}^{(α)}

and

σ_{n}^{(α)}

of approximations (67) and (68) for

n = 3, \dots, N

are computed recursively in

O (N)

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}^{'}

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

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

α

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:

\{\begin{matrix} D_{t}^{α} u (x, t) = D u_{x x} (x, t) + F (x, t), \\ u (x, 0) = u_{0} (x), u (0, t) = u_{1} (t), u (1, t) = u_{2} (t) . \end{matrix}

(76)

where

(x, t) \in [0, 1] \times [0, 1]

and

D_{t}^{α} u (x, t)

denotes the Caputo derivative of the solution

u (x, t)

in time. Let

τ = 1 / M

and

h = 1 / N

, where M and N are positive integers, and let

J

be a rectangular grid with nodes

(x_{n}, y_{n}) = (n h, m h)

. The numerical solution

U_{n}^{m}

on the first three rows of

J

is computed with the approximation (74) of the partial fractional derivative and second-order central difference approximation of the partial derivative in space

\frac{U_{n}^{m} - U_{n}^{0}}{Γ (2 - α) {(m τ)}^{α}} = \frac{D}{h^{2}} (U_{n - 1}^{m} - 2 U_{n}^{m} + U_{n + 1}^{m}) + F_{n}^{m},

U_{n}^{m} - U_{n}^{0} = \frac{D Γ (2 - α) {(m τ)}^{α}}{h^{2}} (U_{n - 1}^{m} - 2 U_{n}^{m} + U_{n + 1}^{m}) + Γ (2 - α) {(m τ)}^{α} F_{n}^{m} .

Denote

η_{m} = D Γ (2 - α) {(m τ)}^{α} / h^{2}

. The numerical solution satisfies the system of linear equations

- η_{m} U_{n - 1}^{m} + (1 + 2 η_{m}) U_{n}^{m} - η_{m} U_{n + 1}^{m} = U_{n}^{0} + Γ (2 - α) {(m τ)}^{α} F_{n}^{m}

and has initial conditions

U_{n}^{0} = u_{0} (n h)

for

n = 0, 1, \dots, N

and boundary conditions

U_{0}^{m} = u_{1} (m τ), U_{N}^{m} = u_{2} (m τ)

for

m = 1, 2, 3

. The numerical solution on row m of the grid

J

, where

m > 3

, which uses approximation (5) of the fractional derivative, satisfies

\frac{1}{2 Γ (1 - α) τ^{α}} \sum_{k = 0}^{m} σ_{k}^{(α)} U_{n}^{m - k} = \frac{D}{h^{2}} (U_{n - 1}^{m} - 2 U_{n}^{m} + U_{n + 1}^{m}) + F_{n}^{m},

σ_{0}^{(α)} U_{n}^{m - k} + \sum_{k = 1}^{m} σ_{k}^{(α)} U_{n}^{m - k} = \frac{2 D Γ (1 - a) τ^{α}}{h^{2}} (U_{n - 1}^{m} - 2 U_{n}^{m} + U_{n + 1}^{m}) + 2 Γ (1 - a) τ^{α} F_{n}^{m} .

Denote

η = 2 D Γ (1 - a) τ^{α} / h^{2}

. The numerical solution on row m satisfies the tridiagonal system of linear equations

- η U_{n - 1}^{m} + (σ_{0}^{(α)} + 2 η) U_{n}^{m} - η U_{n + 1}^{m} = Γ (2 - α) τ^{α} F_{n}^{m} - \sum_{k = 1}^{m} σ_{k}^{(α)} U_{n}^{m - k} .

(77)

Consider the fractional subdiffusion equation

\{\begin{matrix} D_{t}^{α} u (x, t) = D u_{x x} (x, t) + e^{x} (t^{1 - α} E_{1, 2 - α} (t) - D e^{t}), \\ u (x, 0) = e^{x}, u (0, t) = e^{t}, u (1, t) = e^{1 + t} . \end{matrix}

(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 = τ = 0.01

are shown in Figure 2. Experimental results for the error and order of numerical solution (77) of fractional subdiffusion Equation (78) and

α = 0.5, 0.75

are given in Table 16 and Table 17.

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

O (N)

operations. The numerical solutions of partial differential equations are computed with

O (M N)

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 - α

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.

Author Contributions

Conceptualization, Y.D. and S.G.; data curation, V.T.; formal analysis, Y.D. and V.T.; funding acquisition, S.G.; investigation, Y.D.; methodology, V.T. and S.G.; project administration, Y.D. and V.T.; resources, S.G.; software, S.G.; supervision, Y.D.; validation, V.T.; visualization, S.G.; writing—original draft, Y.D. and S.G.; writing—review and editing, Y.D. and V.T. All authors have read and agreed to the published version of the manuscript.

Funding

This study is supported by the Bulgarian National Science Fund under Project KP-06-M62/1 “Numerical deterministic, stochastic, machine and deep learning methods with applications in computational, quantitative, algorithmic finance, biomathematics, ecology and algebra” from 2022.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Akgül, A.; Conejero, J.A. Fractal fractional derivative models for simulating chemical degradation in a bioreactor. Axioms 2024, 13, 151. [Google Scholar] [CrossRef]
Akgül, A.; Khoshnaw, S.H.A. Application of fractional derivative on non-linear biochemical reaction models. Int. J. Intell. Netw. 2020, 1, 2–58. [Google Scholar] [CrossRef]
Baba, I.A.; Rihan, F.A. A fractional–order model with different strains of COVID-19. Phys. A Stat. Mech. Appl. 2022, 603, 127813. [Google Scholar] [CrossRef] [PubMed]
Gu, Y.; Yu, Z.; Lan, P.; Lu, N. fractional derivative viscosity of ANCF cable element. Actuators 2023, 12, 64. [Google Scholar] [CrossRef]
Hilfer, R. Applications of Fractional Calculus in Physics; World Scientific: Singapore, 2003. [Google Scholar]
Kumar, D.; Singh, J. Fractional Calculus in Medical and Health Science; CRC Press: Boca Raton, FL, USA, 2020. [Google Scholar]
Dimitrov, Y.; Georgiev, S.; Miryanov, R.; Todorov, V. Convergence of the L1 two-term equation scheme. J. Phys. Conf. Ser. 2023, 2675, 012027. [Google Scholar] [CrossRef]
Jin, B.; Lazarov, R.; Zhou, Z. An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal. 2016, 36, 197–221. [Google Scholar] [CrossRef]
Li, B.; Xie, X.; Yan, Y. L1 scheme for solving an inverse problem subject to a fractional diffusion equation. Comput. Math. Appl. 2023, 134, 112–123. [Google Scholar] [CrossRef]
Scherer, R.; Kalla, S.L.; Tang, Y.; Huang, J. The Grünwald—Letnikov method for fractional differential equations. Comput. Math. Appl. 2011, 62, 902–917. [Google Scholar] [CrossRef]
Arshad, S.; Defterli, O.; Baleanu, D. A second order accurate approximation for fractional derivatives with singular and non-singular kernel applied to a HIV model. Appl. Math Comput. 2020, 374, 125061. [Google Scholar] [CrossRef]
Nasir, H.M.; Nafa, K. A new second order approximation for fractional derivatives with applications. Sultan Qaboos Univ. J. Sci. 2018, 23, 43–55. [Google Scholar] [CrossRef]
Alikhanov, A.A. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys. 2015, 280, 424–438. [Google Scholar] [CrossRef]
Gao, G.-H.; Sun, Z.-Z.; Zhang, H.-W. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. J. Comput. Phys. 2014, 259, 33–50. [Google Scholar] [CrossRef]
Xing, Y.; Yan, Y. A higher order numerical method for time fractional partial differential equations with nonsmooth data. J. Comput. Phys. 2018, 357, 305–323. [Google Scholar] [CrossRef]
Chen, M.H.; Deng, W.H. Fourth order accurate scheme for the space fractional diffusion equations. SIAM J. Numer. Anal. 2014, 52, 1418–1438. [Google Scholar] [CrossRef]
Gunarathna, W.A.; Nasir, H.M.; Daundasekera, W.B. An explicit form for higher order approximations of fractional derivatives. Appl. Numer. Math. 2019, 143, 51–60. [Google Scholar] [CrossRef]
Hao, Z.-P.; Sun, Z.-Z.; Cao, W.-R. A fourth-order approximation of fractional derivatives with its applications. J. Comput. Phys. 2015, 281, 787–805. [Google Scholar] [CrossRef]
Li, L.; Zhao, D.; She, M.; Chen, X. On high order numerical schemes for fractional differential equations by block-by-block approach. Appl. Math. Comput. 2022, 425, 127098. [Google Scholar] [CrossRef]
Lubich, C. Discretized fractional calculus. SIAM J. Math. Anal. 1986, 17, 704–719. [Google Scholar] [CrossRef]
Nasrollahzadeh, F.; Hosseini, S.M. An implicit difference ADI method for two-dimensional space-time fractional diffusion equation. Iran. J. Math. Sci. Inform. 2016, 11, 71–86. [Google Scholar] [CrossRef]
Wang, Z.; Liang, Y.; Mo, Y. A novel high order compact ADI scheme for two dimensional fractional integro-differential equations. Appl. Numer. Math. 2021, 167, 257–272. [Google Scholar] [CrossRef]
Dimitrov, Y. Approximations for the Caputo derivative (I). J. Fract. Calc. Appl. 2018, 9, 15–44. [Google Scholar]
Ahnert, K.; Abel, M. Numerical differentiation of experimental data: Local versus global methods. J. Comput. Phys. Commun. 2005, 177, 764–774. [Google Scholar] [CrossRef]
Hanke, M.; Scherzer, O. Inverse problems light: Numerical differentiation. Am. Math. Mon. 2001, 108, 512–521. [Google Scholar] [CrossRef]
Apostolov, S.; Dimitrov, Y.; Todorov, V. Constructions of second order approximations of the Caputo fractional derivative. In Lecture Notes in Computer Science; Lirkov, I., Margenov, S., Eds.; Large-Scale Scientific Computing. LSSC 2021; Springer: Berlin/Heidelberg, Germany, 2022; p. 13127. [Google Scholar]
Todorov, V.; Dimitrov, Y.; Dimov, I. Second order shifted approximations for the first derivative. In Studies in Computational Intelligence; Dimov, I., Fidanova, S., Eds.; Advances in High Performance Computing. HPC 2019; Springer: Cham, Switzerland, 2011; p. 902. [Google Scholar]
Dimitrov, Y. A second order approximation for the Caputo fractional derivative. J. Fract. Calc. Appl. 2016, 7, 175–195. [Google Scholar]
Al-Bar, R. On the approximate solution of fractional logistic differential equation using operational matrices of Bernstein polynomials. Appl. Math. 2015, 6, 2096–2103. [Google Scholar] [CrossRef]
Iannelli, M.; Pugliese, A. An Introduction to Mathematical Population Dynamics: Along the Trail of Volterra and Lotka, 4th ed.; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
Kot, M. Elements of Mathematical Ecology; Cambridge University Press: Cambridge, UK, 2001. [Google Scholar]
Hahn, D.W.; Özişik, M.N. Heat Conduction, 3rd ed.; John Wiley & Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
Incorpera, F.P.; DeWitt, D.P.; Bergman, T.L.; Lavine, A.S. Fundamentals of Heat and Mass Transfer, 6th ed.; John Wiley: New York, NY, USA, 2006. [Google Scholar]
Majumdar, P. Computational Methods for Heat and Mass Transfer; CRC Press: Boca Raton, FL, USA, 2006. [Google Scholar]
Morton, K.W.; Mayers, D.F. Numerical Solution of Partial Differential Equations: An Introduction, 2nd ed.; Cambridge University Press: Cambridge, UK, 2005. [Google Scholar]
Zhao, S. A matched alternating direction implicit (ADI) method for solving the heat equation with interfaces. J. Sci. Comput. 2015, 63, 118–137. [Google Scholar] [CrossRef]
Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: Cambridge, UK, 2007. [Google Scholar]
Datta, B.N. Numerical Linear Algebra and Applications, 2nd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2010. [Google Scholar]
Suárez-Carreño, F.; Rosales-Romero, L. Convergency and stability of explicit and implicit schemes in the simulation of the heat equation. Appl. Sci. 2021, 11, 4468. [Google Scholar] [CrossRef]
Sundaram, A. Numerical solution for the heat equation using Crank-Nicolson difference method. Int. J. Res. Appl. Sci. Eng. Technol. 2023, 11, 158–164. [Google Scholar] [CrossRef]
Ames, W. Numerical Methods for Partial Differential Equations; Academic Press: Boston, MA, USA, 1992. [Google Scholar]
Smith, G.D. Numerical Solution of Partial Differential Equations: Finite Difference Method; Oxford University Press: Oxford, UK, 1992. [Google Scholar]
Kolotilina, L.Y. Bounds for the infinity norm of the inverse for certain M- and H-matrices. Linear Algebra Appl. 2009, 430, 692–702. [Google Scholar] [CrossRef]
Varah, J.M. A lower bound for the smallest singular value of a matrix. Linear Algebra Appl. 1975, 11, 3–5. [Google Scholar] [CrossRef]
Kulkarni, D.; Schmidt, D.; Tsui, S.-K. Eigenvalues of tridiagonal pseudo-Toeplitz matrices. Linear Algebra Appl. 1999, 297, 63–80. [Google Scholar] [CrossRef]
Dimitrov, Y.; Georgiev, S.; Todorov, V. Approximation of Caputo fractional derivative and numerical solutions of fractional differential equations. Fractal Fract. 2023, 7, 750. [Google Scholar] [CrossRef]
Li, C.; Zeng, F. Numerical Methods for Fractional Calculus; Chapman and Hall/CRC: Boca Raton, FL, USA, 2015. [Google Scholar]
Avcı, İ.; Mahmudov, N.I. Numerical solutions for multi-term fractional order differential equations with fractional Taylor operational matrix of fractional integration. Mathematics 2020, 8, 96. [Google Scholar] [CrossRef]
Diethelm, K.; Ford, N.J. Analysis of fractional differential equations. J. Math. Anal. Appl. 2002, 265, 229–248. [Google Scholar] [CrossRef]
Edwards, J.T.; Ford, N.J.; Simpson, A.C. The numerical solution of linear multi-term fractional differential equations: Systems of equations. J. Comput. Appl. Math. 2002, 148, 401–418. [Google Scholar] [CrossRef]
El-Mesiry, A.E.M.; El-Sayed, A.M.A.; El-Saka, H.A.A. Numerical methods for multi-term fractional (arbitrary) orders differential equations. Appl. Math. Comput. 2005, 160, 683–699. [Google Scholar] [CrossRef]
Ali, U.; Sohail, M.; Abdullah, F.A. An efficient numerical scheme for variable-order fractional sub-diffusion equation. Symmetry 2020, 12, 1437. [Google Scholar] [CrossRef]
Ashurov, R.; Umarov, S. Determination of the order of fractional derivative for subdiffusion equations. Fract. Calc. Appl. Anal. 2020, 23, 1647–1662. [Google Scholar] [CrossRef]
Langlands, T.A.M.; Henry, B.I. The accuracy and stability of an implicit solution method for the fractional diffusion equation. J. Comput. Phys. 2005, 205, 719–736. [Google Scholar] [CrossRef]
Mainardi, F.; Mura, A.; Pagnini, G.; Gorenflo, R. Sub-diffusion equations of fractional order and their fundamental solutions. In Mathematical Methods in Engineering; Taş, K., Tenreiro Machado, J.A., Baleanu, D., Eds.; Springer: Berlin/Heidelberg, Germany, 2007; pp. 23–55. [Google Scholar]
Podlubny, I. Fractional Differential Equations; Academic Press: Cambridge, MA, USA, 1999. [Google Scholar]

$Fractalfract 08 00608 g001$

Figure 1. Graphs of numerical solution (59) (left) and its error (right) for

h = τ = 0.01

Figure 1. Graphs of numerical solution (59) (left) and its error (right) for

h = τ = 0.01

$Fractalfract 08 00608 g001$

$Fractalfract 08 00608 g002$

Figure 2. Graphs of numerical solution (77) (left) and its error (right) for

h = τ = 0.01

Figure 2. Graphs of numerical solution (77) (left) and its error (right) for

h = τ = 0.01

$Fractalfract 08 00608 g002$

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

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

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

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

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

a = 1 - h^{0.8}

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

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

Table 8. Error and order of numerical method (45) of Equation (50) and

a = (s h^{p} - 1) / (s h^{p} + 1)

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

Table 9. Error and order of numerical method (45) of Equation (50) and

a = 1 - s h^{p}

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

Table 10. Error and order of numerical method (45) of Equation (50) and

a = - 1 + s h^{p}

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

Table 11. Error and order of numerical methods (51) and (53) of Equation (50) and

L = 1, y_{0} = 2

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

Table 12. Error and order of numerical solution (59) and

N = 1000

Table 12. Error and order of numerical solution (59) and

N = 1000

$τ$	$α = - 0.5$		$a = (τ - 1) / (τ + 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$

Table 13. Error and order of numerical solution (59) and

M = 1000

Table 13. Error and order of numerical solution (59) and

M = 1000

h	$α = - 0.5$		$a = (τ - 1) / (τ + 1)$
h	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$

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

h	NS(71), $α = 0.25, L = 1$		NS(72), $α = 0.5, L = 10$		$α = 0.75, L = 5, b = - 0.5$
h	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$

Table 15. Error and order of numerical solution (72) of Equation (75) and

b = α / 2

Table 15. Error and order of numerical solution (72) of Equation (75) and

b = α / 2

h	NS(72), $α = 0.25, L = 1$		NS(72), $α = 0.5, L = 10$		NS(72), $α = 0.75, L = 5$
h	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$

Table 16. Error and order of numerical solution (77) and

N = 1000

Table 16. Error and order of numerical solution (77) and

N = 1000

$τ$	$α = 0.5$		$α = 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$

Table 17. Error and order of numerical solution (77) and

N = 1000

Table 17. Error and order of numerical solution (77) and

N = 1000

h	$α = 0.5$		$α = 0.75$
h	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$

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Dimitrov, Y.; Georgiev, S.; Todorov, V. First Derivative Approximations and Applications. Fractal Fract. 2024, 8, 608. https://doi.org/10.3390/fractalfract8100608

AMA Style

Dimitrov Y, Georgiev S, Todorov V. First Derivative Approximations and Applications. Fractal and Fractional. 2024; 8(10):608. https://doi.org/10.3390/fractalfract8100608

Chicago/Turabian Style

Dimitrov, Yuri, Slavi Georgiev, and Venelin Todorov. 2024. "First Derivative Approximations and Applications" Fractal and Fractional 8, no. 10: 608. https://doi.org/10.3390/fractalfract8100608

APA Style

Dimitrov, Y., Georgiev, S., & Todorov, V. (2024). First Derivative Approximations and Applications. Fractal and Fractional, 8(10), 608. https://doi.org/10.3390/fractalfract8100608

Article Menu

First Derivative Approximations and Applications

Abstract

1. Introduction

2. Asymptotic Formula

3. Error Estimates and Convergence

4. Numerical Solutions of First-Order ODEs

5. Numerical Solution of Heat Equation

6. Approximations of the Fractional Derivative

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI