Next Article in Journal
Stochastic Modeling of Three-Species Prey–Predator Model Driven by Lévy Jump with Mixed Holling-II and Beddington–DeAngelis Functional Responses
Next Article in Special Issue
Some Results on Fractional Boundary Value Problem for Caputo-Hadamard Fractional Impulsive Integro Differential Equations
Previous Article in Journal
Distributed Adaptive Optimization Algorithm for Fractional High-Order Multiagent Systems Based on Event-Triggered Strategy and Input Quantization
Previous Article in Special Issue
On Stability of a Fractional Discrete Reaction–Diffusion Epidemic Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Mathematics, Physics and Informatics, University of Forestry, 1756 Sofia, Bulgaria
2
Department of Informational Modeling, Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria
3
Department of Applied Mathematics and Statistics, University of Ruse, 7017 Ruse, Bulgaria
4
Department of Parallel Algorithms, Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria
*
Author to whom correspondence should be addressed.
Fractal Fract. 2023, 7(10), 750; https://doi.org/10.3390/fractalfract7100750
Submission received: 4 September 2023 / Revised: 28 September 2023 / Accepted: 2 October 2023 / Published: 11 October 2023
(This article belongs to the Special Issue Advances in Fractional Modeling and Computation)

Abstract

:
In this paper, we consider an approximation of the Caputo fractional derivative and its asymptotic expansion formula, whose generating function is the polylogarithm function. We prove the convergence of the approximation and derive an estimate for the error and order. The approximation is applied for the construction of finite difference schemes for 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 and order of the finite difference schemes and to obtain bounds for the error of the numerical methods. The theoretical results for the order and error of the methods are illustrated by the results of the numerical experiments.

1. Introduction

Fractional differential equations are an essential means of modeling complex processes in different fields of science [1,2,3,4,5]. The Riemann–Liouville fractional integral of order α > 0 is defined as as follows:
a I x α f ( x ) = 1 Γ ( α ) 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 follows:
a C D x n + α f ( x ) = I 1 α D n + 1 f ( x ) = 1 Γ ( 1 α ) a x f ( n + 1 ) ( t ) ( x t ) α d t , a R L D x n + α f ( x ) = D n + 1 I 1 α f ( x ) = 1 Γ ( 1 α ) d n + 1 d x n + 1 a x f ( t ) ( x t ) α d t .
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 ) + k = 0 n f ( k ) ( 0 ) x k n α Γ ( k α n + 1 ) .
In this paper, we considered the case when the order of fractional differentiation is between zero and one. In the remainder of the paper, the parameter α ( 0 , 1 ) . Without loss of generality, the lower limit of the integral in the definition of the fractional derivative is zero. The Caputo fractional derivative of order α is defined as follows:
f ( α ) ( x ) = D x α f ( x ) = 1 Γ ( 1 α ) 0 x 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–Leffler function:
E β 1 , β 2 ( x ) = k = 0 x k Γ ( β 1 k + β 2 ) .
Approximations of fractional derivatives involve the use of special functions. The Riemann zeta function is an analytic function that is defined for complex numbers with the real part greater than one as follows:
ζ ( x ) = n = 1 1 n x .
The polylogarithm function of order β is defined as follows:
Li β ( x ) = n = 1 x n n β
and Li β ( e x ) has a series expansion:
Li β ( e x ) = Γ ( 1 β ) x β 1 + n = 0 ζ ( β n ) n ! ( x ) n .
Finite difference schemes for the 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 the L1 approximation [6,7,8,9,10]. The 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 of integer-order derivatives [11,12].
Let n be a positive integer, h = x / n , and y k = y ( k h ) for k = 0 , 1 , , n . The L1 approximation has an order 2 α and is defined as follows:
h α Γ ( 2 α ) k = 0 m σ k ( α ) f n k = f n ( α ) + 𝒪 h 2 α ,
where σ 0 ( α ) = 1 , σ n ( α ) = ( n 1 ) 1 α n 1 α and
σ k ( α ) = ( k 1 ) 1 α k 1 α + ( k + 1 ) 1 α , ( k = 1 , , n 1 ) .
The L1 approximation has a generating function ( x 1 ) 2 Li α 1 ( x ) x Γ ( 2 α ) and a second-order expansion formula [13]:
h α Γ ( 2 α ) k = 0 m σ k ( α ) f n k = f n ( α ) + ζ ( α 1 ) Γ ( 2 α ) f n h 2 α + 𝒪 h 2 ,
The construction of the L1 approximation uses central difference approximations of the first derivative on a uniform net. Over the past few decades, there has been an increasing interest in developing efficient numerical techniques for solving fractional differential equations. The growing interest in this field has resulted in the development of high-order approximations for fractional derivatives and innovative methods for investigating their properties and practical applications. Approximations of the fractional derivative of order 3 α , which are called the L 1 2 and L 1 2 σ formulas, were constructed by Alikhanov [14] and Gao et al. [15]. In [13,16], we used the L1 approximation and approximations of the second derivative for the construction of second-order approximations of the fractional derivative. Finite difference schemes using L2 formula approximations of the fractional derivative of order 3 α and their convergence were studied by Alikhanov [17], Lv and Xu [18], and Wong and Ren [19]. High-order approximations of the fractional derivative and difference schemes for fractional differential equations were constructed in [20,21,22,23,24,25,26,27]. Navot [28] used Taylor polynomials to deal with the singularity of the fractional integral 0 1 t α f ( t ) d t . He derived the asymptotic formula for Riemann sums on a uniform net, called the extended Euler–Maclaurin summation formula. The formula was extended by Navot [29] to functions with a logarithmic singularity. In [30], we derived the asymptotic formula of the Riemann sum of the fractional integral using the series expansion Formula (2) of the generating function Li 1 α x . Many of the techniques, originally formulated for ordinary and partial differential equations, have been modified and extended to address fractional differential equations. Spectral methods employ orthogonally based functions such as the Jacobi and Chebyshev polynomials to approximate solutions of fractional differential equations [31,32]. Fractional multistep methods, which can be viewed as an extension of classical multistep methods, were studied in [33,34]. Methods using spline interpolation were presented in [35,36]. Numerical methods using fractional wavelets were explored in [37,38].
This paper is a continuation of the work in [39]. In the paper, we study the properties of an approximation of order 2 α of the Caputo derivative and the convergence of the corresponding difference schemes for the two-term ordinary fractional differential equation and the time fractional Black–Scholes equation for options pricing. A main approach for constructing approximations of fractional derivatives involves interpolating the function f in the definitions of fractional derivatives over a suitable stencil on a uniform grid. Using this method, formulas for L1 and L2 approximations are derived, as well as high-order approximations of fractional derivatives. The approximations obtained in this way have generating functions that include the polylogarithm function, and their expansion formulas are derived using (2). The direct application of this method leads to the construction of approximations of fractional derivatives with weights whose exponents are greater than 1 , as the fractional integral of a negative order is divergent. Naturally, the question arises of constructing and studying the properties of approximations whose weights contain terms with exponents smaller than 1 . A direct application of (2) for the function Li 1 + α ( e x ) leads to the construction of an approximation of the fractional derivative and its expansion formula, which has weights of 1 / k 1 + α . The construction of this approximation and the corresponding approximations of order 2 α were presented in [39]. The comparison performed with the L1 approximation shows that these approximations have the same order, similar properties of the weights, and comparable performance for the numerical solution of fractional differential equations. In this paper, we derive an estimate for the error of the approximation of order 2 α of the fractional derivative with a generating function Li 1 + α ( x ) and the corresponding difference schemes for the numerical solution of fractional differential equations. The motivation of the study this approximation stems from the simpler formulas of its weights, its mathematical properties and applications, and for constructing efficient methods for numerically solving fractional differential equations.
In [39], we constructed an approximation of the Caputo fractional derivative and its asymptotic formula of order 2 α :
1 h α k = 1 n 1 f n k k 1 + α = Γ ( α ) f n ( α ) + ζ ( 1 + α ) h α f n ζ ( α ) f n h 1 α + 𝒪 h 2 α .
Approximation (5) has a generating function Li 1 + α ( x ) . By substituting the derivative in the right-hand side of (1) using the first-order backward difference ( f n f n 1 ) / h , we find an approximation of the fractional derivative:
L n f ( x ) = h α Γ ( α ) k = 0 m σ k ( α ) f n k = f n ( α ) + 𝒪 h 2 α ,
σ 0 ( α ) = ζ ( α ) ζ ( 1 + α ) , σ 1 ( α ) = 1 ζ ( α ) , σ k ( α ) = 1 k 1 + α , ( 2 k n ) .
Approximations (5) and (6) have an order 2 α when the function f satisfies the condition f ( 0 ) = f ( 0 ) = 0 . In [39], we extended the approximation (6) to all functions in the class C 2 [ 0 , x ] by assigning the values of the last two weights.
A n f ( x ) = h α Γ ( α ) k = 0 n σ k ( α ) f n k = f n ( α ) + 𝒪 h 2 α ,
σ 0 ( α ) = ζ ( α ) ζ ( 1 + α ) , σ 1 ( α ) = 1 ζ ( α ) , σ k ( α ) = 1 k 1 + α , ( 2 k n 2 ) ,
σ n 1 ( α ) = 1 ( n 1 ) 1 + α n S n [ 1 + α ] + S n [ α ] n 1 α α ( 1 α ) ,
σ n ( α ) = ( n 1 ) S n [ 1 + α ] S n [ α ] + n 1 α α ( 1 α ) ,
where S n [ α ] = k = 1 n 1 k α ζ ( α ) . In this paper, we investigate the properties of the approximations (6) and (7) and applications of the approximations for the 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 the approximation (7) and we derive the following inequality:
k = 1 n 1 ( n k ) α k 1 + α < ζ ( 1 + α ) n α 2 .
In Section 3, we derive estimates for the errors of the approximations (6) and (7). In Section 4 and Section 5, we construct finite difference schemes for numerically solving the two-term ordinary fractional differential equation and the time fractional Black–Scholes equation for option pricing. The convergence of the finite difference schemes is proven, and bounds for the errors of the methods are derived using Inequality (11) and the estimate for the truncation error of the approximation.
The approximation (5) is constructed using the series expansion Formula (2) of the function Li 1 + α ( x ) and can be extended to arbitrary orders. Equation (2) holds for any value of β that is not an integer. This indicates that the approximations (5)–(7) also apply to fractional derivatives of order greater than one. One direction for future work related to the paper is to consider applications of the discussed approximations for numerically solving fractional differential equations of orders greater than one. The L1 approximation (4) and the approximation (7) have an order 2 α and similar performance and properties of the weights [39]:
σ 0 ( α ) > 0 , σ 1 ( α ) < σ 2 ( α ) < < σ n 1 ( α ) < 0 .
The properties of the L1 approximation and the approximation (7) presented in Section 2 enable an efficient analysis of the convergence and error of difference schemes for fractional differential equations. One approach for constructing approximations of integer-order and fractional derivatives was presented in [16]. Another direction for future work is constructing high-order approximations for fractional derivatives that have the properties (12) of the weights, investigating the properties of the constructed high-order approximations, and analyzing the convergence of the difference schemes for fractional differential equations. Additionally, it is important to examine the behavior of the numerical solution for values of the parameters where the denominator has a small value.

