Sage Journals: Discover world-class research

Abstract

In this article, we construct a Jacobi spectral collocation scheme to approximate the Caputo fractional derivative based on Jacobi–Gauss quadrature. The convergence analysis is provided in anisotropic Jacobi-weighted Sobolev spaces. Furthermore, the convergence rate is presented for solving Caputo fractional derivative with noisy data by invoking the mollification regularization method. Lastly, numerical examples illustrate the effectiveness and stability of the proposed method.

Keywords

Fractional derivative Jacobi collocation mollification Gaussian quadrature

1. Introduction

Fractional calculus and derivative is a branch of mathematics that deals with the generalization of integration and differentiation to non-integer orders. It has been verified that substantial physical systems in science and engineering such as fractals, anomalous diffusion, viscoelastic materials, porous or fractured media and viscous fluid flows, can be modeled more effectively by invoking fractional-order rather than integer-order derivatives [2,14,22] due to the cumulative memory effect of wall-friction [5,12,21].

There are extensive amount of work to calculate the fractional derivative numerically over the recent decades years. Lin and Xu [13] developed a hybrid scheme for time-fractional diffusion problem, discretizing the time-fractional derivative by using Finite Difference Method and treating the integer-order spatial derivative by a Legendre spectral method, and rigorously established that the global error is dominated by the low-order temporal accuracy. Zayernouri and Karniadakis [23–25] first introduced fractional Lagrange interpolants to construct the exponentially accurate fractional spectral collocation method to solving fractional Sturm–Liouville eigenproblems and linear/nonlinear fractional partial differential equations with field-variable order. And their numerical results well confirmed the exponential convergence of the fractional collocation methods. Esmaeili and Milovanović[7] derived a numerical scheme to calculate the usual fractional derivative by applying nonstandard Gauss–Jacobi–Lobatto quadratures, and implemented the scheme by using a Mathematica package of Orthogonal Polynomials but without any convergence analysis. Chen et al. [3] constructed Petrov–Galerkin methods for a class of fractional differential equations by employing generalized Jacobi functions. Their numerical results and theoretical analysis demonstrated exponential convergence for the exact solution to be non-smooth except specific singularity. Chen et al. [4] proposed a fast algorithm to solve a subdiffusion problem, which contains Caputo time-fractional derivative with order $ν \in (0, 1)$ . They developed a spectral-Galerkin method in time with Log orthogonal functions by applying a log mapping to the Laguerre functions and the Galerkin finite element method in space. They obtained an error estimate with spectral accuracy in time for weakly singular solutions and numerically confirmed the theoretical results. Mao and Karniadakis [15] studied a one-dimensional diffusion equation with general two-sided fractional derivative. They derived a Sturm–Liouville eigenproblem with the two-sided Jacobi polyfractonomials as eigenfunctions by employing a Petrov–Galerkin projection method. They obtained an exponential convergence and numerical examples verified the theoretical results. Dehghan and Mingarelli [6] presented a general solution of three distinct fractional Sturm–Liouville eigenvalue problems, which consist of various compositions of Riemann–Liouville and Caputo fractional derivatives. Nie et al. [18] numerically approximated a fractional sub-diffusion equation with Laplace operator, by employing a modified central finite difference scheme in spatial direction and two kinds of averaged schemes(averaged $L 1$ and averaged second order backward difference scheme) to discretize the Riemann–Liouville fractional derivative in temporal direction. And obtain a better spatial convergence $O (h^{min (σ + \frac{1}{2} - ϵ, 2)})$ instead of $O (h^{2})$ , and an temporal convergence $O (τ^{2})$ for all fractional order $α \in (0, 1)$ . Jin and Zhou [11] considered an inverse problem of reconstructing a spatial-temporal diffusion coefficient in the sub-diffusion model from the distributed observation. They calculated numerically the direct problem by using the backward Euler convolution quadrature in time to approximate the Caputo fractional derivative of order $α \in (0, 1)$ , and the Galerkin Finite Element Method in space. They provided a complete error analysis of the fully discrete algorithm and established two conditional stability results theoretically by invoking a novel test function, and numerical experiments in agreement with the theoretical analysis. Motivated by the link between singular kernel of Caputo derivative and weight functions of Jacobi polynomials, we desire to study a fast and high accuracy numerical scheme for fractional derivative with order $μ \in (0, 1)$ by using Jacobi spectral collocation method, and provide the convergence analysis for calculating numerical fractional derivative from data with noise by mollification regularization method.

The rest of this paper is organized as follows. In Section 2, we introduce some results about Gaussian quadrature and Jacobi orthogonal polynomials. In Section 3, we construct the Jacobi collocation method for Caputo fractional derivatives using the Jacobi weights and polynomials. The convergence analysis for the numerical algorithm are presented in Section 4. In Section 5, the convergence rate is described by virtue of the mollification regularization method when data function f contains noise. Several numerical tests are performed to demonstrate the accuracy and efficiency of the proposed method in Section 6. The paper is concluded in Section 7.

2. Preliminaries

For the purpose of readability, we state some properties and conclusions about the Gaussian quadrature and Jacobi polynomials in this section. We refer the reader to, for example, [1,8,10,20] and [9] for more details.

The mechanism of a Gauss-type quadrature is to seek the best numerical approximation of an integral by selecting optimal nodes at which the integrand is evaluated. It belongs to the family of the numerical quadratures: $\begin{aligned} (1) & \int_{a}^{b} ω (x) f (x) d x = \sum_{k = 0}^{N} f (x_{k}) ω_{k} + E_{N} [f], \end{aligned}$ where $ω (x)$ is the weight function, ${x_{k}, ω_{k}}_{k = 0}^{N}$ are the quadrature nodes and weights, and if $f (x) \in C^{N + 1} [a, b]$ , the truncated error satisfies $\begin{aligned} (2) & E_{N} [f] = \frac{1}{(N + 1)!} \int_{a}^{b} f^{(N + 1)} (ξ (x)) \prod_{k = 0}^{N} (x - x_{k}) d x, ξ (x) \in [a, b] . \end{aligned}$

