A Fractional Gauss-Jacobi quadrature rule for approximating fractional integrals and derivatives

TL;DR

This paper presents a new fractional Gauss-Jacobi quadrature rule that improves the accuracy of approximating fractional integrals and derivatives, facilitating solutions to fractional boundary value problems and Euler-Lagrange equations.

Contribution

It introduces an efficient algorithm based on Hale and Townsend's method for roots and weights, enhancing accuracy over previous Golub-Welsch based approaches.

Findings

The proposed method outperforms previous algorithms in accuracy.

It effectively solves fractional boundary value problems.

The approach is validated through illustrative examples.

Abstract

We introduce an efficient algorithm for computing fractional integrals and derivatives and apply it for solving problems of the calculus of variations of fractional order. The proposed approximations are particularly useful for solving fractional boundary value problems. As an application, we solve a special class of fractional Euler-Lagrange equations. The method is based on Hale and Townsend algorithm for finding the roots and weights of the fractional Gauss-Jacobi quadrature rule and the predictor-corrector method introduced by Diethelm for solving fractional differential equations. Illustrative examples show that the given method is more accurate than the one introduced in [Comput. Math. Appl. 66 (2013), no. 5, 597--607], which uses the Golub-Welsch algorithm for evaluating fractional directional integrals.

Tables

Table 1. Table 1. The errors obtained for Example 4.1 from ( 13 ), with α = 0.25 , 0.5 , 0.75 𝛼 0.25 0.5 0.75 \alpha=0.25,0.5,0.75 and n = 5 , 6 , 7 , 8 , 16 𝑛 5 6 7 8 16 n=5,6,7,8,16 in ( 10 ).

$\alpha$	$n=5$	$n=6$	$n=7$	$n=8$	$n=16$
$0.25$	$3.22\times {10}^{-6}$	$5.14\times {10}^{-8}$	$6.1\times {10}^{-10}$	$5.58\times {10}^{-12}$	$2.81\times {10}^{-15}$
$0.5$	$4.85\times {10}^{-6}$	$7.75\times {10}^{-8}$	$9.18\times {10}^{-10}$	$8.37\times {10}^{-12}$	$7.12\times {10}^{-16}$
$0.75$	$5.35\times {10}^{-6}$	$8.35\times {10}^{-8}$	$9.65\times {10}^{-10}$	$8.6\times {10}^{-12}$	$1.39\times {10}^{-15}$

Table 2. Table 2. The errors obtained for Example  4.2 from ( 13 ) with n = 2 , 3 , … , 7 𝑛 2 3 … 7 n=2,3,\ldots,7 in ( 12 ) for approximating the Liouville–Caputo fractional derivative of order α = 0.5 𝛼 0.5 \alpha=0.5 of function f ​ ( t ) = t 4 𝑓 𝑡 superscript 𝑡 4 f(t)=t^{4} .

$n=2$	$n=3$	$n=4$	$n=5$	$n=6$	$n=7$
$4.0\times {10}^{-17}$	$8.8\times {10}^{-16}$	$1.9\times {10}^{-16}$	$6.4\times {10}^{-16}$	$3.8\times {10}^{-16}$	$3.8\times {10}^{-16}$

Equations

$${}_a{I}_x^{\alpha }f(x)=\frac{1}{\Gamma (\alpha )}{\int }_a^x(x-t{)}^{\alpha -1}f(t)dt$$

$${}_x{I}_b^{\alpha }f(x)=\frac{1}{\Gamma (\alpha )}{\int }_x^b(t-x{)}^{\alpha -1}f(t)dt,$$

$$\Gamma (x)={\int }_0^{\infty }{t}^{x-1}{e}^{-t}dt,$$

$${}_a{I}_x^{\alpha }{}_a{I}_x^{\beta }={{}_a{I}_x}^{\alpha +\beta },{}_a{I}_x^{\alpha }{}_a{I}_x^{\beta }={}_a{I}_x^{\beta }{}_a{I}_x^{\alpha },$$

$${}_a{I}_x^{\alpha }{x}^{\mu }=\frac{\Gamma (\mu +1)}{\Gamma (\alpha +\mu +1)}{x}^{\alpha +\mu },$$

$${}_a{D}_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_a^x(x-t{)}^{n-\alpha -1}f(t)dt$$

$${}_x{D}_b^{\alpha }f(x)=\frac{(-1{)}^n}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_x^b(t-x{)}^{n-\alpha -1}f(t)dt,$$