2. Properties of the Weights

Approximation (7) has the property that the values of A n 1 and A n x are equal to the fractional derivatives of the functions:
A n 1 = D x α 1 = 0 , A n x = D x α x = x 1 α Γ ( 2 α ) .
The values of the weights of σ n 1 ( α ) and σ n ( α ) are computed from the system of equations:
k = 0 n σ k ( α ) = 0 , k = 0 n k σ k ( α ) = n 1 α Γ ( 2 α ) .
In this section, we prove Inequality (11) and we derive an estimate for the weight σ n 1 ( α ) of Approximation (6). The proof uses the inequalities in Claims 1 and 2.
Claim 1.
Let 0 < α < 1 and 1 k n . Then,
( n k ) α < n α α k n 1 α .
Proof. 
From the binomial formula,
( n k ) α = n α 1 k n α = n α i = 0 ( 1 ) i α i k n i .
The numbers ( 1 ) i α i are negative for i 1 . Then,
( n k ) α n α + α k n 1 α = n α i = 2 ( 1 ) i α i k n i < 0 .
Claim 2.
Let 0 < α < 1 . Then,
1 α + α 2 α ( 1 α ) + α ζ ( α ) 2 1 α > 2 .
Proof. 
The Riemann zeta function has a series expansion [40]:
ζ ( α ) = 1 α 1 + γ + k = 1 γ n n ! ( 1 α ) n ,
where γ is an Euler constant and γ n are generalized Euler constants:
γ n = lim N k = 1 N log n k k log n + 1 k n + 1 ,
γ = 0.5772157 , γ 1 = 0.0728158 , γ 2 = 0.0096904 , γ 3 = 0.0020538 .
Then,
ζ ( α ) > 1 α 1 + γ + ( 1 α ) γ 1 + ( 1 α ) 2 2 γ 2 > 1 α 1 + 0.49 + 0.08 α .
The value of ln 2 = 0.693147 < 0.7 . From Taylor’s Theorem,
2 ( 1 α ) = e ( 1 α ) ln 2 < 1 ( 1 α ) ln 2 + ( 1 α ) 2 ln 2 2 2 < 1 ln 2 ln 2 2 2 ( 1 α ) < 1 0.45 ( 1 α ) = 0.55 + 0.45 α .
From (16) and (17),
1 α + α 2 α ( 1 α ) + α ζ ( α ) 2 1 α > 1 α + α 2 α ( 1 α ) + α 1 α 1 + 0.497 + 0.08 α ( 0.55 + 0.45 α ) > 1 α + 0.036 a 3 + 0.26765 a 2 + 0.72335 α .
Inequality (14) holds when
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 ] . □
Lemma 1.
Let n 2 . Then
k = 1 n 1 ( n k ) α k 1 + α < ζ ( 1 + α ) n α 2 .
Proof. 
From Inequality (13),
k = 1 n 1 ( n k ) α k 1 + α < k = 1 n 1 1 k 1 + α n α α k n 1 α < n α k = 1 n 1 1 k 1 + α α n 1 α k = 1 n 1 1 k α .
The sum of the α -th powers of the first n 1 integers has an asymptotic formula [41]:
k = 1 n 1 k α = ζ ( α ) + n 1 + α 1 + α p = 1 1 + α p b p n p ,
where b p are the Bernoulli numbers. Equation (19) is derived using the Euler–Maclaurin formula, and the actual value of k = 1 n 1 k α lies between two consecutive partial sums [41]. Then,
k = 1 n 1 1 k 1 + α < ζ ( 1 + α ) 1 α n α 1 2 n 1 + α 1 + α 12 n 2 + α + ( 1 + α ) ( 2 + α ) ( 3 + α ) 720 n 4 + α ,
k = 1 n 1 1 k α > n 1 α 1 α + ζ ( α ) 1 2 n α α 12 n 1 + α
and
n α k = 1 n 1 1 k 1 + α < ζ ( 1 + α ) n α 1 α 1 2 n 1 + α 12 n 2 + 1 30 n 4 , α n 1 α k = 1 n 1 1 k α > α 1 α + α ζ ( α ) n 1 α α 2 n α 2 12 n 2 .
Therefore,
k = 1 n 1 ( n k ) α k 1 + α < n α k = 1 n 1 1 k 1 + α α n 1 α k = 1 n 1 1 k α < ζ ( 1 + α ) n α 1 α + α 2 α ( 1 α ) α ζ ( α ) n 1 α 1 α 2 n 1 + α α 2 12 n 2 + 1 30 n 4 < ζ ( 1 + α ) n α 1 α + α 2 α ( 1 α ) α ζ ( α ) n 1 α 1 12 n 2 + 1 30 n 4 < ζ ( 1 + α ) n α 1 α + α 2 α ( 1 α ) α ζ ( α ) n 1 α .
From Claim 2, it follows that
k = 1 n 1 ( n k ) α k 1 + α < ζ ( 1 + α ) n α 2 .
Now, we derive an estimate for the weight σ n 1 ( α ) of Approximation (7). Denote
σ ¯ n 1 ( α ) = σ n 1 ( α ) 1 ( n 1 ) 1 + α = S n [ α ] n S n [ 1 + α ] n 1 α α ( 1 α ) .
Claim 3.
Let n 2 . Then
σ ¯ n 1 ( α ) < 1 + 2 α 12 n 1 + α .
Proof. 
From (19),
S n [ α ] = n 1 α 1 α 1 2 n α + R 1 , n S n [ α + 1 ] = n 1 α α 1 2 n α + R 2 ,
where
| R 1 | < 1 + α 12 n 1 + α , | R 2 | < α 12 n 1 + α .
Then,
σ ¯ n 1 ( α ) = S n [ α ] n S n [ 1 + α ] n 1 α α ( 1 α ) < | R 1 | + | R 2 | < 1 + 2 α 12 n 1 + α

