William McLean

Hong-lin Liao, William McLean, Jiwei Zhang

We consider a class of numerical approximations to the Caputo fractional derivative. Our assumptions permit the use of nonuniform time steps, such as is appropriate for accurately resolving the behavior of a solution whose derivatives are singular at~$t=0$. The main result is a type of fractional Grönwall inequality and we illustrate its use by outlining some stability and convergence estimates of schemes for fractional reaction-subdiffusion problems. This approach extends earlier work that used the familiar L1 approximation to the Caputo fractional derivative, and will facilitate the analysis of higher order and linearized fast schemes.

Kassem Mustapha, William McLean

We consider an initial-boundary value problem for $\partial_tu-\partial_t^{-α}\nabla^2u=f(t)$, that is, for a fractional diffusion ($-1<α<0$) or wave ($0<α<1$) equation. A numerical solution is found by applying a piecewise-linear, discontinuous Galerkin method in time combined with a piecewise-linear, conforming finite element method in space. The time mesh is graded appropriately near $t=0$, but the spatial mesh is quasiuniform. Previously, we proved that the error, measured in the spatial $L_2$-norm, is of order $k^{2+α_-}+h^2\ell(k)$, uniformly in $t$, where $k$ is the maximum time step, $h$ is the maximum diameter of the spatial finite elements, $α_-=\min(α,0)\le0$ and $\ell(k)=\max(1,|\log k|)$. Here, we generalize a known result for the classical heat equation (i.e., the case $α=0$) by showing that at each time level $t_n$ the solution is superconvergent with respect to $k$: the error is of order $(k^{3+2α_-}+h^2)\ell(k)$. Moreover, a simple postprocessing step employing Lagrange interpolation yields a superconvergent approximation for any $t$. Numerical experiments indicate that our theoretical error bound is pessimistic if $α<0$. Ignoring logarithmic factors, we observe that the error in the DG solution at $t=t_n$, and after postprocessing at all $t$, is of order $k^{3+α_-}+h^2$.

William McLean

We solve a fractional diffusion equation using a piecewise-constant, discontinuous Galerkin method in time combined with a continuous, piecewise-linear finite element method in space. If there are $N$ time levels and $M$ spatial degrees of freedom, then a direct implementation of this method requires $O(N^2M)$ operations and $O(NM)$ active memory locations, owing to the presence of a memory term: at each time step, the discrete evolution equation involves a sum over \emph{all} previous time levels. We show how the computational cost can be reduced to $O(MN\log N)$ operations and $O(M\log N)$ active memory locations.

Kim Ngan Le, William McLean, Kassem Mustapha

We present a new stability and convergence analysis for the spatial discretization of a time-fractional Fokker--Planck equation in a convex polyhedral domain, using continuous, piecewise-linear, finite elements. The forcing may depend on time as well as on the spatial variables, and the initial data may have low regularity. Our analysis uses a novel sequence of energy arguments in combination with a generalized Gronwall inequality. Although this theory covers only the spatial discretization, we present numerical experiments with a fully discrete scheme employing a very small time step, and observe results consistent with the predicted convergence behavior.

Q. T. Le Gia, William McLean

We propose a method to construct numerical solutions of parabolic equations on the unit sphere. The time discretization uses Laplace transforms and quadrature. The spatial approximation of the solution employs radial basis functions restricted to the sphere. The method allows us to construct high accuracy numerical solutions in parallel. We establish $L_2$ error estimates for smooth and nonsmooth initial data, and describe some numerical experiments.

Kim Ngan Le, William McLean, Bishnu Lamichhane

An initial-boundary value problem for the time-fractional diffusion equation is discretized in space using continuous piecewise-linear finite elements on a polygonal domain with a re-entrant corner. Known error bounds for the case of a convex polygon break down because the associated Poisson equation is no longer $H^2$-regular. In particular, the method is no longer second-order accurate if quasi-uniform triangulations are used. We prove that a suitable local mesh refinement about the re-entrant corner restores second-order convergence. In this way, we generalize known results for the classical heat equation due to Chatzipantelidis, Lazarov, Thomée and Wahlbin.