$${}_a^C{D}_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}{\int }_a^x(x-t{)}^{n-\alpha -1}{f}^{(n)}(t)dt$$

$${}_x^C{D}_b^{\alpha }f(x)=\frac{(-1{)}^n}{\Gamma (n-\alpha )}{\int }_x^b(t-x{)}^{n-\alpha -1}{f}^{(n)}(t)dt,$$

$$({}_a^C{D}_x^{\alpha }{}_a{I}_x^{\alpha }f)(x)=f(x)$$

$$({}_a{I}_x^{\alpha }{}_a^C{D}_x^{\alpha }f)(x)=f(x)-\sum _{k=0}^{n-1}{f}^{(k)}(0^{+})\frac{(x-a{)}^k}{k!},t>0.$$

$$D_{GL}^{\alpha }=\underset{n\to \infty }{\lim }\frac{(\frac{t}{N}{)}^{-\alpha }}{\Gamma (-\alpha )}\sum _{j=0}^{n-1}\frac{\Gamma (j-\alpha )}{\Gamma (j+1)}f\left(t-\frac{tj}{n}\right).$$

$$(1-x{)}^{\lambda }(1+x{)}^{\nu },\lambda ,\nu >-1,$$

$${\int }_{-1}^1(1-x{)}^{\lambda }(1+x{)}^{\nu }{P}_n^{(\lambda ,\nu )}(x)P_m^{(\lambda ,\nu )}(x)dx={\partial }_n^{\lambda ,\nu }{\delta }_{mn},$$

$${\partial }_n^{\lambda ,\nu }=\Vert P_n^{(\lambda ,\nu )}{\Vert }_{w^{\lambda ,\nu }}^2,w^{\lambda ,\nu }(x)=(1-x{)}^{\lambda }(1+x{)}^{\nu },$$

