Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations

Preprint

Article

Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations

Altmetrics

Downloads

392

Views

Comments

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

04 September 2023

Posted:

06 September 2023

You are already at the latest version

Alerts

Abstract

The fractional derivatives are a generalization the derivatives of integer order and find applications in studying memory processes in various scientific fields. Numerical methods are used to solve and analyze fractional models of real world problems. In this paper, we consider an approximation of the Caputo fractional derivative and its asymptotic formula, whose generating function is the polylogarithmic function. In the paper, we prove the convergence of the approximation and derive an estimate for the error and order. We consider an application of the approximation for the construction of finite difference schemes for numerical solution of the two-term ordinary fractional differential equation and the time-fractional Black-Scholes equation for option pricing. The properties of the approximation are used to prove the convergence of the methods used for numerical solution of the fractional differential equations. The theoretical results for the error and order of the methods are illustrated by the results of numerical experiments.

Keywords:

Subject: Computer Science and Mathematics - Computational Mathematics

1. Introduction

Fractional differential equations are an essential means of modeling complex processes in different fields of science [3,10,35,50,51,53]. The Riemann-Liouville fractional integral of order

α > 0

is defined as

_{a} I_{x}^{α} f (x) = \frac{1}{Γ (α)} \int_{a}^{x} {(x - t)}^{α - 1} f (t) d t .

When n is a positive integer and

0 < α < 1

, the Caputo and Riemann-Liouville fractional derivatives of order

α

are defined as

\begin{matrix} _{a}^{C} D_{x}^{n + α} f (x) & = I^{1 - α} D^{n + 1} f (x) = \frac{1}{Γ (1 - α)} \int_{a}^{x} \frac{f^{(n + 1)} (t)}{{(x - t)}^{α}} d t, \\ _{a}^{R L} D_{x}^{n + α} f (x) & = D^{n + 1} I^{1 - α} f (x) = \frac{1}{Γ (1 - α)} \frac{d^{n + 1}}{d x^{n + 1}} \int_{a}^{x} \frac{f (t)}{{(x - t)}^{α}} d t . \end{matrix}

The Caputo and Riemann-Liouville fractional derivatives are related as

_{a}^{R L} D_{x}^{n + α} f (x) =_{a}^{C} D_{x}^{n + α} f (x) + \sum_{k = 0}^{n} \frac{f^{(k)} (0) x^{k - n - α}}{Γ (k - α - n + 1)} .

Without loss of generality the lower limit of the integral in the definition of fractional derivative is zero. When

0 < α < 1

the Caputo derivative is defined as

f^{(α)} (x) = D_{x}^{α} f (x) = \frac{1}{Γ (1 - α)} \int_{0}^{x} \frac{f^{'} (t)}{{(x - t)}^{α}} d t .

Fractional derivatives and integrals are nonlocal and have a singularity at the endpoint. The properties of the fractional derivatives make them an important tool for describing memory processes. The exponential, sine and cosine functions have Caputo derivatives,

D_{x}^{α} e^{x} = x^{1 - α} E_{1, 2 - α} (x), D_{x}^{α} sin x = x^{1 - α} E_{2, 2 - α} (- x^{2}), D_{x}^{α} cos x = - x^{2 - α} E_{2, 3 - α} (- x^{2}),

where

E_{β_{1}, β_{2}} (x)

is the Mittag-Lefler function

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

Finite difference schemes for numerical solution of fractional differential equations use discretizations of the fractional derivative. Two important discretizations of the Caputo fractional derivative are the Grünwald difference approximation and L1 approximation [8,26,30,49]. Grünwald difference approximation has a first order accuracy and a generating function

{(1 - x)}^{α}

. Central difference approximations of integer order derivatives have a second order accuracy. Grünwald difference approximation is a second order shifted approximation of the fractional derivative with a shift parameter

α / 2

, and it is a generalization of finite difference approximations [37,58]. L1 approximation has an order

2 - α

and a generating function

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

. In Apostolov [7] and Dimitrov [16] we use the L1 approximation and approximations of the second derivative for construction of a second order approximations of the fractional derivative. Alikhanov [4] and Gao et al. [22] construct approximations of the fractional derivative of order

3 - α

, which are called

L 1 - 2

and

L 1 - 2_{σ}

formulas. Finite difference schemes using L2 formula approximations of the fractional derivative of order

3 - α

, and their convergence are studied by Alikhanov [5], Lv and Xu [33], Wong and Ren [55]. High order approximations of the fractional derivative and finite difference schemes for fractional differential equations are constructed in [8,9,12,14,29,32,44,48,57]. Navot [38] uses Taylor polynomials to deal with the singularity of the fractional integral

\int_{0}^{1} t^{α} f (t) d t

. He derives the asymptotic formula for Riemann sums on a uniform net, called extended Euler-Maclaurin summation formula. The formula is extended in Navot [39] to functions with a logarithmic singularity. In [18] we derive the asymptotic formula of the Riemann sum of the fractional integral using the series expansion of the generating function

{Li}_{1 - α} (e^{- x})

. Other types of numerical methods for fractional differential equations include fractional multistep methods and methods using spline interpolation and wavelets [1,20,31,34,46,47].

Let n be a positive integer,

h = x / n

and

y_{k} = y (k h)

for

k = 0, 1, \dots, n

. In [17] we derive an approximation of Caputo fractional derivative and its asymptotic formula

\begin{matrix} \frac{1}{h^{α}} \sum_{k = 1}^{n - 1} \frac{f_{n - k}}{k^{1 + α}} = Γ (- α) f_{n}^{(α)} + \frac{ζ (1 + α)}{h^{α}} f_{n} - ζ (α) f_{n}^{'} h^{1 - α} + O (h^{2 - α}) . \end{matrix}

(1)

Approximation (1) has a generating function

{Li}_{1 + α} (x)

, the polylogarithmic function of order

1 + α

. By substituting the derivative in the right-hand side of (1) using first order backward difference

(f_{n} - f_{n - 1}) / h

we find an approximation of the fractional derivative

\begin{matrix} L_{n} f (x) = \frac{h^{- α}}{Γ (- α)} \sum_{k = 0}^{m} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + O (h^{2 - α}), \end{matrix}

(2)

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

Approximations (1) and (2) have an order

2 - α

when the function f satisfies the condition

f (0) = f^{'} (0) = 0

. In paper [17] we extend approximation (2) to all functions in the class

C^{2} [0, x]

by assigning values of the last two weights.

\begin{matrix} A_{n} f (x) = \frac{h^{- α}}{Γ (- α)} \sum_{k = 0}^{n} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + O (h^{2 - α}), \end{matrix}

(3)

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

(4)

\begin{matrix} σ_{n - 1}^{(α)} = \frac{1}{{(n - 1)}^{1 + α}} - n S_{n} [1 + α] + S_{n} [α] - \frac{n^{1 - α}}{α (1 - α)}, \end{matrix}

(5)

\begin{matrix} σ_{n}^{(α)} = (n - 1) S_{n} [1 + α] - S_{n} [α] + \frac{n^{1 - α}}{α (1 - α)}, \end{matrix}

(6)

where

S_{n} [α] = \sum_{k = 1}^{n - 1} k^{- α} - ζ (α)

. In this paper we investigate the properties of approximations (2) and (3) and applications of the approximations for numerical solution of ordinary and partial fractional differential equations. The outline of the paper is as follows. In section 2 we study the properties of the weights of approximation (3) and we derive the inequality

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} < ζ (1 + α) n^{α} - 2 . \end{matrix}

(7)

In section 3 we derive estimates for the errors of approximations (2) and (3). In sections 4 and 5 we construct numerical solutions of the two-term ordinary fractional differential equation and the time-fractional Black-Scholes equation for option pricing and we prove the convergence and order of the methods.

2. Properties of the Weights

L1 approximation and approximation (2) have an order

2 - α

and similar performance and properties of the weights [17]

σ_{0}^{(α)} > 0, σ_{1}^{(α)} < σ_{2}^{(α)} < \dots < σ_{n - 1}^{(α)} < 0, \sum_{k = 0}^{n} σ_{k}^{(α)} = 0 .

The weights of

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

and

σ_{n}^{(α)}

of approximation (2) are chosen such that the values of

A_{n} 1

and

A_{n} x

are equal to the fractional derivatives.

A_{n} 1 = D_{x}^{α} 1 = 0, A_{n} x = D_{x}^{α} x = \frac{x^{1 - α}}{Γ (2 - α)} .

In this section we prove inequality (7) and we derive an estimate for the weight

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

of approximation (2). The proof uses the inequalities in Claim 1 and Claim 2.

Claim 1.

Let

0 < α < 1

and

1 \leq k \leq n

. Then

\begin{matrix} {(n - k)}^{α} < n^{α} - \frac{α k}{n^{1 - α}} . \end{matrix}

(8)

Proof.

From the binomial formula

\begin{matrix} {(n - k)}^{α} = n^{α} {(1 - \frac{k}{n})}^{α} = n^{α} \sum_{i = 0}^{\infty} {(- 1)}^{i} (\binom{α}{k}) {(\frac{k}{n})}^{i} . \end{matrix}

The numbers

{(- 1)}^{i} (\binom{α}{k})

are negative for

k \geq 1

. Then

{(n - k)}^{α} - n^{α} + \frac{α k}{n^{1 - α}} = n^{α} \sum_{i = 2}^{\infty} {(- 1)}^{i} (\binom{α}{k}) {(\frac{k}{n})}^{i} < 0 .

□

Claim 2.

Let

0 < α < 1

. Then

\begin{matrix} \frac{1 - α + α^{2}}{α (1 - α)} + \frac{α ζ (α)}{2^{1 - α}} > 2 . \end{matrix}

(9)

Proof.

The Riemann zeta function has a series expansion [36]

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

(10)

where

γ = 0.5772

is Euler’s constant and

γ_{k}

are Stieltjes constants

γ_{1} = - 0.0728158, γ_{2} = - 0.0096904, γ_{3} = 0.0020538 .

Then

\begin{matrix} ζ (α) & > \frac{1}{α - 1} + γ + (1 - α) γ_{1} + \frac{{(1 - α)}^{2}}{2} γ_{2} > \frac{1}{α - 1} + γ + (1 - α) (γ_{1} + \frac{γ_{2}}{2}) \\ > \frac{1}{α - 1} + γ - (1 - α) 0.08 > \frac{1}{α - 1} + 0.49 + 0.08 α . \end{matrix}

(11)

The value of

ln 2 = 0.693147 < 0.7

. From Taylor’s Theorem

\begin{matrix} 2^{- (1 - α)} = e^{- (1 - α) ln 2} & < 1 - (1 - α) ln 2 + {(1 - α)}^{2} \frac{{ln}^{2} 2}{2} \\ < 1 - (ln 2 - \frac{{ln}^{2} 2}{2}) (1 - α) \\ < 1 - 0.45 (1 - α) = 0.55 + 0.45 α . \end{matrix}

(12)

From (11) and (12)

\begin{matrix} \frac{1 - α + α^{2}}{α (1 - α)} + \frac{α ζ (α)}{2^{1 - α}} & > \frac{1 - α + α^{2}}{α (1 - α)} + α (\frac{1}{α - 1} + 0.497 + 0.08 α) (0.55 + 0.45 α) \\ > \frac{1}{α} + 0.036 a^{3} + 0.26765 a^{2} + 0.72335 α . \end{matrix}

Inequality (9) holds when

\frac{1}{α} + 0.036 α^{3} + 0.26765 α^{2} + 0.72335 α > 2,

0.036 a^{4} + 0.26765 a^{3} + 0.72335 α^{2} - 2 α + 1 > 0 .

The function

f (α) = 0.036 a^{4} + 0.26765 a^{3} + 0.72335 α^{2} - 2 α + 1

has a first derivative

f^{'} (α) = 0.144 α^{3} + 0.80295 α^{2} + 1.4467 α - 2

and a minimum value

f (0.882) = 0.00414

. Therefore

f (α)

is positive on

[0, 1]

. □

The asymptotic formula of

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

is obtained from the asymptotic formula of the fractional integral [18,38]

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} = Γ (1 + α) Γ (- α) & + \sum_{k = 0}^{\infty} {(- 1)}^{k} (\binom{α}{k}) \frac{ζ (1 + α - k)}{n^{k - α}} \\ + \sum_{k = 0}^{\infty} {(- 1)}^{k} (\binom{- 1 - α}{k}) \frac{ζ (- α - k)}{n^{α + k + 1}} . \end{matrix}

(13)

Let

β = min {2 + α, 3 - α}

. From the properties of gamma function

Γ (1 + α) Γ (- α) = - \frac{π}{sin (π α)} .

From (13) we find the asymptotic formula of

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

of order

β

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} = ζ (1 + α) n^{α} - \frac{π}{sin (π α)} - \frac{α ζ (α)}{n^{1 - α}} + \frac{ζ (- α)}{n^{1 + α}} - α (1 - α) \frac{ζ (α - 1)}{2 n^{2 - α}} + O (n^{- β}) . \end{matrix}

Lemma 3.

Let

n \geq 2

. Then

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} < ζ (1 + α) n^{α} - 2 . \end{matrix}

(14)

Proof.

From inequality (8)

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} < \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} (n^{α} - \frac{α k}{n^{1 - α}}) < n^{α} \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} - \frac{α}{n^{1 - α}} \sum_{k = 1}^{n - 1} \frac{1}{k^{α}} . \end{matrix}

