Efficient high order algorithms for fractional integrals and fractional differential equations

TL;DR

This paper introduces an efficient, high-order algorithm for fractional integrals and differential equations using Runge--Kutta convolution quadrature, with rigorous error estimates and demonstrated numerical performance.

Contribution

It presents a novel integral representation and quadrature method for fractional integrals, enabling high-order, efficient solutions with rigorous error analysis.

Findings

Algorithm achieves high accuracy with low memory and computational cost.

Numerical results validate the efficiency and accuracy of the method.

Application to fractional diffusion equations demonstrates practical utility.

Abstract

We propose an efficient algorithm for the approximation of fractional integrals by using Runge--Kutta based convolution quadrature. The algorithm is based on a novel integral representation of the convolution weights and a special quadrature for it. The resulting method is easy to implement, allows for high order, relies on rigorous error estimates and its performance in terms of memory and computational cost is among the best to date. Several numerical results illustrate the method and we describe how to apply the new algorithm to solve fractional diffusion equations. For a class of fractional diffusion equations we give the error analysis of the full space-time discretization obtained by coupling the FEM method in space with Runge--Kutta based convolution quadrature in time.

Tables

Table 1. Table 1: Dependence of the total number of quadrature points on time step h ℎ h and final time T 𝑇 T . The other parameters are fixed at n 0 = 5 subscript 𝑛 0 5 n_{0}=5 , B = 3 𝐵 3 B=3 , tol = 10 − 6 tol superscript 10 6 \mathrm{tol}=10^{-6} , α = 0.5 𝛼 0.5 \alpha=0.5 . On the left the data is for backward Euler and on the right for the 2-stage Radau IIA CQ.

$h\backslash T$	$1$	$10$	$100$	$1000$
${10}^{-1}$	$20$	$30$	$40$	$49$
${10}^{-2}$	$27$	$36$	$44$	$52$
${10}^{-3}$	$31$	$39$	$46$	$50$
${10}^{-4}$	$34$	$40$	$45$	$48$

Table 2. Table 2: Dependence of the total number of quadrature points on the tolerance tol tol \mathrm{tol} and the fractional power α 𝛼 \alpha . The other parameters are fixed at h = 10 − 2 ℎ superscript 10 2 h=10^{-2} , T = 50 𝑇 50 T=50 , n 0 = 5 subscript 𝑛 0 5 n_{0}=5 , B = 3 𝐵 3 B=3 . Again the data on the left is for backward Euler and on the right for 2-stage Radau IIA.

$\mathrm{tol}\backslash \alpha$	$0.1$	$0.3$	$0.5$	$0.7$	$0.9$
${10}^{-2}$	$11$	$11$	$10$	$8$	$6$
${10}^{-4}$	$27$	$27$	$26$	$25$	$21$
${10}^{-6}$	$45$	$44$	$45$	$43$	$36$
${10}^{-8}$	$66$	$65$	$64$	$61$	$55$
${10}^{-10}$	$86$	$87$	$85$	$82$	$74$

Table 3. Table 3: We show the error and the memory requirements for the new method and the standard implementation of CQ.

$N$	error	memory (MB)	standard err.	standard mem. (MB)
32	$2.94\times {10}^{-1}$	39.1	$2.94\times {10}^{-1}$	59.2
64	$3.07\times {10}^{-2}$	40.3	$3.07\times {10}^{-2}$	98.7
128	$2.61\times {10}^{-3}$	42.8	$2.61\times {10}^{-3}$	177.6
256	$2.98\times {10}^{-4}$	44.0	$3.01\times {10}^{-4}$	335.4

Equations

$${\mathcal{I}}^{\alpha }[f](t)=\frac{1}{\Gamma (\alpha )}{\int }_0^t(t-s{)}^{\alpha -1}f(s)ds,$$

$${\mathcal{I}}^{\alpha }[f](t_n)\approx \sum _{j=0}^n{\mathbf{\omega }}_{n-j}{\mathbf{f}}_j,$$

$${\mathbf{\omega }}_n=\frac{h\sin (\pi \alpha )}{\pi }{\int }_0^{\infty }{x}^{-\alpha }{\mathbf{e}}_n(-hx)dx,$$

$$t^{\alpha -1}=\frac{1}{\Gamma (1-\alpha )}{\int }_0^{\infty }{x}^{-\alpha }{e}^{-tx}dx,$$

$${\partial }_t^{\beta }u-\Delta u=f,u^{(j)}(0)=0,j=0,\dots ,m-1,$$

$$t^{\alpha -1}=\frac{1}{\Gamma (1-\alpha )}{\int }_0^L{x}^{-\alpha }{e}^{-tx}dx+\tau (L),$$

$$|\tau (L)|\le \frac{h^{\alpha -1}}{\Gamma (1-\alpha )}{\int }_A^{\infty }{x}^{-\alpha }{e}^{-n_{0}x}dx.$$

$$\begin{array}{rl}|\tau (L)|& =\frac{h^{\alpha -1}}{\Gamma (1-\alpha )}{\int }_A^{\infty }{x}^{-\alpha }{e}^{-\frac{t}{h}x}dx\\ & \le \frac{h^{\alpha -1}{A}^{-\alpha }}{\Gamma (1-\alpha )}{\int }_A^{\infty }{e}^{-n_{0}x}dx\\ & =\frac{A^{-\alpha }{h}^{\alpha -1}{e}^{-n_{0}A}}{n_0\Gamma (1-\alpha )}.\end{array}$$

$$A+\frac{\alpha }{n_0}\log (A)\ge \frac{1}{n_0}\log \left(\frac{h^{\alpha -1}}{n_0\Gamma (1-\alpha )\mathrm{tol}}\right).$$

$$A=\log \left(\frac{1}{n_0\Gamma (1-\alpha )\mathrm{tol}}\right)+(1-\alpha )\log \left(\frac{1}{h}\right).$$

$$I_0=\frac{1}{\Gamma (1-\alpha )}{\int }_0^{L_0}{x}^{-\alpha }{e}^{-tx}dx,$$

$$|If-I_{Q}f|\le 4M\frac{{\varrho }^{-2Q+1}}{\varrho -1}{\int }_{-1}^{1}w(x)dx,$$

$$f(x)=\sum _{k=0}^{\infty }{a}_k{T}_k(x)$$

$$|f-f_K|\le \frac{2M{\varrho }^{-K}}{\varrho -1}.$$

$$\begin{array}{rl}|If-I_{Q}f|& =\left|I(f-f_{2Q-1})-I_Q(f-f_{2Q-1})\right|\\ & \le \frac{2M{\varrho }^{-2Q+1}}{\varrho -1}\left({\int }_{-1}^{1}w(x)dx+\sum _{j=1}^{2Q-1}{w}_j\right)\\ & =\frac{4M{\varrho }^{-2Q+1}}{\varrho -1}{\int }_{-1}^{1}w(x)dx,\end{array}$$

$$I_0=\frac{1}{\Gamma (1-\alpha )}{\left(\frac{L_0}{2}\right)}^{1-\alpha }{\int }_{-1}^1{e}^{-t(y+1)L_0/2}(y+1{)}^{-\alpha }dy.$$

$$f_0(x)=\frac{1}{\Gamma (1-\alpha )}{\left(\frac{L_0}{2}\right)}^{1-\alpha }{e}^{-t(x+1)L_0/2},w(x)=(1+x{)}^{-\alpha }$$

$${\tau }_{\text{GJ}}(Q):=I{f}_0-I_Q{f}_0={\int }_{-1}^1{f}_0(x)w(x)dx-\sum _{j=1}^Q{w}_j{f}_0(x_j).$$

$$|{\tau }_{\mathrm{GJ}}(Q)|\le \frac{4L_0^{1-\alpha }}{\Gamma (2-\alpha )}\left(1+\frac{T{L}_0}{4Q}\right){\left(\frac{eT{L}_0}{8Q}\right)}^{2Q}.$$

$$|{\tau }_{\text{GJ}}(Q)|$$

$$|{\tau }_{\text{GJ}}(Q)|\le \frac{4L_0^{1-\alpha }}{\Gamma (2-\alpha )}\underset{\delta >0}{\min }\frac{e^{\delta }}{e^{\delta }-1}{e}^{-2Q\delta +L_{0}T\left(\cosh \delta -1\right)/2}.$$

$$g(\delta )=2Q\delta -L_{0}T\left(\cosh \delta -1\right)/2.$$

$$g^{\prime }(\delta )=2Q-L_{0}T\sinh \delta /2,g^{\prime \prime }(\delta )=-L_{0}T\cosh \delta /2<0,$$

$$2Q-L_{0}T\sinh \delta /2=0⟹\delta ={\sinh }^{-1}\left(\frac{4Q}{T{L}_0}\right).$$

$${\sinh }^{-1}y=\log \left(y+\sqrt{1+y^2}\right),\cosh x=\sqrt{1+{\sinh }^{2}x},$$

$$\begin{array}{rl}|{\tau }_{\text{GJ}}(Q)|& \le \frac{4L_0^{1-\alpha }}{\Gamma (2-\alpha )}\left(1+\frac{T{L}_0}{4Q}\right)e^{-2Q\delta +L_{0}T(\cosh \delta -1)/2}\\ & \le \frac{4L_0^{1-\alpha }}{\Gamma (2-\alpha )}\left(1+\frac{T{L}_0}{4Q}\right){\left(\frac{T{L}_0}{8Q}\right)}^{2Q}{e}^{L_{0}T\left(-1+\sqrt{1+(4Q/(T{L}_0){)}^2}\right)/2}\\ & \le \frac{4L_0^{1-\alpha }}{\Gamma (2-\alpha )}\left(1+\frac{T{L}_0}{4Q}\right){\left(\frac{T{L}_0}{8Q}\right)}^{2Q}{e}^{2Q},\end{array}$$

$$\frac{1}{\Gamma (1-\alpha )}{\int }_{L_0}^L{x}^{-\alpha }{e}^{-xt}dx=\sum _{j=1}^J{I}_j,$$

$$\begin{array}{rl}{I}_j& =\frac{1}{\Gamma (1-\alpha )}{\int }_{L_{j-1}}^{L_j}{x}^{-\alpha }{e}^{-xt}dx\\ & =\frac{\Delta L_j}{2\Gamma (1-\alpha )}{e}^{-L_{j-1}t}{\int }_{-1}^1{\left(L_{j-1}+\frac{\Delta L_j}{2}(y+1)\right)}^{-\alpha }{e}^{-t(y+1)\Delta L_j/2}dy,\end{array}$$

$${\tau }_j(Q):=I{f}_j-I_Q{f}_j$$

$$f_j(x)=\frac{4\Delta L_j}{\Gamma (1-\alpha )}{e}^{-L_{j-1}t}{\left(L_{j-1}+\frac{\Delta L_j}{2}(x+1)\right)}^{-\alpha }{e}^{-t(x+1)\Delta L_j/2}.$$

Full text

Keywords: fractional integral, fractional differential equations, convolution quadrature, fast and oblivious algorithms.

AMS subject classifications: 65R20, 65L06, 65M15,26A33,35R11.

1 Introduction

Fractional Differential Equations (FDEs) have nowadays become very popular for modeling different physical processes, such as anomalous diffusion [26] or viscoelasticity [1, 25]. In the present paper we develop a fast and memory efficient method to compute the fractional integral

[TABLE]

for a given $\alpha\in(0,1)$. A standard discretization of (1) is obtained by convolution quadrature (CQ) based on a Runge-Kutta scheme [18, 5]

[TABLE]

where the convolution weights $\boldsymbol{\omega}_{n}$ can be expressed as, see Lemma 9,

[TABLE]

with ${\bf e}_{n}(\cdot)$ a function that depends on the Runge-Kutta scheme. For discretizations based on linear multistep methods, see [14].

To compute up to time $T=Nh$ using formula (2) requires $O(N)$ memory and $O(N^{2})$ arithmetic operations. Algorithms based on FFT can reduce the computational complexity to $O(N\log N)$ [16] or $O(N\log^{2}N)$ [9], but not the memory requirements; for an overview of FFT algorithms see [7]. Here we develop algorithms that reduce the memory requirement to $O(|\log\varepsilon|\log N)$ and the computational cost to $O(|\log\varepsilon|N\log N)$, with $\varepsilon$ the accuracy in the computation of the convolution weights. Hence, our algorithm has the same complexity as the fast and oblivious quadratures of [19] and [23], but as we will see, a simpler construction.

The algorithms will depend on an efficient quadrature of (3) for $n\geq n_{0}$, with a very moderate threshold value for $n_{0}$, say $n_{0}=5$. As ${\bf e}_{n}(z)=r(z)^{n}{\bf q}(z)$ and $r(z)=e^{z}+O(z^{p+1})$, where $p$ is the order of the underlying RK method, this is intimately related to the construction of an efficient quadrature for the integral representation of the convolution kernel

[TABLE]

with $t\in[n_{0}h,T]$. Note that as $\Gamma(1-\alpha)\Gamma(\alpha)=\pi/\sin(\pi\alpha)$, $h^{-1}\boldsymbol{\omega}_{n}$ is an approximation of $\tfrac{1}{\Gamma(\alpha)}t^{\alpha-1}$, i.e., the kernel of (1).

Even though we eventually only require the quadrature for (3), we begin with developing a quadrature formula for (4) for a number of reasons: the calculation for (4) is cleaner and easier to follow, such a quadrature allows for efficient algorithms that are not based on CQ, and finally once this is available the analysis for (3) is much shorter. The quadrature we develop for (4) is closely related to the one developed in [13], the main difference being our treatment of the singularity at $x=0$ by Gauss-Jacobi quadrature and the restriction of $t$ to the finite interval rather than semi-infinite as used in [13]. Both these decisions allow us to substantially reduce constants in the above asymptotic estimates of memory and computational costs. Recent references [27, 12, 2] also consider fast computation of (1), but do not address the approximation of the convolution quadrature approximation exploiting (3). Our main contribution here is the development of an efficient quadrature to approximate (3) and its use in a fast and memory efficient scheme for computing the discrete convolution (2).

The stability and convergence properties of RK convolution quadrature are well understood, see [18, 6]. This allows us to apply convolution quadrature not only to the evaluation of fractional integrals, but also to the solution of fractional subdiffusion or diffusion-wave equations of the form

[TABLE]

with $\beta\in(0,2)$. Here, $\partial_{t}^{\beta}=\mathcal{I}^{m-\beta}\partial_{t}^{m}$, with $m=\lceil\beta\rceil$, denotes the Caputo fractional derivative. Solutions of such equations typically have low regularity at $t=0$, but a discussion of adaptive or modified quadratures for this case is beyond the scope of the current paper. For a careful analysis of BDF2 based convolution quadrature of fractional differential equations see [8].

To our knowledge, underlying high order solvers for ODEs have been considered for the approximation of (1) only at experimental level in [2, 3] and in [23], where a fast and oblivious implementation of RK based CQ is considered for more general applications than (1). The fast and oblivious quadratures of [19] and [23] have the same asymptotic complexity as our algorithm, but have a more complicated memory management structure and require the optimization of the shape of the integration contour. Our new algorithm has the advantage of being much easier to implement, as it does not require sophisticated memory management and the optimization of quadrature parameters is much simpler, and furthermore only real arithmetic is required. The new method is also much better suited for the extension to variable steps — this will be investigated in a follow up work. On the other hand, the present algorithm is specially tailored to the application to (1) and related FDEs, whereas the algorithms in [19, 23] allow for a wider range of applications.

The paper is organized as follows. In Section 2 we develop and fully analyze a special quadrature for (1), which uses the same nodes and weights for every $t\in[n_{0}h,T]$. In Section 3, we recall Convolution Quadrature based on Runge–Kutta methods and derive the special representation of the associated weights already stated in (3). In Section 4 we derive a special quadrature for (3), which uses the same nodes and weights for every $n\in[n_{0},N]$, with $T=hN$. In Section 5 we explain how to turn our quadrature for the CQ weights into a fast and memory saving algorithm. In Section 6 we test our algorithm with a scalar problem and in Section 7 we consider the application to a fractional diffusion equation. We provide a complete error analysis of the discretization in space and time of a class of fractional diffusion equations.

2 Efficient quadrature for $t^{\alpha-1}$

In the following we fix an integer $n_{0}>0$, time step $h>0$, and the final computational time $T>0$. Throughout, the parameter $\alpha$ is restricted to the interval $(0,1)$. We develop an efficient quadrature for (4) accurate for $t\in[n_{0}h,T]$.

2.1 Truncation

First of all we truncate the integral

[TABLE]

where $\tau(L)$ denotes the truncation error.

Lemma 1**.**

For $t\geq n_{0}h$ and $L=A/h$ we have that

[TABLE]

Proof.

[TABLE]

∎

Remark 2**.**

Given a tolerance $\mathrm{tol}>0$, $|\tau(L)|\leq\mathrm{tol}$ if

[TABLE]

Assuming $A\geq 1$ we can choose

[TABLE]

However, in practice it is advantageous to use the bound (5) to numerically find the optimal $A$.

2.2 Gauss-Jacobi quadrature for the initial interval

We choose an initial integration interval

[TABLE]

along which we will perform Gauss-Jacobi integration.

Recall the Bernstein ellipse $\mathcal{E}_{\varrho}$, which is given as the image of the circle of radius $\varrho>1$ under the map $z\mapsto(z+z^{-1})/2$. The largest imaginary part on $\mathcal{E}_{\rho}$ is $(\varrho-\varrho^{-1})/2$ and the largest real part is $(\varrho+\varrho^{-1})/2$.

Theorem 3**.**

Let $f$ be analytic inside the Bernstein ellipse $\mathcal{E}_{\varrho}$ with $\varrho>1$ and bounded there by $M$. Then the error of Gauss quadrature with weight $w(x)$ is bounded by

[TABLE]

where $If=\int_{-1}^{1}w(x)f(x)dx$ and $I_{Q}f=\sum_{j=1}^{Q}w_{j}f(x_{j})$ is the corresponding Gauss formula, with weights $w_{j}>0$.

Proof.

A proof of this result for $w(x)\equiv 1$ can be found in [24, Chapter 19]. The same proof works for the weighted Gauss quadrature as well. We give the details next.

First of all note that we can expand $f$ in Chebyshev series

[TABLE]

with $|a_{k}|\leq 2M\varrho^{-k}$ [24, Theorem 8.1]. If we denote by $f_{K}(x)=\sum_{k=0}^{K}a_{k}T_{k}(x)$ the truncated series then

[TABLE]

As $I_{Q}$ is exact for polynomials of degree $2Q-1$, we have that

[TABLE]

where we have used the fact the weights are positive and integrate constants exactly. ∎

Changing variables to the reference interval $[-1,1]$ we obtain

[TABLE]

We apply Theorem 3 to the case

[TABLE]

and denote

[TABLE]

Theorem 4**.**

For $t\in[0,T]$ and any $Q\geq 1$ we have the bound

[TABLE]

Proof.

Since $f_{0}$ in (7) is an entire function, by Theorem 3 we can estimate

[TABLE]

Let $\varrho=e^{\delta}$ with $\delta>0$. Then the error bound can be written as

[TABLE]

We now choose $\delta$ so that it maximises the function

[TABLE]

As

[TABLE]

we have a maximum at

[TABLE]

Using the identities

[TABLE]

we derive an error estimate with the above choice of $\delta$:

[TABLE]

where in the last step above we have used that $-1+\sqrt{1+x^{2}}\leq x$ for $x>0$. This gives the stated result. ∎

2.3 Gauss quadrature on increasing intervals

We next split the remaining integral as

[TABLE]

where

[TABLE]

where $\Delta L_{j}=L_{j}-L_{j-1}$, $j=1,\dots,J$, with $L_{J}=L$. The intervals are chosen so that for some $B\geq 1$, $\Delta L_{j}=BL_{j-1}$, i.e.,