$$\delta_{mn}=\left\{\begin{array}[]{ll}1,&m=n,\\ 0,&m\neq n\end{array}\right.$$

$$\begin{array}{c}{P}_0^{(\lambda ,\nu )}(x)=1,P_1^{(\lambda ,\nu )}(x)=\frac{1}{2}(\lambda +\nu +2)x+\frac{1}{2}(\lambda -\nu ),\\ P_{n+1}^{(\lambda ,\nu )}(x)=\left(a_n^{\lambda ,\nu }x-b_n^{\lambda ,\nu }\right)P_n^{(\lambda ,\nu )}(x)-c_n^{\lambda ,\nu }{P}_{n-1}^{(\lambda ,\nu )}(x),n\ge 2,\end{array}$$

$$a_n^{\lambda ,\nu }$$

$$b_n^{\lambda ,\nu }$$

$$c_n^{\lambda ,\nu }$$

$$P_n^{(\lambda ,\nu )}(x)=\sum _{k=0}^n\frac{2^{n-k}n!(k+\nu +1{)}_{n-k}}{(n-k)!(n+k+\lambda +\nu +1{)}_{n-k}k!}(t-1{)}^k,$$

$$(a{)}_l=a(a+1)(a+2)\cdots (a+l-1)$$

$$P_n^{(\lambda ,\nu )}(-x)=(-1{)}^n{P}_n^{(\lambda ,\nu )}(x),$$

$$\frac{d}{dx}{P}_n^{(\lambda ,\nu )}(-x)=(-1{)}^n\frac{d}{dx}{P}_n^{(\lambda ,\nu )}(x),$$

$$\begin{gathered} (2n+\lambda +\nu )(1-x^2)\frac{d}{dx}{P}_n^{(\lambda ,\nu )}(x) \\ =n\left(\lambda -\nu -(2n+\lambda +\nu )x\right)P_n^{(\lambda ,\nu )}(x)+2(n+\lambda )(n+\nu )P_{n-1}^{(\lambda ,\nu )}(x), \end{gathered}$$

$$P_n^{(\lambda ,\nu )}(1)=\left(\genfrac{}{}{0ex}{}{n+\lambda }{n}\right),P_n^{(\lambda ,\nu )}(-1)=(-1{)}^j\frac{\Gamma (n+\nu +1)}{n!\Gamma (\nu +1)}.$$

$${\int }_{-1}^1(1-x{)}^{\lambda }(1+x{)}^{\nu }f(x)dx=\sum _{k=1}^n{A}_{k}f(x_k)+R_n(f),$$

$$A_k=\frac{2^{\lambda +\nu +1}\Gamma (\lambda +n+1)\Gamma (\nu +n+1)}{n!\Gamma (\lambda +\nu +n+1){\left(\frac{d}{dx}{P}_n^{(\lambda ,\nu )}(x_k)\right)}^2(1-x_k^2)}$$

$$R_n(f)=\frac{f^{(2n)}(\eta )}{(2n)!}\frac{2^{\lambda +\nu +2n+1}n!\Gamma (\lambda +\eta +1)\Gamma (\nu +\eta +1)\Gamma (\nu +\lambda +\eta +1)}{(\lambda +\nu +2n+1){\left(\Gamma (\lambda +\nu +2n+1)\right)}^2}$$

$${\int }_{-1}^1{w}^{\alpha -1}(x)f_n^{\alpha -1}(x)f_m^{\alpha -1}(x)dx={\zeta }_n^{\alpha -1}{\delta }_{mn},$$

$${\zeta }_n^{\alpha -1}=\Vert f_n^{\alpha -1}{\Vert }_{w^{\alpha -1}}^2,w^{\alpha -1}(x)=(1-x{)}^{\alpha -1}.$$

Full text

Part of first author’s PhD project. Partially supported by Islamic Azad University, Tehran, Iran; and CIDMA-FCT, Aveiro, Portugal, within project UID/MAT/04106/2013.

1. Introduction

Finding numerical approximations of fractional integrals or fractional derivatives of a given function is one of the most important problems in theory of numerical fractional calculus. The operators of fractional integration and fractional differentiation are more complicated than the classical ones, so their evaluation is also more difficult than the integer order case. Li et al. use spectral approximations for computing the fractional integral and the Liouville–Caputo derivative [22]. They also developed numerical algorithms to compute fractional integrals and Liouville–Caputo derivatives and for solving fractional differential equations based on piecewise polynomial interpolation [21]. In [29], Pooseh et al. presented two approximations derived from continuous expansions of Riemann–Liouville fractional derivatives into series involving integer order derivatives and they present application of such approximations to fractional differential equations and fractional problems of the calculus of variations. Some other computational algorithms are also introduced in [24, 28, 32]. For increasing the accuracy of the calculation, using the Gauss–Jacobi quadrature rule is appropriate for removing the singularity of the integrand. So considering the nodes and weights of the quadrature rule is an important problem. There are many good papers in the literature addressing the question of how to find the nodes and weights of the Gauss quadrature rule—see [6, 36, 37] and references therein. The more applicable and developed method is the Golub–Welsch (GW) algorithm [14, 15], that is used by many of the mathematicians who work in numerical analysis. This method takes $O(n^{2})$ operations to solve the problem of finding the nodes and weights. Here we use a new method introduced by Hale and Townsend [16], which is based on the Glasier–Liu–Rokhlin (GLR) algorithm [13]. It computes all the nodes and weights of the $n$-point quadrature rule in a total of $O(n)$ operations.

The structure of the paper is as follows. In Section 2, we introduce the definitions of fractional operators and some relations between them. Section 3 discusses the Gauss–Jacobi quadrature rule of fractional order and its application to approximate the fractional operators. In Section 4 we present two methods for finding the nodes and weights of Gauss–Jacobi and discuss their advantages and disadvantages. Two illustrative examples are solved. In Section 5 applications to ordinary fractional differential equations are presented. In Section 6 we investigate problems of the calculus of variations of fractional order and present a new algorithm for solving boundary value problems of fractional order. We end with Section 7 of conclusions and possible directions of future work.

2. Preliminaries and notations about fractional calculus

In this section we give some necessary preliminaries of the fractional calculus theory [18, 27], which will be used throughout the paper.

Definition 2.1**.**

The left and right Riemann–Liouville fractional integrals of order $\alpha$ of a given function $f$ are defined by

[TABLE]

and

[TABLE]

respectively, where $\Gamma$ is Euler’s gamma function, that is,

[TABLE]

$\alpha>0$ with $n-1<\alpha\leq n$, $n\in\mathbb{N}$, and $a<x<b$. The left Riemann–Liouville fractional operator has the following properties:

[TABLE]

[TABLE]

where $\alpha,\beta\geq 0$ and $\mu>-1$. Similar relations hold for the right Riemann–Liouville fractional operator. On the other hand, we have the left and right Riemann–Liouville fractional derivatives of order $\alpha>0$ that are defined by

[TABLE]

and

[TABLE]

respectively.

There are some disadvantages when trying to model real world phenomena with fractional differential equations, when fractional derivatives are taken in Riemann–Liouville sense. One of them is that the Riemann–Liouville derivative of the constant function is not zero. Therefore, a modified definition of the fractional differential operator, which was first considered by Liouville and many decades later proposed by Caputo [7], is considered.

Definition 2.2**.**

The left and right fractional differential operators in Liouville–Caputo sense are given by

[TABLE]

and

[TABLE]

respectively.

The Liouville–Caputo derivative has the following two properties for $n-1<\alpha\leq n$ and $f\in L_{1}[a,b]$:

[TABLE]

and

[TABLE]

Remark 2.3*.*

Using the linearity property of the ordinary integral operator, one deduces that left and right Riemann–Liouville integrals, left and right Riemann–Liouville derivatives and left and right Liouville–Caputo derivatives are linear operators.

Another definition of a fractional differential operator, that is useful for numerical approximations, is the Grünwald–Letnikov derivative, which is a generalization of the ordinary derivative. It is defined as follows:

[TABLE]

3. Fractional Gauss–Jacobi quadrature rule

It is well known that the Jacobi polynomials $\{P_{n}^{(\lambda,\nu)}(x)\}_{n=0}^{\infty}$, $\lambda,\nu>-1$, $x\in[-1,1]$, are the orthogonal system of polynomials with respect to the weight function

[TABLE]

on the segment $[-1,1]$:

[TABLE]

where

[TABLE]

[TABLE]

(see, e.g., [35]). These polynomials satisfy the three-term recurrence relation

[TABLE]

where

[TABLE]

The explicit form of the Jacobi polynomials is

[TABLE]

where we use Pchhammer’s notation:

[TABLE]

(see [20]). Furthermore, the Jacobi polynomials satisfy in the following relations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

On the other hand, the Gauss–Jacobi quadrature formula with remainder term is given by

[TABLE]

where $x_{k}$, $k=1,\ldots,n$, are the zeros of $P_{n}^{(\lambda,\nu)}$, the weights $A_{k}$ are given by

[TABLE]

and the error is given by

[TABLE]

(see [35]). We know that formula (5) is exact for all polynomials of degree up to $2n-1$ and is valid if $f(x)$ possesses no singularity in $[-1,1]$ except at points $\pm 1$.

Now, as a special case, we introduce the fractional Jacobi polynomials $f_{n}^{\alpha-1}(x)$, that is, $\lambda=\alpha-1$ and $\nu=0$ in (1). The set of fractional Jacobi polynomials $\{f_{n}^{\alpha-1}\}_{n=0}^{\infty}$ is an orthogonal system of polynomials that are orthogonal with respect to the weight function $(1-x)^{\alpha-1}$, $\alpha>0$, on the segment $[-1,1]$. It means that

[TABLE]

where

[TABLE]

The three-term recurrence relation for fractional Jacobi polynomial is given by

[TABLE]

where

[TABLE]

Using (4), we have that

[TABLE]

Furthermore, the fractional Gauss–Jacobi quadrature rule is

[TABLE]

where $x_{k}$, $k=1,2,\ldots,n$, are the zeros of $f_{n}^{\alpha-1}$, the weights $l_{k}$ are given by

[TABLE]

and

[TABLE]

Consider the left Riemann–Liouville fractional integral of order $\alpha>0$ for $x>0$:

[TABLE]

For calculating the above integral at the collocation nodes, first we transform the interval $[0,x]$ into the segment $[-1,1]$ using the changing of variable

[TABLE]

Therefore, we can rewrite (8) as

[TABLE]

where

[TABLE]

So, we have

[TABLE]

and using the $n$ point fractional Gauss–Jacobi quadrature rule (7) for (9) we get

[TABLE]

where

[TABLE]

[TABLE]

and $x_{k}$, $k=1,\ldots,n$, are quadrature nodes. In this way we just need to compute the nodes and weights of the fractional Gauss–Jacobi quadrature rule using a fast and accurate algorithm. Similarly, an approximation formula for computing the left Liouville–Caputo fractional derivative of a smooth function $f$, using the suggested method, is given by

[TABLE]

where $x_{k}$, $k=1,2,\ldots,n$, are the zeros of $f_{n}^{-\alpha}$,

[TABLE]

and

[TABLE]

Analogous formulas hold for the right Liouville–Caputo fractional derivative. Next we present two methods for computing the nodes and weights of the fractional Gauss–Jacobi quadrature rule (7). The first one is based on Golub–Welsch algorithm, while the second is a recent method introduced by Hale and Townsend [16], which is based on Newton’s iteration for finding roots and a good asymptotic formula for the weights.

4. Two methods for calculating nodes and weights

Pang et al. [26] use Gauss–Jacobi and Gauss–Jacobi–Lobatto quadrature rules for computing fractional directional integrals, together with the Golub–Welsch (GW) algorithm for computing nodes and weights of quadrature rules. The GW algorithm exploits the three-term recurrence relations (2) satisfied by all real orthogonal polynomials. This relation gives rise to a symmetric tridiagonal matrix,

[TABLE]

where

[TABLE]

in which

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It is easy to prove that $x_{i}$ is the zero of $P_{n}^{(\lambda,\nu)}$ if and only if $x_{i}$ is the eigenvalue of the tridiagonal matrix [26]. Moreover, Golub and Welcsh proved that the weights of quadrature are the first component of corresponding eigenvectors [15]. This algorithm takes $O(n^{2})$ operations to solve the eigenvalue problem by taking advantage of the structure of the matrix and noting that only the first component of the normalized eigenvector need to be computed. An alternative approach for computing the nodes and weights of the Gauss–Jacobi rule is to use the same three-term recurrence relation in order to compute Newton’s iterates, which converge to the zeros of the orthogonal polynomial [30]. Since the recurrence requires $O(n)$ operations for each evaluation of the polynomial and its derivative, we have the total complexity of order $O(n^{2})$ for all nodes and weights. Furthermore, the relative maximum error in the weights is of order $O(n)$ and for the nodes is independent of $n$. Here we develop the new technique introduced by Hale and Townsend, which utilizes asymptotic formulas for both accurate initial guesses of the roots and efficient evaluation of the fractional Jacobi polynomial $f_{n}^{\alpha-1}$ inside Newton’s method [16]. With this method it is possible to compute the nodes and weights of the fractional Gauss–Jacobi rule in just $O(n)$ operations and almost full double precision of order $O(1)$. For performance of the method, first put $\theta_{k}=\cos^{-1}x_{k}$ to avoid the existing clustering in computing $w_{k}$, because of presence of the term $(1-x_{k}^{2})$ in the denominator of (11) that lead to cancellation error for $x_{k}=\pm 1$. Then, using a simple method, like the bisection one, for finding $\theta_{k}^{[0]}$ as the initial guess of the $k$th root, we construct the Newton iterates for finding the nodes as follows:

[TABLE]

Once the iterates have converged, the nodes are given by $x_{k}=\cos\theta_{k}$ and using (11) the weights are given by

[TABLE]

Since the zeros of the orthogonal fractional Jacobi polynomial are simple [35], we conclude that the Newton iterates converge quadratically [19]. Furthermore, since all fractional Jacobi polynomials satisfy the relation (3), we just need to consider all calculations for $x\in[0,1]$, i.e., $\theta\in[0,\frac{\pi}{2}]$. There are three asymptotic formulas for finding the nodes.

One is given by Gatteschi and Pittaluga [12] for Jacobi polynomials and here to the case of fractional Jacobi polynomials:

[TABLE]

$|\alpha|\leq\frac{1}{2}$, where $\rho=n+\frac{\alpha}{2}$ and $f_{k}=(k+\frac{\alpha}{2}-\frac{3}{4})\frac{\pi}{\rho}$. 2. 2)

Let $j_{\alpha-1,k}$ denote the $k$th root of Bessel’s function $J_{\alpha-1}(z)$. Then the approximation given by Gatteschi and Pittaluga [12] for the nodes near $x=1$ becomes

[TABLE]

$\alpha\in[-\frac{1}{2},\frac{1}{2}]$, where $\gamma=\sqrt{\rho^{2}+(2-\alpha^{2}-2\alpha)/12}$. 3. 3)

Other asymptotic formulas for the nodes near $x=1$ are given by Olver et al. [25]:

[TABLE]

$\alpha>-\frac{1}{2}$, where $\psi_{k}=\dfrac{j_{\alpha-1,k}}{\rho}$.

For computing the weights $l_{k}$ we just need to evaluate $\dfrac{d}{d\theta}f_{n}^{\alpha-1}(\cos\theta)$ in the $\theta_{k}$. To this end, we use relation (6) between $f_{n}^{\alpha-1}$ and $(f_{n}^{\alpha-1})^{\prime}$. So we just need to compute the value of $f_{n}^{\alpha-1}$ at $x=\cos\theta$. This work is done by using an asymptotic formula introduced in [16], which takes the following form for fractional Jacobi polynomials:

[TABLE]

where $\rho=n+\alpha/2$, $B(\alpha,\beta)$ is the Beta function,

[TABLE]

[TABLE]

and

[TABLE]

where for $\alpha\in(-1/2,1/2)$ the error term $V_{M,n}$ is less than twice the magnitude of the first neglected term. We give two examples illustrating the usefulness of the method.

Example 4.1**.**

The fractional integral of the function $f(t)=\sin(t)$ is defined as

[TABLE]

The explicit expression of the integral, presented in Table 9.1 of [34], is

[TABLE]

where $i=\sqrt{-1}$ and ${}_{1}F_{1}$ is the generalized hyper-geometric function defined by

[TABLE]

Choosing the test points in the set $\{s_{k}=\frac{k\pi}{8},k=0,1,\ldots,16\}$, the error is given by

[TABLE]

where $E_{k}$ is the exact result at $s_{k}$ and $A_{k}$ is the approximated result by our method. The values of (13) for $\alpha=0.25$, $0.5$ and $0.75$ are listed in the Table 1. We used the command jacpts of the chebfun software [16] for finding the nodes and weights of the fractional Gauss–Jacobi quadrature (fgjq) rule. The errors show that the method is more accurate than the method introduced in [26].

Example 4.2**.**

In this example we consider the function $f(t)=t^{4}$. The Liouville–Caputo fractional derivative of $f$ of order $0<\alpha<1$ is defined as

[TABLE]

The explicit form of the result is

[TABLE]

(see [29]). Using (12), we can write

[TABLE]

Here the nodes and weights coincide with the nodes and weights of Example 4.1 for $\alpha=0.5$. With the test point set $\{s_{k}=\frac{k}{10},k=0,1,2,\ldots,10\}$ the error defined by (13) is listed in Table 2. Looking to Table 2, one can see that the best approximation is obtained for $n=2$ and that by increasing $n$ we do not get better results.

5. Application to fractional ordinary differential equations

Consider the fractional initial value problem

[TABLE]

where $y(t)$ is the unknown function, $g:\mathbb{R}^{2}\to\mathbb{R}$ is given and $b_{k}$, $k=0,1,\ldots,n-1$, are real constants. The following result establishes existence and uniqueness of solution for (14).

Theorem 5.1** (See [18, 27]).**

Consider the fractional differential equation (14). Let $K>0$, $h^{*}>0$, $b_{0},b_{1},\ldots,b_{n-1}\in\mathbb{R}$, $G:=[0,h^{*}]\times[b_{0}-K,b_{0}+K]$, and function $g:G\to\mathbb{R}$ be continuous. Then, there exists $h>0$ and a function $y\in C[0,h]$ solving the Liouville–Caputo fractional initial value problem (14). If $0<\alpha<1$, then the parameter $h$ is given by

[TABLE]

Furthermore, if $g$ fulfils a Lipschitz condition with respect to the second variable, that is,

[TABLE]

for some constant $L>0$ independent of $t,y_{1}$ and $y_{2}$, then the function $y\in C[0,h]$ is unique.

In order to solve problem (14) using our method, we need the following theorem.

Theorem 5.2**.**

Under the assumptions of Theorem 5.1, function $y\in C[0,h]$ is a solution to the Liouville–Caputo fractional differential equation (14) if and only if it is a solution to the Volterra integral equation of second kind

[TABLE]

Proof.

Using the Laplace transform formula for the Liouville–Caputo fractional derivative,

[TABLE]

(see [18, 27]). Thus, using (14), we can write that

[TABLE]

or

[TABLE]

where $G(s,y(s))=\mathcal{L}\{g(t,y(t))\}$. Applying the inverse Laplace transform gives

[TABLE]

where we used the relations

[TABLE]

and

[TABLE]

The proof is complete. ∎

Now, using (9), (10) and Theorem 5.2, we have

[TABLE]

where

[TABLE]

[TABLE]

So, using the fractional Gauss–Jacobi quadrature rule introduced in Section 2, we can approximate $y(t)$ by the following formula:

[TABLE]

where $l_{k}$ and $x_{k}$ are the weights and nodes of our quadrature rule, respectively. Also, let $t\in[0,1]$ and $t_{j}=\frac{j}{n}$ for $j=0,1,\ldots,n$. Furthermore, let the numerical values of $y(t)$ at $t_{0},t_{1},\ldots,t_{n}$ have been obtained, and let these values be $y_{0},y_{1},\ldots,y_{n}$, respectively, with $y_{0}=b_{0}$. Now we are going to compute the value of $y(t)$ at $t_{n+1}$, i.e., $y_{n+1}$. To this end we use the predictor-corrector method introduced by Diethelm et al. in [10], but where for the prediction part we use interpolation at equispaced nodes instead of piecewise linear interpolation. Using first (15), we know that

[TABLE]

where $l_{k}$ and $x_{k}$, $k=1,2,\ldots,N$, are the weights and nodes of the fractional Gauss–Jacobi quadrature rule introduced in Section 3. Looking to (16) we see that for calculating the value of $y(t_{n+1})$ we need to know the value of $g$ at point $\left(\frac{t_{n+1}}{2}(1+x_{k}),y\left(\frac{t_{n+1}}{2}(1+x_{k})\right)\right)$. For obtaining these values we act as follows. First, using the predictor-corrector method introduced in [9], we calculate the values of $y_{1}(t_{i})$, $i=0,1,\ldots,n+1$, using the following formulas:

[TABLE]

where

[TABLE]

and

[TABLE]