3. Estimate for the Error

In this section, we use the method from [28] to derive estimates for the errors of Approximations (6) and (7). Let 0 < t < x and
g ( t ) = f ( t ) f ( x ) + ( x t ) f ( x ) ( x t ) 2 2 f ( x ) .
The function g ( t ) satisfies g ( x ) = g ( x ) = g ( x ) = 0 . From Taylor’s Theorem,
g ( t ) = ( t x ) 3 6 f ( ξ 0 , t ) , g ( t ) = ( 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 ) | < M 3 6 ( x t ) 3 , | g ( t ) | < 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 is not defined at the point x.
Claim 4.
Let 0 < t < x . Then,
| G ( t ) | < 5 6 M 3 ( x t ) 1 α , | G ( t ) | < 4 M 3 ( x t ) α .
Proof. 
G ( t ) = ( x t ) g ( x ) + ( 1 + α ) g ( t ) ( x t ) 2 + α ,
G ( t ) = ( x t ) 2 g ( t ) + 2 ( α + 1 ) ( x t ) g ( t ) + ( α + 1 ) ( α + 2 ) g ( t ) ( x t ) 3 + α .
Therefore,
| G ( t ) | < 1 2 M 3 ( x t ) 1 α + 1 + α 6 M 3 ( x t ) 1 α < 5 6 M 3 ( x t ) 1 α , | G ( t ) | < ( 2 + α ) ( 1 + α ) M 3 6 + ( 1 + α ) M 3 + M 3 ( x t ) α < 4 M 3 ( x t ) α .
Lemma 2.
Let f ( 0 ) = f ( 0 ) = 0 . Then,
0 x G ( t ) d t = Γ ( α ) f ( α ) ( x ) + f ( x ) α x α + f ( x ) x 1 α 1 α f ( x ) x 2 α 2 ( 2 α ) .
Proof. 
0 x G ( t ) d t = 0 x g ( t ) ( x t ) 1 + α d t = 0 x g ( t ) d ( x t ) α α = g ( t ) ( x t ) α α 0 x 1 a 0 x g ( t ) ( x t ) α d t = g ( 0 ) α x α + Γ ( α ) Γ ( 1 α ) 0 x g ( t ) ( x t ) α d t .
The function g ( t ) satisfies
g ( 0 ) = f ( x ) + f ( x ) x f ( x ) 2 x 2 .
Hence,
0 x g ( t ) ( x t ) 1 + α d t = Γ ( α ) g ( α ) ( x ) + 1 α f ( x ) x α x 1 α f ( x ) + x 2 α 2 f ( x ) .
In view of Equation (20), the fractional derivative of the function g ( x ) satisfies
g ( α ) ( x ) = f ( α ) ( x ) x 1 α f ( x ) ( 1 α ) Γ ( 1 α ) + x 2 α f ( x ) ( 2 α ) Γ ( 1 α ) ,
which can be written in the form:
Γ ( α ) g ( α ) ( x ) = Γ ( α ) f ( α ) ( x ) + x 1 α f ( x ) α ( 1 α ) x 2 α f ( x ) α ( 2 α ) .
From (21) and (22), it follows that
0 x g ( t ) ( x t ) 1 + α d t = Γ ( α ) f ( α ) ( x ) + f ( x ) α x α + f ( x ) x 1 α 1 α f ( x ) x 2 α 2 ( 2 α ) .
The trapezoidal rule of a function F ( t ) has a second-order accuracy when F C 2 [ 0 , x ] . The error of the trapezoidal rule E T F satisfies [42]
| E T F | < 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 be the error of the trapezoidal rule of G ( t ) in the interval [ 0 , x ] .
Lemma 3.
| T n | 1 3 M 3 x ( x 2 + 1 ) h 2 α .
Proof. 
Let T n 1 be the of trapezoidal rule of G ( t ) in the interval [ 0 , x h ] .
T n 1 1 12 h 2 ( x h ) 3 max 0 < t < x h | G ( t ) | 1 12 h 2 x 3 4 M 3 h α 1 3 M 3 x 3 h 2 α .
Let T n 2 be the error of the trapezoidal rule of G ( t ) in the interval [ x h , x ] .
x h x G ( t ) d t = x h x G ( t ) d t + h 2 x = h 2 G ( x h ) x h x t + h 2 x G ( t ) d t .
The error T n 2 satisfies
T n 2 = x h x t + h 2 x G ( t ) d t ,
| T n 2 | < x h x t + h 2 x | G ( t ) | d t < 5 6 M 3 h 1 α x h x t + h 2 x d t < 5 24 M 3 h 3 α .
From Equations (23) and (24), we obtain a bound for the error T n = T n 1 + T n 2 of the function G ( t ) in the interval [ 0 , x ] :
| T n | | T n 1 | + | T n 2 | < 1 3 M 3 x 3 h 2 α + 5 24 M 3 h 3 α < 1 3 M 3 x 3 + h h 2 α < 1 3 M 3 x ( x 2 + 1 ) h 2 α .
The trapezoidal rule of the function G ( t ) satisfies
1 h α k = 1 n 1 g n k k 1 + α + g 0 2 n 1 + α = 0 x g ( t ) ( x t ) 1 + α d t + T n
The function g ( t ) has a value at zero:
g ( 0 ) 2 n 1 + α h α = f ( x ) 2 n 1 + α h α + f ( x ) x 2 n 1 + α h α f ( x ) x 2 4 n 1 + α h α = h 2 f ( x ) x 1 + α + f ( x ) x α f ( x ) x 1 α 2 .
Therefore,
1 h α k = 1 n 1 g n k k 1 + α = Γ ( α ) f ( α ) ( x ) + f ( x ) α x α + f ( x ) h 2 x 1 + α + x 1 α f ( x ) 1 α f ( x ) h 2 x α + x 1 α f ( x ) h 4 x 2 α f ( x ) 2 ( 2 α ) + T n .
Let E n be the error of Approximation (5).
1 h α k = 1 n 1 f n k k 1 + α = Γ ( α ) f n ( α ) + ζ ( 1 + α ) h α f n ζ ( α ) f n h 1 α + E n .
Theorem 1.
Let f ( 0 ) = f ( 0 ) = 0 and f C 3 [ 0 , x n ] . Then,
| E n | < 2 M 3 x ( x 2 + 1 ) + M 2 6 h 2 α .
Proof. 
From the definition (20) of the function g ( t ) ,
1 h α k = 1 n 1 g n k k 1 + α = 1 h α k = 1 n 1 f n k k 1 + α f ( x ) h α k = 1 n 1 1 k 1 + α + f ( x ) h 1 α k = 1 n 1 1 k α f ( x ) 2 h 2 α k = 1 n 1 k 1 α .
From Equation (19), we find
k = 1 n 1 1 k 1 + α = ζ ( 1 + α ) 1 α n α 1 2 n 1 + α + R 1 , k = 1 n 1 1 k α = n 1 α 1 α + ζ ( α ) 1 2 n α + R 2 ,
k = 1 n 1 k 1 α = n 2 α 2 α n 1 α 2 + ζ ( α 1 ) + R 3 ,
where
| R 1 | < 1 + α 12 n 2 + α , | R 2 | < α 12 n 1 + α , | R 3 | < 1 α 12 n α .
Then,
f ( x ) h α k = 1 n 1 1 k 1 + α = ζ ( 1 + α ) f ( x ) h α f ( x ) α n α h α f ( x ) 2 n 1 + α h α + R 1 f ( x ) h α = ζ ( 1 + α ) f ( x ) h α f ( x ) α x α f ( x ) h 2 x 1 + α + R 1 f ( x ) h α .
Similarly,
f ( x ) h 1 α k = 1 n 1 1 k α = f ( x ) h 1 α ζ ( α ) + f ( x ) x 1 α 1 α f ( x ) h 2 x α + f ( x ) h 1 α R 2 ,
f ( x ) 2 h 2 α k = 1 n 1 k 1 α = ζ ( α 1 ) f ( x ) 2 h 2 α + f ( x ) x 2 α 2 α f ( x ) x 1 α h 4 + R 3 f ( x ) 2 h 2 α .
Hence,
1 h α k = 1 n 1 f n k k 1 + α = Γ ( α ) f ( α ) ( x ) + f ( x ) h α ζ ( 1 + α ) f ( x ) h 1 α ζ ( α ) + f ( x ) 2 h 2 α ζ ( α 1 ) + E n ,
where
E n = T n + R 1 f ( x ) h α f ( x ) h 1 α R 2 + R 3 f ( x ) 2 h 2 α ,
| E n | < | T n | + | R 1 f ( x ) h α | + | f ( x ) h 1 α R 2 | + | R 3 f ( x ) 2 h 2 α | .
From Taylor’s Theorem,
f ( x ) = 1 2 f ( ξ 1 ) x 2 , f ( x ) = f ( ξ 2 ) x ,
where 0 < ξ 2 < x for i = 1 , 2 . Therefore,
M 0 1 2 M 2 x 2 , M 1 M 2 x .
The terms on the right-hand side of (25) satisfy the estimates:
| R 1 f ( x ) h α | < ( 1 + α ) M 0 h α 12 n 2 + α = ( 1 + α ) M 0 h 2 12 x 2 + α < ( 1 + α ) M 2 h 2 α 24 , | R 2 f ( x ) h 1 α | < α M 1 h 1 α 12 m 1 + α = α M 1 h 2 12 x 1 + α < α M 2 h 2 12 x α < α M 2 h 2 α 12 , | R 3 f ( x ) 2 h 2 α | < ( 1 α ) M 2 h 2 α 24 n α = ( 1 α ) M 2 h 2 24 x α < ( 1 α ) M 2 h 2 α 24 .
Hence,
| E n | < 1 3 M 3 x ( x 2 + 1 ) h 2 α + ( 1 + α ) M 2 h 2 α 12 < 2 M 3 x ( x 2 + 1 ) + M 2 6 h 2 α .
Let L n be the error of Approximation (6).
L n f ( x ) = h α Γ ( α ) k = 1 m σ k ( α ) f n k = f n ( α ) + L n ,
σ 0 ( α ) = ζ ( α ) ζ ( 1 + α ) , σ 1 ( α ) = 1 ζ ( α ) , σ k ( α ) = 1 k 1 + α , ( 2 k n ) .
Approximation (6) is constructed from Approximation (5) by substituting f n with first-order backward difference approximation ( f n f n 1 ) / h .
Claim 5.
Let f C 3 [ 0 , x ] and f ( 0 ) = f ( 0 ) = 0 . Then,
| L n | < 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 + h 2 2 f ( ξ ) ,
f n f n f n 1 h < h 2 f ( ξ ) ,
where ξ ( x n 1 , x n ) . The error L n of Approximation (6) satisfies
L n = 1 Γ ( α ) E n 1 2 ζ ( α ) f ( ξ ) h 2 α ,
| L n | < 1 | Γ ( α ) | 2 M 3 x ( x 2 + 1 ) + M 2 6 h 2 α + | ζ ( α ) | M 2 2 h 2 α < 2 M 3 x ( x 2 + 1 ) + ( 1 + 3 | ζ ( α ) | ) M 2 6 | Γ ( α ) | h 2 α .
Let A n be the error of Approximation (7).
A n f ( x ) = h α Γ ( α ) k = 1 n σ k ( α ) f n k = f n ( α ) + A n .
The weights of Approximation (7) are defined with (8)–(10). Denote
f ˜ ( t ) = f ( t ) f ( 0 ) f ( 0 ) t .
The function f ˜ satisfies f ˜ ( 0 ) = f ˜ ( 0 ) = 0 , and the second and third derivatives of f and f ˜ are equal. Therefore,
max 0 t x | f ˜ ( t ) | = M 2 , max 0 t x | f ˜ ( t ) | = M 3 .
Lemma 4.
Let f C 3 [ 0 , x ] . Then,
| A n | < 4 M 3 x ( x 2 + 1 ) + 3 ( 1 + 2 | ζ ( α ) | ) M 2 12 | Γ ( α ) | h 2 α .
Proof. 
A n f ( x ) = A n f ˜ ( x ) + A n [ f ( x ) f ˜ ( x ) ] = L n f ˜ ( x ) + σ ¯ n 1 ( α ) f ˜ 1 Γ ( α ) h α + D x α [ f ( x ) f ˜ ( x ) ] .
From Claim 3,
σ ¯ n 1 ( α ) < 1 + 2 α 12 n 1 + α < 1 4 n < 1 8 .
From Taylor’s Theorem,
| f ˜ 1 | = h 2 2 | f ˜ ( ξ ) | h 2 M 2 2 .
Then,
| A n | < | L n | + 1 | Γ ( α ) | h α σ ¯ n 1 ( α ) | f ˜ 1 | < 2 M 3 x ( x 2 + 1 ) + ( 1 + 3 | ζ ( α ) | ) M 2 6 | Γ ( α ) | + M 2 16 Γ ( α ) h 2 α < 4 M 3 x ( x 2 + 1 ) + 3 ( 1 + 2 | ζ ( α ) | ) M 2 12 | Γ ( α ) | h 2 α .

