Arbitrary order of convergence for Riesz fractional derivative via central difference method

Riesz derivative is the symmetrical version of the fractional Riemann–Liouville derivative. Fractional calculus is an important topic because fractional derivatives can account for non-local effect, which is present in many physics problems. For example, the Riesz derivative has been applied to model wave propagation in viscoelastic material and fractional electrical impedance in human skin [1, 2]. One of the most efficient methods to approximate spatial fractional derivatives is based on Grünwald-Letnikov derivatives. The development began with the derivation of the central difference scheme to approximate the Riesz potential in [3]. The order of accuracy was found to be of second order later in [4]. In [5, 6], the order of convergence was increased to 4 by applying another finite difference on the operator such that the second order term in the Maclaurin series is eliminated. Because the matrix system resulting from the fractional difference operation is not sparse, there is practically no additional cost in using a highest order approximation in a dynamic problem aside from the potential overhead in the evaluation of the coefficients for the operator. This motivates us to study methods to increase the order of approximation to an arbitrary order so that the convergence is faster, leading to reduced number of nodes required to achieve accurate results. In this paper, we propose a prefilter approach to generalise the canceling of error terms in [5, 6]. The result is a filter which maximizes the flatness of the fractional power of sinc function at the origin.

In order to understand our approach, we first discuss the development of higher order finite difference approximation using the Grünwald-Letnikov derivatives for Riesz derivative. Starting from the very basic, Grünwald-Letnikov derivatives are obtained by applying fundamental theorem of calculus recursively as

and so on. Shifting the difference to the centre, we have

Without the limits, the first Grünwald-Letnikov derivative is equivalent to first order accurate finite difference. Because of this, the subsequent derivatives are also first order accurate while the central difference version has second order of accuracy.

Grünwald-Letnikov derivatives of a function can be written as a discrete convolution of the function evaluated at the discrete points with a stencil. For the central difference, the derivatives can be expressed as

where the discrete convolution operator is defined as $\left(f\ast g\right)\left[n\right]\triangleq\sum_{m}f\left[n-m\right]g\left[m\right]$, $f\left[n\right]=f\left(nh\right)$, $\mathrm{S}_{t}$ is the shift operator defined as $\mathrm{S}_{t}\left\{f\right\}\left[n\right]=f\left[n+t\right]$, and the array is 0 indexed. Now, let $H_{k}$ be the stencil of $k$-th derivative. Then the $k$-th derivative can be written as

with the recursive relationship

Applying discrete time Fourier transform to $H_{k}$, we have

This leads to the solution $\{H_{k}\}_{m}=\left(-1\right)^{m}\binom{k}{m}$. To understand the error convergence, we replace $m$ with $m\frac{h}{2}$ in (12) to obtain the Fourier transform of the $k$-th derivative of $f$

Since $\operatorname{\mathrm{sinc}}\left(x\right)$ can be expanded as a Maclaurin series as $1+O\left(x^{2}\right)$ and it would converge for $x$ from $-\pi$ to $\pi$, this approximation is second order accurate, assuming $f$ is sufficiently smooth such that the magnitude of its frequency response drops significantly below the difference between 1 and the sinc function in its main lobe.

The binomial series (12) can be generalised to non-integer $k$ for fractional derivatives, as seen in [3, 7], because the magnitude of the variable in the series is not bigger than 1. Consider the Riesz derivative defined as

where $f\left(x\right)=0$ everywhere except $x\in\left[a,b\right]$, $C_{\alpha}=-\frac{1}{2\cos\frac{\pi\alpha}{2}}$, and the left and right Riemann-Liouville fractional derivatives are defined as

where $m=\left\lceil\alpha\right\rceil$. Since $f$ can be scaled by introducing dimensionless variables, we consider $a=0$, $b=1$ for the Riesz derivative throughout the paper. The Riesz derivative has the frequency response $-\left|\omega\right|^{\alpha}$. Therefore, for a stencil which approximates the Riesz derivative to second order of accuracy, we seek a solution of the inverse transform of $\left|\omega\operatorname{\mathrm{sinc}}\left(\frac{\omega h}{2}\right)\right|^{\alpha}$. Let the fractional convolution kernel be $K_{\alpha}$ such that

Then the solution is the Fourier series coefficient of its transform given by

The coefficients are symmetrical because the frequency response is real and symmetric about the origin. Therefore, it does not matter whether the exponent is negative or not but the negative choice will be apparent later. Owning to the binomial series and Cauchy residue theorem, each biomial coefficent can be expressed in the integral form as

for any real $a$ and $b>-1$. Eq. (18) can be rewritten in the same form to arrive at the solution as

This method to solve binomial related problems is known as the Egorychev method[8]. Therefore, to increase the convergence order, we propose a filter with a frequency response $G_{\alpha}$ such that $G_{\alpha}\cdot K_{\alpha}$ now has a response of $\left|\omega\right|^{\alpha}\left(1+O\left(h^{N}\right)\right)$, where $N$ is the order of convergence.

This paper is organized as follows. In Section 2, we start with presenting our method, which is then followed by the details of the algorithms. In Section 3, we discuss the convergence behaviour of the method, which shows that the convergence is indeed $O\left(h^{N}\right)$ subject to the smoothness of the function, and we comment on the filter’s positivity and the eigenvalue problem. Numerical results are presented in Section 4 to validate the error prediction in the previous section before the conclusion in the final section.

In this section, we derive the method to reduce the error of approximation of the Riesz derivative and algorithms to compute the stencil for this purpose. As discussed previously, we seek a filter to improve the flatness of the sinc function in the neighbourhood of the origin such that resulting frequency response of the finite difference operator is closer to $-\left|\omega\right|^{\alpha}$. We expand the fractional power of the sinc function in (18) as the Maclaurin series as

where $N$ is the choice of convergence order, $N_{h}=\left\lceil\frac{N}{2}\right\rceil-1$, $a_{n}=\frac{1}{\left(2n\right)!}\frac{\mathrm{d}\operatorname{\mathrm{sinc}}\left(\frac{x}{2}\right)^{\alpha}}{\mathrm{d}x}|_{x=0}$. The odd coefficients are ignored because it is obvious that they are 0 for an even function. Because this series converges up to $\left|x\right|<\pi$, we know that each higher order term is always smaller in magnitude than the respective 1 lower order term below the Nyquist frequency. Then it is clear that by cancelling out the lower order ones, we will improve the flatness of this function, provided that the filter does not raise the magnitude of the coefficients of the higher order terms. Now, let us consider only the problem of eliminating the lower order terms first. Let the filter be $G_{\alpha}$ such that

The transform is a cosine series because the filter is obviously real symmetric, and there are no odd terms in the sinc expansion to cancel. Expanding each $\cos\left(nx\right)$ as Maclaurin series, we have

where $c_{n,m}=\frac{\left(-1\right)^{k}n^{2k}}{\left(2k\right)!}.$ Let $b_{n}$ be the series coefficients of $\tilde{G}_{\alpha}$. They can be expressed in matrix form as

The coefficients of the series from the product of the two polynomials can be obtained through discrete convolution, and we set the coefficients for $n>0$ up to $2N_{h}$ be 0 to eliminate the lower order terms. This can be expressed in matrix form as

The matrix involving $c_{n,m}$ is not convenient to solve but it gives us insight of the actual problem size. We rewrite it as a larger problem, in terms of the expansion of the exponential function, leading to

where $\mathbf{V}\left(\left\{x_{0},x_{1},\ldots,x_{N}\right\}\right)$ is the Vandermonde matrix defined as

Let $v_{i,j}^{N}=\left\{\mathbf{V}^{-1}\left(\left\{x_{0},x_{1},\ldots,x_{N}\right\}\right)\right\}_{i,j}$, the inverse of the Vandermonde matrix can be solved using the following recurrence equations [9]

and for $0\leq j<k$

where $k$ is looped over $k=0,\ldots,N$. All $v_{i,j}^{0}$ are 0 initialised except for $v_{0,0}^{0}=1$. However, for our special case, there are some properties which can be exploited for further optimisation. In this special case, $k$ is even, and $x_{\frac{k}{2}+n}=x_{\frac{k}{2}-n}=n$ for $n=0,\,\ldots,\,\frac{k}{2}$. Since we are only concerned about even $k$, we map $k\mapsto\frac{k}{2}$ and let the size grow 2 per step. To find the new recurrence relationship for each step, we write out the Lagrange polynomial basis as

Each basis has the following symmetry

which leads to $v_{2i,N+j}^{k}=v_{2i,N-j}^{k}$. This relationship also arises from the fact that $g_{j}$ is symmetric. Therefore, we need not consider the left half of the Vandermonde matrix inverse, and we can shift the indexing to the left by $N$.

The second property is

so the odd terms are not required for the evaluation of future coefficients. To prove this, we write out the recurrence relationship for this special case first. For $j=k$, we have

which leads to the recurrence relationship of the coefficient after substituting in (32),

and for $j<k$,

leads to

One can verify that the property (34) is true for $k=1$ using the recurrence relationships. We assume this is also true for other $k$ up to $N-1$. Then for $k=N$, $j=N$, we substitute (36) into (34), we get

This is similarly proven by induction for the case $j<N$ as such

For the case $i=1$, one can see it is also true because from (37) and (40), we have $v_{0,j}^{k}=0$ for all $j>0$ and $v_{0,0}^{k}=1$. Therefore, by induction, (34) is true for all $k>0$. All together, these $v_{i,j}^{N}$ at even $i$ and right half of the Vandermonde matrix inverse constitute the inverse of the matrix product involving $c_{n,m}$ and the diagonal matrix on the left hand side of (24).

Next, we discuss the evaluation of $a_{j}$, the derivatives of the fractional power of sinc function at 0. This function is a composite function and its derivatives are given by Faà di Bruno’s formula. However, it is not convenient for use because of the involvement of combinatorics. Instead, one can simply find the corresponding Maclaurin series coefficients. They may be efficiently recovered from the inverse transform of fractional power of the fast Fourier transform of the Maclaurin series coefficients of the sinc function given by

while the odd cases are all 0. To prove this, one can simply apply Cauchy residue theorem again similar to (19) for the infinite case. And because the Maclaurin series converges about the unit circle for the sinc function, the result is valid. Then to prove that this is also true for the finite case, one just has to realise that the fractional number can be split up into an integer numerator $p$ and an integer denominator $q$. Then there exists some transform coefficients such that their $q$ power gives the original transform coefficients because the integer power of transform coefficients is just convolution of the coefficients with themselves an integer number of times. Therefore, the same coefficients from the infinite case can be recovered using only the finite transform coefficients. Since we would already have the transform of those coefficients, (25) can be solved by directly applying the inverse transform to the reciprocal of the transform coefficients. All in all, let $\check{a}_{n}$ be the $2n$-th Maclaurin series coefficients of the sinc function given by (43), $b_{n}$ in (25) are given by

where $\omega_{m}=\frac{2\pi m}{2N+1}$.

To summarise, we provide the pseudo code to obtain the stencil for higher order central difference scheme for the Riesz derivative here. Algorithm 1 describes the procedure to invert the transpose of the matrix product on the left hand side of (24). The columns are scaled to better balance the scaling of the sinc Maclaurin series coefficients. Algorithm 2 outputs the right half of the stencil for $N_{x}$ number of elements between 0 and 1. There are $N_{x}-1$ number of elements in the output because the derivatives at the end points are not needed in dynamic problems as the boundary condition is set there. Let $\mathbf{k}$ be the vector returned from Algorithm 2, a $\left(N_{x}-2\right)\times\left(N_{x}-2\right)$ Toeplitz matrix $\mathbf{D}$ can be constructed from $\mathbf{k}$ with $\left\{\mathbf{D}\right\}_{n,m}=$$\left\{\mathbf{k}\right\}_{\left|n+1-m\right|}$ to approximate the derivative as

If adaptive grid size is required, one can save the last $2N_{h}$ elements of $\mathbf{k}$ before the convolution as well as the filter coefficients of $\mathbf{g}$ for the generation of elements in higher indices. The generation of the coefficients can then be resumed, after scaling the coefficients for the larger stencil, from line 14 of Algorithm 2.

First, we study the error convergence. We follow a similar approach given in [5, 6]. A lemma for the smoothness requirement of the function is given before the error analysis. Because of this lemma, which can be considered an extension to the Riemann-Lebesgue lemma, the order requirement for derivative continuity is lower than the ones given in [5, 6].

In Theorem 1, we have limited $\alpha$ to be a positive number smaller than 2, which is typical for physics problems, but it can also be generalised to higher order. The smoothness requirement means that this method does require 0 boundary condition. Because the algorithm and analysis is based on the expansion around the origin, the convergence order drops off as we approach Nyquist frequency. For demonstration, we plot the frequency response of the central difference operator divided by the response of Riesz operator given by $\operatorname{\mathrm{sinc}}^{\alpha}\left(\frac{x}{2}\right)\tilde{G}_{\alpha}\left(x\right)$ and the rate of convergence with respect to frequency given by

where $\tilde{G}_{\alpha}$ is given by (22), for $x\in\left[0,\pi\right]$ and $\alpha=1.3$. The plots are similar for various $\alpha$, therefore we only plot them for a single value. We see in Figure 1(a) that the flatness of the curve around the origin is improved with increasing $N$. From Figure 1(b), we observe that the higher the convergence order parameter $N$ is, the greater the drop in the convergence order with respect to the normalized angular frequency. This decrease will limit the convergence rate depending on the bandwidth of the function. Despite this, a higher $N$ still results in a greater convergence order as long as the bandwidth of the function is less than 2 times the Nyquist frequency.

Because the Vandermonde matrix can be decomposed into lower and upper triangular matrices without any zero elements in the diagonal entries [10], solutions exist for any $N$. However, because the matrix system (25) does not give us control over derivatives beyond $N$-th order, there is no guarantee that this method ensures that there is no overshoot in the resulting response and the response is always ‘flatter’ than the result of a lower order one. How close the bound (54) is depends on the constant $C_{1}$, related to the coefficient of $\left(\omega h\right)^{N}$ in the expansion (52), which the formulation does not have explicit control. Nonetheless, numerical results show that there is improvement on the frequency response up to $N=46$ without extra constraints imposed on the magnitude of higher order Maclaurin series coefficients. The proposed algorithm does become unstable at higher $N$ using double precision floating point.

The Riesz derivative is a negative definite operator, so $\mathbf{g}$ has to be positive definite for the correct results. We do not provide a proof for this but numerical results suggest $\left|\mathbf{g}_{0}\right|>\sum_{n>0}\left|\mathbf{g}_{n}\right|$, which is sufficient for the filter to be positive definite because of Gershgorin circle theorem. Therefore, positivity can be checked with minimal overhead during application of the method. Moreover, because the filter response is symmetric about the point $\pi$, which also means that point is an extremum, and also because the filter should be increasing from the origin to cancel out the downward slope of the sinc function, we should expect that the filter has a response greater than 1 with a maximum at $\pi$. As for the magnitude of this maximum, which leads to the upper bound of the eigenvalue of the Toeplitz matrix in (46), if we assume that $0<\left|\operatorname{\mathrm{sinc}}\left(\frac{x}{2}\right)\right|^{\alpha}\tilde{G}_{\alpha}\left(x\right)\leq 1$ for all $N$ then the upper bound for the maxima is $\pi^{\alpha}$ and so the magnitude of the eigenvalue of the Toeplitz matrix is bounded by $\left(\pi N_{x}\right)^{\alpha}$.

To verify the error analysis, we apply the derivative operator to