First Derivative Approximations and Applications

Preprint

Article

First Derivative Approximations and Applications

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

06 September 2024

Posted:

09 September 2024

You are already at the latest version

Alerts

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 first derivative for 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:

Subject: Computer Science and Mathematics - Computational Mathematics

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 have been used to describe phenomena that cannot be adequately explained by integer-order derivatives. Fractional differential equations have found applications in various fields, including physics, engineering, chemistry and biology [1,2,3,4,5,6]. Numerical methods are used to analyze models that involve fractional differential equations. The finite difference schemes for solving numerically fractional differential equations use approximations of the fractional derivative. Consider the interval

[l, x]

and a uniform grid with a step size of

h = (x - l) / N

, where N is a positive integer. Denote

x_{n} = l + 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)}^{α}

. Grünwald approximation is a second-order shifted approximation of he 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 are the finite-difference approximations which have first-order accuracy and second-order at the midpoint. Finite-difference approximations are local numerical differentiation methods and are a special case of the Grünwald difference approximation, where the order is an integer. Methods for global numerical differentiation are studied in [24,25]. The methods for constructing approximations of fractional derivatives using a generating function also apply to the construction of integer-order derivative approximations. In [26,27] we construct approximations of the first and second derivatives whose generating functions are transformations of the exponential and logarithmic functions. In [26,28] we construct second-order approximations of the fractional derivative using the asymptotic formula of the L1 approximation and second derivative approximations. In this paper we apply the method from [26,27] for constructing an approximation of first derivative and its second-order asymptotic formula

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

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 get

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

(NS1)

The value of

u_{n}

is computed with one addition and one multiplication, assuming that the values

F_{n}

of the function F at the nodes

x_{n}

of the grid are known. The total number of arithmetic operations of Euler method is

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

(11)

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

(NS2)

Numerical solution NS1 is a special case of NS2 when the parameter a is zero,

a = 0

. Direct computation of numerical solution NS2 using the above formula requires

O (N^{2})

operations. The weights of

A_{2}

consist of powers of the parameter a, which allows to compute NS2 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 NS2 is given below.

\begin{matrix} Step 1 : & A = 1 - a; H = h / A; S = u_{0} = y_{0}; u_{1} = u_{0} + h * F_{1}; \\ Step 2 : & 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 on Step 1 and five operations for calculating each value of the sequence

u_{n}

for

n = 2, \dots, N

in Step 2. Total number of operations for computation of numerical solution NS2 is

5 N - 1

. Consider the ordinary differential equation

y^{'} (x) = e^{x}, y (0) = 1,

(12)

where

x \in [0, 1]

. The numerical results for the error and order of numerical solutions 11 and NS1 of equation (NS2) are given in Table 1. The calculation of the numerical solutions was performed with Mathematica software. Although formula (9) is valid only asymptotically, the experimental results in Table 1 suggest that numerical method NS1 has a first-order accuracy. This fact is proven in the next section. The time for computation of the values of

F_{n} = e^{n h}

is greater than the computational time of the numerical solutions. We observed a small difference of the computational time of the numerical solutions in Table 1, which is a fraction of a second. The finite geometric series satisfies

\sum_{k = 0}^{n - 1} a^{k} = \frac{1 - a^{n}}{1 - a} .

(13)

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

(14)

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

(15)

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 .

(16)

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 (13) and (14) we get

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

(17)

From (16) and (17) 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)

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 get

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

(NS3)

Numerical solution NS3 has initial conditions

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

. The algorithm and pseudocode for computing the numerical solution NS3 are given below:

u_{n} = \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} Step 1 : & 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; \\ Step 2 : & 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 on 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 NS3 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}),

(18)

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

(19)

Equation (19) has a solution

y (x) = sin x - cos x

. The experimental results of numerical solution NS3 for equation (19) and values of the parameter

a = 0.5, a = - 1

and

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

are given in Table 2.

3. Error Estimates and Convergence

In this section we derive estimates for the errors of approximations

A_{2}

and

A_{3}

of the first derivative and the errors of numerical solutions NS2 and NS3 of equation (10). Denote

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

Lemma 2.

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,