Sum of the

α

-th powers of the first

n - 1

integers has an asymptotic formula [21]

\begin{matrix} \sum_{k = 1}^{n - 1} k^{α} = ζ (- α) + \frac{n^{1 + α}}{1 + α} \sum_{p = 1}^{\infty} (\binom{1 + α}{p}) \frac{b_{p}}{n^{p}}, \end{matrix}

(15)

where

b_{p}

are the Bernoulli numbers. Formula (15) is derived using Euler-Maclaurin formula and the actual value of

\sum_{k = 1}^{n - 1} k^{α}

lies between two consecutive partial sums [21]. Then

\sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} < ζ (1 + α) - \frac{1}{α n^{α}} - \frac{1}{2 n^{1 + α}} - \frac{1 + α}{12 n^{2 + α}} + \frac{(1 + α) (2 + α) (3 + α)}{720 n^{4 + α}},

\sum_{k = 1}^{n - 1} \frac{1}{k^{α}} > \frac{n^{1 - α}}{1 - α} + ζ (α) - \frac{1}{2 n^{α}} - \frac{α}{12 n^{1 + α}}

and

\begin{matrix} n^{α} \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} < ζ (1 + α) n^{α} - \frac{1}{α} - \frac{1}{2 n} - \frac{1 + α}{12 n^{2}} + \frac{1}{30 n^{4}}, \\ \frac{α}{n^{1 - α}} \sum_{k = 1}^{n - 1} \frac{1}{k^{α}} > \frac{α}{1 - α} + \frac{α ζ (α)}{n^{1 - α}} - \frac{α}{2 n} - \frac{α^{2}}{12 n^{2}} . \end{matrix}

Therefore

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} & < n^{α} \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} - \frac{α}{n^{1 - α}} \sum_{k = 1}^{n - 1} \frac{1}{k^{α}} \\ < ζ (1 & + α) n^{α} - \frac{1 - α + α^{2}}{α (1 - α)} - \frac{α ζ (α)}{n^{1 - α}} - \frac{1 - α}{2 n} - \frac{1 + α - α^{2}}{12 n^{2}} + \frac{1}{30 n^{4}} \\ < ζ (1 & + α) n^{α} - \frac{1 - α + α^{2}}{α (1 - α)} - \frac{α ζ (α)}{n^{1 - α}} - \frac{1}{12 n^{2}} + \frac{1}{30 n^{4}} \\ < ζ (1 & + α) n^{α} - \frac{1 - α + α^{2}}{α (1 - α)} - \frac{α ζ (α)}{n^{1 - α}} . \end{matrix}

From Claim 2 it follows that

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{{(n - k)}^{α}}{k^{1 + α}} < ζ (1 + α) n^{α} - 2 . \end{matrix}

□

Now we derive an estimate for the weight

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

of approximation (3). Denote

{\bar{σ}}_{n - 1}^{(α)} = σ_{n - 1}^{(α)} - \frac{1}{{(n - 1)}^{1 + α}} = S_{n} [α] - n S_{n} [1 + α] - \frac{n^{1 - α}}{α (1 - α)} .

Claim 4.

Let

n \geq 2

. Then

|{\bar{σ}}_{n - 1}^{(α)}| < \frac{1 + 2 α}{12 n^{1 + α}} .

Proof.

From (15)

\begin{matrix} S_{n} [α] = \frac{n^{1 - α}}{1 - α} - \frac{1}{2 n^{α}} + R_{1}, n S_{n} [α + 1] = - \frac{n^{1 - α}}{α} - \frac{1}{2 n^{α}} + R_{2}, \end{matrix}

where

| R_{1} | < \frac{1 + α}{12 n^{1 + α}}, | R_{2} | < \frac{α}{12 n^{1 + α}} .

Then

\begin{matrix} |{\bar{σ}}_{n - 1}^{(α)}| & = |S_{n} [α] - n S_{n} [1 + α] - \frac{n^{1 - α}}{α (1 - α)}| \\ < | R_{1} | + | R_{2} | < \frac{1 + 2 α}{12 n^{1 + α}} \end{matrix}

□

3. Estimate for the Error

In this section we use the method from [38] to derive estimates for the errors of approximations (2) and (3). Let

0 < t < x

and