Shev MacNamara, William McLean, Kevin Burrage

Contour integrals in the complex plane are the basis of effective numerical methods for computing matrix functions, such as the matrix exponential and the Mittag-Leffler function. These methods provide successful ways to solve partial differential equations, such as convection--diffusion models. Part of the success of these methods comes from exploiting the freedom to choose the contour, by appealing to Cauchy's theorem. However, the pseudospectra of non-normal matrices or operators present a challenge for these methods: if the contour is too close to regions where the norm of the resolvent matrix is large, then the accuracy suffers. Important applications that involve non-normal matrices or operators include the Black--Scholes equation of finance, and Fokker--Planck equations for stochastic models arising in biology. Consequently, it is crucial to choose the contour carefully. As a remedy, we discuss choosing a contour that is wider than it might otherwise have been for a normal matrix or operator. We also suggest a semi-analytic approach to adapting the contour, in the form of a parabolic bound that is derived by estimating the field of values. To demonstrate the utility of the approaches that we advocate, we study three models in biology: a monomolecular reaction, a bimolecular reaction and a trimolecular reaction. Modelling and simulation of these reactions is done within the framework of Markov processes. We also consider non-Markov generalisations that have Mittag-Leffler waiting times instead of the usual exponential waiting times of a Markov process.

William McLean, Vidar Thomée

In earlier work we have studied a method for discretization in time of a parabolic problem which consists in representing the exact solution as an integral in the complex plane and then applying a quadrature formula to this integral. In application to a spatially semidiscrete finite element version of the parabolic problem, at each quadrature point one then needs to solve a linear algebraic system having a positive definite matrix with a complex shift, and in this paper we study iterative methods for such systems. We first consider the basic and a preconditioned version of the Richardson algorithm, and then a conjugate gradient method as well as a preconditioned version thereof.

Hong-lin Liao, William McLean, Jiwei Zhang

Stability and convergence of a time-weighted discrete scheme with nonuniform time steps are established for linear reaction-subdiffusion equations. The Caupto derivative is approximated at an offset point by using linear and quadratic polynomial interpolation. Our analysis relies on two tools: a discrete fractional Grönwall inequality and the global consistency analysis. The new consistency analysis makes use of an interpolation error formula for quadratic polynomials, which leads to a convolution-type bound for the local truncation error. To exploit these two tools, some theoretical properties of the discrete kernels in the numerical Caputo formula are crucial and we investigate them intensively in the nonuniform setting. Taking the initial singularity of the solution into account, we obtain a sharp error estimate on nonuniform time meshes. The fully discrete scheme generates a second-order accurate solution on the graded mesh provided a proper grading parameter is employed. An example is presented to show the sharpness of our analysis.

William McLean

Given $β>0$ and $δ>0$, the function $t^{-β}$ may be approximated for $t$ in a compact interval $[δ,T]$ by a sum of terms of the form $we^{-at}$, with parameters $w>0$ and $a>0$. One such an approximation, studied by Beylkin and Monzón, is obtained by applying the trapezoidal rule to an integral representation of $t^{-β}$, after which Prony's method is applied to reduce the number of terms in the sum with essentially no loss of accuracy. We review this method, and then describe a similar approach based on an alternative integral representation. The main difference is that the new approach achieves much better results before the application of Prony's method; after applying Prony's method the performance of both is much the same.

Kim Ngan Le, William McLean, Kassem Mustapha

We study two schemes for a time-fractional Fokker-Planck equation with space- and time-dependent forcing in one space dimension. The first scheme is continuous in time and is discretized in space using a piecewise-linear Galerkin finite element method. The second is continuous in space and employs a time-stepping procedure similar to the classical implicit Euler method. We show that the space discretization is second-order accurate in the spatial $L_2$-norm, uniformly in time, whereas the corresponding error for the time-stepping scheme is $O(k^α)$ for a uniform time step $k$, where $α\in(1/2,1)$ is the fractional diffusion parameter. In numerical experiments using a combined, fully-discrete method, we observe convergence behaviour consistent with these results.