(20)

where the coefficient

E_{n}^{1}

satisfies the estimate

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

(21)

Proof.

From Taylor theorem

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

(22)

for

k = 1, \dots, n

and

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

. By substituting

f_{n - k}

with (22) in formula (20) for

A_{2}

we get

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 (14) and (15) we get

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

| E_{n} | \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)} .

□

Lemma 3.

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

(23)

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 Lemma 2 and Lemma 3 depends on the step size h and approximations

A_{2}

and

A_{3}

can be defined for any interval

[l, X]

. Asymptotic formula

A_{2}

is a first-order approximation of the first derivative for

n > ln N

while

A_{3}

is an approximation of the first derivative for all

n \geq 2

. Now we use the bound for the errors of approximations

A_{2}

and

A_{3}

to prove the convergence of numerical methods NS2 and NS3. Let

e_{k} = y_{k} - u_{k}

be the error of numerical solution NS3 at the point

x_{k}

for

k = 1, \dots, N

. First we estimate the errors

e_{1}

and

e_{2}

. From Taylor Theorem the solution of equation (10) satisfies

y_{1} = y_{0} + h F_{1} - \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 Theorem 4 and Theorem 5 we derive bounds for the errors of numerical solutions NS2 and NS3 of equation (10). The theorems are proved by induction.

Theorem 4.

The errors of numerical solution NS3 satisfies

| e_{n} | < \frac{M X h}{{(1 - a)}^{2}}, (1 \leq n \leq N) .

(24)

Proof.

The estimate (24) for the error of numerical method NS3 holds for

n = 1

and

n = 2

. Suppose that (24) holds for

n = 1, \dots, m - 1

. The solution of equation (10) and the error of NS3 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} .

(25)

By applying (25) successively

m - 3

times, we obtain

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

Denote

E = 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 NS2.

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

The error of numerical solution NS2 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) .

(26)

Proof.

The estimate (26) for the error of numerical method NS2 holds for

n = 1

and

n = 2

. Suppose that (26) holds for

n = 1, \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}

(27)

By applying (27) 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 6.

Consider the first-order linear ODE with constant coefficients

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

(28)

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

(NS4)

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 (28) 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 (28) 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}) .

(29)

The numerical solution (29) is calculated with

O (N)

arithmetic operations using an algorithm that is similar to the algorithms for calculating NS2 and NS3. From formula (29), the value of

u_{n}

is expressed in terms of

u_{n - 1}

as follows.

(1 + 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 (28) 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})) .

(NS5)

Consider the following first-order linear ODE

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

(30)

Experimental results for the error and order of numerical solutions NS4 and NS5 of equation (30) and positive values of L are are given in Table 3. Now we prove the convergence of numerical method NS5 when the value of L is positive. Denote by

e_{n} = y_{n} - u_{n}

the error of NS5 at the point

x_{n}

. First, we estimate the errors

e_{1}

and

e_{2}

. The solution of equation (28) satisfies

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

(31)

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

(32)

From 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 NS5 of equation (30) 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 7.

Let

L > 0

. Then

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

(33)

Proof.

Inequality (33) is satisfied for

n = 1

and

n = 2

. Assume that (33) 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 3:

| 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 (33) holds for all

n = 1, \dots, N

. □

We will compare the performance of numerical methods NS4 and NS5 and negative values of the parameter L of equation (30). Experimental results for the error and order of numerical solutions NS4 and NS5 and

L = - 20, - 30, - 50

are given in Tables 4-6. The results of the numerical experiments presented in this section refer to ODSs that are defined in the interval

[0, 1]

. The numerical results in Table 4 show that method NS4 has a first-order accuracy and its error can be very large for small values of the step size h of the method. In Table 5 the values of the parameter a are fixed numbers close to one. In Table 6 the parameter

a = 1 - h^{0.8}

and numerical solution NS5 has an order

0.2

. Although the order of NS5 in Tables 5-6 is smaller than one, the errors of the method are significantly smaller than the corresponding errors of numerical solution NS4 in Table 4.

Example 8.

Consider the first-order linear ODE

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

(34)

The numerical solutions NS6 and NS7 for equation (34) 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},

(NS6)

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

(35)

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

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

(NS7)

The first-order liner ODE

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

(36)

has a solution

y (x) = sin x - cos x

. The experimental results of numerical solutions NS6 and NS7 of equation (36) are given in Table 7.

Numerical solution NS7 depends on a parameter a. The method can have any order in the interval

(0, 2]

with a suitable choice of the parameter. The coefficient

C_{a}

in an asymptotic formula (3) has a value

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

. The order of NS7 is

1 + p

when the coefficient

C_{a}

is equal to

s h^{p}

\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 NS7 has an order

1 + p

. When

a = 1 - s h^{p}

the coefficient

C_{a} = (2 - s h^{p}) h^{- p} / s

and NS7 has an order

1 - p

. A summary of these results is given below.

• Numerical method NS7 has an order

1 + p

when

s > 0

and

- 1 < p \leq 1

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

(37)

Values of the parameter p greater than one,

p > 1

lead to a second-order method. The experimental results for

p = - 0.3, p = 0.7

and

p = 1

are given in Table 8.

•NS7 has an order

1 - p

, when

s > 0

and

0 < p < 1

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

(38)

The experimental results for

p = 0.2, p = 0.5

and

p = 0.7

are given in Table 9.

•NS7 has an order

1 + p

, when

s > 0

and

0 < p \leq 1

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

(39)

When

p > 1

, numerical solution NS7 has a second-order accuracy. The experimental results for

p = 0.25, p = 0.75

and

p = 1.5

are given in Table 10.

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

Example 9.

Consider the logistic equation

y^{'} (x) = L y (x) (1 - y (x)), y (0) = y_{0} .

(40)

Equation (40) has a solution

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

The explicit schemes of equation (40) 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},

(NS8)

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

(41)

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

(NS9)

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

5. Numerical Solutions of Heat Equation

The heat equation is a partial differential equation that models the distribution of temperature in a one-dimensional object. Heat conduction a classical problem in physics and engineering that can be generalized to higher dimensions and has a wide range of applications. The analytical and numerical solutions of the heat equation are extensively studied in [32,33,34,35,36]. In this section we construct a finite difference scheme for the heat equation which uses approximation

A_{3}

for the partial derivative in time and second order difference approximation for the partial derivative in space. The stability and convergence of the method is proved and its performance is compared with the standard method using first order difference approximation for the partial derivative in time.

Example 10.

Consider the following heat equation

\{\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}

(42)

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 by

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

. By approximating the partial derivative in time by first-order backward difference approximation and the partial derivative in space by second-order central difference approximation we obtain

\frac{u_{n}^{m} - u_{n}^{m - 1}}{τ} = D \frac{u_{n - 1}^{m} - 2 u_{n}^{m} + u_{n + 1}^{m}}{h^{2}} + F_{n}^{m} + C_{n}^{m} τ + K_{n}^{m} h^{2},

- η_{0} u_{n - 1}^{m} + (1 + 2 η_{0}) u_{n}^{m} - η_{0} u_{n + 1}^{m} = u_{n}^{m - 1} + τ F_{n}^{m} + τ (C_{n}^{m} τ + K_{n}^{m} h^{2}) .

(43)

where

η_{0} = D τ / h^{2}

and

C_{n}^{m} τ + K_{n}^{m} h^{2}

is the truncation error of the approximations of the partial derivatives at the point

(x_{n}, y_{m})

of the grid

J

. The coefficients

C_{n}^{m}

and

K_{n}^{m}

satisfy

| C_{n}^{m} | < M_{2} / 2, | K_{n}^{m} | < M_{4} / 12

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

The numerical solution

{U_{n}^{m}}

of heat equation (42) on layer m of the grid

J

is a solution of the system of linear equations

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

(44)

and has initial conditions

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

for

n = 0, \dots, N

and boundary conditions

U_{0}^{m} = u_{1} (m τ), U_{N}^{m} = u_{2} (m τ), (1 \leq m \leq M) .

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

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

(NS10)

where

U_{m}

and

F_{m}

are vectors of dimension

N - 1

which have entries

