1. Introduction
Fractional calculus [
1,
2,
3] has attracted increased interest over the last decade and has been applied in several fields including finance, control theory, electronic circuit theory, mechanics, physics, and signal processing [
4,
5,
6,
7,
8,
9,
10,
11]. There are two popular definitions of the fractional differentiation: the Riemann–Liouville derivative and the Caputo derivative. Let
,
n be a positive integer with
, and
.
Riemann–Liouville derivative: The Riemann–Liouville derivative of a function
starting at the point
a is
Caputo derivative: The Caputo derivative of a function
starting at the point
a is
The comparison of these two definitions can be found in [
12] and the definitions of fractional derivatives are also revised in some studies [
11,
12,
13,
14].
The trapezoidal rule was used for integration or differential equations in the following papers [
15,
16,
17]. However, the functions of the integrand are assumed to be regular. This paper is devoted to the computation of the Caputo fractional derivative on financial derivatives [
18,
19,
20,
21]. In some of them, the functions of the stock or option prices are only of Lipschitz continuity. Our goal is to calculate the Caputo fractional integral for non-smooth functions. This calculation will also encounter the difficulty induced by the singular kernel. In [
18], an implicit numerical discretization is used for the Riemann–Liouville integral to calculate the chaotic behavior for financial models. In [
22], the treatment for a singular kernel involves the linear expansion of the smooth functions and direct integration of the product of the linear polynomial and the singular kernel. In our approach, we consider the function non-smooth. The function could be also singular, and the impact of the function for the integral is similar to the kernel.
Let
n be a positive integer and
be an interval. Define
and
, where
. To explore the niche of this research, let us explain the following examples. The set of
represents the collection of all functions whose domain on
and they are of a continuous
k-th derivative. If
, it is well known in the textbook of numerical analysis, and the approximation is
where
is in
. For the particular case,
, the order of accuracy of the trapezoidal rule method is reduced because the function
belongs exclusively to
.
Definition 1. Let I be an interval and the setwhere is a continuous function for each . For example,
,
, and
. Then,
with
If , then and is continuous on for each . Hence, . Moreover, for a fixed x, the function may not be continuous on c since and for all .
This paper is organized as follows. The order of accuracy for the trapezoidal method on the set
is derived in
Section 2. The proposed method for calculation of Caputo fractional derivative is described in
Section 3, using three examples. Smooth, regular and non-regular functions are used in numerical simulations in
Section 4.
Section 5 shows the analysis of the method to explain the obtained results and
Section 6 demonstrates two applications of the proposed method. The conclusion is given in the last section.
2. Order of Accuracy for Trapezoidal Method on
In this section, we extend the analysis of the order of accuracy for the trapezoidal method on the set . Let us begin to consider the interpolation on the set .
Lemma 1. Let . The linear interpolation of f on has the propertywhere is a continuous function and . Proof. Here, . □
Lemma 2. Let and . Then,where is a continuous function. Proof. This lemma holds. It is followed by Lemma 1 and
where
in
, and the second equality is followed by the weighted mean value theorem. □
Lemma 3. Let and . Then, Proof. From the following,
and taking the subtraction of the above two equations, it yields
□
Moreover, for and . Since h is continuous and bounded by the extremum theorem of continuous functions on a closed interval, Lemma 2 can be re-estimated to be Theorem 1 below.
Theorem 1. Let and . Then, Remark. If then , where between c and x, then .
Theorem 2. Let and . Ifand is uniformly bounded for all , then Proof. Using (
3) as
taking the integration of the above equation on
, we have
Since
is uniformly bounded for all
and
is continuous on
, it implies that
is uniformly bounded for
and
The last equality is followed by
and (
4). Therefore, this theorem holds. □
3. Method
For the sake of simplicity and without loss generality, the case of
is considered in the whole paper. Equation (
1) is equal to
or
here,
.
Let the interval
and
N be a positive integer. The interval
I is divided into
N-subintervals
with the sample points
,
.
Since
is monotonic whenever
, the inverse of
exists. Using the substitution rule,
for fixed
t, the integral
can be rewritten into
where
and
. The linear interpolation of
on the interval
I with the endpoints
and
is
Substituting (
8) into (
7), it yields
The approximation in the last equation listed above represents the trapezoidal method but uses the Riemann–Stieltjes integral. The roles of
f and
g may be interchanged. Equation (
9) is modified to
where
is the Heaviside step function. We refer to the approach in (
10) as the TRSI method. If the function
f is smooth and
is non-smooth, then TRSI in (10) may only use
. On the other hand, the function
is smooth and
f is non-smooth, then TRSI in (10) may only use
. For Caputo fractional derivatives,
is described as the form
and its derivative is singular at its origin. Therefore, if the function
f is smooth, then
only occurs at the singularity of
.
The stability of the TRSI method to use Equation (
6) is to estimate the following:
If
is uniformly bounded for
,
and
are bounded,
, then
and it follows that TRSI is stable. The condition
is uniformly bounded. It also indicates the existence of the Riemann–Stieltjes integral. It is identical to the existence of the Riemann–Stieltjes integral
requires the condition that the discontinuity of
and
cannot occur coincidentally, and vice versa. Therefore, the stability theorem of the TRSI method is stated in the following theorem.
Theorem 3. The TRSI method is stable if the condition that the discontinuity of and φ cannot occur coincidentally is held.
4. Simulations
Let us consider the interval
and there are
N uniform cells; that is, each subinterval
has the length
with the sample points
. We will vary
from
to
. To probe the behavior of the TRSI method, let us define the 1-norm,2-norm and ∞-norm in vectors of numerical solutions by
Furthermore, the order of accuracy is defined as
where
and
is the error between the numerical and exact solutions at the size of zones
. In the following subsection, we adopt three examples as model examples which represent the smooth, regular and non-smooth functions from Example 1 to Example 3 below, respectively.
4.1. Model Examples
Example 1. Let us consider
and
. The polynomial is smooth because
exists for any
n, which is a non-negative integer. The Caputo fractional derivative of
for
is
The analytic solution is
. The errors between the exact and numerical solutions are shown in
Table 1, which demonstrates that the order of accuracy is near 1.5 for 1-norm, 2-norm and
∞-norm.
Example 2. Let us consider
and
. The power function
only can take the first derivate because
is singular at the origin. The Caputo fractional derivative of
for
is
The analytic solution is
. The errors are shown in
Table 2. The results demonstrate that the order of accuracy is near 1.5, 1.45 and 1 for 1-norm, 2-norm and
∞-norm, respectively.
Example 3. Let us consider
and
. The power function
f does not have the first derivative because
does not exist. The Caputo fractional derivative of
for
is
The analytic solution is
. The errors are shown in
Table 3. The order of accuracy is near 0.52, 0.54 and 0.16 for 1-norm, 2-norm and
∞-norm, respectively. In
Figure 1, the top-left panel shows the exact solution (red dot line) and the numerical solution (blue solid line). The errors between the numerical and exact solutions are shown in the top-right panel. The zoom-in profiles on
are shown in the corresponding panels below.
The approximation of the non-smooth or continuous function may improve the accuracy by refining the meshes. However, it is not equivalent to a finer mesh refinement in this case, as the kernel function
is not only non-smooth, but it is singular for fixed
. Therefore, we divide the subinterval by
-zones again. More precisely,
where
,
, with
. The results of fixed
for
,
are shown in
Table 4 and the corresponding profiles are shown in
Figure 2. The errors were reduced from
to
; see
Table 3 and
Table 4, respectively.
4.2. A Comparison Study
The modified trapezoidal rule (MTR) [
22] uses the linear interpolation on
rather than
in the traditional sense for the following integral, and we rewrite it as shown below. The integral can be approximated by
where
The errors are shown in
Table 5 and
Table 6 for model example 1 and 2, respectively. However, Example 3 cannot be simulated by the MTR method because the derivative of the exact function does not exist at the origin.
5. Error Analysis
Let us start to observe the approximation of the function
by the linear interpolation
on
,
The error
on
has the maximum error
where
. Let
,
; the error for
is
This explains that the reason for Example 2 using the trapezoidal method is only of first-order accuracy.
Theorem 4. Let the function be the linear interpolation of the function on each subinterval , and is uniformly bounded for . The modified trapezoidal rule for calculationhas the error bounded by and Proof. It follows that the error is less than
□
Theorem 4 can be applied to explain the results (
Table 5 and
Table 6) for Example 1 and Example 2 obtained using MTR. Next, we will analyze the TRSI method. Let us first recall the error analysis for smooth functions as Theorem 5 below for the trapezoidal method in comparison with the estimation of the errors for the functions in
shown in Theorem 6 below.
Let
and
. Then,
if
has the third continuous derivative. Furthermore, if
and
are continuous whenever
, then
The above approximation leads to the theorem below.
Theorem 5. If exists and is continuous and , thenfor . Furthermore,and for , the following approximation is reduced to From (
7) and Theorem 2, if
, then we have
for
and for
, it is reduced to
Furthermore, the term
. On the other hand,
Theorem 6. If and , thenfor . Furthermore,and for , the following approximation is reduced to Let
k be a positive integer and
. If an integration scheme has the order of accuracy
,
where
and
for some numerical method, then the refinable approach using the mesh
is read as
This implies that the order of accuracy is .
7. Conclusions
The analysis of the trapezoidal method was extended from to and, for each , has the order of accuracy . The trapezoidal method using the Riemann–Stieltjes integral on Caputo fractional derivatives for non-smooth functions was proposed, and the approximation ability was also investigated using three models of examples of smoothness, regularity and non-smoothness. The product of the integrand reveals that, if and the integration is approximated by using the differential , then the trapezoidal method has the second order of accuracy compared to the traditional one. On the other hand, if the integration is approximated by using the differential , , then the order of accuracy for the trapezoidal method is of the fractional order of accuracy. The novelty of this method can be addressed to automatically choose the non-smooth functions or the singular kernel for linear interpolation.
The errors in
Table 3 show that increasing the number of zones cannot significantly improve the accuracy, and the order of accuracy is 0.16 for
∞-norm. Therefore, a refining mesh shown in
Table 4 demonstrated that the order of accuracy is 1.59 for the
∞-norm. To confirm this point, we further apply the refinable approach to MTR. The result for the MTR method using a refinable approach is shown in
Table 9; the order of accuracy improves from 1.0 to 1.50 for the
∞-norm, see
Table 6 and
Table 9.