1. Introduction
In recent decades, fractional calculus, due to its ability for modeling memory and hereditary properties of various materials and processes, has been applied to many fields of science and engineering, including relaxation and oscillation, viscoelasticity, anomalous diffusion, control problems, etc. [
1,
2,
3,
4,
5,
6,
7,
8]. At present, theories and applications of fractional calculus have attracted much interest and have become a vibrant research area. Fractional differential equations, including the existence, uniqueness and stability of solutions, were studied by some scholars [
2,
3,
4,
6,
9,
10]. In particular, new analytical and numerical methods were proposed [
3,
4,
6,
7,
11,
12,
13]. Lie symmetry analysis and conservation laws were investigated for fractional evolution equations [
14].
In [
15], the difference method and its convergence for the space-time fractional advection–diffusion equation were investigated. In [
16], a matrix representation of discrete analogues of fractional differentiation and integration was suggested and used to the numerical solution of fractional integral and differential equations. In [
17], a local discontinuous Galerkin finite element method was suggested for Caputo-type fractional partial differential equations. In [
18], a numerical Laplace transform technique was used to solve the irrational fractional-order systems. In [
19,
20,
21], numerical methods based on spline functions were presented for the fractional differential equations. In [
22], an Adams-type predictor–corrector method for the numerical solution of fractional differential equations was proposed. In [
23], a numerical differentiation formula for the Caputo fractional derivative was developed by means of the quadratic interpolation approximation on three nodes. In [
24], a survey and a MATLAB software tutorial for numerical methods were presented. In [
25], Wang et al. proposed an asymptotic approximation method for a class of linear weakly singular Volterra integral equations based on the Laplace transform.
We recall the basic concepts in the fractional calculus. Let
be piecewise continuous on
and integrable on any finite subinterval of
. The Riemann–Liouville fractional integral of
has the definition
where
is a positive real constant, and
is Euler’s gamma function.
If
exists and
is integrable on any finite subinterval of
, then applying the integration by parts in Equation (
1), we have
It follows that
So, we rationally define
for complementarity. The fractional integral satisfies the following equalities,
The Riemann–Liouville fractional derivative of
of order
is defined, when it exists, as
while the Caputo fractional derivative of
of order
is
From the definitions, for the Caputo fractional derivative of polynomial functions, the following equality holds
and for the power function
,
, the Caputo fractional derivative is
where
denotes the least integer greater than or equal to
. We use the Caputo fractional derivative in this article in view of its convenience for formation of initial value condition. We denote the operators
as
and
as
for short.
In this article, we consider approximate solution by quadratic spline interpolation function for the initial value problem of the fractional differential equation with two Caputo fractional derivatives
where
are constants, and
is a given continuous function on the interval
. In next section, we derive the quadratic spline interpolation function. In
Section 3, the numerical schemes of the two fractional derivatives based on the quadratic spline interpolation are devised. In
Section 4, the recursion scheme of numerical solutions for the fractional differential equation is generated. Two numerical examples are used to check the proposed method. Additionally, comparisons with the
–
numerical solutions are conducted.
2. Quadratic Spline Interpolation Function
Suppose
, for
, are known, where
and
. Additionally, suppose
is known. Then, the quadratic spline interpolation function
with the nodes
,
, satisfying
for
and
, exists uniquely [
26,
27]. This means on each subinterval,
,
is a quadratic polynomial,
and
is continuous on the whole interval
.
First, we introduce the notations
for
which serve as the interim parameters, to derive the quadratic spline interpolation function. Due to
being a spline function of degree 1, interpolating the values
,
so
has the form on the subinterval
,
Operating indefinite integration leads to
where
are the integral constants. By setting
we obtain
, and so Equation (
5) becomes
The parameters
,
may be determined by the continuity of the function
on the interval
as
Applying the condition to Equation (
6), we have
By the iterations in Equation (
7), we give expressions to each
,
in terms of
,
as
The general form is
where the sum
vanishes if
.
Substituting Equation (
8) into Equation (
6), we obtain the quadratic spline interpolation function
on the subinterval
expressed through
,
,
We indicate that due to
,
, are constants, Equations (
4) and (
6) will be used in the numerical computing process of fractional derivatives, and the expression of
in Equation (
8) will be substituted at the final procedure to avoid large expressions by using Equation (
9).
For estimation of interpolation remainder
, it was proved that if
with
of bounded variation, then there exists
such that [
26,
27]
Li and Huang [
28] proved the result under the assumption
.
3. Numerical Computation of Fractional Derivatives
We calculate numerically the fractional derivatives
and
at each nodes
,
, using the quadratic spline interpolation function
. First, the
-th order derivative is approximated as
Form Equation (
4),
is piecewise constants on the interval
, and has the form on the interval
as
Integrating piecewise in Equation (
11) and applying Equation (
12) yield
Substituting
leads to
Regrouping the right hand side leads to the following equation
where
Substituting the derivatives
,
, in Equation (
8) into Equation (
14), we obtain the fractional derivative of order
at
in terms of
as
where
For the
-th order fractional derivative
at
, in a similar manner, we have
Instead, here
is a linear function as in Equation (
4). The sub-domain integration is calculated as
Substituting it into Equation (
18) and regrouping according to
lead to
where
Substituting the derivatives
,
, in Equation (
8) into Equation (
19), we obtain the fractional derivative of order
at
in terms of
,
as
where
We remark for the two fractional derivatives that the integral in Equation (
18) is a little more tactical than that in Equation (
13), and
in Equation (
20) has a different form compared with
in Equation (
15). Nevertheless, except the expressions of
and
, the expressions in Equations (
19)–(
22) present the same layouts as in Equations (
14)–(
17).
For the error estimation, from Equation (
10) we have
and
for
4. Solution of Fractional Differential Equation
At
, the fractional differential Equation (
2) becomes
where
are known values. Approximating the two fractional derivatives and the first order derivative by the counterparts of the quadratic spline interpolation function
yields
where the truncation error is estimated from Equations (
10), (
23) and (
24) as
where
M is a constant related to
.
Inserting the results about the derivative
in Equation (
8) and the fractional derivatives in Equations (
16) and (
21) into Equation (
25), we have
Leaving out the error term
and replacing
by their numerical approximations
, we obtain the recursion scheme of the numerical approximations
,
, from
as
Thus the recursion scheme (
27) gives the numerical solutions
, derived from quadratic splines, and Equation (
9) gives the quadratic spline approximate solution by replacing
by
,
.
We will compare the present algorithm with the usual
–
algorithm to approximate the fractional derivatives
and
[
1,
23]. We note that the
method utilizes the quadratic interpolation polynomials on three nodes to approximate the function
, while the
method approximate the function
by using piecewise linear interpolation. So for the present problem, Equations (
2) and (3), the first two node values need to be given as the iterative initial values in the
–
numerical solutions.
On the interval
,
, approximate
by the quadratic interpolation polynomials
on the nodes
, while on the first interval
,
is approximated by
. Thus, the
method derives the following approximation for the fractional derivative:
The first-order derivative is also approximated by using the quadratic interpolation polynomials as
for
. For the fractional derivative of order
, using the
method we have
Thus, by discretizing Equation (
2) at
,
, and using Equations (
28) and (
29), the
–
numerical solutions are obtained as
Here, we use , and , , to denote the quadratic spline numerical solutions derived from the quadratic spline interpolation and the – numerical solutions derived from the and methods, respectively.
Next, we consider two numerical examples, one has a monotonically increasing excitation and another has a sinusoidal excitation.
Example 1. Consider the initial value problem for the fractional differential equation The exact solution can be expressed in terms of the generalized Mittag–Leffler functions [
3,
4].
where the generalized Mittag–Leffler function is defined as
Calculating the convolution in Equation (
30) yields the exact solution in the following form
We take
,
,
to compute the quadratic spline numerical solutions
on the interval
from Equation (
27) and compare them with the exact solution in Equation (
31) and the
–
numerical solutions
. In
Figure 1a, the black dash line is depicted from the exact solution in Equation (
31), while the blue circles denote the numerical solutions
and the red crosses denote the numerical solutions
by using
. In
Figure 1b–d, the errors
of the numerical solutions
are plotted for the different step-sizes
h = 0.1, 0.05 and 0.025, respectively.
We examined the numerical solutions
derived from the
–
methods with the same step-sizes and found that the errors also decrease oscillatorily as
t increases, but the maximal error is about five times of that of the numerical solutions
. For example, in
Figure 2, we depict the plot of the errors
of the numerical solutions
derived from the
–
methods with the step-size
on the interval
. Compared with the errors of the numerical solutions
with the identical step-size in
Figure 1c, the maximal errors of the numerical solutions
increases to about five times.
Example 2. Consider the initial value problem for the fractional differential equation For this example, we use the high-precision numerical inverse of the Laplace transform proposed by Wang et al. [
29] as a reference to the exact solution. The Laplace transform of the solution
is
The numerical solutions obtained by the high-precision numerical Laplace inverse transform are denoted by .
In
Figure 3a, the three numerical solutions on the interval
with the step-size
obtained from the high-precision numerical Laplace inverse transform, the present quadratic spline interpolation and the
–
methods are displayed. In the initial stage, we can readily see that the quadratic spline numerical solutions
are closer to
than the
–
numerical solutions.
The differences of the two numerical solutions
and
are shown in
Figure 3b–d for
, 0.05 and 0.025, respectively. We also examined the differences of the numerical solutions
and
and found that for a same step-size, the maximum value of
is about three times of that of
. In
Figure 4, the differences
of the numerical solutions from the high-precision inverse Laplace transform and the
–
methods with the step-size
are shown. Compared with
Figure 3c, the maximum value in
Figure 4 enlarges to about three times.
The considered equation belongs to the forced fractional oscillation equation. In Example 1, the transient oscillation evolves into steady-state increasing, while in Example 2, the transient oscillation grows into steady-state oscillation. Two examples show that the decreasing of the step-size h can effectively enhance the accuracy of the numerical solutions, and the present quadratic spline numerical solutions have higher precision than the – numerical solutions. It is worth noting that for the considered problem, in the initial stage the numerical solutions emerge slightly large errors, but with the process evolution errors can fall off in an oscillatory manner. In Example 1, the second-order derivative of the exact solution does not exist at . This responds the unusual errors in the initial stage.
We note that if a differentiable nonlinearity
is added in the left hand side of Equation (
2), then we can approximate the nonlinearity at
,
, as
to derive an explicit scheme of numerical solutions in the nonlinear case. Inasmuch as we focus on the numerical schemes for the fractional derivatives by using the quadratic splines, examples in the nonlinear case are not involved.