U_{n}^{m}

and

F_{n}^{m}

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 NS10 has an accuracy

O (τ + h^{2})

and it is computed with

O (M N)

operations because

I + η Δ

is a tridiagonal M-matrix [37,38]. Now we construct the finite difference scheme of equation (42) 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 (42) 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}

(45)

From Lemma 3 the coefficient

S_{n}^{m}

satisfies the estimate

| S_{n}^{m} | < \frac{(1 + a) M_{2}}{2 (1 - a)}

. The formula (45) 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}

(46)

From formulas (45) and (46) 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 (51) 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}

(47)

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 (42) on row m of the grid

J

which corresponds to (47) 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}

(48)

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

(I + η Δ) U_{m} = (I + a η Δ) U_{m - 1} + τ {\bar{F}}_{m} + η I_{m},

(NS11)

where

{\bar{F}}_{m}

is the vector of dimension

N - 1

with elements

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

. Consider the following heat equation

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

(49)

where

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

. Experimental results for error and order of numerical methods 44 and NS11 of equation (49) and

h = τ

are given in Table 12. The orders of the methods with respect to

τ

and h are computed with formulas

{log}_{2} | e_{τ} / e_{2 τ} |

and

{log}_{2} | (e_{h} - 2 e_{2 h}) / (e_{2 h} - 2 e_{4 h}) |

, where

e_{τ} = e_{h}

is he error of the method for the given values of

τ = h

in Table 12. The two methods, 44 and NS11, have similar performance with respect to computational time and accuracy. Numerical solution NS11 has an order

1 + p

when the parameter

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

. The numerical results for

p = - 0.5, 0.5, 1

and

h = 2 τ

are given in Table 13.

A Crank-Nicolson scheme for heat equation has second order accuracy [39,40]. The convergence of finite difference schemes for heat equation are studied in [41,42]. Now we derive a bound for the error of finite difference scheme NS11. Let

e_{n}^{m} = U_{n}^{m} - u_{n}^{m}

be the error of 59 at the point

(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 11 we estimate the errors of difference scheme NS11 on the first row of the grid

J

Lemma 11.

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}

(50)

Proof.

The solution of equation (42) 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}

(51)

for

n = 1, \dots, N - 1

and

e_{0}^{1} = e_{N}^{1} = 0

. The errors of numerical method NS11 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},

(52)

where

W_{1}

is the vector with elements

ε_{n}^{1}

. From (62) 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},

(53)

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

(54)

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 12 we prove that finite difference scheme NS11 is unconditionally stable and has an accuracy

O (τ + h^{2})

Theorem 12.

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 (54) 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 11 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}) .

□

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

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

(55)

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 Caputo derivative using the asymptotic formula of 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 - α}),

(56)

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 - α}),

(57)

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 13 we show that the coefficient

σ_{1}^{(α)}

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

b = α / 2

Lemma 13.

Let

b = a / 2

. Then

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

(58)

Proof.

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

Inequality (58) holds when

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

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

Taking the logarithm of both sides, we get

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 (58) 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 next examples, we will consider an application of approximation (57) for construction of numerical schemes for the two-term ordinary fractional differential equation and the fractional subdiffusion equation.

Example 14.

Consider the following two-term equation

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

(59)

The numerical solution

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

of equation (59) which uses L1 approximation (56) 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}) .

(NS12)

The numerical solution of equation (59) which uses approximations (56) and (57) 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})

(NS13)

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

(60)

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

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

(61)

Consider the two-term equation

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

(62)

Equation (62) 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 NS12 and NS13 of equation (62) on the interval

[0, 1]

are given in Table 14 and Table 15. The second column of Table 14 pertains to NS13 using approximation (57), while the results in the third column relate to NS13 using approximation (58) with a parameter

b = - 0.5

. The numerical results in Table 15 were obtained using approximation (5) for the fractional derivative. The sequence

s_{n}

and the coefficients

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

and

σ_{n}^{(α)}

of approximations (56) and (57) for

n = 3, \dots, N

are computed recursively in

O (N)

arithmetic operations. The NS12 and NS13 methods have similar computational time and accuracy. The approximation (56) is obtained from the formula (2) by replacing the first derivative