The Jacobi polynomials with parameters α and β, denoted by $J_{n}^{α, β} (x)$ , are orthogonal with respect to the Jacobi weight function $ω^{α, β} (x) : = {(1 - x)}^{α} {(1 + x)}^{β}$ over the interval $I : = (- 1, 1)$ , namely, $\begin{aligned} (3) & \int_{- 1}^{1} ω^{α, β} (x) J_{n}^{α, β} (x) J_{m}^{α, β} (x) d x = \{\begin{array}{ll} 0, & m \neq n, \\ ‖ J_{n}^{α, β} ‖_{ω^{α, β}}^{2}, & m = n, \end{array} \end{aligned}$ where the weight function $ω^{α, β}$ belongs to $L^{1} (I)$ if and only if $α, β > - 1$ .

In fact, the Jacobi polynomials are the eigenfunctions of a singular Sturm–Liouville operator defined by $\begin{aligned} (4) & L_{α, β} u : = - {(1 - x)}^{- α} {(1 + x)}^{- β} \frac{d}{d x} ({(1 - x)}^{α + 1} {(1 + x)}^{β + 1} \frac{d}{d x} u (x)) . \end{aligned}$ More precisely, the following properties of Jacobi polynomials are hold.

Lemma 1 ([20, Theorem 3.16]).

The Jacobi polynomials are the eigenfunctions of the singular Sturm–Liouville problem: $\begin{aligned} (5) & L_{α, β} J_{n}^{α, β} (x) = λ_{n}^{α, β} J_{n}^{α, β} (x), \end{aligned}$ and the corresponding eigenvalues are $\begin{aligned} (6) & λ_{n}^{α, β} = n (n + α + β + 1) . \end{aligned}$

Lemma 2 ([20, Recurrence Formulas]).

The Jacobi polynomials are generated by the three-term recurrence relation: $\begin{aligned} (7) & \begin{aligned} J_{n + 1}^{α, β} (x) = (a_{n}^{α, β} x - b_{n}^{α, β}) J_{n}^{α, β} (x) - c_{n}^{α, β} J_{n - 1}^{α, β} (x), n ⩾ 1, \\ J_{1}^{α, β} (x) = 1, J_{1}^{α, β} (x) = \frac{1}{2} (α + β + 2) x + \frac{1}{2} (α - β) . \end{aligned} \end{aligned}$ where $\begin{aligned} (8) & a_{n}^{α, β} = \frac{(2 n + α + β + 1) (2 n + α + β + 2)}{2 (n + 1) (n + α + β + 1)}, \\ (9) & b_{n}^{α, β} = \frac{(β^{2} - α^{2}) (2 n + α + β + 1)}{2 (n + 1) (n + α + β + 1) (2 n + α + β)}, \\ (10) & c_{n}^{α, β} = \frac{(n + α) (n + β) (2 n + α + β + 2)}{(n + 1) (n + α + β + 1) (2 n + α + β)} . \end{aligned}$

Lemma 3 ([20, Theorem 3.25]).

Let Jacobi–Gauss nodes ${ξ_{k}}_{k = 0}^{N}$ being the zeros of the Jacobi polynomials $J_{N + 1}^{α, β} (x)$ and the corresponding weights given by $\begin{aligned} ω_{k} & = \frac{G_{N}^{α, β}}{J_{N}^{α, β} (ξ_{k}) \partial_{x} J_{N + 1}^{α, β} (ξ_{k})} \\ (11) & = \frac{{\tilde{G}}_{N}^{α, β}}{(1 - ξ_{k}^{2}) {[\partial_{x} J_{N + 1}^{α, β} (ξ_{k})]}^{2}}, \end{aligned}$ where $\begin{aligned} (12) & G_{N}^{α, β} = \frac{2^{α + β} (2 N + α + β + 2) Γ (N + α + 1) Γ (N + β + 1)}{(N + 1)! Γ (N + α + β + 2)}, \\ (13) & {\tilde{G}}_{N}^{α, β} = \frac{2^{α + β + 1} Γ (N + α + 2) Γ (N + β + 2)}{(N + 1)! Γ (N + α + β + 2)}, \end{aligned}$ then the following Jacobi–Gauss quadrature formula $\begin{aligned} (14) & \int_{- 1}^{1} ω^{α, β} (x) f (x) d x = \sum_{k = 0}^{N} ω_{k} f (ξ_{k}) + E_{N} [f], \end{aligned}$ is exact (namely, $E_{N} [f] = 0$ ) for any $f (x) \in P_{2 N + 1}$ with Jacobi–Gauss quadrature pairs ${ξ_{k}, ω_{k}}_{k = 0}^{N}$ .

Lemma 4 ([20, Differentiation Matrix in the Physical Space]).

Let ${ξ_{k}}_{k = 0}^{N}$ be a set of Jacobi–Gauss-type points, and ${h_{k}}_{k = 0}^{N}$ be the associated Lagrange basis polynomials. Assume that $f (x) \in P_{N}$ is an approximation of the unknown solution then $\begin{aligned} (15) & f (x) = \sum_{j = 0}^{N} h_{j} (x) f (ξ_{j}), f^{(m)} (ξ_{k}) = \sum_{j = 0}^{N} h_{j}^{(m)} (ξ_{k}) f (ξ_{j}), 0 ⩽ k ⩽ N . \end{aligned}$ Denote $\begin{aligned} (16) & \begin{aligned} f^{(m)} : = {(f^{(m)} (ξ_{0}), f^{(m)} (ξ_{1}), \dots, f^{(m)} (ξ_{N}))}^{T}, f : = f^{(0)}, \\ D^{(m)} : = {(d_{k j}^{(m)})}_{k, j = 0}^{N}, D : = D^{(1)}, \end{aligned} \end{aligned}$ then the following properties holds $(\forall m ⩾ 1)$ $\begin{aligned} (17) & D^{(m)} = D D \dots D = D^{m}, f^{(m)} = D^{(m)} f, \end{aligned}$ where $\begin{aligned} (18) & h_{j} (x) = \prod_{i = 0, i \neq j}^{N} \frac{x - ξ_{i}}{ξ_{j} - ξ_{i}}, h_{j}^{'} (ξ_{k}) = d_{k j}, \end{aligned}$ $\begin{array}{c} (19) & d_{k j} = \{\begin{array}{ll} \frac{\partial_{x} J_{N + 1}^{α, β} (ξ_{k})}{\partial_{x} J_{N + 1}^{α, β} (ξ_{j})} \frac{1}{ξ_{k} - ξ_{j}} & 0 ⩽ k \neq j ⩽ N, \\ \frac{α - β + (α + β + 2) ξ_{k}}{2 (1 - ξ_{k}^{2})} & 1 ⩽ k = j ⩽ N, \end{array} \end{array}$ and $D = {(d_{k j})}_{k, j = 0}^{N}$ is called the first-order differentiation matrix in the physical space.

Remark 1.
Notice that all the nodes ${ξ_{k}}_{k = 0}^{N}$ of Gauss formula lie in the interior of the interval $(- 1, 1)$ . Below, in order to impose potential boundary conditions in partial differential equations, we consider the Gauss–Lobatto quadrature which include both endpoints as nodes, and the similar conclusions for Gauss-Radau quadrature can be derived straightforwardly through the same idea.

3. Jacobi spectral method for Caputo fractional derivatives

Without loss of generality, we just focus on computations of the Caputo fractional derivative. The left-side Caputo fractional derivative with order $μ \in (0, 1)$ , $t \in [0, T]$ is defined by [19] $\begin{aligned} (20) & _{0}^{C} D_{t}^{μ} f (t) = \frac{1}{Γ (1 - μ)} \int_{0}^{t} \frac{1}{{(t - τ)}^{μ}} \frac{d f (τ)}{d τ} d τ . \end{aligned}$

In view of the weakly singular kernel, we consider the spectral approximation based on the Jacobi–Gauss–Lobatto points ${t_{j}}_{j = 0}^{M}$ to construct the Jacobi-collocation method.

Let ${ξ_{k}, ω_{k}}_{k = 0}^{N}$ be a set of the Jacobi–Gauss–Lobatto points and weights with $α = - μ$ , $β = 0$ , and the linear transform $\begin{aligned} (21) & τ = \frac{t (ξ + 1)}{2}, \end{aligned}$ maps points ${ξ_{k}}_{k = 0}^{N} \in [- 1, 1]$ into ${τ_{k}}_{k = 0}^{N} \in [0, t]$ , then we have $\begin{aligned} _{0}^{C} D_{t}^{μ} f (t) & = \frac{1}{Γ (1 - μ)} \int_{0}^{t} \frac{f^{'} (τ)}{{(t - τ)}^{μ}} d τ \\ = \frac{1}{Γ (1 - μ)} {(\frac{t}{2})}^{1 - μ} \int_{- 1}^{1} {(1 - ξ)}^{- μ} \partial_{τ} f (τ (t, ξ)) d ξ \\ (22) & = \frac{1}{Γ (1 - μ)} {(\frac{t}{2})}^{- μ} \int_{- 1}^{1} {(1 - ξ)}^{- μ} \partial_{ξ} f (τ (t, ξ)) d ξ, \end{aligned}$ and the corresponding Jacobi-collocation method for (22) is: $\begin{array}{c} (23) & \{\begin{array}{l} Find f_{N} \in P_{N} such that for 0 ⩽ j ⩽ M, \\ _{0}^{C} D_{t}^{μ} f_{N} (t_{j}) = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} \sum_{k = 0}^{N} ω_{k} \partial_{ξ} f_{N} (τ (t_{j}, ξ_{k})) . \end{array} \end{array}$

In the light of the first-order differentiation matrix $\begin{aligned} (24) & \partial_{ξ} f_{N} (ξ_{k}) = \sum_{j = 0}^{N} d_{k j} f_{N} (ξ_{j}), \end{aligned}$ the scheme (23) becomes $\begin{aligned} (25) & _{0}^{C} D_{t}^{μ} f_{N} (t_{j}) = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} \sum_{k = 0}^{N} ω_{k} \sum_{i = 0}^{N} d_{k i} f_{N} (τ (t_{j}, ξ_{i})), \end{aligned}$ or the matrix form $\begin{aligned} (26) & _{0}^{C} D_{t}^{μ} f_{N} (t_{j}) = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} W D f_{j}, \end{aligned}$ where $\begin{aligned} (27) & W = {[\begin{array}{cccc} ω_{0}, & ω_{1}, & \dots, & ω_{N} \end{array}]}_{(1, N + 1)}, f_{j} = {[\begin{array}{c} f_{N} (τ (t_{j}, ξ_{0})) \\ f_{N} (τ (t_{j}, ξ_{1})) \\ \dots \\ f_{N} (τ (t_{j}, ξ_{N})) \end{array}]}_{(N + 1, 1)}, \end{aligned}$ and $D = {(d_{k j})}_{k, j = 0}^{N}$ is the first-order differentiation matrix defined by (19).

Remark 2.
The numerical scheme (26) for left-side Caputo fractional derivative can be easily extend to calculating left-side Riemann–Liouville fractional derivative by taking advantage of the well-known relation between Caputo and Riemann–Liouville derivative, $\forall μ \in (0, 1)$ : $\begin{aligned} (28) & _{0}^{C} D_{t}^{μ} f (t) = \frac{f (0)}{t^{μ} Γ (1 - μ)} +_{0}^{RL} D_{t}^{μ} f (t), \end{aligned}$ where $_{0}^{RL} D_{t}^{μ} f (t) = \frac{1}{Γ (1 - μ)} \frac{d}{d t} \int_{0}^{t} \frac{f (τ)}{{(t - τ)}^{μ}} d τ$ is left-side Riemann–Liouville fractional derivative.

Next, we will analyze the error of scheme (25), which approximate the Caputo fractional derivative (20) in anisotropic Jacobi-weighted Sobolev spaces.
4. Convergence analysis

For simplicity, we assume that the points ${t_{j}}$ and ${ξ_{j}}$ are the same (Jacobi–Gauss–Lobatto), although they can be chosen differently in type and in number. For clarity of presentation, we introduce the non-uniformly Jacobi-weighted Sobolev space: $\begin{aligned} (29) & B_{α, β}^{m} (I) : = {u : \partial_{x}^{k} u \in L_{ω^{α + k, β + k}}^{2} (I), 0 ⩽ k ⩽ m}, m \in N, \end{aligned}$ equipped with the inner product, norm and semi-norm $\begin{aligned} {(u, v)}_{B_{α, β}^{m}} = \sum_{k = 0}^{m} {(\partial_{x}^{k} u, \partial_{x}^{k} v)}_{ω^{α + k, β + k}}, \\ (30) & ‖ u ‖_{B_{α, β}^{m}} = {(u, v)}_{B_{α, β}^{m}}^{1 / 2}, | u |_{B_{α, β}^{m}} = {‖ \partial_{x}^{m} u ‖}_{ω^{α + m, β + m}} . \end{aligned}$ Let ${(\cdot, \cdot)}_{ω^{- μ, 0}}$ and ${⟨ \cdot, \cdot ⟩}_{N, ω^{- μ, 0}}$ be the weighted continuous and discrete inner product, respectively, and $I_{N}^{- μ, 0}$ be the corresponding interpolation operator. Then we introduce some lemmas as follows.

Lemma 5 ([20, Lemma 5.1]).

For any non-negative integer r and $0 < κ < 1$ , there exists a linear transform $T_{N} : C^{r, κ} (I) \to P_{N}$ such that $\begin{aligned} (31) & ‖ v - T_{N} v ‖_{\infty} ⩽ c_{r, κ} N^{- (r + κ)} ‖ v ‖_{C^{r, κ}}, \forall v \in C^{r, κ} (I), \end{aligned}$ where $c_{r, κ}$ is a positive constant, $C^{r, κ} (I)$ is the space of functions whose r-th derivatives are Hölder continuous with exponent κ, endowed with the usual norm $\begin{array}{r} ‖ v ‖_{C^{r, κ}} = max_{0 ⩽ l ⩽ r} max_{x \in I} | \partial_{x}^{l} v (x) | + max_{0 ⩽ l ⩽ r} sup_{x \neq y} \frac{| \partial_{x}^{l} v (x) - \partial_{x}^{l} v (y) |}{| x - y |^{κ}} . \end{array}$

Lemma 6 ([16, Theorem 3.1]).

Let ${h_{i}}_{i = 0}^{N}$ be the Lagrange basis polynomials associated with the Jacobi–Gauss–Lobatto interpolations with the parameter pair ${- μ, 0}$ . Then, for $- 1 / 2 ⩽ μ < 3 / 2$ , we have $\begin{aligned} (32) & Λ_{N} : = max_{| x | ⩽ 1} \sum_{i = 0}^{N} | h_{i} (x) | \sim ln N . \end{aligned}$

Next, we reveal the error estimate between the exact value and Jacobi collocation approximation for Caputo fractional derivative as Theorem 1.

Theorem 1.
Assume $f \in L^{\infty} (I) \cap B_{- 1, - 1}^{m} (I)$ with $f^{(s)} (0) = 0$ , $s = 1, 2, \dots, m$ and $\begin{aligned} (33) & f_{m}^{} : = max_{0 ⩽ j ⩽ N} {(\int_{0}^{t_{j}} {| \partial_{τ}^{m + 1} f (τ) |}^{2} {(t_{j} - τ)}^{m - 1 - μ} τ^{m - 1} d τ)}^{1 / 2} < \infty, \end{aligned}$ for certain integer $1 ⩽ m ⩽ N + 1$ . Then we have the estimate: $\begin{aligned} (34) & {‖_{0}^{C} D_{t}^{μ} f -_{0}^{C} D_{t}^{μ} f_{N} ‖}_{\infty} ⩽ c f_{m}^{} N^{- m + 1} ln N, \end{aligned}$ where c is a positive constant independent of N, m and f.
Proof.
We rewrite (22) and (25) respectively as follows $\begin{aligned} (35) & _{0}^{C} D_{t}^{μ} f (t_{j}) = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} {(\partial_{ξ} f (τ (t_{j}, \cdot)), 1)}_{ω^{- μ, 0}}, \\ (36) & _{0}^{C} D_{t}^{μ} f_{N} (t_{j}) = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} {⟨ \partial_{ξ} f_{N} (τ (t_{j}, \cdot)), 1 ⟩}_{N, ω^{- μ, 0}} . \end{aligned}$ Denote $e =_{0}^{C} D_{t}^{μ} f -_{0}^{C} D_{t}^{μ} f_{N}$ , then we have $\begin{aligned} e (t_{j}) & = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} {{(\partial_{ξ} f (τ (t_{j}, \cdot)), 1)}_{ω^{- μ, 0}} - {⟨ \partial_{ξ} f_{N} (τ (t_{j}, \cdot)), 1 ⟩}_{N, ω^{- μ, 0}}} \\ (37) & = \frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} {(\partial_{ξ} (f - f_{N}) (τ (t_{j}, \cdot)), 1)}_{ω^{- μ, 0}} + G (t_{j}), \end{aligned}$ where $\begin{aligned} (38) & G (t) = \frac{1}{Γ (1 - μ)} {(\frac{t}{2})}^{- μ} {{(\partial_{ξ} f_{N} (τ (t, \cdot)), 1)}_{ω^{- μ, 0}} - {⟨ \partial_{ξ} f_{N} (τ (t, \cdot)), 1 ⟩}_{N, ω^{- μ, 0}}}, \end{aligned}$ equivalently, we write (37) as $\begin{aligned} (39) & I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} f) -_{0}^{C} D_{t}^{μ} f_{N} = I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} (f - f_{N})) + I_{N}^{- μ, 0} G . \end{aligned}$ Consequently, $\begin{aligned} (40) & e =_{0}^{C} D_{t}^{μ} (f - f_{N}) + G_{1} + G_{2} + I_{N}^{- μ, 0} G, \end{aligned}$ where $\begin{aligned} G_{1} =_{0}^{C} D_{t}^{μ} f - I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} f), \\ G_{2} = I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} (f - f_{N})) -_{0}^{C} D_{t}^{μ} (f - f_{N}) . \end{aligned}$ Now, we estimate the four terms on the right hand side of (40). Firstly, by using [20, Theorem 3.43] and Cauchy–Schwarz inequality, we have $\begin{aligned} {‖_{0}^{C} D_{t}^{μ} (f - f_{N}) ‖}_{\infty} & = {‖ \frac{1}{Γ (1 - μ)} \int_{0}^{t} {(t - τ)}^{- μ} \partial_{τ} (f - f_{N}) (τ) d τ ‖}_{\infty} \\ = {‖ \frac{1}{Γ (1 - μ)} {(\frac{t}{2})}^{- μ} \int_{- 1}^{1} {(1 - ξ)}^{- μ} \partial_{ξ} (f - f_{N}) (τ (t, ξ)) d ξ ‖}_{\infty} \\ ⩽ | \frac{1}{Γ (1 - μ)} {(\frac{t}{2})}^{- μ} {(\frac{2^{1 - μ}}{1 - μ})}^{1 / 2} | {‖ \partial_{ξ} (f - f_{N}) (τ (t, \cdot)) ‖}_{ω^{- μ, 0}} \\ (41) & ⩽ c_{1} N^{- m + 1} {‖ \partial_{ξ}^{m} f ‖}_{ω^{m - μ, m}} . \end{aligned}$ To estimate $G_{1}$ , let $I_{N}$ be the Legendre-Gauss–Lobatto polynomial interpolation operator. The Sobolev inequality gives $\begin{aligned} ‖ G_{1} ‖_{\infty} & = {‖_{0}^{C} D_{t}^{μ} f - I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} f) ‖}_{\infty} \\ = {‖_{0}^{C} D_{t}^{μ} f - I_{N} (_{0}^{C} D_{t}^{μ} f) + I_{N}^{- μ, 0} (I_{N} (_{0}^{C} D_{t}^{μ} f) -_{0}^{C} D_{t}^{μ} f) ‖}_{\infty} \\ ⩽ (1 + c_{2} ln N) {‖_{0}^{C} D_{t}^{μ} f - I_{N} (_{0}^{C} D_{t}^{μ} f) ‖}_{\infty} \\ ⩽ c_{3} ln N {‖_{0}^{C} D_{t}^{μ} f - I_{N} (_{0}^{C} D_{t}^{μ} f) ‖}^{1 / 2} {‖_{0}^{C} D_{t}^{μ} f - I_{N} (_{0}^{C} D_{t}^{μ} f) ‖}_{1}^{1 / 2} \\ ⩽ c_{4} \sqrt{\frac{(N - m + 1)!}{N!}} {(N + m)}^{- \frac{m}{2}} ln N {‖_{0}^{C} D_{t}^{m + μ} f ‖}_{\infty} \\ (42) & ⩽ c_{5} N^{- m + \frac{1}{2}} ln N ‖ \partial_{ξ} f ‖_{ω^{m - μ, 0}} . \end{aligned}$ From Lemmas 5-6 and [20, Lemma 5.2], we obtain $\begin{aligned} ‖ G_{2} ‖_{\infty} & = {‖ I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} (f - f_{N})) -_{0}^{C} D_{t}^{μ} (f - f_{N}) ‖}_{\infty} \\ ⩽ {‖ I_{N}^{- μ, 0} (_{0}^{C} D_{t}^{μ} (f - f_{N})) - T_{N} (_{0}^{C} D_{t}^{μ} (f - f_{N})) ‖}_{\infty} \\ + {‖ T_{N} (_{0}^{C} D_{t}^{μ} (f - f_{N})) -_{0}^{C} D_{t}^{μ} (f - f_{N}) ‖}_{\infty} \\ ⩽ (1 + c_{6} ln N) {‖ T_{N} (_{0}^{C} D_{t}^{μ} (f - f_{N})) -_{0}^{C} D_{t}^{μ} (f - f_{N}) ‖}_{\infty} \\ (43) & ⩽ c_{7} N^{- κ} ln N {‖_{0}^{C} D_{t}^{μ} (f - f_{N}) ‖}_{C^{0, κ}} ⩽ c_{8} N^{- κ - m + 1} ln N ‖ f ‖_{ω^{- μ, 0}} . \end{aligned}$ By Lemma 6, we yield $\begin{aligned} {‖ I_{N}^{- μ, 0} G ‖}_{\infty} & ⩽ max_{0 ⩽ j ⩽ N} | G (t_{j}) | \sum_{j = 0}^{N} | h_{j} (x) | ⩽ c_{9} ln N max_{0 ⩽ j ⩽ N} | G (t_{j}) | \\ ⩽ c_{10} \sqrt{\frac{(N - m + 1)!}{N!}} {(N + m)}^{- \frac{m + 1}{2}} ln N \\ (44) & max_{0 ⩽ j ⩽ N} {\frac{1}{Γ (1 - μ)} {(\frac{t_{j}}{2})}^{- μ} {‖ \partial_{ξ}^{m + 1} f ‖}_{ω^{m - 1 - μ, m - 1}}} . \end{aligned}$ The definition shows that $\begin{aligned} (45) & {‖ \partial_{ξ}^{m + 1} f (τ (t_{j}, \cdot)) ‖}_{ω^{m - 1 - μ, m - 1}} = {(\frac{t_{j}}{2})}^{\frac{3 + μ}{2}} {(\int_{0}^{t_{j}} {| \partial_{τ}^{m + 1} f (τ) |}^{2} {(t_{j} - τ)}^{m - 1 - μ} τ^{m - 1} d τ)}^{1 / 2} . \end{aligned}$ Therefore, we have $\begin{aligned} (46) & {‖ I_{N}^{- μ, 0} G ‖}_{\infty} & ⩽ c_{11} f_{m}^{*} N^{- m} ln N, \end{aligned}$ where the above $c_{k} (1 ⩽ k ⩽ 11)$ is a positive constant independent of N, m and f.

Finally, a combination of (40), (41), (42), (43) and (46) leads to the estimate (34). □

5. Mollification regularization and its convergence analysis

Thereafter, we consider the mollification regularization of the numerical scheme (26) when function f is contaminated by noise. First, we introduce the mollification method [17].

Assume that $g (t) \in L^{1} [0, 1]$ and $t \in [p ε, 1 - p ε]$ , the ε-mollification of $g (t)$ is defined by the convolution $\begin{aligned} (47) & J_{ε} g (t) = (ρ_{ε} * g) (t) = \int_{t - p ε}^{t + p ε} ρ_{ε} (t - s) g (s) d s, p > 0, ε > 0, \end{aligned}$ where the nonnegative kernel satisfies $\begin{matrix} ρ_{ε} (t) = \{\begin{array}{ll} A_{p} ε^{- 1} exp (- \frac{t^{2}}{ε^{2}}), & | t | ⩽ p ε, \\ 0, & | t | > p ε, \end{array} \int_{- p ε}^{p ε} ρ_{ε} (t) d t = 1, A_{p} = \frac{1}{\int_{- p}^{p} exp (- s^{2}) d s}, \end{matrix}$ $g^{δ} (t)$ denotes the noisy data with noise level δ of exact data function $g (t)$ (i.e. $‖ g (t) - g^{δ} (t) ‖_{\infty} ⩽ δ$ ), ε is the radius of mollification, and the mollified solution in the sense of convolution implies that $\begin{aligned} (48) & J_{ε} (_{0}^{C} D_{t}^{μ} g^{δ}) (t) =_{0}^{C} D_{t}^{μ} (J_{ε} g^{δ}) (t) . \end{aligned}$

Next, we express two Lemmas as follows and finally obtain a convergence Theorem about the mollification regularization.

Lemma 7 ([17, Theorem 2.1 and 4.1]).

Let ε is the radius of mollification, $v^{δ}$ is the measurement data with noise level δ of the exact v (i.e. $‖ v - v^{δ} ‖_{\infty} ⩽ δ$ ), if $v^{'}$ and $v^{δ}$ are uniformly Lipschitz on the interval $[0, 1]$ and then there exists constants $c_{12}$ , $c_{13}$ , independent of ε, such that $\begin{array}{r} {‖ v - J_{ε} v^{δ} ‖}_{\infty} ⩽ c_{12} ε + δ, \end{array}$ and $\begin{array}{r} {‖_{0}^{C} D_{t}^{μ} v -_{0}^{C} D_{t}^{μ} (J_{ε} v^{δ}) ‖}_{\infty} ⩽ \frac{c_{13}}{Γ (2 - μ)} (ε + \frac{δ}{ε}) . \end{array}$

Lemma 8.
Assume that v satisfy the prerequisites of both Lemma 5 and Lemma 7 at the same time, we have $\begin{aligned} (49) & {‖ v - T_{N} (J_{ε} v^{δ}) ‖}_{\infty} ⩽ (c_{12} ε + δ) + c_{r, κ} N^{- (r + κ)} {‖ J_{ε} v^{δ} ‖}_{C^{r, κ}}, \end{aligned}$ where $c_{r, κ}$ and $c_{12}$ are constants stated in Lemma 5 and Lemma 7 , respectively.
Proof.
From Lemma 5, Lemma 7 and the triangle inequality, the conclusion of Lemma 8 is easily to be established. More precisely, $\begin{aligned} {‖ v - T_{N} (J_{ε} v^{δ}) ‖}_{\infty} & = {‖ v - J_{ε} v^{δ} + J_{ε} v^{δ} - T_{N} (J_{ε} v^{δ}) ‖}_{\infty} \\ ⩽ {‖ v - J_{ε} v^{δ} ‖}_{\infty} + {‖ J_{ε} v^{δ} - T_{N} (J_{ε} v^{δ}) ‖}_{\infty} \\ ⩽ (c_{12} ε + δ) + c_{r, κ} N^{- (r + κ)} {‖ J_{ε} v^{δ} ‖}_{C^{r, κ}} . \end{aligned}$ □
Theorem 2.
Assume $f \in L^{\infty} (I) \cap B_{- 1, - 1}^{m} (I)$ with $f^{(s)} (0) = 0$ , $s = 1, 2, \dots, m$ and $\begin{aligned} (50) & {(J_{ε} f^{δ})}_{m}^{} : = max_{0 ⩽ j ⩽ N} {(\int_{0}^{t_{j}} {| \partial_{τ}^{m + 1} (J_{ε} f^{δ}) (τ) |}^{2} {(t_{j} - τ)}^{m - 1 - μ} τ^{m - 1} d τ)}^{1 / 2} < \infty, \end{aligned}$ for certain integer $1 ⩽ m ⩽ N + 1$ . Then we have the estimate: $\begin{aligned} (51) & {‖_{0}^{C} D_{t}^{μ} f -_{0}^{C} D_{t}^{μ} {(J_{ε} f^{δ})}_{N} ‖}_{\infty} ⩽ \frac{c_{13}}{Γ (2 - μ)} (ε + \frac{δ}{ε}) + c {(J_{ε} f^{δ})}_{m}^{} N^{- m + 1} ln N, \end{aligned}$ where c and $c_{13}$ are constants stated in Lemma 1 and Lemma 7 , respectively. In other words, the convergence rate could reach $O (δ^{1 / 2})$ when $ε = δ^{1 / 2}$ and $N = δ^{\frac{1}{4 - 2 m}}$ for $m ⩾ 3$ .
Proof.
By virtue of Theorem 1 and Lemma 7, it reveals $\begin{aligned} {‖_{0}^{C} D_{t}^{μ} f -_{0}^{C} D_{t}^{μ} {(J_{ε} f^{δ})}_{N} ‖}_{\infty} & = {‖_{0}^{C} D_{t}^{μ} f -_{0}^{C} D_{t}^{μ} (J_{ε} f^{δ}) +_{0}^{C} D_{t}^{μ} (J_{ε} f^{δ}) -_{0}^{C} D_{t}^{μ} {(J_{ε} f^{δ})}_{N} ‖}_{\infty} \\ ⩽ {‖_{0}^{C} D_{t}^{μ} f -_{0}^{C} D_{t}^{μ} (J_{ε} f^{δ}) ‖}_{\infty} + {‖_{0}^{C} D_{t}^{μ} (J_{ε} f^{δ}) -_{0}^{C} D_{t}^{μ} {(J_{ε} f^{δ})}_{N} ‖}_{\infty} \\ ⩽ \frac{c_{13}}{Γ (2 - μ)} (ε + \frac{δ}{ε}) + c {(J_{ε} f^{δ})}_{m}^{*} N^{- m + 1} ln N . \end{aligned}$ □

6. Numerical examples for fractional derivative

To test the computational effect of the numerical scheme (26), we take the following two functions for examples.

Example 1.
Take the function $f_{1} (t) = t^{6}$ , $t \in [0, T]$ , then the exact Caputo fractional derivative is $\begin{aligned} (52) & _{0}^{C} D_{t}^{μ} f_{1} (t) = \frac{1}{Γ (1 - μ)} \int_{0}^{t} \frac{6 τ^{5}}{{(t - τ)}^{μ}} d τ = \frac{720 t^{6 - μ}}{Γ (7 - μ)} . \end{aligned}$
Example 2.
Take the function $f_{2} (t) = e^{2 t} - 2 t - 2 t^{2} - \frac{4}{3} t^{3} - \frac{2}{3} t^{4}$ , $t \in [0, T]$ , then the exact Caputo fractional derivative reads $\begin{aligned} _{0}^{C} D_{t}^{μ} f_{2} (t) & = \frac{1}{Γ (1 - μ)} \int_{0}^{t} \frac{f_{2}^{'} (τ)}{{(t - τ)}^{μ}} d τ \\ = \sum_{k = 0}^{\infty} \frac{t^{- μ} {(2 t)}^{k + 1}}{Γ (k + 2 - μ)} - 2 (\frac{t^{1 - μ}}{Γ (2 - μ)} + \frac{2 t^{2 - μ}}{Γ (3 - μ)} + \frac{4 t^{3 - μ}}{Γ (4 - μ)} + \frac{8 t^{4 - μ}}{Γ (5 - μ)}) \\ (53) & = 2 t^{1 - μ} E_{1, 2 - μ} (2 t) - \sum_{k = 1}^{4} \frac{2^{k} t^{k - μ}}{Γ (k + 1 - μ)}, \end{aligned}$ where $E_{η, ζ} (z)$ is the Mittag-Leffler function $\begin{aligned} (54) & E_{η, ζ} (z) = \sum_{k = 0}^{\infty} \frac{z^{k}}{Γ (η k + ζ)}, η > 0, ζ \in R . \end{aligned}$
Remark 3.
In the calculation of above two examples, the following formulas would to be involved. $\begin{aligned} (55) & e^{t} = \sum_{k = 0}^{\infty} \frac{t^{k}}{k!}, \frac{1}{Γ (1 - μ)} \int_{0}^{t} \frac{τ^{k}}{{(t - τ)}^{μ}} d τ = \frac{k! t^{k + 1 - μ}}{Γ (k + 2 - μ)} . \end{aligned}$

Next, we will compare the computational fractional derivative approximated by numerical scheme (25) from the measurement $f_{1}^{δ} (t)$ , $f_{2}^{δ} (t)$ of exact value (52) and (53) with the noise level δ. Take $T = 1$ , $M = 10$ , $N = 16$ , and $L^{1} errors$ as $\begin{aligned} (56) & L^{1} errors = \frac{1}{M + 1} \sum_{i = 0}^{M} |_{0}^{C} D_{t}^{μ} f (t_{i}) -_{0}^{C} D_{t}^{μ} f_{N} (t_{i}) | . \end{aligned}$

Figure 1 illustrate the numerical errors of scheme (25) for $f_{1} (t)$ and $f_{2} (t)$ with various order μ (from 0.01 to 0.99 with step 0.01) when $δ = 0$ .

It can be seen from Fig. 1 that the numerical scheme (25) of Caputo fractional derivative is of high resolution for data without any noise.

Fig. 1.
Numerical error under $μ \in [0.01, 0.99]$ with $δ = 0$ . Solid line is for $f_{1} (t)$ and the dashed line for $f_{2} (t)$ .

Fig. 2.
Numerical error under $μ \in [0.01, 0.99]$ with various noise level $δ = 0.01, 0.05$ . Solid line is for $f_{1} (t)$ and the dashed line for $f_{2} (t)$ .

Fig. 3.
Numerical results for various order μ with noise level $δ = 0$ . The first row is for Example 1 and the second row for Example 2.

Fig. 4.
Numerical results for various order μ with noise level $δ = 0.01$ without mollification. The first row is for Example 1 and the second row for Example 2.

Fig. 5.
Numerical results for various order μ with noise level $δ = 0.01$ with mollification. The first row is for Example 1 and the second row for Example 2.

Fig. 6.
Numerical results for various order μ with noise level $δ = 0.05$ without mollification. The first row is for Example 1 and the second row for Example 2.

Figures 1 and 3 show that the scheme (26) performs excellent to calculate the fractional derivative numerically for various order $μ \in (0, 1)$ under the data without any noise. Figures 2 exhibit that the effect of numerical scheme with mollification regularization under various fractional order (from 0.01 to 0.99) is significantly better than the case without mollification when data contains noise. Figures 4 and 5 illustrate a comparison of numerical fractional derivative for diverse order $μ \in (0, 1)$ without- and with- mollification method under the noise level $δ = 0.01$ , and the corresponding results for noise level $δ = 0.05$ are elaborated in Fig. 6 and Fig. 7. The above figures state the truth that the numerical algorithm performs overwhelmingly well as data is exactly, and the mollification regularization method helps the numerical scheme achieve an effective solution when data contains noise. The experimental results are in good agreement with the theoretical results.
7. Conclusion

We propose a numerical scheme to approximate Caputo fractional derivative with order $μ \in (0, 1)$ by using Jacobi–Gauss–Lobatto collocation method on the sense of Gaussian quadrature, and the convergence analysis is presented as well. Some numerical experiments well verified the theoretical results, and confirmed the efficiency to approximate fractional derivative by invoking mollification as function data comprises noise. The proposed method provides a versatile and effective approach to calculate the fractional derivative. The same technique can be extend to other quadrature rules which include one endpoint or all interior points.

Footnotes

Acknowledgements

The authors would like to thank editors and reviewers for their valuable suggestions and comments. This work is partially supported by National Natural Science Foundation of China (11861007, 11961002, 12061008, 12261004), Natural Science Foundation of Jiangxi Province of China (20232BAB201019).

Conflict of interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

References

Canuto,

M.Y.

Hussaini,

Quarteroni and

T.A.

Zang, Spectral Methods. Fundamentals in Single Domains, Springer-Verlag, Berlin, 2006.

Carpinteri and

Mainardi, Fractals and Fractional Calculus in Continuum Mechanics, Springer-Verlag, Vienna, 1998.

Chen,

Shen and

L.L.

Wang, Generalized Jacobi functions and their applications to fractional differential equations, Mathematics of Computation 85(300) (2016), 1603–1638. doi:10.1090/mcom3035.

Chen,

Shen,

Z.M.

Zhang and

Zhou, A spectrally accurate approximation to subdiffusion equations using the log orthogonal functions, SIAM Journal on Scientific Computing 42(2) (2020), A849–A877. doi:10.1137/19M1281927.

Chester, Resonant oscillations in closed tubes, Journal of Fluid Mechanics 18(1) (1964), 44–64. doi:10.1017/S0022112064000040.

Dehghan and

A.B.

Mingarelli, Fractional Sturm–Liouville eigenvalue problems (I), RACSAM 114 (2020), 46. doi:10.1007/s13398-019-00756-8.

Esmaeili and

G.V.

Milovanović, Nonstandard Gauss–Lobatto quadrature approximation to fractional derivatives, Fractional Calculus and Applied Analysis 17(4) (2014), 1075–1099. doi:10.2478/s13540-014-0215-z.

B.Y.

Guo and

L.L.

Wang, Jacobi interpolation approximations and their applications to singular differential equations, Advances in Computational Mathematics 14 (2001), 227–276. doi:10.1023/A:1016681018268.

B.Y.

Guo and

L.L.

Wang, Jacobi approximations in non-uniformly Jacobi-weighted Sobolev spaces, Journal of Approximation Theory 128(1) (2004), 1–41. doi:10.1016/j.jat.2004.03.008.

10.

Hesthaven,

Gottlieb and

Gottlieb, Spectral Methods for Time-Dependent Problems, Cambridge University Press, Cambridge, 2007.

11.

B.T.

Jin and

Zhou, Recovery of a space-time-dependent diffusion coefficient in subdiffusion: Stability, approximation and error analysis, IMA Journal of Numerical Analysis drac051 (2022).

12.

J.J.

Keller, Propagation of simple non-linear waves in gas filled tubes with friction, Z. angew. Math. Phys. 32 (1981), 170–181. doi:10.1007/BF00946746.

13.

Y.M.

Lin and

C.J.

Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, Journal of Computational Physics 225(2) (2007), 1533–1552. doi:10.1016/j.jcp.2007.02.001.

14.

R.L.

Magin, Fractional Calculus in Bioengineering, Begell House Inc., Redding, 2006.

15.

Z.P.

Mao and

G.E.

Karniadakis, A spectral method (of exponential convergence) for singular solutions of the diffusion equation with general two-sided fractional derivative, SIAM Journal on Numerical Analysis 56 (2018), 24–49. doi:10.1137/16M1103622.

16.

Mastroianni and

Occorsio, Optimal systems of nodes for Lagrange interpolation on bounded intervals. A survey, Journal of Computational and Applied Mathematics 134 (2001), 325–341. doi:10.1016/S0377-0427(00)00557-4.

17.

D.A.

Murio, On the stable numerical evaluation of Caputo fractional derivatives, Computers & Mathematics with Applications 51(9–10) (2006), 1539–1550. doi:10.1016/j.camwa.2005.11.037.

18.

D.X.

Nie,

Sun and

W.H.

Deng, Sharp error estimates for spatial-temporal finite difference approximations to fractional sub-diffusion equation without regularity assumption on the exact solution, Fractional Calculus and Applied Analysis 26 (2023), 1421–1464. doi:10.1007/s13540-023-00162-3.

19.

Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999.

20.

Shen,

Tang and

L.L.

Wang, Spectral Methods: Algorithms, Analysis and Applications, Springer-Verlag, Berlin Heidelberg, 2011.

21.

Sugimoto and

Kakutani, Generalized Burgers’ equation for nonlinear viscoelastic waves, Wave Motion 7(5) (1985), 447–458. doi:10.1016/0165-2125(85)90019-8.

22.

B.J.

West,

Bologna and

Grigolini, Physics of Fractal Operators, Springer-Verlag, New York, 2003.

23.

Zayernouri and

G.E.

Karniadakis, Fractional Sturm–Liouville eigen-problems: Theory and numerical approximation, Journal of Computational Physics 252 (2013), 495–517. doi:10.1016/j.jcp.2013.06.031.

24.

Zayernouri and

G.E.

Karniadakis, Fractional spectral collocation method, SIAM Journal on Scientific Computing 36(1) (2014), A40–A62. doi:10.1137/130933216.

25.

Zayernouri and

G.E.

Karniadakis, Fractional spectral collocation methods for linear and nonlinear variable order FPDEs, Journal of Computational Physics 293 (2015), 312–338. doi:10.1016/j.jcp.2014.12.001.

A Jacobi spectral method for calculating fractional derivative based on mollification regularization

Abstract

Keywords

1. Introduction

2. Preliminaries

Lemma 1 ([20, Theorem 3.16]).

Lemma 2 ([20, Recurrence Formulas]).

Lemma 3 ([20, Theorem 3.25]).

Lemma 4 ([20, Differentiation Matrix in the Physical Space]).

Lemma 5 ([20, Lemma 5.1]).

Lemma 6 ([16, Theorem 3.1]).

Lemma 7 ([17, Theorem 2.1 and 4.1]).

Footnotes

Acknowledgements

Conflict of interest

References