\begin{matrix} g (t) = f (t) - f (x) + (x - t) f^{'} (x) - \frac{{(x - t)}^{2}}{2} f^{''} (x) . \end{matrix}

(16)

The function

g (t)

satisfies

g (x) = g^{'} (x) = g^{''} (x) = 0

. From Taylor’s Theorem

g (t) = \frac{{(t - x)}^{3}}{6} f^{'''} (ξ_{0, t}), g^{'} (t) = \frac{{(t - x)}^{2}}{2} f^{'''} (ξ_{1, t}), g^{''} (t) = (t - x) f^{'''} (ξ_{2, t})

where

t < ξ_{i, t} < x

for

i = 0, 1, 2

. Let

M_{i} = max_{0 < t < x} | f^{(i)} (t) |, (i = 0, 1, 2, 3) .

Then

| g (t) | < \frac{M_{3}}{6} {(x - t)}^{3}, | g^{'} (t) | < \frac{M_{3}}{2} {(x - t)}^{2}, | g^{''} (t) | < M_{3} (x - t) .

Denote

G (t) = g (t) / {(x - t)}^{1 + α}

. The function

G (t)

satisfies

G (x) = G^{'} (x) = 0

and its second derivative of is not defined at the point x.

Claim 5.

Let

0 < t < x

. Then

| G^{'} (t) | < \frac{5}{6} M_{3} {(x - t)}^{1 - α}, | G^{''} (t) | < \frac{4 M_{3}}{{(x - t)}^{α}} .

Proof.

G^{'} (t) = \frac{(x - t) g^{'} (x) + (1 + α) g (t)}{{(x - t)}^{2 + α}},

G^{''} (t) = \frac{{(x - t)}^{2} g^{''} (t) + 2 (α + 1) (x - t) g^{'} (t) + (α + 1) (α + 2) g (t)}{{(x - t)}^{3 + α}} .

Therefore

\begin{matrix} | G^{'} (t) | & < \frac{1}{2} M_{3} {(x - t)}^{1 - α} + \frac{1 + α}{6} M_{3} {(x - t)}^{1 - α} < \frac{5}{6} M_{3} {(x - t)}^{1 - α}, \\ | G^{''} (t) | & < ((2 + α) (1 + α) \frac{M_{3}}{6} + (1 + α) M_{3} + M_{3}) {(x - t)}^{- α} < \frac{4 M_{3}}{{(x - t)}^{α}} . \end{matrix}

□

Lemma 6.

Let

f (0) = f^{'} (0) = 0

. Then

\int_{0}^{x} G (t) d t = Γ (- α) f^{(α)} (x) + \frac{f (x)}{α x^{α}} + \frac{f^{'} (x) x^{1 - α}}{1 - α} - \frac{f^{''} (x) x^{2 - α}}{2 (2 - α)} .

Proof.

\begin{matrix} \int_{0}^{x} G (t) d t & = \int_{0}^{x} \frac{g (t)}{{(x - t)}^{1 + α}} d t = \int_{0}^{x} g (t) d \frac{{(x - t)}^{- α}}{α} \\ = g (t) {\frac{{(x - t)}^{- α}}{α}|}_{0}^{x} - \frac{1}{a} \int_{0}^{x} \frac{g^{'} (t)}{{(x - t)}^{α}} d t \\ = - \frac{g (0)}{α x^{α}} + \frac{Γ (- α)}{Γ (1 - α)} \int_{0}^{x} \frac{g^{'} (t)}{{(x - t)}^{α}} d t . \end{matrix}

The function

g (t)

satisfies

g (0) = - f (x) + f^{'} (x) x - \frac{f^{''} (x)}{2} x^{2} .

Hence

\begin{matrix} \int_{0}^{x} \frac{g (t)}{{(x - t)}^{1 + α}} d t = Γ (- α) g^{(α)} (x) + \frac{1}{α} (\frac{f (x)}{x^{α}} - x^{1 - α} f^{'} (x) + \frac{x^{2 - α}}{2} f^{''} (x)) . \end{matrix}

(17)

The fractional derivative of the function

g (x)

satisfies

\begin{matrix} g^{(α)} (x) & = f^{(α)} (x) - \frac{f^{'} (x)}{Γ (1 - α)} \int_{0}^{x} {(x - t)}^{- α} d t + \frac{f^{''} (x)}{2 Γ (1 - α)} \int_{0}^{x} 2 {(x - t)}^{1 - α} d t \\ = f^{(α)} (x) - \frac{x^{1 - α} f^{'} (x)}{(1 - α) Γ (1 - α)} + \frac{x^{2 - α} f^{''} (x)}{(2 - α) Γ (1 - α)}, \end{matrix}

\begin{matrix} Γ (- α) g^{(α)} (x) = Γ (- α) f^{(α)} (x) + \frac{x^{1 - α} f^{'} (x)}{α (1 - α)} - \frac{x^{2 - α} f^{''} (x)}{α (2 - α)} . \end{matrix}

(18)

From (17) and (18) it follows that

\int_{0}^{x} \frac{g (t)}{{(x - t)}^{1 + α}} d t = Γ (- α) f^{(α)} (x) + \frac{f (x)}{α x^{α}} + \frac{f^{'} (x) x^{1 - α}}{1 - α} - \frac{f^{''} (x) x^{2 - α}}{2 (2 - α)} .

□

The trapezoidal rule of a function

F (t)

has a second order accuracy when

F \in C^{2} [0, x]

. The error of the trapezoidal rule

E_{T}^{F}

satisfies [56]

| E_{T}^{F} | < \frac{1}{12} h^{2} x^{3} max_{0 < t < x} | F^{''} (t) | .

The second derivative of the function

G (t)

is undefined at the point x, which leads to a lower order of accuracy of the trapezoidal rule in the interval

[0, x]

. Now we estimate the error of the trapezoidal rule of the function

G (t)

. Let

h = x / n

and

T_{n}^{1}

be the error of trapezoidal rule of

G (t)

in the interval

[0, x - h]

\begin{matrix} h (\frac{G (0)}{2} + \sum_{k = 1}^{n - 2} G (k h) + \frac{G (x - h)}{2}) = \int_{0}^{x - h} G (t) d t + T_{n}^{1} . \end{matrix}

(19)

Claim 7.

| T_{n}^{1} | < \frac{1}{3} M_{3} x^{3} h^{2 - α} .

Proof.

The error

T_{n}^{1}

of trapezoidal rule (19) satisfies

\begin{matrix} |T_{n}^{1}| & \leq \frac{1}{12} h^{2} {(x - h)}^{3} max_{0 < t < x - h} | G^{''} (t) | \\ \leq \frac{1}{12} h^{2} x^{3} \frac{4 M_{3}}{h^{α}} \leq \frac{1}{3} M_{3} x^{3} h^{2 - α} . \end{matrix}

□

Let

T_{n}^{2}

be the error of the trapezoidal rule of

G (t)

in the interval

[x - h, x]

\begin{matrix} \frac{h G (x - h)}{2} = \int_{x - h}^{x} G (t) d t + T_{n}^{2} . \end{matrix}

(20)

Claim 8.

| T_{n}^{2} | < \frac{5}{24} M_{3} h^{3 - α} .

Proof.

\int_{x - h}^{x} G (t) d t = \int_{x - h}^{x} G (t) d (t + \frac{h}{2} - x) = \frac{h}{2} G (x - h) - \int_{x - h}^{x} (t + \frac{h}{2} - x) G^{'} (t) d t .

The error

T_{n}^{2}

satisfies

T_{n}^{2} = \int_{x - h}^{x} (t + \frac{h}{2} - x) G^{'} (t) d t,

| T_{n}^{2} | < \int_{x - h}^{x} |t + \frac{h}{2} - x| | G^{'} (t) | d t < \frac{5}{6} M_{3} h^{1 - α} \int_{x - h}^{x} |t + \frac{h}{2} - x| d t < \frac{5}{24} M_{3} h^{3 - α} .

□

Let

T_{n} = T_{n}^{1} + T_{n}^{2}

be the error of the trapezoidal rule of

G (t)

in the interval

[0, x]

Lemma 9.

| T_{n} | \leq \frac{1}{3} M_{3} x (x^{2} + 1) h^{2 - α} .

Proof.

From Claim 7 and Claim 8

\begin{matrix} | T_{n} | \leq | T_{n}^{1} | + | T_{n}^{2} | & < \frac{1}{3} M_{3} x^{3} h^{2 - α} + \frac{5}{24} M_{3} h^{3 - α} \\ < \frac{1}{3} M_{3} (x^{3} + h) h^{2 - α} < \frac{1}{3} M_{3} x (x^{2} + 1) h^{2 - α} . \end{matrix}

□

Let

E_{n}

be the error of approximation (1).

\frac{1}{h^{α}} \sum_{k = 1}^{n - 1} \frac{f_{n - k}}{k^{1 + α}} = Γ (- α) f_{n}^{(α)} + \frac{ζ (1 + α)}{h^{α}} f_{n} - ζ (α) f_{n}^{'} h^{1 - α} + E_{n} .

Theorem 10.

Let

f (0) = f^{'} (0) = 0

. Then

| E_{n} | < \frac{2 M_{3} x (x^{2} + 1) + M_{2}}{6} h^{2 - α} .

Proof.

The trapezoidal rule of the function

G (t)

on the interval

[0, x]

satisfies

\frac{1}{h^{α}} (\sum_{k = 1}^{n - 1} \frac{g_{n - k}}{k^{1 + α}} + \frac{g_{0}}{2 n^{1 + α}}) = \int_{0}^{x} \frac{g (t)}{{(x - t)}^{1 + α}} d t + T_{n}

and the function

g (t)

has a value at zero

\begin{matrix} \frac{g (0)}{2 n^{1 + α} h^{α}} & = - \frac{f (x)}{2 n^{1 + α} h^{α}} + \frac{f^{'} (x) x}{2 n^{1 + α} h^{α}} - \frac{f^{''} (x) x^{2}}{4 n^{1 + α} h^{α}} = \frac{h}{2} (- \frac{f (x)}{x^{1 + α}} + \frac{f^{'} (x)}{x^{α}} - \frac{f^{''} (x) x^{1 - α}}{2}) . \end{matrix}

Then

\begin{matrix} \frac{1}{h^{α}} \sum_{k = 1}^{n - 1} \frac{g_{n - k}}{k^{1 + α}} = \int_{0}^{x} \frac{g (t)}{{(x - t)}^{1 + α}} d & t + (\frac{f (x)}{x^{1 + α}} - \frac{f^{'} (x)}{x^{α}} + \frac{f^{''} (x) x^{1 - α}}{2}) \frac{h}{2} + T_{n} \\ = Γ (- α) f^{(α)} (x) + & \frac{f (x)}{α x^{α}} + \frac{f (x) h}{2 x^{1 + α}} + \frac{x^{1 - α} f^{'} (x)}{1 - α} - \frac{f^{'} (x) h}{2 x^{α}} \\ + & \frac{x^{1 - α} f^{''} (x) h}{4} - \frac{x^{2 - α} f^{''} (x)}{2 (2 - α)} + T_{n} . \end{matrix}

From the definition (16) of the function

g (t)

\frac{1}{h^{α}} \sum_{k = 1}^{n - 1} \frac{g_{n - k}}{k^{1 + α}} = \frac{1}{h^{α}} \sum_{k = 1}^{n - 1} \frac{f_{n - k}}{k^{1 + α}} - \frac{f (x)}{h^{α}} \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} + f^{'} (x) h^{1 - α} \sum_{k = 1}^{n - 1} \frac{1}{k^{α}} - \frac{f^{''} (x)}{2} h^{2 - α} \sum_{k = 1}^{n - 1} k^{1 - α} .

From (15)

\begin{matrix} \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} = ζ (1 + α) - \frac{1}{α n^{α}} - \frac{1}{2 n^{1 + α}} + R_{1}, \\ \sum_{k = 1}^{n - 1} \frac{1}{k^{α}} = \frac{n^{1 - α}}{1 - α} + ζ (α) - \frac{1}{2 n^{α}} + R_{2}, \\ \sum_{k = 1}^{n - 1} k^{1 - α} = \frac{n^{2 - α}}{2 - α} - \frac{n^{1 - α}}{2} + ζ (α - 1) + R_{3}, \end{matrix}

where

| R_{1} | < \frac{1 + α}{12 n^{2 + α}}, | R_{2} | < \frac{α}{12 n^{1 + α}}, | R_{3} | < \frac{1 - α}{12 n^{α}} .

Then

\begin{matrix} \frac{f (x)}{h^{α}} \sum_{k = 1}^{n - 1} \frac{1}{k^{1 + α}} & = ζ (1 + α) \frac{f (x)}{h^{α}} - \frac{f (x)}{α n^{α} h^{α}} - \frac{f (x)}{2 n^{1 + α} h^{α}} + R_{1} \frac{f (x)}{h^{α}} \\ = ζ (1 + α) \frac{f (x)}{h^{α}} - \frac{f (x)}{α x^{α}} - \frac{f (x) h}{2 x^{1 + α}} + R_{1} \frac{f (x)}{h^{α}} . \end{matrix}

Similarly

\begin{matrix} f^{'} (x) h^{1 - α} \sum_{k = 1}^{n - 1} \frac{1}{k^{α}} = f^{'} (x) h^{1 - α} ζ (α) + \frac{f^{'} (x) x^{1 - α}}{1 - α} - \frac{f^{'} (x) h}{2 x^{α}} + f^{'} (x) h^{1 - α} R_{2}, \end{matrix}

\begin{matrix} \frac{f^{''} (x)}{2} h^{2 - α} \sum_{k = 1}^{n - 1} k^{1 - α} = ζ (α - 1) \frac{f^{''} (x)}{2} h^{2 - α} + \frac{f^{''} (x) x^{2 - α}}{2 - α} - \frac{f^{''} (x) x^{1 - α} h}{4} + R_{3} \frac{f^{''} (x)}{2} h^{2 - α} . \end{matrix}

Hence

\begin{matrix} \frac{1}{h^{α}} \sum_{k = 1}^{n - 1} \frac{f_{n - k}}{k^{1 + α}} = Γ (- α) f^{(α)} (x) + \frac{f (x)}{h^{α}} & ζ (1 + α) - f^{'} (x) h^{1 - α} ζ (α) \\ + \frac{f^{''} (x)}{2} h^{2 - α} ζ (α - 1) + E_{n}, \end{matrix}

where

E_{n} = T_{n} + R_{1} \frac{f (x)}{h^{α}} - f^{'} (x) h^{1 - α} R_{2} + R_{3} \frac{f^{''} (x)}{2} h^{2 - α},

\begin{matrix} | E_{n} | < | T_{n} | + | R_{1} \frac{f (x)}{h^{α}} | + | f^{'} (x) h^{1 - α} R_{2} | + | R_{3} \frac{f^{''} (x)}{2} h^{2 - α} | . \end{matrix}

(21)

From Taylor’s Theorem

f (x) = \frac{1}{2} f^{''} (ξ_{1}) x^{2}, f^{'} (x) = f^{''} (ξ_{2}) x,

where

0 < ξ_{2} < x

for

i = 1, 2

. Therefore

M_{0} \leq \frac{1}{2} M_{2} x^{2}, M_{1} \leq M_{2} x .

The terms on the right hand side of (21) satisfy the estimates

\begin{matrix} | R_{1} f (x) h^{- α} | < \frac{(1 + α) M_{0} h^{- α}}{12 n^{2 + α}} = \frac{(1 + α) M_{0} h^{2}}{12 x^{2 + α}} < \frac{(1 + α) M_{2} h^{2 - α}}{24}, \\ | R_{2} f^{'} (x) h^{1 - α} | < \frac{α M_{1} h^{1 - α}}{12 m^{1 + α}} = \frac{α M_{1} h^{2}}{12 x^{1 + α}} < \frac{α M_{2} h^{2}}{12 x^{α}} < \frac{α M_{2} h^{2 - α}}{12}, \\ | R_{3} \frac{f^{''} (x)}{2} h^{2 - α} | < \frac{(1 - α) M_{2} h^{2 - α}}{24 n^{α}} = \frac{(1 - α) M_{2} h^{2}}{24 x^{α}} < \frac{(1 - α) M_{2} h^{2 - α}}{24} . \end{matrix}

Hence

\begin{matrix} | E_{n} | & < \frac{1}{3} M_{3} x (x^{2} + 1) h^{2 - α} + \frac{(1 + α) M_{2} h^{2 - α}}{12} \\ < \frac{1}{3} M_{3} x (x^{2} + 1) h^{2 - α} + \frac{M_{2} h^{2 - α}}{6} = \frac{2 M_{3} x (x^{2} + 1) + M_{2}}{6} h^{2 - α} . \end{matrix}

□

Let

L_{n}

be the error of approximation (2).

L_{n} f (x) = \frac{h^{- α}}{Γ (- α)} \sum_{k = 1}^{m} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + L_{n} h^{2 - α},

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

Approximation (2) is constructed from approximation (1) by substituting

f_{n}^{'}

with first order backward difference approximation

(f_{n} - f_{n - 1}) / h

Claim 11.

Let

f \in C^{3} [0, x]

and

f (0) = f^{'} (0) = 0

. Then

| L_{n} | < \frac{2 M_{3} x (x^{2} + 1) + (1 + 3 | ζ (α) |) M_{2}}{6 | Γ (- α) |} h^{2 - α} .

Proof.

From Taylor’s Theorem

f_{n - 1} = f_{n} - h f_{n}^{'} + \frac{h^{2}}{2} f^{''} (ξ),

|f_{n}^{'} - \frac{f_{n} - f_{n - 1}}{h}| < \frac{h}{2} |f^{''} (ξ)|,

where

ξ \in (x_{n - 1}, x_{n})

. The error

L_{n}

of approximation (2) satisfies

L_{n} = \frac{1}{Γ (- α)} (E_{n} - \frac{1}{2} ζ (α) f^{''} (ξ) h^{2 - α}),

\begin{matrix} | L_{n} | & < \frac{1}{| Γ (- α) |} (\frac{2 M_{3} x (x^{2} + 1) + M_{2}}{6} h^{2 - α} + \frac{| ζ (α) | M_{2}}{2} h^{2 - α}) \\ < \frac{2 M_{3} x (x^{2} + 1) + (1 + 3 | ζ (α) |) M_{2}}{6 | Γ (- α) |} h^{2 - α} . \end{matrix}

□

Let

A_{n}

be the error of approximation (3).

\begin{matrix} A_{n} f (x) = \frac{h^{- α}}{Γ (- α)} \sum_{k = 1}^{n} σ_{k}^{(α)} f_{n - k} = f_{n}^{(α)} + A_{n} . \end{matrix}

The weights of approximation (3) are defined with (4), (5) and (6). Denote

\tilde{f} (t) = f (t) - f (0) - f^{'} (0) t .

The function

\tilde{f}

satisfies

\tilde{f} (0) = {\tilde{f}}^{'} (0) = 0

and the second and third derivatives of f and

\tilde{f}

are equal. Therefore

max_{0 \leq t \leq x} | {\tilde{f}}^{''} (t) | = M_{2}, max_{0 \leq t \leq x} | {\tilde{f}}^{'''} (t) | = M_{3} .

Lemma 12.

Let

f \in C^{3} [0, x]

. Then

| A_{n} | < \frac{4 M_{3} x (x^{2} + 1) + 3 (1 + 2 | ζ (α) |) M_{2}}{12 | Γ (- α) |} h^{2 - α} .

Proof.

\begin{matrix} A_{n} f (x) & = A_{n} \tilde{f} (x) + A_{n} [f (x) - \tilde{f} (x)] \\ = L_{n} \tilde{f} (x) + \frac{{\bar{σ}}_{n - 1}^{(α)} {\tilde{f}}_{1}}{Γ (- α) h^{α}} + D_{x}^{α} [f (x) - \tilde{f} (x)] . \end{matrix}

From Claim 4

|{\bar{σ}}_{n - 1}^{(α)}| < \frac{1 + 2 α}{12 n^{1 + α}} < \frac{1}{4 n} < \frac{1}{8} .

From Taylor’s Theorem

| {\tilde{f}}_{1} | = \frac{h^{2}}{2} | {\tilde{f}}^{''} (ξ) | \leq \frac{h^{2} M_{2}}{2} .

Then

\begin{matrix} | A_{n} | & < | L_{n} | + \frac{1}{| Γ (- α) | h^{α}} |{\bar{σ}}_{n - 1}^{(α)}| | {\tilde{f}}_{1} | \\ < (\frac{2 M_{3} x (x^{2} + 1) + (1 + 3 | ζ (α) |) M_{2}}{6 | Γ (- α) |} + \frac{M_{2}}{16 Γ (- α)}) h^{2 - α} \\ < \frac{4 M_{3} x (x^{2} + 1) + 3 (1 + 2 | ζ (α) |) M_{2}}{12 | Γ (- α) |} h^{2 - α} . \end{matrix}

□

4. Numerical Solution of Two-Term Equation

The two-term equation is an ordinary fractional differential equation in the form

\begin{matrix} y^{(α)} (t) + D y (t) = F (t), y (0) = y_{0} . \end{matrix}

(22)

Numerical and analytical solutions of ordinary fractional differential equations are studied in [6,24,25,41,43,45]. In this section we construct the numerical solution of the two-term equation (22) which uses approximation (3) of the fractional derivative and prove its convergence and order. Let

h = x / N

and

x_{n} = n h

for

n = 0, \dots, N

. By approximating the fractional derivative at the point

x_{n}

we find

\begin{matrix} \frac{1}{Γ (- α) h^{α}} \sum_{k = 0}^{n} σ_{k}^{(α)} y_{n - k} + D y_{n} = F_{n} + A_{n}, \end{matrix}

\begin{matrix} (σ_{0}^{(α)} + Γ (- α) D h^{α}) y_{n} + \sum_{k = 1}^{n} σ_{k}^{(α)} y_{n - k} + D y_{n} = Γ (- α) h^{α} F_{n} + B_{n} h^{2}, \end{matrix}

where

B_{n} = Γ (- α) A_{n} h^{α - 2}

\begin{matrix} (ζ (α) - ζ (1 + α) + Γ (- α) D h^{α}) y_{n} - ζ (α) y_{n - 1} + \sum_{k = 1}^{n} \frac{y_{n - k}}{k^{1 + α}} - {\bar{σ}}_{n - 1}^{(α)} u_{1} = Γ (- α) h^{α} F_{n} + B_{n} h^{2} . \end{matrix}

The numerical solution

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

of equation (22) is computed as

\begin{matrix} (ζ (α) - ζ (1 + α) + Γ (- α) D h^{α}) u_{n} - ζ (α) u_{n - 1} + \sum_{k = 1}^{n} \frac{u_{n - k}}{k^{1 + α}} + {\bar{σ}}_{n - 1}^{(α)} u_{1} + D u_{n} = Γ (- α) h^{α} F_{n}, \end{matrix}

\begin{matrix} u_{n} = \frac{1}{ζ (α) - ζ (1 + α) + Γ (- α) D h^{α}} (Γ (- α) h^{α} F_{n} + ζ (α) u_{n - 1} - \sum_{k = 1}^{n} \frac{u_{n - k}}{k^{1 + α}} - {\bar{σ}}_{n - 1}^{(α)} u_{1}) . \end{matrix}

(23)

The value of the numerical solution

u_{1}

on the first step is computed with the following approximation [17]

\begin{matrix} y_{1}^{(α)} = \frac{y_{1} - y_{0}}{Γ (2 - α) h^{α}} + a_{1} h^{2 - α} . \end{matrix}

(24)

From the estimate for the error of L1 approximation [8]

\begin{matrix} | a_{1} | < \frac{1}{Γ (2 - α)} (\frac{2 + α}{2^{α} (2 - α)} - \frac{11 + α}{12}) max_{0 \leq t \leq h} |y^{''} (t)| < \frac{2 M_{2}}{Γ (2 - α)}, \end{matrix}

(25)

where

M_{2} = {max}_{0 \leq t \leq x} |y^{''} (t)|

. Numerical solution (23) has initial conditions [17]

u_{0} = y (0), u_{1} = \frac{y (0) + Γ (2 - α) F (h) h^{α}}{1 + D Γ (2 - α) h^{α}} .

When the solution of two-term equation (22) satisfies the condition

y (0) = y^{'} (0) = 0

both approximations (2) and (3) can be used for computation of numerical solution (23).

Example 1: Consider the following two-term equation

y^{(α)} + D y = - α^{2} t^{2 - α} E_{2, 3 - α} (- {(α t)}^{2}) + D (cos (α t) - 1), y_{0} = 0 .

(26)

Two-term equation (26) has a solution

y (x) = cos (α t) - 1

which satisfies

y (0) = y^{'} (0) = 0

. Experimental results for the error and order of numerical solution (23) of equation (26) which uses approximation (2) of the fractional derivative are given in Table 1 and Table 2.

Example 2:

y^{(α)} + D y = α t^{1 - α} E_{1, 2 - α} (α t) + D e^{α t}, y_{0} = 1 .

(27)

Equation (27) has a solution

y (t) = e^{α t}

. Experimental results for the error and order of numerical solution (23) of equation (27) are given in Table 3 and Table 4. The experimental results in Table 1, Table 2, Table 3 and Table 4 suggest that numerical solution (23) of two-term equation (22) converges and has an order

2 - α

for all positive values of the parameter D and negative values with small and large modulus.

Denote by

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

the error of numerical solution (23) at the point

x_{n} = n h

, where

n = 0, 1, \dots, N

. The errors

e_{n}

satisfy

\begin{matrix} e_{n} = \frac{1}{ζ (α) - ζ (1 + α) + Γ (- α) D h^{α}} (ζ (α) e_{n - 1} - \sum_{k = 1}^{n - 1} \frac{e_{n - k}}{k^{1 + α}} - {\bar{σ}}_{n - 1}^{(α)} e_{1} + B_{n} h^{2}), \end{matrix}

(28)

e_{0} = 0, e_{1} = \frac{Γ (2 - α) a_{1} h^{2}}{1 + D Γ (2 - α) h^{α}} .

In Theorem 13 and Theorem 14 we prove the convergence of numerical solution (23) of two-term equation (22) and derive estimates for the error, depending on the value of the parameter D. The proofs use induction on n, where

n = 1, 2, \dots, N

Theorem 13.

Let

D \in (- \infty, - \frac{2}{h^{α}} \frac{ζ (α) - ζ (1 + α)}{Γ (- α)}) \cup (0, \infty)

and

y \in C^{2} [0, 1]

. Then

\begin{matrix} | e_{n} | < C h^{2} n^{α}, (1 \leq n \leq N), \end{matrix}

(29)

where

C = 2 M_{2} + \frac{1}{2} max_{1 \leq n \leq N} | B_{n} |

Proof.

The value of the error

e_{1}

on the first step satisfies

\begin{matrix} | e_{1} | = \frac{Γ (2 - α) | a_{1} | h^{2}}{| 1 + D Γ (2 - α) h^{α} |} < \frac{2 M_{2} h^{2}}{| 1 + D Γ (2 - α) h^{α} |} . \end{matrix}

(30)

When

D > 0

the denominator

| 1 + D Γ (2 - α) h^{α} | > 1

and

| e_{1} | < 2 M_{2} h^{2} < C h^{2} .

Now we estimate the error

e_{1}

when the parameter

D < - \frac{2}{h^{α}} \frac{ζ (α) - ζ (1 + α)}{Γ (- α)}

D Γ (- α) h^{α} > 2 (| ζ (α) | + ζ (1 + α)),

| D | Γ (2 - α) h^{α} > 2 α (1 - α) (| ζ (α) | + ζ (1 + α)) .

From (10)

ζ (1 + α) > \frac{1}{α} + γ, | ζ (α) | > \frac{1}{1 - α} - γ,

| ζ (α) | + ζ (1 + α) > \frac{1}{1 - α} + \frac{1}{α} = \frac{1}{α (1 - α)} .

Hence

| D | Γ (2 - α) h^{α} > 2 .

In this case the denominator of (30) is again greater than one,

| 1 + D Γ (2 - α) h^{α} | > 1

and the error

e_{1}

satisfies

| e_{1} | < 2 M_{2} h^{2} < C h^{2} .

Suppose that inequality (29) holds for all

n = 1, \dots, m - 1

. The values of

ζ (α)

and

Γ (- α)

are negative. In both cases

D > 0

and

D < - \frac{2}{h^{α}} \frac{ζ (α) - ζ (1 + α)}{Γ (- α)}

the following estimate holds

| ζ (α) - ζ (1 + α) + Γ (- α) D h^{α} | > | ζ (α) | + ζ (1 + α) .

From (28) and Claim 4

\begin{matrix} | e_{m} | & < \frac{1}{| ζ (α) | + ζ (1 + α)} (| ζ (α) | e_{m - 1} + \sum_{k = 1}^{m - 1} \frac{|e_{m - k}|}{k^{1 + α}} + |{\bar{σ}}_{m - 1}^{(α)} e_{1}| h^{2} + | B_{m} | h^{2}) \\ < \frac{1}{| ζ (α) | + ζ (1 + α)} (| ζ (α) | e_{m - 1} + \sum_{k = 1}^{m - 1} \frac{|e_{m - k}|}{k^{1 + α}} + \frac{M_{2}}{4} h^{2} + | B_{m} | h^{2}) \\ < \frac{1}{| ζ (α) | + ζ (1 + α)} (| ζ (α) | e_{m - 1} + \sum_{k = 1}^{m - 1} \frac{|e_{m - k}|}{k^{1 + α}} + 2 C h^{2}) . \end{matrix}

By induction hypoyhesis

\begin{matrix} | e_{m} | & < \frac{C h^{2}}{| ζ (α) | + ζ (1 + α)} ({| ζ (α) | (m - 1)}^{α} + \sum_{k = 1}^{m - 1} \frac{{(m - k)}^{α}}{k^{1 + α}} + 2) . \end{matrix}

From inequality (7)

\begin{matrix} | e_{m} | & < \frac{C h^{2}}{| ζ (α) | + ζ (1 + α)} (| ζ (α) | m^{α} + (ζ (1 + α) m^{α} - 2) + 2) = C h^{2} m^{α} . \end{matrix}

□

Theorem 14.

Let

\frac{1}{Γ (- α)} < D < 0

and

y \in C^{2} [0, 1]

. Then

\begin{matrix} | e_{n} | < K h^{2} n^{α}, (1 \leq n \leq N), \end{matrix}

(31)

where

K = 3 M_{2} + max_{1 \leq n \leq N} | B_{n} |

Proof.

The error

e_{1}

on the first step satisfies the bound

\begin{matrix} | e_{1} | < \frac{2 M_{2} h^{2}}{1 - α (1 - α) D Γ (- α) h^{α}} < \frac{2 M_{2} h^{2}}{1 - α (1 - α)} \leq \frac{8}{3} M_{2} h^{2} . \end{matrix}

Assume that (31) holds for all

n = 1, \dots, m - 1

. Then

\begin{matrix} | e_{m} | & < \frac{1}{| ζ (α) | + ζ (1 + α) - Γ (- α) D h^{α}} (| ζ (α) | e_{m - 1} + \sum_{k = 1}^{m - 1} \frac{|e_{m - k}|}{k^{1 + α}} + |{\bar{σ}}_{m - 1}^{(α)} e_{1}| h^{2} + | B_{m} | h^{2}) \\ < \frac{1}{| ζ (α) | + ζ (1 + α) - h^{α}} (| ζ (α) | K h^{2} {(m - 1)}^{α} + K h^{2} \sum_{k = 1}^{m - 1} \frac{{(m - k)}^{α}}{k^{1 + α}} + K h^{2}) \\ < \frac{K h^{2}}{| ζ (α) | + ζ (1 + α) - h^{α}} (| ζ (α) | m^{α} + ζ (1 + α) m^{α} - 1) \\ < \frac{K h^{2} m^{α}}{| ζ (α) | + ζ (1 + α) - h^{α}} (| ζ (α) | + ζ (1 + α) - \frac{1}{m^{α}}) \leq K h^{2} m^{α} . \end{matrix}

□

From Theorem 13 and Theorem 14 the errors

e_{n}

of numerical solution (23) of the two-term equation in the interval

[0, 1]

, satisfy the estimate

| e_{n} | < K h^{2} n^{α} \leq K h^{2 - α} .

(32)

for all

n = 1, \dots, N

. Now we generalize the results of Theorem 13 and Theorem 14 to an arbitrary interval

[a, b]

. Consider the two-term equation

y^{(α)} (t) + D y (t) = F (t), y (a) = y_{0}, t \in [a, b] .

(33)

Substitute

t = a + (b - a) s

and

z (s) = y (t) = y (a + (b - a) s) .

The function

z (s)

satisfies a two-term equation

z^{(α)} (s) + {(b - a)}^{α} D z (s) = G (s) = {(b - a)}^{α} F (a + (b - a) s), z (0) = y_{0}, s \in [0, 1]

(34)

Let

E = {e_{n}}_{n = 1}^{N}

be an N-dimensional vector with the errors of numerical solution (23) of equation (33). The

L_{\infty}

(maximum) norm of a vector is the maximum of the absolute values of its elements. From Theorem 13 and Theorem 14 we get conditions for the convergence of numerical solution (23) of two-term equations (33) and (34).

Corollary 15.

The maximum error of numerical solution (23) of two-term equation (33), satisfies

∥ E ∥ < K {(b - a)}^{α} h^{2 - α}

when the solution

y \in C^{2} [a, b]

and

h = (b - a) / N

D \in (- \infty, - \frac{2 (ζ (α) - ζ (1 + α))}{Γ (- α) h^{α}}) \cup (\frac{1}{Γ (- α) {(b - a)}^{α}}, \infty)

□

5. Time Fractional Black–Scholes Equation

The time fractional Black–Scholes equation is a fractional partial differential equation which is used for modeling the prices of the options [1,11,13,28,42,54,59].

\begin{matrix} \frac{\partial^{α} C (S, t)}{\partial t^{α}} + \frac{1}{2} σ^{2} S^{2} \frac{\partial^{2} C (S, t)}{\partial S^{2}} + r S \frac{\partial C (S, t)}{\partial S} - r C (S, t) = 0, \end{matrix}

(35)

where C is the option price, T is the expiry time, r is the risk-free rate,

σ

is the volatility and S is the underlying stock price. Fractional Black-Scholes equation (35) has terminal and boundary conditions

C (S, T) = C_{0} (S), C (0, t) = C_{1} (t), C (\infty, t) = C_{2} (t)

and the fractional derivative is

\begin{matrix} \frac{\partial^{α} C (S, t)}{\partial t^{α}} = \frac{1}{Γ (1 - α)} \frac{d}{d t} \int_{t}^{T} \frac{C (S, ξ) - C (S, T)}{{(ξ - z)}^{α}} d ξ . \end{matrix}

(36)

The terminal condition is transformed to an initial condition with the substitution [59]

t = T - z, η = T - ξ, v (S, μ) = C (S, ξ) = C (S, t - μ) .

By applying the substitution to (36) we find [59]

\begin{matrix} \frac{\partial^{α} C (S, t)}{\partial t^{α}} & = - \frac{1}{Γ (1 - α)} \frac{d}{d z} \int_{z}^{T} \frac{C (S, ξ) - C (S, T)}{{(ξ - z)}^{α}} d ξ \\ = - \frac{1}{Γ (1 - α)} \frac{d}{d z} \int_{0}^{z} \frac{v (S, μ) - v (S, 0)}{{(z - μ)}^{α}} d μ = D_{z}^{α} v (S, z) . \end{matrix}

The function

v (S, z)

satisfies the partial fractional differential equation

\begin{matrix} D_{z}^{α} v (S, z) = \frac{1}{2} σ^{2} S^{2} \frac{\partial^{2} v (S, z)}{\partial S^{2}} + r S \frac{\partial v (S, z)}{\partial S} - r C (S, z) . \end{matrix}

(37)

The substitution

u (x, t) = v (e^{x}, z)

transforms (37) into a linear partial fractional differential equation [59]

D_{t}^{α} u (x, t) = \frac{1}{2} σ^{2} \frac{\partial^{2} u (x, t)}{\partial x^{2}} + (r - \frac{σ^{2}}{2}) \frac{\partial u (x, t)}{\partial x} - r u (x, t) .

(38)

Finite difference schemes for the time fractional Black-Scholes equation are constructed in [19,23,52,59]. Now we construct an implicit finite difference scheme which uses approximation (3) of the fractional derivative and central difference approximations of the partial derivatives for the following time fractional Black-Scholes equation

D_{t}^{α} u (x, t) = \frac{1}{2} σ^{2} \frac{\partial^{2} u (x, t)}{\partial x^{2}} + (r - \frac{σ^{2}}{2}) \frac{\partial u (x, t)}{\partial x} - r u (x, t) + F (x, t),

(39)

where

(x, t) \in R = [B_{l}, B_{r}] \times [0, T]

. Equation (39) has initial and boundary conditions

u (x, 0) = u_{0} (x), u (B_{l}, t) = u_{1} (t), u (B_{r}, t) = u_{2} (t) .

Let M and N be positive integers and

G

be a rectangular grid on

R

which has a step size

h = (B_{r} - B_{l}) / N

in space and

τ = T / M

in time,

G = {(n h, m τ) | 0 \leq n \leq N, 0 \leq m \leq M} .

Denote

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

. The central difference approximations of the partial derivatives of equation (39) have second order accuracy

\frac{\partial u (n h, m τ)}{\partial x} = \frac{u_{n + 1}^{m} - u_{n - 1}^{m}}{2 h} + e_{n, m}^{1} h^{2},

\frac{\partial^{2} u (n h, m τ)}{\partial x^{2}} = \frac{u_{n + 1}^{m} - 2 u_{n}^{m} + u_{n - 1}^{m}}{h^{2}} + e_{n, m}^{2} h^{2},

where

|e_{n, m}^{1}| < {\bar{M}}_{3} / 3

and

|e_{n, m}^{2}| < {\bar{M}}_{4} / 12

. The numbers

{\bar{M}}_{3}

and

{\bar{M}}_{4}

are bounds for the third and fourth order partial derivatives

{\bar{M}}_{3} = max_{(x, t) \in R} |\frac{\partial^{3} u (x, t)}{\partial x^{3}}|, {\bar{M}}_{4} = max_{(x, t) \in R} |\frac{\partial^{4} u (x, t)}{\partial x^{4}}| .

The numerical solution

{U_{n}^{m}}_{n = 1}^{N - 1}

of equation (39) on layer m of the grid

G

satisfies

\begin{matrix} \frac{τ^{- α}}{Γ (- α)} \sum_{k = 0}^{m} σ_{k}^{(α)} U_{n}^{m - k} = & \frac{σ^{2}}{2 h^{2}} (U_{n + 1}^{m} - 2 U_{n}^{m} + U_{n - 1}^{m}) \\ + \frac{1}{2 h} (r - \frac{σ^{2}}{2}) (U_{n + 1}^{m} - U_{n - 1}^{m}) - r U_{n}^{m} + F_{n}^{m} \end{matrix}

for

m = 2, \dots, M

and has has initial conditions

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

Denote

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

and

K_{1} = - (\frac{σ^{2}}{2} - (r - \frac{σ^{2}}{2}) h) η, K_{2} = σ_{0}^{(α)} + σ^{2} η + r h^{2} η, K_{3} = - (\frac{σ^{2}}{2} + (r - \frac{σ^{2}}{2}) h) η .

The numerical solution of time-fractional Black-Scholes equation (39) is a solution of the system of linear equations

\begin{matrix} K_{1} U_{n - 1}^{m} - K_{2} U_{n}^{m} + K_{3} U_{n + 1}^{m} = - \sum_{k = 1}^{m} σ_{k}^{(α)} U_{n}^{m - k} + h^{2} η F_{n}^{m} . \end{matrix}

(40)

Denote by

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

the error of finite difference scheme (40) at the point

(n h, m τ)

and by

A_{n}^{m} τ^{2 - α}

and

B_{n}^{m} h^{2}

be the the truncation errors of the approximations of the left-hand and right-hand sides of (39). The errors

E_{n}^{m}

on row

m \geq 2

of the grid

G

satisfy

K_{1} E_{n - 1}^{m} + K_{2} E_{n}^{m} + K_{3} E_{n + 1}^{m} = - \sum_{k = 1}^{m} σ_{k}^{(α)} E_{n}^{m - k} + Γ (- α) τ^{α} (A_{n}^{m} τ^{2 - α} + B_{n}^{m} h^{2}) .

The errors are equal to zero,

E_{n}^{0} = E_{0}^{m} = E_{N}^{m} = 0

at the boundary points of

G

. Let

K

be the

(N - 1)

-dimensional tridiagonal matrix with entries

K_{2}

on the main diagonal and

K_{1}

and

K_{3}

bellow and above the main diagonal. The error vector

E^{m}

of row m is a solution of the matrix equation

K E^{m} = Q^{m}; E^{m} = K^{- 1} Q^{m},

where

\begin{matrix} Q_{n}^{m} = ζ (α) E_{n}^{m - 1} - \sum_{k = 1}^{m} \frac{E_{n}^{m - k}}{k^{1 + α}} - {\bar{σ}}_{m - 1}^{(α)} E_{n}^{1} + Γ (- α) τ^{α} (A_{n}^{m} τ^{2 - α} + B_{n}^{m} h^{2}) . \end{matrix}

(41)

The numerical solution on the first row of

G

is computed with approximation (24)

\begin{matrix} \frac{U_{n}^{1} - U_{n}^{0}}{Γ (2 - α) τ^{α}} = & \frac{σ^{2}}{2 h^{2}} (U_{n + 1}^{1} - 2 U_{n}^{1} + U_{n - 1}^{1}) + \frac{1}{2 h} (r - \frac{σ^{2}}{2}) (U_{n + 1}^{1} - U_{n - 1}^{1}) - r U_{n}^{1} + F_{n}^{1} . \end{matrix}

Denote

\tilde{η} = Γ (2 - α) τ^{α} / h^{2}

. Then

{\tilde{K}}_{1} U_{n - 1}^{1} + {\tilde{K}}_{2} U_{n}^{1} + {\tilde{K}}_{3} U_{n + 1}^{1} = U_{n}^{0} + h^{2} η F_{n}^{1},

where

1 \leq n \leq N

and

{\tilde{K}}_{1} = - (\frac{σ^{2}}{2} - (r - \frac{σ^{2}}{2}) h) \tilde{η}, {\tilde{K}}_{2} = 1 + σ^{2} \tilde{η} + r h^{2} \tilde{η}, {\tilde{K}}_{3} = - (\frac{σ^{2}}{2} + (r - \frac{σ^{2}}{2}) h) \tilde{η} .

The errors of the numerical solution of the fractional Black-Scholes equation (39) on the first row of

G

are the solutions of the system of linear equations

{\tilde{K}}_{1} E_{n - 1}^{1} + {\tilde{K}}_{2} E_{n}^{1} + {\tilde{K}}_{3} E_{n + 1}^{1} = Γ (2 - α) τ^{α} (A_{n}^{1} τ^{2 - α} + B_{n}^{1} h^{2}) .

Let

\tilde{K}

be the

(N - 1)

-dimensional tridiagonal matrix with elements

{\tilde{K}}_{1}, {\tilde{K}}_{2}

and

{\tilde{K}}_{3}

. The error vector for the first row satisfies

\tilde{K} E^{1} = Q^{1}; E^{1} = {\tilde{K}}^{- 1} Q^{1},

where

Q_{n}^{1} = Γ (2 - α) τ^{α} (A_{n}^{1} τ^{2 - α} + B_{n}^{1} h^{2}) .

Let

C = {max}_{n, m} {| A_{n}^{m} |, | B_{n}^{m} |}

for all

n = 1, \dots, N; m = 1, \dots, M

and

K = | Γ (- α) C |

. The

L_{\infty}

norm of a matrix is the maximum of the absolute row sums.

Theorem 16.

Let

u \in C^{4} [B_{l}, B_{r}] \times C^{2} [0, T]

and

r < \frac{σ^{2}}{2 h}

. Then

\begin{matrix} ∥ E^{m} ∥ < K m^{α} τ^{α} (h^{2} + τ^{2 - α}) \end{matrix}

(42)

for all

m = 1, \dots M

Proof.

When

0 < α < 1

, the gamma function satisfies [15]

0.88 < Γ (2 - α) < 1,

| Γ (- α) | = \frac{Γ (2 - α)}{α (1 - α)} > 3 .

Then

∥ Q^{1} ∥ < C Γ (2 - α) τ^{α} (τ^{2 - α} + h^{2}) < C τ^{α} (τ^{2 - α} + h^{2}) .

When

r < σ^{2} / (2 h)

we have that

|K_{1}| + | K_{3} | = |\frac{σ^{2}}{2} - (r - \frac{σ^{2}}{2}) h| | η | + |\frac{σ^{2}}{2} + (r - \frac{σ^{2}}{2}) h| | η | \leq σ^{2} | η |,

|K_{2}| - |K_{1}| - | K_{3} | > | σ_{0}^{(α)} | .

Similarly

|{\tilde{K}}_{2}| - |{\tilde{K}}_{1}| - | {\tilde{K}}_{3} | > 1 .

The matrices

K

and

\tilde{K}

are M-matrices. From the Ahlberg-Nilson-Varah bound [27,40]

∥ {\tilde{K}}^{- 1} ∥ < 1, ∥ K^{- 1} ∥ < \frac{1}{| σ_{0}^{(α)} |} = \frac{1}{| ζ (α) | + ζ (1 + α)} .

The

L_{\infty}

norm of the error vector of the first row of

G

satisfies the bound

∥ E^{1} ∥ \leq ∥ {\tilde{K}}^{- 1} ∥ ∥ Q^{1} ∥ < C τ^{α} (τ^{2 - α} + h^{2}) < K τ^{α} (τ^{2 - α} + h^{2}) .

Suppose that inequality (42) holds for all

m = 1, \dots, p - 1

. From (41)

\begin{matrix} |Q_{n}^{p}| & < | ζ (α) | ∥ E^{p - 1} ∥ + \sum_{k = 1}^{p} \frac{∥ E^{p - k} ∥}{k^{1 + α}} + |{\bar{σ}}_{p - 1}^{(α)} E_{n}^{1}| + Γ (- α) C τ^{α} (h^{2} + τ^{2 - α}) \\ < K τ^{α} (h^{2} + τ^{2 - α}) (| ζ (α) | p^{α} + \sum_{k = 1}^{p} \frac{{(p - k)}^{α}}{k^{1 + α}} + \frac{1}{8 | Γ (- α) |} + 1) \\ < K τ^{α} (h^{2} + τ^{2 - α}) (| ζ (α) | p^{α} + ζ (1 + α) p^{α} - 1 + \frac{1}{8 | Γ (- α) |}) . \end{matrix}

Therefore

∥ Q^{p} ∥ < (| ζ (α) | + ζ (1 + α)) K p^{α} τ^{α} (h^{2} + τ^{2 - α}),

∥ E^{p} ∥ < ∥ K^{- 1} ∥ ∥ Q^{p} ∥ < K p^{α} τ^{α} (h^{2} + τ^{2 - α}) .

□

Corollary 17.

Let

r < \frac{σ^{2}}{2 h}

. Then

∥ E^{m} ∥ < K T^{α} (h^{2} + τ^{2 - α})

for all

m = 1, 2, \dots, M

Proof.

From Theorem 16

∥ E^{m} ∥ < K m^{α} τ^{α} (h^{2} + τ^{2 - α}) \leq K M^{α} τ^{α} (h^{2} + τ^{2 - α}) = K T^{α} (h^{2} + τ^{2 - α}) .

□

Example 3: Consider the following time fractional Black-Scholes equation

\{\begin{matrix} D_{t}^{α} u (x, z) = \frac{1}{2} σ^{2} \frac{\partial^{2} u (x, t)}{\partial S^{2}} + (r - \frac{σ^{2}}{2}) \frac{\partial u (x, t)}{\partial S} - r u (x, t) + t^{1 - α} e^{x} E_{1, 2 - α} (t), \\ u (x, 0) = e^{x}, u (0, t) = e^{t}, u (1, t) = e^{1 + t} . \end{matrix}

(43)

Equation (43) has a solution

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

. Experimental results of the numerical solution of equation (43) with parameters

σ = 0.35, r = 0.05

are given in Table 5 and Table 6. The orders of convergence of the finite difference scheme are

{log}_{\frac{τ_{1}}{τ_{2}}} \frac{{error}_{1}}{{error}_{2}}

in time and

{log}_{\frac{h_{1}}{h_{2}}} \frac{{error}_{1}}{{error}_{2}}

in space.

Example 4: Consider the time fractional Black-Scholes equation (37) for pricing European call options, with a source term

F (x, t) = 0

and initial and boudary conditions

u_{0} (x) = {(e^{x} - K)}^{+}, u_{1} (t) = 0, u_{2} (t) = e^{B_{r}} - K e^{- r t} .

The European call option premium curves for different values of

α

are given in Figure 1, for values of the parameters

r = 2 %

σ = 30 %

B_{l} = 0

B_{r} = 200

T = 1

(year) and strike price

K = 100

. Regarding near-the-money options, lower values of

α

prices them lower, and evaluates higher out-of-the-money and in-the-money options, compared to the classical Black-Scholes dynamics.

The graphs of call option and put option premiums for

α = 0.7

and all t are given in Figure 2 and Figure 3.

6. Conclusions

In this paper we study the convergence and order of numerical solutions of ordinary and partial fractional differential equations which use approximation (3) of the fractional derivative. The results of the numerical experiments given in the paper illustrate the theoretical results. In future work we will use the method from [7] for construction of secon-order and high-order approximations of the fractional derivative and we will study the convergence of the finite difference schemes for numerical solution of fractional differential equations.

Author Contributions

Conceptualization, Y.D and V.T.; Methodology, Y.D. and S.A.; Software, S.G.; Validation, V.T.; Formal analysis, Y.D.; Writing–original draft preparation, Y.D. and S.G.; writing–review and editing, Y.D. and V.T.; Visualization, S.G.; Project administration, Y.D. and V.T.; Funding acquisition, S.G. 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 and by the National Program “Young Scientists and Postdoctoral Researchers - 2” – Bulgarian Academy of Sciences.

Conflicts of Interest

The authors declare no conflict of interest.

References

Abbas, M.; Bibi, A.; Alzaidi, A.S.M.; Nazir, T.; Majeed, A.; Akram, G. Numerical solutions of third-order time-fractional differential equations using cubic B-spline functions. Fractal Fract. 2022, 6, 528. [Google Scholar] [CrossRef]
Abdi, N.; Aminikhah, H.; Sheikhani, A.H.R. High-order compact finite difference schemes for the time-fractional Black-Scholes model governing European options. Chaos Solit. Fractals 2022, 162, 112423. [Google Scholar] [CrossRef]
Acay, B.; Bas, E.; Abdeljawad, T. Fractional economic models based on market equilibrium in the frame of different type kernels. Chaos Solit. Fractals. 2020, 130, 109438. [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]
Alikhanov, A.A.; Huang, C. A high-order L2 type difference scheme for the time-fractional diffusion equation. Appl. Math. Comput. 2021, 411, 1–19. [Google Scholar] [CrossRef]
Anjara, F.; Solofoniaina, J. Solution of general fractional oscillation relaxation equation by Adomian’s method. Gen. math. notes 2014, 20(2), 1–11. [Google Scholar]
Apostolov, S.; Dimitrov, Y.; Todorov, V. Constructions of second order approximations of the Caputo fractional derivative. Lect. Notes Comput. Sci. 2022, 13127, 31–39. [Google Scholar]
Cai, M.; Li, C. Numerical approaches to fractional integrals and derivatives: a review. Mathematics 2020, 8(1), 43. [Google Scholar] [CrossRef]
Cao, J.; Cai, Z. Numerical analysis of a high-order scheme for nonlinear fractional differential equations with uniform accuracy. Numerical Mathematics: Theory, Methods and Applications 2020, 14, 71–112. [Google Scholar] [CrossRef]
Cardone, A.; Donatelli, M.; Durastante, F.; Garrappa, R.; Mazza, M.; Popolizio, M. Fractional Differential Equations "Modeling, Discretization, and Numerical Solvers". Springer, Singapore, 2023;
Cen, Z.; Huang, J.; Xu, A.; Le, A. Numerical approximation of a time-fractional Black–Scholes equation. Comput. Math. Appl. 2018, 75(8), 2874–2887. [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]
Chen, X.; Xu, X.; Zhu, S.-P. Analytically pricing double barrier options based on a time-fractional Black-Scholes equation. Comput. Math. Appl. 2015, 69(12), 1407–1419. [Google Scholar] [CrossRef]
Dehghan, M.; Safarpoor, M.; Abbaszadeh, M. Two high-order numerical algorithms for solving the multi-term time fractional diffusion-wave equations. J. Comput. Appl. Math. 2015, 290, 174–195. [Google Scholar] [CrossRef]
Deming, W.; Colcord, C. The minimum in the gamma function. Nature 1935, 135, 917. [Google Scholar] [CrossRef]
Dimitrov, Y. A second order approximation for the Caputo fractional derivative. J. Fractional Calc. & Appl. 2016, 7, 175–195. [Google Scholar]
Dimitrov, Y. Approximations for the Caputo derivative (I). J. Fractional Calc. & Appl. 2018, 9, 15–44. [Google Scholar]
Dimitrov, Y. Approximations of the fractional integral and numerical solutions of fractional integral equations. Commun. Appl. Math. Comput. 2021, 2021 3, 545–569. [Google Scholar] [CrossRef]
Dimitrov, Y.; Vulkov, L. Three-point compact finite difference scheme on non-uniform meshes for the time-fractional Black-Scholes equation. AIP Conf Proc . 2015, 1690, 040022. [Google Scholar]
Duan, J.-S.; Li, M.; Wang, Y.; An, Y.-L. Approximate solution of fractional differential equation by quadratic splines. Fractal Fract. 2022, 6, 369. [Google Scholar] [CrossRef]
Edwards, H.M. Riemann’s Zeta Function. Academic Press, New York, USA, 1974;
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]
Grzegorz Krzyżanowski, G.; Magdziarz, M.; Plociniczak, L. A weighted finite difference method for subdiffusive Black–Scholes model, Comput. Math. Appl. 2020, 80(5), 653–670. [Google Scholar]
Güulsu, M.; Öztürk, Y.; Anapalı, A. Numerical approach for solving fractional relaxation-oscillation equation. Appl. Math. Model. 2013, 37(8), 5927–5937. [Google Scholar] [CrossRef]
Hu, Y.; Luo, Y.; Lu, Z. Analytical solution of the linear fractional differential equation by Adomian decomposition method. J. Comput. Appl. Math. 2008, 215(1), 220–229. [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(1), 197–221. [Google Scholar] [CrossRef]
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]
Kumar, S.; Kumar, D.; Singh, J. Numerical computation of fractional Black–Scholes equation arising in financial market. Egyptian Journal of Basic and Applied Sciences 2014, 1, 177–183. [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]
Lin, Y.; Xu, C. Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput. Phys. 2007, 225(2), 1533–1552. [Google Scholar] [CrossRef]
Lubich, C. Fractional linear multistep methods for Abel-Volterra integral equations of the second kind. Math. Comp. 1985, 45, 463–469. [Google Scholar] [CrossRef]
Lubich, C. Discretized fractional calculus. SIAM J. Math. Anal. 1986, 17(3), 704–719. [Google Scholar] [CrossRef]
Lv, C.; Xu, C. Error analysis of a high order method for time-fractional diffusion equations. SIAM J. Sci. Comput. 2016, 38, A2699–A2724. [Google Scholar] [CrossRef]
Marasi, H.; Derakhshan, M.H.; Joujehi, A.S.; Kumar, P. Higher-order fractional linear multi-step methods. Phys. Scr. 2023, 98, 024004. [Google Scholar] [CrossRef]
Mascarenhas, P.V.S.; Moraes, R.M.; Cavalcante, A.L.B. Using a shifted Grünwald–Letnikov scheme for the Caputo derivative to study anomalous solute transport in porous medium. Int. J. Numer. Anal. Methods Geomech. 2019, 43, 1956–1977. [Google Scholar] [CrossRef]
Matsuoka, Y. On the power series coefficients of the Riemann zeta function. Tokyo J. Math. 1989, 12(1), 49–58. [Google Scholar] [CrossRef]
Nasir, H.M.; Nafa, K. A new second order approximation for fractional derivatives with applications. SQUJS 2018, 23, 43–55. [Google Scholar]
Navot, I. An extension of the Euler-Maclaurin summation formula to functions with a branch singularity. J. Math. Phys. 1961, 40, 271–276. [Google Scholar] [CrossRef]
Navot, I. A further extension of the Euler-Maclaurin summation formula. J. Math. Phys. 1962, 41, 155–163. [Google Scholar] [CrossRef]
Nilson, E.N.; Ahlberg, J.H. Convergence properties of the spline fit. J. SIAM 1963, 11, 95–104. [Google Scholar]
Odibat, Z.M.; Momani, S.M. An algorithm for the numerical solution of differential equations of fractional order. J. Appl. Math. Inform. 2008, 26, 15–27. [Google Scholar]
Nuugulu, S.M.; Gideon, F.; Patidar, K.C. A robust numerical solution to a time-fractional Black–Scholes equation. Adv. Differ. Equ. 2021, 2021, 123. [Google Scholar] [CrossRef]
Podlubny, I. Fractional Differential Equation. Academic Press, New Jersey, USA, 1999;
Ramezani, M.; Mokhtari, R.; Haase, G. Some high order formulae for approximating Caputo fractional derivatives. Appl. Numer. Math 2020, 153, 300–318. [Google Scholar] [CrossRef]
Ray, S.S.; Bera, R.K. Analytical solution of the Bagley Torvik equation by Adomian decomposition method. Appl. Math. Comput 2005, 168, 398–410. [Google Scholar] [CrossRef]
Ray, S.S.; Gupta, A.K. Wavelet Methods for Solving Partial Differential Equations and Fractional Differential Equations, Chapman and Hall/CRC; 2018.
ur Rehman, M.; Baleanu, D.; Alzabut, J.; Ismail, .; Saeed, . Green–Haar wavelets method for generalized fractional differential equations. Adv Differ Equ 2020, 515 (2020). 2020.
Roul, P.; Rohil, V. A novel high-order numerical scheme and its analysis for the two-dimensional time-fractional reaction-subdiffusion equation. Numer. Algor. 2022, 90, 1357–1387. [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]
Shams, M.; Kausar, N.; Samaniego, C.; Agarwal, P.; Ahmed, S.F.; Momani, S. On efficient fractional Caputo-type simultaneous scheme for finding all roots of polynomial equations with biomedical engineering applications. Fractals 2023, 31(4), 2340075. [Google Scholar] [CrossRef]
Singh, H.; Kumar, D.; Baleanu, D. Methods of Mathematical Modelling "Fractional Differential Equations", CRC Press, Taylor & Francis, 2019;
Song, L.; Wang, W. ; Solution of the fractional Black-Scholes option pricing model by finite difference method. Abstr. Appl. Anal. 2013, 194286. [Google Scholar] [CrossRef]
Sun, D.; Liu, J.; Su, X.; Pei, G. Fractional differential equation modeling of the HBV infection with time delay and logistic proliferation. Front. Public Health 2022, 10, 1036901. [Google Scholar] [CrossRef] [PubMed]
Wang, X.-T. Scaling and long-range dependence in option pricing I: Pricing European option with transaction costs under the fractional Black–Scholes model. Phys. A: Stat. Mech. 2010, 389, 438–444. [Google Scholar] [CrossRef]
Wang, Y.-M.; Ren, L. A high-order L2-compact difference method for Caputo-type time-fractional sub-diffusion equations with variable coefficients, Appl. Math. Comput. 2019, 342, 71–93. [Google Scholar]
Weideman, J.A.C. Numerical integration of periodic functions: a few examples. Am. Math. Mon. 2002, 109(1), 21–36. [Google Scholar] [CrossRef]
Wu, R.F.; Ding, H.F.; Li, C.P. Determination of coefficients of high-order schemes for Riemann-Liouville derivative. Sci. World J. 2014, 402373. [Google Scholar] [CrossRef] [PubMed]
Zeng, F.; Li, C.; Liu, F.; Turner, I. Numerical algorithms for time-fractional subdiffusion equation with second-order accuracy. SIAM J. Sci. Comput. 2015, A55–A78. [Google Scholar] [CrossRef]
Zhang, H.; Liu, F.; Turner, I.; Yang, Q. Numerical solution of the time fractional Black–Scholes model governing European options, Comput. Math. Appl. 2016, 71(9), 1772–1783. [Google Scholar]

Figure 1. Option premia for various

α

(left) and difference between them and the option premium for integer-order model (right).

Figure 1. Option premia for various

α

(left) and difference between them and the option premium for integer-order model (right).

Figure 2. Call option premium for all t.

Figure 3. Put option premium for all t.

Table 1. Error and order of numerical solution (23) of equation (26) in

[0, 1]

Table 1. Error and order of numerical solution (23) of equation (26) in

[0, 1]

$h$	$α = 0.1, D = 1$		$α = 0.3, D = 100$		$α = 0.5, D = - 1$
$h$	$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.005$	$5.1055 \times 10^{- 9}$	$1.7046$	$1.0160 \times 10^{- 8}$	$1.2454$	$6.0036 \times 10^{- 5}$	$1.5741$
$0.0025$	$1.3680 \times 10^{- 9}$	$1.7130$	$3.0907 \times 10^{- 9}$	$1.2473$	$3.4169 \times 10^{- 5}$	$1.5808$
$0.00125$	$3.6656 \times 10^{- 10}$	$1.7196$	$9.3780 \times 10^{- 10}$	$1.2484$	$2.1249 \times 10^{- 6}$	$1.5857$

Table 2. Error and order of numerical solution (23) of equation (26) in

[0, 2]

Table 2. Error and order of numerical solution (23) of equation (26) in

[0, 2]

$h$	$α = 0.7, D = - 10$		$α = 0.9, D = - 20$		$α = 0.5, D = - 1000$
$h$	$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.005$	$2.0798 \times 10^{7}$	$3.1743$	$3.6733 \times 10^{8}$	$4.2249$	$1.6892 \times 10^{- 8}$	$1.4931$
$0.0025$	$5.1807 \times 10^{6}$	$2.0052$	$7.0672 \times 10^{7}$	$2.3778$	$6.0126 \times 10^{- 9}$	$1.4902$
$0.00125$	$1.7366 \times 10^{6}$	$1.5768$	$2.2313 \times 10^{7}$	$1.6632$	$2.1463 \times 10^{- 9}$	$1.4861$

Table 3. Error and order of numerical solution (23) of equation (27) in

[0, 1]

Table 3. Error and order of numerical solution (23) of equation (27) in

[0, 1]

$h$	$α = 0.2, D = 1$		$α = 0.4, D = 100$		$α = 0.6, D = - 0.1$
$h$	$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.005$	$9.6187 \times 10^{- 8}$	$1.7998$	$6.3486 \times 10^{- 8}$	$1.5989$	$2.1642 \times 10^{- 4}$	$1.3975$
$0.0025$	$2.7624 \times 10^{- 8}$	$1.7999$	$2.0950 \times 10^{- 8}$	$1.5994$	$8.2091 \times 10^{- 5}$	$1.3985$
$0.00125$	$7.9333 \times 10^{- 9}$	$1.7999$	$6.9123 \times 10^{- 9}$	$1.5997$	$3.1124 \times 10^{- 5}$	$1.3992$

Table 4. Error and order of numerical solution (23) of equation (27) in

[0, 2]

Table 4. Error and order of numerical solution (23) of equation (27) in

[0, 2]

$h$	$α = 0.8, D = - 20$		$α = 0.75, D = - 50$		$α = 0.75, D = - 1000$
$h$	$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.005$	$7.0524 \times 10^{34}$	$15.733$	$4.6487 \times 10^{46}$	$27.426$	$4.2435 \times 10^{- 8}$	$1.4986$
$0.0025$	$7.9788 \times 10^{32}$	$6.4657$	$4.9857 \times 10^{43}$	$9.8648$	$1.5010 \times 10^{- 8}$	$1.4993$
$0.00125$	$7.8671 \times 10^{31}$	$3.3422$	$2.0753 \times 10^{42}$	$4.5863$	$5.3082 \times 10^{- 9}$	$1.4996$

Table 5. Error and order of the numerical solution of equation (43) and

N = 100

Table 5. Error and order of the numerical solution of equation (43) and

N = 100

$τ$	$α = 0.25$		$α = 0.5$		$α = 0.75$
$τ$	$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 / 10$	$8.1169 \times 10^{- 5}$		$2.0202 \times 10^{- 3}$		$1.9422 \times 10^{- 2}$
$1 / 20$	$2.4655 \times 10^{- 5}$	$1.7190$	$7.4716 \times 10^{- 4}$	$1.4350$	$8.6786 \times 10^{- 3}$	$1.1622$
$1 / 40$	$7.3993 \times 10^{- 6}$	$1.7364$	$2.6940 \times 10^{- 4}$	$1.4716$	$3.7472 \times 10^{- 3}$	$1.2116$
$1 / 80$	$2.2091 \times 10^{- 6}$	$1.7439$	$9.6139 \times 10^{- 5}$	$1.4865$	$1.5955 \times 10^{- 3}$	$1.2318$
$1 / 160$	$6.5748 \times 10^{- 7}$	$1.7483$	$3.4152 \times 10^{- 5}$	$1.4931$	$6.7505 \times 10^{- 4}$	$1.2409$

Table 6. Error and order of the numerical solution of equation (43) and

M = 1000

Table 6. Error and order of the numerical solution of equation (43) and

M = 1000

h	$α = 0.2$		$α = 0.5$		$α = 0.8$
h	$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 / 4$	$5.4801 \times 10^{- 4}$		$5.3179 \times 10^{- 4}$		$6.5937 \times 10^{- 4}$
$1 / 8$	$1.4906 \times 10^{- 4}$	$1.8782$	$1.5812 \times 10^{- 4}$	$1.7497$	$2.1445 \times 10^{- 4}$	$1.6204$
$1 / 16$	$3.8174 \times 10^{- 5}$	$1.9653$	$4.4648 \times 10^{- 5}$	$1.8244$	$6.6887 \times 10^{- 5}$	$1.6809$
$1 / 32$	$9.6285 \times 10^{- 6}$	$1.9872$	$1.2065 \times 10^{- 5}$	$1.8877$	$1.9485 \times 10^{- 5}$	$1.7793$
$1 / 64$	$2.4210 \times 10^{- 6}$	$1.9917$	$3.0974 \times 10^{- 6}$	$1.9617$	$5.3211 \times 10^{- 6}$	$1.8726$

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.

© 2023 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

Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations

Abstract

1. Introduction

2. Properties of the Weights

3. Estimate for the Error

4. Numerical Solution of Two-Term Equation

5. Time Fractional Black–Scholes Equation

6. Conclusions

Author Contributions

Funding

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe