A collocation methods based on the quadratic quadrature technique for fractional differential equations

1.   Introduction

In this paper, we discuss the numerical solution of the fractional differential equations (FDEs) initial value problem

where $D^{\alpha}_{*}$ is fractional Caputo derivative of order $\alpha$ and defined as

where $n =\left\lceil \alpha \right\rceil$ is the smallest integer greater than or equal to $\alpha$ and the right hand side function $f(t, y)$ is assumed to be continuous with respect to two variables $t$ and $y$. Here, we note that there are several definitions of fractional derivatives such as Caputo, Grunwald-Letnikov, Marchaud, Riemann-Liouville, Weyl, etc. In this study, we just focus on the Caputo type of fractional derivative.

Fractional calculus is a research topic in many areas of science and engineering, such as signal processing, control engineering, electromagnetism, bioscience, fluid mechanics, electrochemistry, diffusion processes, continuum and statistical mechanics and propagation of spherical flames. Due to this reason, during the past few decades, mathematical theories and numerical analysis of fractional differential equations have received lots of attention and several numerical methods have been developed. For example, Deithem et al. ^[8,11] introduced the theory and numerical schemes for the predictor-corrector type of Adams methods. In ^[15], a non-polynomial collocation method was proposed for fractional equations with having non-smooth solutions. Yan et al. ^[30] developed a higher order predictor-corrector methods using quadratic quadrature techniques based on fractional Adam-type methods ^[11] for both linear and nonlinear cases. In ^[28], authors proposed the usage of a suitable truncated series expressed in terms of fractional powers of the independent variable for ordinary fractional differential equation. Besides, there are relevant researches ^{[8,9,10,13,14,22,23,25,26,27,28,32]}.

The main contribution of this paper is, unlike other numerical schemes described above, to split two different time intervals-one is the sum of the whole previous time interval and the other one is the current time interval. Then, different numerical schemes are provided to the two time intervals to obtain desired solutions at the current time interval-piecewise quadratic quadrature techniques are applied to the whole previous time intervals and Lagrangian collocation methods are cast for the current time step. There are few attempts to use the mixed numerical schemes. For example, in ^[31], a new product integration scheme is introduced by using the idea of local Fourier expansion and several types of quadrature rules. Inspired by this idea, in this paper, a new version of the mixed numerical scheme using the piecewise quadratic quadrature techniques and the Lagrangian collocation methods are introduced.

The procedure of this mixed numerical method can be described as follows: Since the IVP (1.1) is equivalent to the following Volterra integral equation of the second kind, Eq (1.1) is rewritten as the Volterra integral equation :

First of all, a given time interval $[t_0, t_{final}]$ is divided into several sub-intervals. In each sub-interval $[t_i, t_{i+1}]$, we estimate new solutions by using solutions calculated in all previous intervals $[t_0, t_i]$. However, it is not easy to calculate the integral equations directly because FDEs originally have a singularity at the end point of integral interval.

To hurdle this drawback, the Lagrangian interpolating formulation is applied so that the beta function property can be used for the fractional integration, and it enables to eliminate the singularity by removing the fractional integral. For this, we introduce Chebyshev node points in the sub-interval $[t_i, t_{i+1}]$ and calculate solutions at the collocation points. Moreover, solutions at all collocation points through the whole previous intervals are accumulated. Based on all these accumulated solutions, we can easily calculate the fractional integral by using any quadrature rule in the whole previous time interval $[t_0, t_i]$. Here, for the higher accuracy, a piecewise quadratic quadrature rule is introduced. Note that it has been already known ^[6,20] that multi-stage methods such as collocation methods have lots of good properties in the sense of stability and accuracy, compared with multi-step methods.

This paper is organized as follows. In Section 2, we briefly review the basic background such as Lagrangian interpolations. In Section 3, a Lagrangian interpolation technique is applied to calculate the fractional integral in a certain time sub-interval. Moreover, a piecewise quadratic interpolation polynomial is introduced to approximate an integral with all known solutions calculated in all previous time intervals. In Section 4, to examine the effectiveness and efficiency of the proposed scheme, several numerical results are presented and shown numerically a convergence order of the propose method. Finally in Section 5, we summarize our results and discuss possibilities to increase the efficiency of the propose scheme for solving other types of fractional differential equations, such as multi-term FDEs or fractional partial differential equations, etc.

2.   Preliminaries

In this section, we briefly review the numerical techniques required to develop the numerical methods to solve the FDE initial value problems (IVPs).

2.1.   Lagrangian scheme

Suppose that $s_0 < s_1 < \cdots < s_n$ are $n+1$ Chebyshev-Gauss-Lobatto (CGL) node points on $[0, 1]$, where $s_0 = 0$ and $s_n = 1$. We usually discretize the whole time interval into several sub-intervals to solve IVP (1.1) and solve the IVP in each sub-interval. Suppose $y_m$ is approximated for the exact value $y(t_m)$ at time $t_m$. Based on the approximated solutions $y_m$, we need to approximate the solution at time $t_{m+1}$. At first, $n+1$ CGL-points on $[t_m, t_{m+1}]$ are required through the following linear transformation

for $\quad k = 0, 1, \cdots, n,$ where $h = t_{m+1}-t_m$. Once the numerical solutions $y_{m, 1}, y_{m, 2}, ..., y_{m, n}$ at the nodes $t_{m, 1}, t_{m, 2}$, $..., t_{m, n}$ are obtained, we can write the solution $y(t)$ and the function function $f(t, y(t))$ in Lagrange interpolation form as follows:

where $L_k(t)$ us the Lagrange interpolation polynomial of order $n$, given as

Also, $L_k(t)$ should be transformed into following expressions :

where the coefficients $c_{j, k}$ can be computed by it matrix form,

where $I_{n+1}$ is the identity matrix of order $n+1$. Moreover, the fractional integral containing $L_k(t)$ can be easily rewritten by using the Beta function property as follows:

As seen above, it turns out that the Lagrangian polynomial can remove the fractional integral, so the singularity in FDE can be resolved.

3.   Development of the proposed algorithm

In this section, we introduce a mixed numerical method to solve the FDEs for $0 < \alpha < 2$ since the case for $\alpha > 2$ is not our primary practical concern ^[12]. Note that we formulate all equations in terms of the Caputo sense.

3.1.   Chebyshev interpolation

The fractional IVP (1.1) is equivalent to the following Volterra integral equation

with $y(0) = y_0$. The other condition $y'(0) = y'_0$ is needed only for $1 < \alpha < 2$, so, for $0 < \alpha < 1$, $y'(0) = y'_0$ is not necessary, so we set $y'_0 = 0$ for $0 < \alpha < 1$. Here, we suppose the the function $f(\tau)$ satisfies the Lipschitz condition.

First of all, we discretize the given time interval $[0, T]$ into N sub-intervals equally and the step size $h = \frac{T}{N}$. Based on all the solutions calculated in all previous intervals $[0, t_m]$, we approximate the solution at $t_{m+1}$ as follows:

Note that the integral $\int_0^{t_{m}} (t_{m+1}-\tau)^{\alpha-1} f(\tau, y(\tau))d\tau$, the third term of right hand side of Eq (3.2) can be calculated numerically since the approximated values of $y(t_i)$ in $t_i \in [0, t_m]$ are calculated in the previous intervals $[t_i, t_{i+1}]$ for $i = 0, \cdots, m-1$. The details of calculation for the integration will be described in the following subsection. Instead, we discuss the calculation of the last term $\int_{t_m}^{t_{m+1}} (t_{m+1}-\tau)^{\alpha-1} f(\tau, y(\tau))d\tau$ of right hand side of Eq (3.2). Since the first three terms of the right hand side of Eq (3.2) are known values, we just let the three terms a constant $C_m$, so the Eq (3.2) is rewritten as

where

To discretize the integral equation Eq (3.3) in a time interval $[t_m, t_{m+1}]$, we introduce Lagrangian interpolation described in Section 2.1. With $h = t_{m+1}-t_m$, we let $\tau = hs+ t_m, \quad 0 \leq s \leq 1$. Here, we use Gauss-Legendre-Lobatto nodes for $s$ described in Section 2.1.

Eq (3.3) leads to

Using the notations of Lagrangian interpolation

where

Eq (3.4) can be rewritten as

where $\binom ji$ denotes a combination of $j$ objects taken $i$. Using the definition of the beta distribution described in Eq (2.4) and the derivation of Eq (3.5), the following formula can be derived,

where $t$ is a collocation point contained in $[t_m, t_{m+1}]$ defined above. Using the notation of the Lagrangian interpolation for $y(t)$, it produces a nonlinear system for $y_k$,

By solving the nonlinear system Eq (3.7) for $y_k, k = 0, \cdots, n$, we can finally approximate the solution $y(t_{m+1})$.

Theorem 3.1. Let $p_N = \sum_{k = 0}^{n} f_k L_k(x)$ defined above. Then, for $f \in C^N([0, T])$,

for some positive constant $\bar{c}$.

Proof. The details of the proof can be found in ^[1,16,18].

3.2.   Quadrature rule

In this subsection, we explain how to calculate $C_m$ in Eq (3.7) using a piecewise quadratic quadrature rules. There are several well-developed quadrature rules ^{[3,4,5,17,19,24]}. For example, ^[4] presented explicit quadrature rules for spaces of quintic splines with uniform knot sequences over finite domains by using only 2 quadrature points per element. In ^[17,19], an efficient ruless for NURBS-based isogeometric analysis was presented for spaces arising in the calculation, for which the number of quadrature points in an optimal rule is almost equal to half the number of degrees-of-freedom. This idea was extended to the practical computation of quadrature rules for univariate non-uniform splines up to any precision. Despite various choices of the quadrature rules, in this work, we simply use a piecewise quadratic quadrature for calculation of $C_m$ in Eq (3.7), since the quadrature efficiency is not the main focus of this work. Later, we will apply various higher order quadrature rules to increase accuracy of the proposed scheme and report it in the future.

Remind that the $C_m$ is

To calculate the integration in $C_m$, we equally discretize the given whole integration interval $[0, t_{m}]$ into N subintervals. As described in subsection 3.1, in each subtime interval $[t_i, t_{i+1}]$, we use p-Gaussian node points to approximate the solutions. That is, there are p-points in the subinterval $[t_i, t_{i+1}]$:

By Eq (3.7) with a notation of Lagrangian interpolation, solutions at $p$ node points in the interval $[t_i, t_{i+1}]$ can be calculated. Finally, at the time interval $[t_m, t_{m+1}]$, we can have estimated solutions in $mp$ node points in whole previous time interval $[0, t_m]$. To denote the points explicitly, we define $\hat{t}$ by

With having all calculated solutions in the whole previous time interval set $[0, t_m]$ and a notation in Eq (3.9), we consider a sub-interval $[ \hat{t}_{2i}, \hat{t}_{2i+2}]$ having 3 time points ($\hat{t}_{2i}, \hat{t}_{2i+1}, \hat{t}_{2i+2}$), and we can write the integral as the 3-point quadrature polynomial to calculate the fractional integration. That is, the given function $f$ is replaced by the quadratic polynomial:

Therefore the fractional integration in the sub-interval $[ \hat{t}_{2i}, \hat{t}_{2i+2}]$ can be calculated as

where $\omega_1, \omega_2$ and $\omega_3$ are easily obtained by calculation of the integration (3.10). Based on the calculation in sub-interval, we can extend this calculation to the whole time interval $[0, t_m]$. We summarize the calculation as the following remark.

Remark 3.2. For $0 < \alpha < 1$, the piecewise quadratic quadrature for fractional integration is

where

with

Note that to calculate Eq (3.12), we apply the quadrature rule in $\frac{mp}{2}$ sub-intervals. In each sub-interval, 3 $w_i$ should be calculated and each $w_i$ is obtained from one of $w_{k, m+1}^i$ defined in Eqs (3.13)–(3.15). Thus, a computational cost in the interval $[0, t_{m}]$ is $\frac{3mp}{2}$. Overall, since we discretize the given interval $[0, T]$ into $N$ sub-intervals, the whole computational costs for quadrature in the proposed scheme can be $\sum_{m = 1}^N \frac{3mp}{2} = O(pN^2)$.

Theorem 3.3. If $f(t) \in C^3[0, T]$, $y(t_k)$ and $y_k$, $k = 0, 1, \cdots, 2m$ and $T = t_{2m}$ be the solutions, then there exists a constant $C_0$ such that

Proof.

where $P(\tau)$ is the piecewise quadratic polynomial defined in Eq (3.10). Therefore,

where $C$ is a constant depending on $\alpha$. Since $\int_{0}^{t_m} (t_{m+1}-\tau)^{\alpha-1} d\tau = \frac{t_{m+1}^{\alpha}}{\alpha}-\frac{h^{\alpha}}{\alpha}$, Eq (3.17) is summarized as

where $C$ is a constant depending on $\alpha$.

4.   Numerical results

In this section, we test several examples to examine the effectiveness of the proposed scheme. The numerical results are compared with exact solutions. For showing the superiority of the methods, the results are compared with those obtained by other methods ^{[11,21,30,31]}. To investigate numerically the convergence order in each example, experimentally determined order of convergence (EOC) is calculated as follows:

where $Error(h)$ denotes the absolute error between the analytic solutions and the numerical solutions simulated with a step size $h$. In addition, if the right hand side of the FDEs is nonlinear, a nonlinear system is derived from Eq (3.7). For calculation of the nonlinear system, the matlab-builtin routine "fsolve" is used. Also, for the initial guess of the nonlinear solver, we use the fractional explicit Euler method, which is the most basic and economical method. Details of each problem will be explained in each subsection.

4.1.   Linear case

As the first example, we consider a linear fractional differential equation described by

with the initial condition $y(0) = 0$ and $y'(0) = 0$. The exact solutions is $y(t) = t^2$. For the experiment, it is marched from $t = 0.0$ to $t = 1.0$ with a step size $h = 0.1$ and 4 Chebyshev node points are used in each sub interval. To check the effectiveness of the proposed scheme, we plot the absolute errors between the proposed scheme and the analytic solution for different the value $\alpha = 0.25, 0.5$, and $0.75$ in Figure 1.

It can be seen that the proposed scheme seems to work well for this problem.

To examine the numerical convergence order, we calculate numerical errors at $T = 1$ and the experimentally determined order of convergence (EOC) by varying the step size $h = 1/10, 1/20, 1/40, 1/80, 1/160$ and $1/320$ for different $\alpha = 0.5$, $0.75$ and $1.5$ with 2 Chebyshev node points. All results are reported in Table 1 and Figure 2.

Table 1 shows that for $\alpha = 0.5$, EOC converges to $2.5$, for $\alpha = 0.75$ EOC to $2.25$, and for $\alpha = 1.5$ EOC to $1.5$. One can guess that for linear problems, the proposed algorithm has $3-\alpha$ convergence order. For more details, we plot the numerical results and straight line $y = (3-\alpha)x$ for $\alpha = 0.75$ in Figure 2. It can be seen that the two lines in the figure are parallel, so the convergence order of the numerical results for $\alpha = 0.75$ is 2.25. For $\alpha = 1.5$, one can convincingly conclude the convergence order, but for $\alpha = 0.5$ it is not so clear.

Note that for the linear problems of $1 < \alpha < 2$, this scheme seems not appropriate since the convergence order is quite low. Also, to be precise, we need to theoretically analyze the convergence orders.

4.2.   Nonlinear case 1

In this subsection, the following nonlinear fractional differential equation is considered

with an initial condition $y(0) = 0$ and $y'(0) = 0$. The exact solution of this problem is given as

For the experiment, we march from $t = 0.0$ to $t = 1.0$ with step size $h = 0.2$. We calculate the numerical errors between the proposed scheme and the analytic solution is calculated by varying the number of Chebyshev node points $n = 2, 4$, and $6$ for the value $\alpha = 0.25, 0.5$, and $0.75$ and all results are plotted in Figure 3. The results show that the proposed scheme has higher accuracy as the number of Chebyshev node points is increasing, which implies why the collocation methods are useful.

For investigating the numerical convergence order for the nonlinear problem, numerical errors and the experimentally determined order of convergence (EOC) are computed by varying the step size $h$ for different $\alpha = 0.25, 0.5$ and $0.75$ with 2 Chebyshev node points in time interval $[0, 1]$ and results are reported in Table 2.

In Table 2, we observe that for $\alpha = 0.25, 0.5$ and $0.75$, the EOCs converge to $3.2, 3.5$ and $3.75$ as the step size is decreasing. Similar to the previous result, for $\alpha = 0.5$ and $0.75$, one can convincingly conclude the convergence order, but for $\alpha = 0.25$ it is not so clear. It is almost $3+\alpha$, unlike the case of the linear problem. It is almost identical to the theoretic convergence order described in Theorem 3.3.

Additionally, to investigate the efficiency of the proposed scheme, we compare the results from the proposed scheme with those from the existing higher order methods introduced in ^[7]. For the experiment, we march from $t = 0$ to $t = 1$ with various time steps $h = 1/10, 1/20, 1/40, 1/80, 1/160$ and $1/320$ to check the EOCs for both schemes at the same time. For the proposed scheme, 2 Chebyshev nodes are used. We calculate the absolute errors and the corresponding EOCs for both methods by difference between the results and analytic solutions. All results are reported in Table 3.

As mentioned, the method introduced in ^[7] is also higher order numerical scheme, so as seen in Table 3, the method in ^[7] have similar convergence order to the proposed scheme. However, it can be easily seen that the proposed scheme has more accurate results in this comparison. Therefore, we can conclude that the proposed scheme is quite efficient for this problem.

4.3.   Nonlinear case 2

As the last example, we consider the following nonlinear fractional differential equation described by

with an initial condition $y(0) = 0$ and $y'(0) = 0$. The exact solution of this problem is

To investigate the accuracy of the proposed scheme, we compare the numerical results of the proposed scheme with those obtained from existing methods-one is a numerical method introduced in ^[21] and the other one is predictor-correction methods ^[11]. For the experiment, it is marched from $t = 0.0$ to $t = 1.0$ with step size $h = 0.025$ for $\alpha = 0.75$. The numerical results introduced in ^[21] have been already represented in ^[21], so the data are directly excerpted from the reference. The predictor-corrector methods ^[11] is implemented for the same setting and the proposed scheme uses 4 Chebyshev node points. Figure 4 represents numerical errors generated from three methods over the given time interval.

The figure shows that the proposed scheme generates more accurate solution for the problem, compared with other existing methods. As well as the numerical methods in references ^[11,21], there are several references to represent the numerical results of this example and the results can be found in ^[2,29,30,31], etc. Even compared with the results in ^[31] where is most recently developed and has higher accuracy, the results of the proposed scheme is quite competitive.

To check the numerical convergence order, we investigate numerical errors at $T = 1$ and the experimentally determined order of convergence (EOC) by varying the step size $h$ for different $\alpha = 0.25, 0.5$ and $0.75$ with 2 Chebyshev node points in Table 4.

The results show that for $\alpha = 0.25, 0.5$ and $0.75$, the EOCs are $3.25, 3.5$ and $3.75$, respectively. That is, as similar to the previous nonlinear example, a numerical convergence order for this problem is about $3+\alpha$. For checking in detail, we plot the EOC for $\alpha = 0.5$ and a straight line with having a slope $3.5$ in Figure 5. The figure shows the two lines are exactly parallel.

To examine the numerical convergence order for $1 < \alpha < 2$, Numerical errors at $T = 1$ and the experimentally determined order of convergence (EOC) are computed by varying the step size $h$ for different $\alpha = 1.5$ and $1.75$ with 2 Chebyshev node points and reported in Table 5.

The results show that for $\alpha = 1.5$ and $1.75$, the EOC is over $4$. It is quite competitive convergence order compared with the existing techniques.

5.   Conclusions

In this paper, a mixed numerical technique is developed for solving fractional differential equations (FDEs) by splitting into two time intervals-the whole previous time interval and current time integration step. At a current time integration, we compute the solution by using Chebyshev collocation methods which can be rewritten as a Lagrangian interpolating form. By the Lagrangian interpolating form, we can remove a singularity which belongs originally to the FDEs. For calculating the fraction integral in the previous time interval, we use a piecewise quadratic quadrature technique to get higher accuracy by using all solutions including solutions at the collocation points. Several numerical examples are presented to show the efficiency of the proposed method and compare it with several existing methods. The numerical results present that the proposed techniques can get competitively better accuracy and and numerically convergence orders $3-\alpha$ for linear cases. Also, the numerical convergence orders are $3+\alpha$ and over $4$ for nonlinear cases, when $0 < \alpha < 1$ and $1 < \alpha < 2$, respectively.

In order to fully explore the efficiency of the proposed scheme, several extended issues are currently being pursued. First, we apply the proposed scheme to other types of fractional differential equations such as several Bagley-Torvik equations, a popular FDE with $\alpha > 1$. Secondly, we are doing theoretically convergence analysis for linear which will be hopefully consistent with the numerical convergence orders. Lastly, an adaptive time stepping method should be considered by calculating the fractional integral for long time simulations. Preliminary results are quite promising. Results along these issues will be reported soon.

Acknowledgments

This work was supported by basic science research program through the National Research Foundation of Korea (NRF) funded by the Korea government (MSIT) (grant number: NRF-2019R1F1A1058378).

Conflict of interest

The author declares no conflict of interest in this paper.