f_{n}^{'}

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

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}

(63)

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

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 (61) 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} .

(NS14)

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}

(64)

Equation (64) has a solution

u (x, t) = e^{x + t}

. Experimental results for the error and order of numerical solution NS14 of equation (64) are given in Table 16.

7. Conclusions

In this paper, we construct approximations of the first derivative whose weights contain powers of a parameter a, where

| a | < 1

. The approximations are applied for the numerical solution of ordinary differential equations and the heat equation. The properties of the weights allow the computation of the numerical solution of ordinary differential equations with

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 the time and N is the number of nodes in space. The numerical methods are compared with the Euler method using a first-order backward approximation of the first derivative. In the examples given in the paper, we observe similar computational time and accuracy for the numerical methods, with the difference in computational time being less than a second. The numerical methods for ODEs and PDEs using approximation (3) of the first derivative can achieve an arbitrary order in

(0, 2]

using values of the parameter (37), (38), and (39). The numerical methods constructed in the paper are easy to implement and have a better performance compared to the Euler method for equation (30) with negative values of L and values of the parameter a close to one, as shown in Tables 4-6. The main goal of constructing the approximations of the first derivative (3) and (4) is to apply them to the construction of approximations of the fractional derivative, whose weights have property (1). An example of a construction of an approximation of the fractional derivative of order

2 - α

which uses asymptotic formula (2) and approximation (3) of the first derivative is given in section 6. In future work, we will extend the method for constructing approximations of the second and higher-order derivatives. We will study constructions of high-order approximations of the fractional derivative that satisfy (1), and we will investigate the properties and convergence of the finite difference schemes for fractional differential equations.

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 conflict 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, IJIN 2020, 1, 2–58. [CrossRef]
Baba, I.A.; Rihan, F.A. A fractional–order model with different strains of COVID-19. Phys. A: Stat. Mech. Appl. 1278. [Google Scholar] [CrossRef]
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, 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(3), 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. SQUJS 2018, 23(1), 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. IJMSI 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(1), 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(6), 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; Large-Scale Scientific Computing. LSSC 2021; Lirkov, I., Margenov, S., Eds.; Springer, 2022; p. 13127. [Google Scholar]
Todorov, V.; Dimitrov, Y.; Dimov, I. Second order shifted approximations for the first derivative. In Studies in Computational Intelligence; Advances in High Performance Computing. HPC 2019; Dimov, I., Fidanova, S., Eds.; Springer: Cham, 2011; p. 902. [Google Scholar]
Dimitrov, Y. A second order approximation for the Caputo fractional derivative, J. Fract. Calc. Appl. 2016, 7(2), 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, 4 ed.; Springer, 2014. [Google Scholar]
Kot, M. Elements of Mathematical Ecology; Cambridge University Press, 2001. [Google Scholar]
Hahn, D.W.; Özişik, M.N. Heat Conduction, 3 ed.; John Wiley & Sons, 2012. [Google Scholar]
Incorpera, F.P.; DeWitt, D.P.; Bergman, T.L.; Lavine, A.S. Fundamentals of Heat and Mass Transfer, 6 ed.; John Wiley: New York, 2006. [Google Scholar]
Majumdar, P. Computational Methods for Heat and Mass Transfer; CRC Press: Boca Raton, 2006. [Google Scholar]
Morton, K.W.; Mayers, D.F. Numerical Solution of Partial Differential Equations: An Introduction, 2 ed.; Cambridge University Press, 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, 2007. [Google Scholar]
Datta, B.N. Numerical Linear Algebra and Applications, 2 ed.; Society for Industrial and Applied Mathematics: 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. IJRASET 2023, 11. [Google Scholar] [CrossRef]
Ames, W. Numerical Methods for Partial Differential Equations; Academic Press: Boston, USA, 1992. [Google Scholar]
Smith, G.D. Numerical Solution of Partial Differential Equations: Finite Difference Method; Oxford, 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(1), 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, 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(2), 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(2), 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(3), 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(2), 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, 2007; pp. 23–55. [Google Scholar]
Podlubny, I. Fractional Differential Equations; Academic Press: USA, 1999. [Google Scholar]

Table 1. Error and order of numerical methods NS1 and NS2 of equation (12).

$h$	NS1		NS2, $a = - 0.1$		NS2, $a = 0.9$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 methods NS3 of equation (19).

$h$	NS3, $a = 0.5$		NS3, $a = - 1$		NS3, $a = (h - 1) / (h + 1)$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS4 and NS5 of equation (30).

$h$	NS4, $L = 2$		NS5, $L = 2, a = - 0.9$		NS5, $L = 20, a = 0.5$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS4 of equation (30).

$h$	NS4, $L = - 20$		NS4, $L = - 30$		NS4, $L = - 50$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$1.56250 \times 10^{- 5}$	$200.1160$	$1.0045$	$2.8992 \times 10^{6}$	$1.0101$	$8.4295 \times 10^{14}$	$1.0282$
$7.81250 \times 10^{- 6}$	$99.9028$	$1.0022$	$1.4445 \times 10^{6}$	$1.0051$	$4.1738 \times 10^{14}$	$1.0141$
$3.90625 \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 NS5 of equation (30).

$h$	NS5, $L = - 20$		NS5, $L = - 30$		NS5, $L = - 50$
	$a = 0.99995$		$a = 0.99999$		$a = 0.99997$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$1.56250 \times 10^{- 5}$	$0.0333$	$0.5963$	$0.0453$	$0.1657$	$0.0179$	$0.4295$
$7.81250 \times 10^{- 6}$	$0.0192$	$0.7976$	$0.0366$	$0.3072$	$0.0113$	$0.6594$
$3.90625 \times 10^{- 6}$	$0.0103$	$0.8943$	$0.0255$	$0.5211$	$0.0064$	$0.8317$

Table 6. Error and order of numerical method NS5 of equation (30) and

a = 1 - h^{0.8}

Table 6. Error and order of numerical method NS5 of equation (30) and

a = 1 - h^{0.8}

$h$	NS5, $L = - 20$		NS5, $L = - 30$		NS5, $L = - 50$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$1.56250 \times 10^{- 5}$	$0.0140$	$0.1779$	$0.0092$	$0.1783$	$0.0054$	$0.1785$
$7.81250 \times 10^{- 6}$	$0.0124$	$0.1805$	$0.0081$	$0.1808$	$0.0048$	$0.1811$
$3.90625 \times 10^{- 6}$	$0.0109$	$0.1828$	$0.0071$	$0.1831$	$0.0042$	$0.1833$

Table 7. Error and order of numerical methods NS6 and NS7 of equation (47).

$h$	NS6		NS7, $a = - 0.3$		NS7, $a = - 0.9$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS7 of equation (40) and

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

Table 8. Error and order of numerical method NS7 of equation (40) and

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

$h$	NS7, $s = 1, p = - 0.3$		NS7, $s = 2, p = 0.7$		NS7, $s = 3, p = 1$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS7 of equation (40) and

a = 1 - s h^{p}

Table 9. Error and order of numerical method NS7 of equation (40) and

a = 1 - s h^{p}

$h$	NS7, $s = 1, p = 0.2$		NS7, $s = 2, p = 0.5$		NS7, $s = 3, p = 0.7$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS7 of equation (40) and

a = - 1 + s h^{p}

Table 10. Error and order of numerical method NS7 of equation (40) and

a = - 1 + s h^{p}

$h$	NS7, $s = 1, p = 0.25$		NS7, $s = 2, p = 0.75$		NS7, $s = 3, p = 1.5$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS8 and NS9 of equation (40) and

L = 1, y_{0} = 2

Table 11. Error and order of numerical methods NS8 and NS9 of equation (40) and

L = 1, y_{0} = 2

$h$	NS8		NS9, $a = 0.3$		NS9, $a = 0.5$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 solutions NS10 and NS11 of equation (49) and

τ = h

Table 12. Error and order of numerical solutions NS10 and NS11 of equation (49) and

τ = h

$h$	NS10			NS11, $α = - 0.5$
	$E r r o r$	$O r d e r (τ)$	$O r d e r (h)$	$E r r o r$	$O r d e r (τ)$	$O r d e r (h)$
$0.02$	$1.8625 \times 10^{- 3}$	$-$	$-$	$6.4127 \times 10^{- 4}$	$-$	$-$
$0.01$	$9.2957 \times 10^{- 4}$	$1.0026$	$-$	$3.1497 \times 10^{- 4}$	$1.0257$	$-$
$0.005$	$4.6436 \times 10^{- 4}$	$1.0013$	$2.0040$	$1.5606 \times 10^{- 4}$	$1.0131$	$1.9979$
$0.0025$	$2.3208 \times 10^{- 4}$	$1.0006$	$2.0644$	$7.7680 \times 10^{- 5}$	$1.0065$	$2.0049$
$0.00125$	$1.1601 \times 10^{- 4}$	$1.0003$	$1.9244$	$3.8751 \times 10^{- 5}$	$1.0033$	$1.9929$

Table 13. Error and order of numerical solutions NS11 of equation (49) and

τ = h / 2

Table 13. Error and order of numerical solutions NS11 of equation (49) and

τ = h / 2

$h$	$a = (1 - \sqrt{τ}) / (1 + \sqrt{τ})$		$a = (\sqrt{τ} - 1) / (\sqrt{τ} + 1)$		$a = (τ - 1) / (τ + 1)$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$0.02$	$8.8706 \times 10^{- 3}$	$-$	$1.1284 \times 10^{- 4}$	$-$	$2.9381 \times 10^{- 5}$	$-$
$0.01$	$6.3492 \times 10^{- 3}$	$0.4824$	$3.7825 \times 10^{- 5}$	$1.5769$	$7.3456 \times 10^{- 6}$	$1.9999$
$0.005$	$4.5312 \times 10^{- 3}$	$0.4867$	$1.2854 \times 10^{- 5}$	$1.5571$	$1.8364 \times 10^{- 6}$	$2.0000$
$0.0025$	$3.2258 \times 10^{- 3}$	$0.4902$	$4.4148 \times 10^{- 6}$	$1.5418$	$4.5912 \times 10^{- 7}$	$2.0000$
$0.00125$	$2.2921 \times 10^{- 3}$	$0.4929$	$1.5283 \times 10^{- 6}$	$1.5304$	$1.1477 \times 10^{- 7}$	$2.0001$

Table 14. Error and order of numerical solutions NS12 and NS13 of equation (62).

$h$	NS12, $α = 0.25, L = 1$		NS13, $α = 0.5, L = 10$		$α = 0.75, L = 5, b = - 0.5$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 NS13 of equation (62) and

b = α / 2

Table 15. Error and order of numerical solution NS13 of equation (62) and

b = α / 2

$h$	$α = 0.25, L = 1$		$α = 0.5, L = 10$		$α = 0.75, L = 5$
	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$	$E r r o r$	$O r d e r$
$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 solutions NS14 of equation (64).

$h = τ$	$α = 0.5, D = 1$			$α = 0.75, D = 2$
	$E r r o r$	$O r d e r (τ)$	$O r d e r (h)$	$E r r o r$	$O r d e r (τ)$	$O r d e r (h)$
$0.04$	$2.2260 \times 10^{- 3}$	$-$	$-$	$4.4963 \times 10^{- 3}$	$-$	$-$
$0.02$	$8.0575 \times 10^{- 4}$	$1.4661$	$-$	$1.9274 \times 10^{- 4}$	$1.2221$	$-$
$0.01$	$2.8873 \times 10^{- 4}$	$1.4806$	$2.2819$	$8.1797 \times 10^{- 4}$	$1.2365$	$2.2811$
$0.005$	$1.0292 \times 10^{- 4}$	$1.4882$	$2.1992$	$3.4556 \times 10^{- 4}$	$1.2431$	$2.2092$
$0.0025$	$3.6571 \times 10^{- 5}$	$1.4927$	$2.1887$	$1.4564 \times 10^{- 4}$	$1.2465$	$2.2361$

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 (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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 Solutions of Heat Equation

6. Approximations of the Fractional Derivative

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe