Exponential sum approximation for Mittag-Leffler function and its application to fractional Zener wave equation

Highlights

• Uniform convergence for $t>0$ of exponential sum approximation of Mittag-Leffler.

• Finite propagation speed for internally damped wave equation.

• Hybrid generalized-alpha method for fractional Zener wave equation.

• Stability analysis of the time stepping scheme based on eigenanalysis.

• Error analysis relates approximation error to second-order convergence of algorithm.

Abstract

The Mittag-Leffler function (MLF) is fundamental to solving fractional calculus problems. An exponential sum approximation for the single-parameter MLF with negative input is proposed. Analysis shows that the approximation, which is based on the Gauss-Legendre quadrature, converges uniformly for all non-positive input. The application to modelling of wave propagation in viscoelastic material with the finite element method is also presented. The propagation speed of the solution to the approximated wave equation is proved to have the same upper bound as the original. The discretized scheme, based on the generalised alpha method, involves only a single matrix inverse whose size is the degree of freedom of the geometric model per time step. It is proved that the scheme is unconditionally stable. Furthermore, the solution converges to the true solution at $O(T^2)$, where T is the interval each time step, provided that the approximation error of the MLF is sufficiently small.

Introduction

Fractional derivatives have applications in modelling viscoelastic materials and control theory [1], [2], [3], [4], [5], [6]. Particularly, the Caputo variant is the most common as its initial conditions involve only integer-order derivatives, which do have physical meaning, in contrast to the fractional derivative initial conditions of the Riemann-Liouville [7]. Because fractional derivatives comprise of convolution with a kernel function with a long tail, specifically $t^{-\beta }$, evaluation is non-local and involves storage of and operations on all history (see methods based on Grünwald–Letnikov fractional derivative for example [5], [8]). Approaches to reduce the cost can be categorized by either compressing the function on which the fractional derivative is applied or compressing the convolution kernel of the fractional derivative operator. In the first case, the to-be-convolved function is approximated by basis functions which have closed-form fractional derivatives [9], [10], [11]. For dynamic problems involving convolution of the kernel with the solution function, a fixed set of basis functions may not be able to approximate the solution uniformly as it evolves with time. Therefore, there is no guarantee that the method converges to the correct solution of the dynamic system. The coefficients may also need to be re-evaluated as time grows, and is thus inefficient. In case of compressing the convolution kernel, the exponential sum approximation was proposed in 2002 by different authors independently to solve the impedance boundary problem and fractional derivatives efficiently [12], [13] (although the work of Yuan and Agrawal lacks convergence analysis). In this paper, the method of exponential sum refers to approximating the kernel function in the form

$$f\left(t\right)\sim \sum _{i=1}^{N_{\text{ex}}}{w}_i{\mathrm{e}}^{-s_{i}t},$$

where $w_i$ and $s_i$ are generally complex numbers as the method was generalised to convolution kernels other than the one in fractional differentiation based on a numerical Laplace inverse transform scheme [14]. Because the kernel may not be uniformly approximated by a fixed set of coefficients for the desirable time interval, an adaptive algorithm was proposed [15].

In the case of the fractional derivative kernel $t^{-\beta }$, the convergence of the exponential sum approximation, with various formulations to the approximation coefficients, has been extensively studied for $\beta \in (0,2)$, which corresponds to the range of practical values in typical physical problems. Direct application of Gaussian quadrature to the Laplace transform identity

$$\mathrm{\Gamma }(1-\alpha )t^{\alpha -1}=\mathcal{L}\left\{s^{-\alpha }\right\}$$

has been analysed in [16], [17], [18]. The multipole expansion method, which is based on the following Laplace transform identity

$$\mathcal{L}\left\{{\mathrm{e}}^{-\lambda t}\right\}=\frac{1}{s+\lambda },$$

has been studied in [19], [20]. Last but not least, double exponential substitution with the trapezoid rule has been explored in [21], [22]. Since the fractional derivative kernel can be scaled simply to $(0,1]$ for any $t>0$, these studies approximate the kernel in the range $[\delta ,1]$ with $0<\delta <1$. Moreover, because the Laplace transform of this kernel has the same form as the kernel itself, the resulting approximation coefficients are all real, whereas, in the general case, the coefficients are complex by applying numerical Laplace inverse transform. δ is typically chosen as the smallest step size required to achieve a certain accuracy in the application of a time-stepping scheme. Due to the finite domain of approximation and finite step size, the Prony's method can be applied to optimise the number of coefficients required to meet an error threshold [21].

The generalized three-parameter Mittag-Leffler function (MLF) is defined in the form of series expansion as [23]

$$E_{\alpha ,\beta }^{\gamma }\left(z\right)=\frac{1}{\mathrm{\Gamma }\left(\gamma \right)}\sum _{n=0}^{\infty }\frac{\mathrm{\Gamma }(\gamma +n)z^n}{\mathit{\text{n}}\text{!}\mathrm{\Gamma }(\alpha n+\beta )},$$

for $\alpha ,\beta ,\gamma \in \mathbb{C}$, with $\Re \left(\alpha \right)>0$. The two-parameter and single-parameter MLFs are the special cases of (1.4) when $\gamma =1$, and $\gamma ,\beta =1$ respectively. It can be shown that, for $\alpha >0$, the single-parameter MLF is the eigenfunction of the Caputo's fractional differential operator, that is [24],

$$D_t^{\alpha }{E}_{\alpha }\left(\lambda t^{\alpha }\right)=\lambda E_{\alpha }\left(\lambda t^{\alpha }\right),\lambda \in \mathbb{R},$$

where $D_t^{\alpha }$, the differential operator, is defined as [7]

$$D_t^{\alpha }f\left(t\right)\triangleq \frac{1}{\mathrm{\Gamma }(m-\alpha )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{m-\alpha -1}\frac{{\mathrm{d}}^{m}f\left(\tau \right)}{\mathrm{d}{\tau }^m}\mathrm{d}\tau ,$$

with $m=\lceil \alpha \rceil$. As a result, the single-parameter MLF can also be applied to solve fractional differential equations. In particular, for viscoelastic problems, the MLF of the form

$$E_{\alpha }(-t^{\alpha })\triangleq \sum _{n=0}^{\infty }\frac{{(-1)}^n{t}^{n\alpha }}{\mathrm{\Gamma }(1+\alpha n)},\alpha \in (0,1),t\ge 0$$

is involved as a convolution kernel for $\alpha \in (0,1)$. For example, the adaptive fast convolution algorithm in [15] with the Laplace transform of (1.7) has been applied to model pressure on artery wall as a result of blood flow in [5] using the finite element method. However, the use of the adaptive scheme results in an algorithm with growing memory consumption per time step, which may not be appropriate for hardware implementation. Moreover, stability and error analysis, particularly for the error introduced by the approximation in the inverse Laplace transform, have not been studied in details. Furthermore, we believe there is merit for choosing the MLF for solving fractional differentiation equations because it is an entire function, in contrast to the singularity at 0 of $t^{-\beta }$, which limits its evaluation to δ as we have seen earlier. Therefore, we consider the study of exponential sum approximation of single-parameter MLF in the form of (1.7) and aim to justify its application to finite element modelling.

In this paper, we propose a scheme based on Gauss-Legendre quadrature to find the suitable coefficients for (1.1) to evaluate (1.7). There are two reasons for the choice of Gaussian quadrature. First, it yields mostly real coefficients except for one possible complex coefficient depending on α. Second, it may achieve uniform convergence for all $t>0$, which is partly enabled by the property that the MLF is an entire function, as we will see later. There were earlier papers which discuss the evaluation of the MLF in terms of exponential functions, but to the best of our knowledge, they cannot achieve uniform convergence for all $t>0$ [25], [26].

The fractional wave equation, which involves a higher integer-order time derivative than the diffusion problem in [5], is presented as an application example of the exponential sum approximation of the MLF. This wave equation is regarded as the most accurate model for sound propagation in a viscoelastic medium such as body tissue [27]. This wave model also finds application in modelling seismic wave propagation [2], and earthquake localisation [28].

This paper is organised as follows. Section 2 reviews the application of the method of exponential sum approximation to fast convolution, which is followed by the connection of the MLF to the viscoelastic problems. The MLF is also arranged into a form suitable for use with Gaussian quadrature and exponential sum approximation. In Section 3, we briefly review the standard approach to Gaussian quadrature error analysis based on complex analysis, prior to proving that the formulation introduced in Section 2 converges uniformly, as well as providing a rough estimate on the number of exponential terms required. Then in Section 4, the numerical scheme for the fractional wave equation is derived. The proofs of stability and error convergence behaviour are then presented. Section 5 gives numerical experiments and an example solution to the fractional wave equation to demonstrate the convergence. Finally, concluding remarks are given in Section 6.