4. Numerical Solution of Two-Term Equation

The two-term equation is an ordinary fractional differential equation in the form:
y ( α ) ( t ) + D y ( t ) = F ( t ) , y ( 0 ) = y 0 .
Numerical and analytical solutions of ordinary fractional differential equations were studied in [43,44,45,46,47,48]. In this section, we construct the numerical solution of the two-term Equation (26), which uses the approximation (7) of the fractional derivative, and prove its convergence and order. Let h = x / N and x n = n h for n = 0 , , N . By approximating the fractional derivative at the point x n , we find
1 Γ ( α ) h α k = 0 n σ k ( α ) y n k + D y n = F n + A n ,
σ 0 ( α ) + Γ ( α ) D h α y n + k = 1 n σ k ( α ) y n k + D y n = Γ ( α ) h α F n + B n h 2 ,
where B n = Γ ( α ) A n h α 2 .
ζ ( α ) ζ ( 1 + α ) + Γ ( α ) D h α y n ζ ( α ) y n 1 + k = 1 n y n k k 1 + α σ ¯ n 1 ( α ) u 1 = Γ ( α ) h α F n + B n h 2 .
The numerical solution { u n } n = 2 N of Equation (26) is computed as
ζ ( α ) ζ ( 1 + α ) + Γ ( α ) D h α u n ζ ( α ) u n 1 + k = 1 n u n k k 1 + α + σ ¯ n 1 ( α ) u 1 + D u n = Γ ( α ) h α F n ,
u n = 1 ζ ( α ) ζ ( 1 + α ) + Γ ( α ) D h α Γ ( α ) h α F n + ζ ( α ) u n 1 k = 1 n u n k k 1 + α σ ¯ n 1 ( α ) u 1 .
The value of the numerical solution u 1 in the first step is computed with the following approximation [39]:
y 1 ( α ) = y 1 y 0 Γ ( 2 α ) h α + a 1 h 2 α .
From the estimate for the error of the L1 approximation [6],
| a 1 | < 1 Γ ( 2 α ) 2 + α 2 α ( 2 α ) 11 + α 12 max 0 t h y ( t ) < 2 M 2 Γ ( 2 α ) ,
where M 2 = max 0 t x y ( t ) . The numerical solution (27) has initial conditions [39]:
u 0 = y ( 0 ) , u 1 = y ( 0 ) + Γ ( 2 α ) F ( h ) h α 1 + D Γ ( 2 α ) h α .
When the solution of two-term Equation (26) satisfies the condition y ( 0 ) = y ( 0 ) = 0 , both Approximations (6) and (7) can be used for the computation of the numerical solution (27).

