Exponentially accurate spectral and spectral element methods for fractional ODEs

Abstract

Current discretizations of fractional differential equations lead to numerical solutions of low order of accuracy. Here, we present different methods for fractional ODEs that lead to exponentially fast decay of the error. First, we develop a Petrov–Galerkin (PG) spectral method for Fractional Initial-Value Problems (FIVPs) of the form ${}_0\mathcal{D}_t^{\nu }u(t)=f(t)$ and Fractional Final-Value Problems (FFVPs) ${}_t\mathcal{D}_T^{\nu }u(t)=g(t)$, where $\nu \in (0,1)$, subject to Dirichlet initial/final conditions. These schemes are developed based on a new spectral theory for fractional Sturm–Liouville problems (FSLPs), which has been recently developed in [1]. Specifically, we obtain solutions to FIVPs and FFVPs in terms of the new fractional (non-polynomial) basis functions, called Jacobi polyfractonomials, which are the eigenfunctions of the FSLP of first kind (FSLP-I). Correspondingly, we employ another space of test functions as the span of polyfractonomial eigenfunctions of the FSLP of second kind (FSLP-II). Subsequently, we develop a Discontinuous Spectral Method (DSM) of Petrov–Galerkin sense for the aforementioned FIVPs and FFVPs, where the basis functions do not satisfy the initial/final conditions. Finally, we extend the DSM scheme to a Discontinuous Spectral Element Method (DSEM) for efficient longer time-integration and adaptive refinement. In these discontinuous schemes, we employ the asymptotic eigensolutions to FSLP-I & -II, which are of Jacobi polynomial forms, as basis and test functions. Our numerical tests confirm the exponential/algebraic convergence, respectively, in p- and h-refinements, for various test cases with integer- and fractional-order solutions.

Introduction

Fractional differential operators of form ${\mathcal{D}}_t^{\nu }\equiv d^{\nu }/d{t}^{\nu }$, where $\nu \in \mathbb{R}$, appear in many systems in science and engineering such as electrochemical processes [2], porous or fractured media [3], viscoelastic materials [4], [5], bioengineering applications [6]. For instance, it has been found that the transport dynamics in complex and/or disordered systems is governed by non-exponential relaxation patterns and anomalous diffusion [7], [8], [9]. For such non-Markovian processes, a time-fractional diffusion equation, in which the time-derivative emerges as ${\mathcal{D}}_t^{\nu }u(t)$, governs the evolution for the Probability Density Function (PDF). Another interesting example occurring in viscous fluid flows is the cumulative memory effect of the wall-friction through the boundary layer, which gives rise to fractional derivatives in equations of fluid motion [10], [11], [12].

Over the last two decades, the notion of fractional derivative has been extended to many ordinary fractional differential equations (FODEs) such as fractional Cauchy equation, fractional Gauss equations [13], [14] and fractional Sturm–Liouville equation [15], in addition to a variety of fractional partial differential equations (FPDEs) such as fractional Fokker–Planck equation [16], fractional Burgersʼ equation [17], and fractional advection–diffusion equation [18]. In these problems, the corresponding differential operators can be defined based on different but closely related ways. The extension of existing numerical methods for integer-order differential equations ([19], [20], [21], [22], [23] and references therein) to their corresponding fractional differential equations (FDEs) is not trivial. While the approximation of these models is computationally demanding due to their long-range history-dependence, the development of numerical schemes in this area does not have a long history, and has undergone a fast evolution. Depending on how (temporal ${\mathcal{D}}_t^{\nu }$ or spatial ${\mathcal{D}}_x^{\nu }$) fractional derivatives are discretized and according to their order of accuracy, different classes of numerical methods have been developed in the literature.

To our knowledge, Lubich [24], [25] is the pioneer of the idea of discretized fractional calculus within the spirit of finite difference method (FDM). Later, Sanz-Serna [26] adopted the idea of Lubich and presented a temporal semi-discrete algorithm for partial integro-differential equations, which was first order accurate. Sugimoto [17] also employed an FDM for approximating the fractional derivative emerging in Burgersʼ equation. Later on, the paper of Metzler and Klafter [8] opened a new horizon toward FPDEs by introducing a fractional dynamics approach to time-fractional diffusion. Subsequently, Gorenflo et al. [27] adopted a finite difference scheme by which they could generate discrete models of random walk in such anomalous diffusion. Diethelm et al. proposed a predictor–corrector scheme in addition to a fractional Adams method [13], [28]. After that, Langlands and Henry [29] considered the fractional diffusion equation, and analyzed the $L^1$ scheme for the time-fractional derivative. Sun and Wu [30] also constructed a difference scheme with $L^{\infty }$ approximation of time-fractional derivative. In order to develop and analyze higher-order FDM schemes Lin and Xu [31] analyzed an FDM for the discretization of the time-fractional diffusion equation with order ($2-\alpha$). Kumar and Agrawal [32] proposed a block-by-block method for a class of fractional initial-value problems which later Huang et al. [33] proved that it enjoys a rate of convergence of at least 3. Recently, Cao and Xu [34] rigorously developed the scheme to $(3+\alpha )$-th order, $\alpha \in (0,1)$. To the best of knowledge, this is the highest-order and the most recent FDM scheme for discretization of fractional derivatives.

Although implementation of such FDM approaches is relatively easy, their biggest issue is that the accuracy is limited. Moreover, these approaches suffer from heavy cost of computing the long-ranged memory in discretization of the fractional derivatives at each point. In fact, FDM is inherently a local approach whereas fractional derivatives are essentially global (nonlocal) differential operators. This property would suggest that global schemes such as Spectral Methods (SMs) are more appropriate tools for discretizing fractional differential equations.

Unlike the attention put on developing FDM schemes, very little effort has been put on developing rigorous high-order spectral methods. A Fourier SM was utilized by Sugimoto [17] in a fractional Burgersʼ equation, linearized by Taylor expansion, and a spline-based collocation method was employed by Blank [35] for numerical treatment of a class of FODEs. This approach was later employed by Rawashdeh [36] for solving fractional integro-differential equations. In these works, the expected high convergence rate was not observed and no error/stability analysis was carried out. Lin and Xu [31] developed a hybrid scheme for time-fractional diffusion problem, treating the time-fractional derivative using FDM and discretizing the integer-order spatial derivative by a Legendre SM. In such mixed approaches, the error associated with the low-order temporal accuracy can easily dominate the global error, for instance when the time-dependent portion of the exact solution is discontinuous, or if is a monomial of form $t^n$, where n is sufficiently large, or is a smooth function e.g., $\sin (n\pi t)$. The idea of collocation was later adopted by Khader [37], where he proposed a Chebyshev collocation method for the discretization of the space-fractional diffusion equation. More recently, Khader and Hendy [38] developed a Legendre pseudospectral method for fractional-order delay differential equations.

The collocation and pseudospectral schemes for fractional equations are relatively easy to implement but their performance has not been tested thoroughly. For instance, when the exact solution is of polynomial form it is claimed that a fast convergence is observed. However, for other test cases no such exponential-like convergence is achieved. The first fundamental work on spectral methods for FPDEs was done by Li and Xu [39], [40] who developed a space–time SM for time-fractional diffusion equation. To the best of our knowledge, they were the first who achieved exponential convergence in their numerical tests in agreement with their error analysis. However, in this scheme, the corresponding stiffness and mass matrices are dense and gradually become ill-conditioned when the fractional order α tends to small values. Moreover, this approach is not effective e.g., when the forcing term exhibits discontinuity in the time-domain. This, in turn, motivates the use of domain decomposition and Finite Element Methods (FEM) and Spectral Element Methods (SEM) in the context of fractional calculus.

A theoretical framework for the least-square finite element approximation of a fractional-order differential equation was developed by Fix and Roop [41], where optimal error estimates are proven for piecewise linear trial elements. The main hurdle to overcome in FEM is the nonlocal nature of the fractional operator which leads to large dense matrices; even construction of such matrices presents difficulties [42]. There are, however, a number of recent works already employed in this area using FEM to obtain more efficient schemes. McLean and Mustapha [43] developed a piecewise-constant discontinuous Galerkin method for the time discretization of a sub-diffusion equation. Hanert [44] also considered the use of a Chebyshev spectral element method for the numerical solution of the fractional-order transport. Recently, the idea of the least-square FEM [41] was extended to the spectral element method by Carella [45]. Despite the spectral expansion, these schemes are not properly formulated and fail to achieve exponential convergence.

In this paper, we develop exponentially accurate numerical schemes of Petrov–Galerkin type for the FODEs of form ${}_0\mathcal{D}_t^{\nu }u(t)=f(t)$ and ${}_t\mathcal{D}_T^{\nu }u(t)=f(t)$, introduced, respectively, as Fractional Initial-Value Problem (FIVP) and Fractional Final-Value Problem (FFVP) subject to Dirichlet initial/final conditions. To this end, we first develop a Petrov–Galerkin (PG) spectral method whose corresponding stiffness matrix is diagonal. Subsequently, we develop a Discontinuous Spectral Method (DSM) of Petrov–Galerkin sense with exact quadrature rules for the aforementioned FIVPs and FFVPs. This scheme is also extended to a discontinuous spectral element method (DSEM) for efficient longer time-integrations and adaptive refinement. These schemes are developed based on a new spectral theory for fractional Sturm–Liouville problems (FSLPs), which has been recently developed in [1]. In order to test the performance of our schemes, p-refinement and h-refinement tests are performed for a range of test cases, where the exact solution is a monomial $t^n$, $n\in \mathbb{N}$, fractonomial $t^{n+\mu }$, $\mu \in (0,1)$, (see [1]), smooth functions of form $t^p\sin (n\pi t)$, $p\in \mathbb{N}$, fractional functions $t^{n_1+{\mu }_1}\sin (n\pi t^{n_2+{\mu }_2})$, $n_1,n_2\in \mathbb{N}$ and ${\mu }_1,{\mu }_2\in (0,1)$, or any combinations of these functions.