and $k$ is the step length, i.e., $k=t_{j+1}-t_{j}$ (for more details, see [9]). Finally, if $P_{n}$ is the interpolating polynomial passing through $\{(t_{i},y1(t_{i})\}_{i=0}^{n+1}$, then we have the following formula for calculating $y(t_{n+1})$:

[TABLE]

We solve two examples illustrating the applicability of our method.

Example 5.3**.**

Consider the nonlinear ordinary differential equation

[TABLE]

subject to $y(0)=0$ and $y^{\prime}(0)=0$, where

[TABLE]

Following [17], we take $k=1$ and $\alpha=\frac{3}{2}$. The exact solution is then given by

[TABLE]

Figure 1 plots the results obtained for $N=16$ in (17).

Example 5.4**.**

Consider the fractional oscillation equation

[TABLE]

subject to initial conditions $y(0)=0$ and $y^{\prime}(0)=0$. The exact solution is

[TABLE]

where $E_{\alpha,\beta}(t)$ is the generalized Mittag–Lefller function [17]. Results are shown in Figure 2 for $N=16$ in (17) and $\alpha=\frac{3}{2}$.

6. Application to fractional variational problems

In this section we apply our method to solve, numerically, some problems of the calculus of variations of fractional order. The calculus of variations is a rich branch of classical mathematics dealing with the optimization of physical quantities (such as time, area, or distance). It has applications in many diverse fields, including aeronautics (maximizing the lift of an aircraft wing), sporting equipment design (minimizing air resistance on a bicycle helmet or optimizing the shape of a ski), mechanical engineering (maximizing the strength of a column, a dam, or an arch), boat design (optimizing the shape of a boat hull) and physics (calculating trajectories and geodesics, in both classical mechanics and general relativity) [33]. The problem of the calculus of variations of fractional order is more recent (see [2, 23, 31] and references therein) and consists, typically, to find a function $y$ that is a minimizer of a functional

[TABLE]

subject to boundary conditions

[TABLE]

One can deduce fractional necessary optimality equations to problem (18) of Euler–Lagrange type: if function $y$ is a solution to problem (18), then it satisfies the fractional Euler–Lagrange differential equation

[TABLE]

It is often hard to find the analytic solution to problems of the calculus of variations of fractional order by solving (19). Therefore, it is natural to use some numerical method to find approximations of the solution of (19). Here we are going to apply our method for solving some examples of such fractional variational problems. To this end, at first we find a class of fractional Lagrangians

[TABLE]

such that the corresponding Euler–Lagrange equation (19) takes the form

[TABLE]

We assume a solution to this problem as follows:

[TABLE]

Then, we evaluate the left-hand side of (20) for (21) and we obtain

[TABLE]

So, if the Lagrangian of the variational problem (18) is of form (21), then we just need to solve the following boundary value problem of fractional order:

[TABLE]

subject to the boundary conditions $y(a)=u_{a}$ and $y(b)=u_{b}$, where

[TABLE]

Now, we apply the operator ${{}_{t}I_{b}^{\alpha}}$ to both sides of equation (22). We get the following boundary value problem of fractional order:

[TABLE]

Theorem 6.1**.**

Consider the Lagrangian (21) be such that $g(t,y(t))$ given by (23) is continuous and let $-{{}_{t}I_{b}^{\alpha}}g(t,y(t))=h(t,y(t))$. Then the function $h:[a,b]\times\mathbb{R}\to\mathbb{R}$ is continuous and the problem (24) is equivalent to the integral equation

[TABLE]

Proof.

The continuity of $h$ is one of the properties of the Riemann–Liuoville fractional integral operator ${{}_{t}I_{b}^{\alpha}}$ (see, for example, [18, 27]). The second part of the theorem is the main result of [38]. ∎

Now, let $t_{j}=jk$, $j=0,1,\ldots,n$, and $k=\frac{b-a}{n}$. Furthermore, let the numerical values of $y(t)$ at $t_{0},t_{1},\ldots,t_{n}$ have been obtained and let these values be $y_{0},y_{1},\ldots,y_{n}$, respectively, and $y_{0}=u_{a}$. Now we are going to compute the value of $y(t)$ at $t_{n+1}$, i.e., $y_{n+1}$. To this end, first we use the same idea of [11] to predict the values of $y_{0},y_{1},\ldots,y_{n}$. Then, we correct them using our method. Let

[TABLE]

and

[TABLE]

Then, using (25), we have

[TABLE]

Therefore,

[TABLE]

Finally, if $P_{n}$ is the interpolation polynomial of the nodes $\{(t_{i},y_{i})\}_{i=0}^{n+1}$, then the following correction formula is an approximation of the solution to problem (24) at the point $t_{n+1}$:

[TABLE]

where $x_{k}$ and $l_{k}$, $k=1,2,\ldots,N$, are the nodes and the weights of the fractional Gauss–Jacobi quadrature rule. Now we are going to apply this method to solve a concrete fractional problem of the calculus of the variations.

Example 6.2**.**

Consider the following Fractional Variational Problem (FVP): to find the minimizer of the functional

[TABLE]

subject to

[TABLE]

When $\alpha\rightarrow 1$, this FVP tends to the classical problem of the calculus of variations

[TABLE]

which has the solution

[TABLE]

On the other hand, the fractional Euler–Lagrange equation associated with (26) is

[TABLE]

Looking to Figure 3 and Figure 4, we can see that the numerical solutions to the fractional Euler–Lagrange equation (30) subject to boundary conditions (27) are approximating the solution (29) to the classical variational problem (28). Note that the exact solution to the FVP (26)–(27) is unknown.

7. Conclusions

We presented a new numerical method for approximating Riemann–Liouvile fractional integrals and Liouville–Caputo fractional derivatives. We applied the method for solving fractional linear and nonlinear initial value problems. Furthermore, we introduced a new method for solving fractional boundary value problems and applied this new method for solving fractional Euler–Lagrange equations. Finally, by considering a concrete problem of the calculus of variations, we showed that the method is useful to solve fractional variational problems.

The results of the paper can be extended in several directions. As an example of further possible developments in the field, we can mention the new fractional derivative with non-local and non-singular kernel proposed by Atangana and Baleanu [4]. Such derivative seems to answer some outstanding questions that were posed within the field of fractional calculus, giving origin to chaotic behaviors [5]. The importance to develop new numerical methods to deal with such fractional calculus is well presented in [1]. We claim that our results can be generalized to cover the Atangana–Baleanu calculus. Similarly can be said about numerical simulations with the Caputo–Fabrizio fractional order derivative [3, 8].

Acknowledgements

This work is part of first author’s PhD project. Partially supported by Islamic Azad University (Science and Research branch of Tehran), Iran; and CIDMA–FCT, Portugal, within project UID/MAT/04106/2013. Jahanshahi was supported by a scholarship from the Ministry of Science, Research and Technology of the Islamic Republic of Iran, to visit the University of Aveiro, Portugal. The hospitality and the excellent working conditions at the University of Aveiro are here gratefully acknowledged. The authors are grateful to two anonymous referees, for several comments and suggestions, which helped them to improve the manuscript.