4.1. Convergence and Error Estimates

Denote by e n = y n u n the error of the numerical solution (27) at the point x n = n h , where n = 0 , 1 , , N . The errors e n satisfy
e n = 1 ζ ( α ) ζ ( 1 + α ) + Γ ( α ) D h α ζ ( α ) e n 1 k = 1 n 1 e n k k 1 + α σ ¯ n 1 ( α ) e 1 + B n h 2 ,
e 0 = 0 , e 1 = Γ ( 2 α ) a 1 h 2 1 + D Γ ( 2 α ) h α .
In Theorems 2 and 3, we prove the convergence of the numerical solution (27) of two-term Equation (26) and derive estimates for the error, depending on the value of the parameter D. The proofs use induction on n, where n = 1 , 2 , , N .
Theorem 2.
Let D , 2 h α ζ ( α ) ζ ( 1 + α ) Γ ( α ) 0 , and y C 2 [ 0 , 1 ] .
Then,
| e n | < C h 2 n α , ( 1 n N ) ,
where C = 2 M 2 + 1 2 max 1 n N | B n | .
Proof. 
The value of the error e 1 in the first step satisfies
| e 1 | = Γ ( 2 α ) | a 1 | h 2 | 1 + D Γ ( 2 α ) h α | < 2 M 2 h 2 | 1 + D Γ ( 2 α ) h α | .
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 < 2 h α ζ ( α ) ζ ( 1 + α ) Γ ( α ) .
D Γ ( α ) h α > 2 ( | ζ ( α ) | + ζ ( 1 + α ) ) ,
| D | Γ ( 2 α ) h α > 2 α ( 1 α ) ( | ζ ( α ) | + ζ ( 1 + α ) ) .
From (15),
ζ ( 1 + α ) > 1 α + γ , | ζ ( α ) | > 1 1 α γ ,
| ζ ( α ) | + ζ ( 1 + α ) > 1 1 α + 1 α = 1 α ( 1 α ) .
Hence, | D | Γ ( 2 α ) h α > 2 . In this case, the denominator of (32) 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 (31) holds for all n = 1 , , m 1 . The values of ζ ( α ) and Γ ( α ) are negative. In both cases, D > 0 and D < 2 h α ζ ( α ) ζ ( 1 + α ) Γ ( α ) , the following estimate holds:
| ζ ( α ) ζ ( 1 + α ) + Γ ( α ) D h α | > | ζ ( α ) | + ζ ( 1 + α ) .
From (30) and Claim 3,
| e m | < 1 | ζ ( α ) | + ζ ( 1 + α ) | ζ ( α ) | e m 1 + k = 1 m 1 e m k k 1 + α + σ ¯ m 1 ( α ) e 1 h 2 + | B m | h 2 < 1 | ζ ( α ) | + ζ ( 1 + α ) | ζ ( α ) | e m 1 + k = 1 m 1 e m k k 1 + α + M 2 4 h 2 + | B m | h 2 < 1 | ζ ( α ) | + ζ ( 1 + α ) | ζ ( α ) | e m 1 + k = 1 m 1 e m k k 1 + α + 2 C h 2 .
By the induction hypothesis,
| e m | < C h 2 | ζ ( α ) | + ζ ( 1 + α ) | ζ ( α ) | ( m 1 ) α + k = 1 m 1 ( m k ) α k 1 + α + 2 .
From Inequality (11),
| e m | < C h 2 | ζ ( α ) | + ζ ( 1 + α ) | ζ ( α ) | m α + ( ζ ( 1 + α ) m α 2 ) + 2 = C h 2 m α .
Theorem 3.
Let 1 Γ ( α ) < D < 0 and y C 2 [ 0 , 1 ] . Then,
| e n | < K h 2 n α , ( 1 n N ) ,
where K = 3 M 2 + max 1 n N | B n | .
Proof. 
The error e 1 in the first step satisfies the bound:
| e 1 | < 2 M 2 h 2 1 α ( 1 α ) D Γ ( α ) h α < 2 M 2 h 2 1 α ( 1 α ) 8 3 M 2 h 2 .
Assume that (33) holds for all n = 1 , , m 1 . Then,
| e m | < 1 | ζ ( α ) | + ζ ( 1 + α ) Γ ( α ) D h α | ζ ( α ) | e m 1 + k = 1 m 1 e m k k 1 + α + σ ¯ m 1 ( α ) e 1 h 2 + | B m | h 2 < 1 | ζ ( α ) | + ζ ( 1 + α ) h α | ζ ( α ) | K h 2 ( m 1 ) α + K h 2 k = 1 m 1 ( m k ) α k 1 + α + K h 2 < K h 2 | ζ ( α ) | + ζ ( 1 + α ) h α | ζ ( α ) | m α + ζ ( 1 + α ) m α 1 < K h 2 m α | ζ ( α ) | + ζ ( 1 + α ) h α | ζ ( α ) | + ζ ( 1 + α ) 1 m α K h 2 m α .
From Theorems 2 and 3, the errors e n of the numerical solution (27) of the two-term equation in the interval [ 0 , 1 ] satisfy the estimate:
| e n | < K h 2 n α K h 2 α
for all n = 1 , , N . Now, we generalize the results of Theorems 2 and 3 to an arbitrary interval [ a , b ] . Consider the two-term equation:
y ( α ) ( t ) + D y ( t ) = F ( t ) , y ( a ) = y 0 , t [ a , b ] .
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 [ 0 , 1 ]
Let E = { e n } n = 1 N be an N-dimensional vector with the errors of the numerical solution (27) of Equation (35). The L (maximum) norm of a vector is the maximum of the absolute values of its elements. From Theorems 2 and 3, we obtain the conditions for the convergence and a bound of the error of the numerical solution (27) of two-term Equations (35) and (36).
Corollary 1.
The maximum error of the numerical solution (27) of two-term Equation (35) satisfies
E < K ( b a ) α h 2 α
when the solution y C 2 [ a , b ] and h = ( b a ) / N :
D , 2 ( ζ ( α ) ζ ( 1 + α ) ) Γ ( α ) h α 1 Γ ( α ) ( b a ) α ,

4.2. Numerical Examples

The numerical results for the error and order of the numerical solution (27) of two-term Equation (26) are given below.
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 .
Two-term Equation (37) has a solution y ( x ) = cos ( α t ) 1 , which satisfies y ( 0 ) = y ( 0 ) = 0 . The experimental results for the error and order of the numerical solution (27) of Equation (37), which uses the approximation (6) 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 .
Equation (38) has a solution, y ( t ) = e α t . The experimental results for the maximum error and order of the numerical solution (27) of Equation (38) are given in Table 3 and Table 4. The experimental results in Table 1, Table 2, Table 3 and Table 4 are consistent with the theoretical results from Theorems 2 and 3.

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 [49,50,51,52,53,54,55].
α C ( S , t ) t α + 1 2 σ 2 S 2 2 C ( S , t ) S 2 + r S C ( S , t ) S r C ( S , t ) = 0 ,
where C is the options 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 (39) has terminal and boundary conditions:
C ( S , T ) = C 0 ( S ) , C ( 0 , t ) = C 1 ( t ) , C ( , t ) = C 2 ( t )
and the fractional derivative is
α C ( S , t ) t α = 1 Γ ( 1 α ) d d t t T C ( S , ξ ) C ( S , T ) ( ξ z ) α d ξ .
The terminal condition is transformed to an initial condition with the substitution [55]:
t = T z , η = T ξ , v ( S , μ ) = C ( S , ξ ) = C ( S , t μ ) .
By applying the substitution to (40), we find [55]
α C ( S , t ) t α = 1 Γ ( 1 α ) d d z z T C ( S , ξ ) C ( S , T ) ( ξ z ) α d ξ = 1 Γ ( 1 α ) d d z 0 z v ( S , μ ) v ( S , 0 ) ( z μ ) α d μ
In view of the definitions of the Riemann–Liouville and Caputo fractional derivatives and (1)
α C ( S , t ) t α = 0 R L D z α v ( S , z ) + v ( S , 0 ) z 1 α Γ ( 2 α ) = D z α v ( S , z ) .
The function v ( S , z ) satisfies the partial fractional differential equation:
D z α v ( S , z ) = 1 2 σ 2 S 2 2 v ( S , z ) S 2 + r S v ( S , z ) S r v ( S , z ) .
The substitution u ( x , t ) = v ( e x , z ) transforms (41) into a linear partial fractional differential equation [55]:
D t α u ( x , t ) = 1 2 σ 2 2 u ( x , t ) x 2 + r σ 2 2 u ( x , t ) x r u ( x , t ) .
Finite difference schemes for the time fractional Black–Scholes equation were constructed in [55,56,57,58]. Now, we construct an implicit finite difference scheme, which uses the approximation (7) 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 ) = 1 2 σ 2 2 u ( x , t ) x 2 + r σ 2 2 u ( x , t ) x r u ( x , t ) + F ( x , t ) ,
where ( x , t ) R = [ B l , B r ] × [ 0 , T ] . Equation (43) 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 n N , 0 m M } .
Denote u n m = u ( n h , m τ ) . The central difference approximations of the partial derivatives of Equation (43) have second-order accuracy:
u ( n h , m τ ) x = u n + 1 m u n 1 m 2 h + e n , m 1 h 2 ,
2 u ( n h , m τ ) x 2 = 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 < M ¯ 3 / 3 and e n , m 2 < M ¯ 4 / 12 . The numbers M ¯ 3 and M ¯ 4 are the bounds for the third- and fourth-order partial derivatives:
M ¯ 3 = max ( x , t ) R 3 u ( x , t ) x 3 , M ¯ 4 = max ( x , t ) R 4 u ( x , t ) x 4 .
The numerical solution { U n m } n = 1 N 1 of Equation (43) on layer m of the grid G satisfies
τ α Γ ( α ) k = 0 m σ k ( α ) U n m k = σ 2 2 h 2 U n + 1 m 2 U n m + U n 1 m + 1 2 h r σ 2 2 U n + 1 m U n 1 m r U n m + F n m
for m = 2 , , 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 = σ 2 2 r σ 2 2 h η , K 2 = σ 0 ( α ) + σ 2 η + r h 2 η , K 3 = σ 2 2 + r σ 2 2 h η .
The numerical solution of the time fractional Black–Scholes Equation (43) is a solution of the system of linear equations:
K 1 U n 1 m K 2 U n m + K 3 U n + 1 m = k = 1 m σ k ( α ) U n m k + h 2 η F n m .
Denote by E n m = U n m u n m the error of the finite difference scheme (44) at the point ( n h , m τ ) and by A n m τ 2 α and B n m h 2 the the truncation errors of the approximations of the left-hand and right-hand sides of (43). The errors E n m on row m 2 of the grid G satisfy
K 1 E n 1 m + K 2 E n m + K 3 E n + 1 m = 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 below 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
Q n m = ζ ( α ) E n m 1 k = 1 m E n m k k 1 + α σ ¯ m 1 ( α ) E n 1 + Γ ( α ) τ α A n m τ 2 α + B n m h 2 .
The numerical solution on the first row of G is computed with Approximation (28):
U n 1 U n 0 Γ ( 2 α ) τ α = σ 2 2 h 2 U n + 1 1 2 U n 1 + U n 1 1 + 1 2 h r σ 2 2 U n + 1 1 U n 1 1 r U n 1 + F n 1 .
Denote η ˜ = Γ ( 2 α ) τ α / h 2 . Then,
K ˜ 1 U n 1 1 + K ˜ 2 U n 1 + K ˜ 3 U n + 1 1 = U n 0 + h 2 η F n 1 ,
where 1 n N and
K ˜ 1 = σ 2 2 r σ 2 2 h η ˜ , K ˜ 2 = 1 + σ 2 η ˜ + r h 2 η ˜ , K ˜ 3 = σ 2 2 + r σ 2 2 h η ˜ .
The errors of the numerical solution of the fractional Black–Scholes Equation (43) on the first row of G are the solutions of the system of linear equations:
K ˜ 1 E n 1 1 + K ˜ 2 E n 1 + K ˜ 3 E n + 1 1 = Γ ( 2 α ) τ α A n 1 τ 2 α + B n 1 h 2 .
Let K ˜ be the ( N 1 ) -dimensional tridiagonal matrix with elements K ˜ 1 , K ˜ 2 , and K ˜ 3 . The error vector for the first row satisfies
K ˜ E 1 = Q 1 ; E 1 = K ˜ 1 Q 1 ,
where
Q n 1 = Γ ( 2 α ) τ α A n 1 τ 2 α + B n 1 h 2 .

5.1. Convergence and Error Estimate

Let C = max n , m { | A n m | , | B n m | } for all n = 1 , , N ; m = 1 , , M and K = | Γ ( α ) C | . The L norm of a matrix is the maximum of the absolute row sums.
Theorem 4.
Let u C 4 [ B l , B r ] × C 2 [ 0 , T ] and r < σ 2 2 h . Then,
E m < K m α τ α h 2 + τ 2 α
for all m = 1 , M .
Proof. 
When 0 < α < 1 , the gamma function satisfies [59]
0.88 < Γ ( 2 α ) < 1 ,
| Γ ( α ) | = Γ ( 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 | = σ 2 2 r σ 2 2 h | η | + σ 2 2 + r σ 2 2 h | η | σ 2 | η | ,
K 2 K 1 | K 3 | > | σ 0 ( α ) | .
Similarly,
K ˜ 2 K ˜ 1 | K ˜ 3 | > 1 .
The matrices K and K ˜ are M-matrices. From the Ahlberg–Nilson–Varah bound [60,61],
K ˜ 1 < 1 , K 1 < 1 | σ 0 ( α ) | = 1 | ζ ( α ) | + ζ ( 1 + α ) .
The L norm of the error vector of the first row of G satisfies the bound:
E 1 K ˜ 1 Q 1 < C τ α τ 2 α + h 2 < K τ α τ 2 α + h 2 .
Suppose that Inequality (46) holds for all m = 1 , , p 1 . From (45),
Q n p < | ζ ( α ) | E p 1 + k = 1 p E p k k 1 + α + σ ¯ p 1 ( α ) E n 1 + Γ ( α ) C τ α h 2 + τ 2 α < K τ α h 2 + τ 2 α | ζ ( α ) | p α + k = 1 p ( p k ) α k 1 + α + 1 8 | Γ ( α ) | + 1 < K τ α h 2 + τ 2 α | ζ ( α ) | p α + ζ ( 1 + α ) p α 1 + 1 8 | Γ ( α ) | .
Therefore,
Q p < ( | ζ ( α ) | + ζ ( 1 + α ) ) K p α τ α h 2 + τ 2 α ,
E p < K 1 Q p < K p α τ α h 2 + τ 2 α .
Corollary 2.
Let r < σ 2 2 h . Then,
E m < K T α h 2 + τ 2 α
for all m = 1 , 2 , , M .
Proof. 
From Theorem 4,
E m < K m α τ α h 2 + τ 2 α K M α τ α h 2 + τ 2 α = K T α h 2 + τ 2 α .

5.2. Numerical Examples

Numerical results for the error and order of the numerical solution of time fractional Black–Scholes Equation (43) are given below.
Example 3.
Consider the following time fractional Black–Scholes equation:
D t α u ( x , z ) = 1 2 σ 2 2 u ( x , t ) S 2 + r σ 2 2 u ( x , t ) 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 .
Equation (47) has a solution u ( x , t ) = e t + x . The experimental results of the numerical solution of Equation (47) 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 τ 1 τ 2 error 1 error 2 in time and log h 1 h 2 error 1 error 2 in space.
The order of the error of the finite difference scheme for the fractional Black-Scholes equation (47) in space is greater than the order in time. When the step size h is small, the errors in the temporal direction dominate the errors in space. The errors of the numerical solution to (47) for different values of N and M are plotted in Figure 1.
Example 4.
Consider the time fractional Black–Scholes Equation (41) for pricing European call options, with a source term F ( x , t ) = 0 and initial and boundary 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 2, 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 α price them lower and evaluate higher out-of-the-money and in-the-money options, compared to the classical Black–Scholes dynamics.
The graphs of the call options and put options premiums for α = 0.7 and all t are given in Figure 3 and Figure 4. To conclude the section with computational simulations, we compare the numerical solutions to (41), using the scheme explained in the paper (7)–(10) and the L1 scheme, as defined, for example, in [55]. We plot the difference as the solution for α = 0.7 , obtained from the L1 scheme, subtracted from the solution, obtained by the A n scheme (7). The results for t = T and for all t are given on Figure 5.

6. Conclusions

In this paper, we studied the convergence and order of the numerical solutions of ordinary and partial fractional differential equations, which use the approximation (7) of the fractional derivative. The results of the numerical experiments given in the paper illustrated the theoretical results. In future work, we will use the method from [16] for the construction of second-order and high-order approximations of the fractional derivative. We will consider applications of approximations for numerical solutions of linear and nonlinear fractional differential equations and analyze their convergence.

Author Contributions

Conceptualization, Y.D. and V.T.; data curation, S.G.; formal analysis, Y.D.; funding acquisition, S.G.; investigation, Y.D.; methodology, Y.D. and S.G.; project administration, Y.D. and V.T.; resources, S.G.; software, S.G.; supervision, V.T.; 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 Project BG05M2OP001-1.001-0004 UNITe: “Virtualization, Digitization and Prototyping”, funded by the Operational Programme “Science and Education for Smart Growth”, co-funded by the European Union trough the European Structural and Investment Funds, 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.

Data Availability Statement

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

Acknowledgments

The authors are grateful to the anonymous referees for the useful suggestions and comments.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Acay, B.; Bas, E.; Abdeljawad, T. Fractional economic models based on market equilibrium in the frame of different type kernels. Chaos Solit. Fract. 2020, 130, 109438. [Google Scholar] [CrossRef]
  2. Cardone, A.; Donatelli, M.; Durastante, F.; Garrappa, R.; Mazza, M.; Popolizio, M. Fractional Differential Equations “Modeling, Discretization, and Numerical Solvers”; Springer: Singapore, 2023. [Google Scholar]
  3. Singh, H.; Kumar, D.; Baleanu, D. Methods of Mathematical Modelling “Fractional Differential Equations”; CRC Press, Taylor & Francis: Boca Raton, FL, USA, 2019. [Google Scholar]
  4. 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, 2340075. [Google Scholar] [CrossRef]
  5. 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]
  6. Cai, M.; Li, C. Numerical approaches to fractional integrals and derivatives: A review. Mathematics 2020, 8, 43. [Google Scholar] [CrossRef]
  7. 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]
  8. Lin, Y.; Xu, C. Finite difference/spectral approximations for the time fractional diffusion equation. J. Comput. Phys. 2007, 225, 1533–1552. [Google Scholar] [CrossRef]
  9. 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]
  10. Scherer, R.; Kalla, S.L.; Tang, Y.; Huang, J. The Grünwald–Letnikov method for fractional differential equations. Comput. Math. Appl. 2011, 62, 902–917. [Google Scholar] [CrossRef]
  11. Nasir, H.M.; Nafa, K. A new second order approximation for fractional derivatives with applications. SQUJS 2018, 23, 43–55. [Google Scholar]
  12. Zeng, F.; Li, C.; Liu, F.; Turner, I. Numerical algorithms for time fractional subdiffusion equation with second-order accuracy. SIAM J. Sci. Comput. 2015, 37, A55–A78. [Google Scholar] [CrossRef]
  13. Dimitrov, Y. A second order approximation for the Caputo fractional derivative. J. Fract. Calc. Appl. 2016, 7, 175–195. [Google Scholar]
  14. Alikhanov, A.A. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys. 2015, 280, 424–438. [Google Scholar] [CrossRef]
  15. 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]
  16. 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]
  17. 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]
  18. 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]
  19. 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] [CrossRef]
  20. 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]
  21. Cao, J.; Cai, Z. Numerical analysis of a high-order scheme for nonlinear fractional differential equations with uniform accuracy. Numer. Math. Theory Methods Appl. 2020, 14, 71–112. [Google Scholar] [CrossRef]
  22. 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]
  23. 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]
  24. Lubich, C. Discretized fractional calculus. SIAM J. Math. Anal. 1986, 17, 704–719. [Google Scholar] [CrossRef]
  25. 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]
  26. 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]
  27. Wu, R.F.; Ding, H.F.; Li, C.P. Determination of coefficients of high-order schemes for Riemann–Liouville derivative. Sci. World J. 2014, 2014, 402373. [Google Scholar] [CrossRef] [PubMed]
  28. 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]
  29. Navot, I. A further extension of the Euler–Maclaurin summation formula. J. Math. Phys. 1962, 41, 155–163. [Google Scholar] [CrossRef]
  30. Dimitrov, Y. Approximations of the fractional integral and numerical solutions of fractional integral equations. Commun. Appl. Math. Comput. 2021, 3, 545–569. [Google Scholar] [CrossRef]
  31. Doha, E.H.; Bhrawy, A.H.; Baleanu, D.; Ezz-Eldien, S.S. On shifted Jacobi spectral approximations for solving fractional differential equations. Appl. Math. Comput. 2013, 219, 8042–8056. [Google Scholar] [CrossRef]
  32. Youssri, Y.H.; Abd-Elhameed, W.M.; Ahmed, H.M. New fractional derivative expression of the shifted third-kind Chebyshev polynomials: Application to a type of nonlinear fractional pantograph differential equations. J. Funct. Spaces 2022, 2022, 3966135. [Google Scholar] [CrossRef]
  33. Lubich, C. Fractional linear multistep methods for Abel-Volterra integral equations of the second kind. Math. Comp. 1985, 45, 463–469. [Google Scholar] [CrossRef]
  34. 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]
  35. 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]
  36. 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]
  37. Ray, S.S.; Gupta, A.K. Wavelet Methods for Solving Partial Differential Equations and Fractional Differential Equations; Chapman and Hall/CRC: Boca Raton, FL, USA, 2018. [Google Scholar]
  38. Ur Rehman, M.; Baleanu, D.; Alzabut, J.; Saeed, I. Green–Haar wavelets method for generalized fractional differential equations. Adv. Differ. Equ. 2020, 515, 2020. [Google Scholar] [CrossRef]
  39. Dimitrov, Y. Approximations for the Caputo derivative (I). J. Fract. Calc. Appl. 2018, 9, 15–44. [Google Scholar]
  40. Matsuoka, Y. On the power series coefficients of the Riemann zeta function. Tokyo J. Math. 1989, 12, 49–58. [Google Scholar] [CrossRef]
  41. Edwards, H.M. Riemann’s Zeta Function; Academic Press: New York, NY, USA, 1974. [Google Scholar]
  42. Weideman, J.A.C. Numerical integration of periodic functions: A few examples. Am. Math. Mon. 2002, 109, 21–36. [Google Scholar] [CrossRef]
  43. Anjara, F.; Solofoniaina, J. Solution of general fractional oscillation relaxation equation by Adomian’s method. Gen. Math. Notes 2014, 20, 1–11. [Google Scholar]
  44. Güulsu, M.; Öztürk, Y.; Anapalı, A. Numerical approach for solving fractional relaxation-oscillation equation. Appl. Math. Model. 2013, 37, 5927–5937. [Google Scholar] [CrossRef]
  45. Hu, Y.; Luo, Y.; Lu, Z. Analytical solution of the linear fractional differential equation by Adomian decomposition method. J. Comput. Appl. Math. 2008, 215, 220–229. [Google Scholar] [CrossRef]
  46. 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]
  47. Podlubny, I. Fractional Differential Equation; Academic Press: Princeton, NJ, USA, 1999. [Google Scholar]
  48. 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]
  49. 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]
  50. Cen, Z.; Huang, J.; Xu, A.; Le, A. Numerical approximation of a time fractional Black–Scholes equation. Comput. Math. Appl. 2018, 75, 2874–2887. [Google Scholar] [CrossRef]
  51. 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. Fract. 2022, 162, 112423. [Google Scholar] [CrossRef]
  52. 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, 1407–1419. [Google Scholar] [CrossRef]
  53. 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]
  54. Kumar, S.; Kumar, D.; Singh, J. Numerical computation of fractional Black–Scholes equation arising in financial market. Egypt. J. Basic Appl. Sci. 2014, 1, 177–183. [Google Scholar] [CrossRef]
  55. 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, 1772–1783. [Google Scholar] [CrossRef]
  56. Grzegorz Krzyżanowski, G.; Magdziarz, M.; Plociniczak, L. A weighted finite difference method for subdiffusive Black–Scholes model. Comput. Math. Appl. 2020, 80, 653–670. [Google Scholar] [CrossRef]
  57. 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]
  58. Song, L.; Wang, W. Solution of the fractional Black–Scholes option pricing model by finite difference method. Abstr. Appl. Anal. 2013, 2013, 194286. [Google Scholar] [CrossRef]
  59. Deming, W.; Colcord, C. The minimum in the gamma function. Nature 1935, 135, 917. [Google Scholar] [CrossRef]
  60. 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]
  61. Nilson, E.N.; Ahlberg, J.H. Convergence properties of the spline fit. J. SIAM 1963, 11, 95–104. [Google Scholar]
Figure 1. Numerical solution to (47) for α = 0.5 .
Figure 1. Numerical solution to (47) for α = 0.5 .
Fractalfract 07 00750 g001
Figure 2. Option premia for various α (left) and the difference between them and the option premium for the integer-order model (right).
Figure 2. Option premia for various α (left) and the difference between them and the option premium for the integer-order model (right).
Fractalfract 07 00750 g002
Figure 3. Call options premiums for all t.
Figure 3. Call options premiums for all t.
Fractalfract 07 00750 g003
Figure 4. Put options premiums for all t.
Figure 4. Put options premiums for all t.
Fractalfract 07 00750 g004
Figure 5. Difference in the numerical solution at t = T (left) and for all t (right).
Figure 5. Difference in the numerical solution at t = T (left) and for all t (right).
Fractalfract 07 00750 g005
Table 1. Error and order of the numerical solution (27) of Equation (37) in [ 0 , 1 ] .
Table 1. Error and order of the numerical solution (27) of Equation (37) in [ 0 , 1 ] .
h α = 0.25 , D = 1 α = 0.75 , D = 100 α = 0.4 , D = 1
Error Order Error Order Error Order
0.005 5.10 × 10 9 1.704 1.02 × 10 8 1.245 6.00 × 10 5 1.574
0.0025 1.37 × 10 9 1.713 3.09 × 10 9 1.2473 3.42 × 10 5 1.581
0.00125 3.7 × 10 10 1.720 9.4 × 10 10 1.2484 2.12 × 10 6 1.586
Table 2. Error and order of the numerical solution (27) of Equation (37) in [ 0 , 2 ] .
Table 2. Error and order of the numerical solution (27) of Equation (37) in [ 0 , 2 ] .
h α = 0.7 , D = 10 α = 0.9 , D = 20 α = 0.5 , D = 1000
Error Order Error Order Error Order
0.005 2.08 × 10 7 3.174 3.67 × 10 8 4.225 1.69 × 10 8 1.493
0.0025 5.18 × 10 6 2.005 7.07 × 10 7 2.378 6.01 × 10 9 1.490
0.00125 1.74 × 10 6 1.577 2.23 × 10 7 1.663 2.15 × 10 9 1.486
Table 3. Error and order of the numerical solution (27) of Equation (38) in [ 0 , 1 ] .
Table 3. Error and order of the numerical solution (27) of Equation (38) in [ 0 , 1 ] .
h α = 0.2 , D = 1 α = 0.4 , D = 100 α = 0.6 , D = 0.1
Error Order Error Order Error Order
0.005 9.62 × 10 8 1.799 6.35 × 10 8 1.599 2.16 × 10 4 1.398
0.0025 2.76 × 10 8 1.799 2.09 × 10 8 1.599 8.21 × 10 5 1.399
0.00125 7.93 × 10 9 1.799 6.91 × 10 9 1.599 3.11 × 10 5 1.399
Table 4. Error and order of the numerical solution (27) of Equation (38) in [ 0 , 2 ] .
Table 4. Error and order of the numerical solution (27) of Equation (38) in [ 0 , 2 ] .
h α = 0.8 , D = 20 α = 0.5 , D = 50 α = 0.5 , D = 1000
Error Order Error Order Error Order
0.005 7.05 × 10 34 15.733 4.65 × 10 46 27.426 4.24 × 10 8 1.498
0.0025 7.98 × 10 32 6.466 4.98 × 10 43 9.864 1.50 × 10 8 1.499
0.00125 7.87 × 10 31 3.342 2.07 × 10 42 4.586 5.31 × 10 9 1.499
Table 5. Error and order of the numerical solution of Equation (47) and N = 100 .
Table 5. Error and order of the numerical solution of Equation (47) and N = 100 .
τ α = 0.25 α = 0.5 α = 0.75
Error Order Error Order Error Order
1 / 10 8.12 × 10 5 2.02 × 10 3 1.94 × 10 2
1 / 20 2.47 × 10 5 1.719 7.47 × 10 4 1.435 8.68 × 10 3 1.162
1 / 40 7.40 × 10 6 1.736 2.69 × 10 4 1.472 3.75 × 10 3 1.212
1 / 80 2.21 × 10 6 1.744 9.61 × 10 5 1.487 1.59 × 10 3 1.232
1 / 160 6.58 × 10 7 1.748 3.42 × 10 5 1.493 6.75 × 10 4 1.241
Table 6. Error and order of the numerical solution of Equation (47) and M = 1000 .
Table 6. Error and order of the numerical solution of Equation (47) and M = 1000 .
h α = 0.2 α = 0.5 α = 0.8
Error Order Error Order Error Order
1 / 4 5.48 × 10 4 5.32 × 10 4 6.59 × 10 4
1 / 8 1.49 × 10 4 1.878 1.58 × 10 4 1.750 2.14 × 10 4 1.620
1 / 16 3.82 × 10 5 1.965 4.46 × 10 5 1.824 6.69 × 10 5 1.681
1 / 32 9.63 × 10 6 1.987 1.21 × 10 5 1.888 1.95 × 10 5 1.780
1 / 64 2.42 × 10 6 1.992 3.10 × 10 6 1.962 5.32 × 10 6 1.873
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.

Share and Cite

MDPI and ACS Style

Dimitrov, Y.; Georgiev, S.; Todorov, V. Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations. Fractal Fract. 2023, 7, 750. https://doi.org/10.3390/fractalfract7100750

AMA Style

Dimitrov Y, Georgiev S, Todorov V. Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations. Fractal and Fractional. 2023; 7(10):750. https://doi.org/10.3390/fractalfract7100750

Chicago/Turabian Style

Dimitrov, Yuri, Slavi Georgiev, and Venelin Todorov. 2023. "Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations" Fractal and Fractional 7, no. 10: 750. https://doi.org/10.3390/fractalfract7100750

APA Style

Dimitrov, Y., Georgiev, S., & Todorov, V. (2023). Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations. Fractal and Fractional, 7(10), 750. https://doi.org/10.3390/fractalfract7100750

Article Metrics

Back to TopTop