Kevin Burrage

Ling Guo, Fanhai Zeng, Ian Turner et al.

In this work, we extend the fractional linear multistep methods in [C. Lubich, SIAM J. Math. Anal., 17 (1986), pp.704--719] to the tempered fractional integral and derivative operators in the sense that the tempered fractional derivative operator is interpreted in terms of the Hadamard finite-part integral. We develop two fast methods, Fast Method I and Fast Method II, with linear complexity to calculate the discrete convolution for the approximation of the (tempered) fractional operator. Fast Method I is based on a local approximation for the contour integral that represents the convolution weight. Fast Method II is based on a globally uniform approximation of the trapezoidal rule for the integral on the real line. Both methods are efficient, but numerical experimentation reveals that Fast Method II outperforms Fast Method I in terms of accuracy, efficiency, and coding simplicity. The memory requirement and computational cost of Fast Method II are $O(Q)$ and $O(Qn_T)$, respectively, where $n_T$ is the number of the final time steps and $Q$ is the number of quadrature points used in the trapezoidal rule. The effectiveness of the fast methods is verified through a series of numerical examples for long-time integration, including a numerical study of a fractional reaction-diffusion model.

Fanhai Zeng, Ian Turner, Kevin Burrage

A unified fast time-stepping method for both fractional integral and derivative operators is proposed. The fractional operator is decomposed into a local part with memory length $ΔT$ and a history part, where the local part is approximated by the direct convolution method and the history part is approximated by a fast memory-saving method. The fast method has $O(n_0+\sum_{\ell}^L{q}_α(N_{\ell}))$ active memory and $O(n_0n_T+ (n_T-n_0)\sum_{\ell}^L{q}_α(N_{\ell}))$ operations, where $L=\log(n_T-n_0)$, $n_0={ΔT}/τ,n_T=T/τ$, $τ$ is the stepsize, $T$ is the final time, and ${q}_α{(N_{\ell})}$ is the number of quadrature points used in the truncated Laguerre--Gauss (LG) quadrature. The error bound of the present fast method is analyzed. It is shown that the error from the truncated LG quadrature is independent of the stepsize, and can be made arbitrarily small by choosing suitable parameters that are given explicitly. Numerical examples are presented to verify the effectiveness of the current fast method.

Fanhai Zeng, Ian Turner, Kevin Burrage et al.

We propose a new class of semi-implicit methods for solving nonlinear fractional differential equations and study their stability. Several versions of our new schemes are proved to be unconditionally stable by choosing suitable parameters. Subsequently, we develop an efficient strategy to calculate the discrete convolution for the approximation of the fractional operator in the semi-implicit method and we derive an error bound of the fast convolution. The memory requirement and computational cost of the present semi-implicit methods with a fast convolution are about $O(N\log n_T)$ and $O(Nn_T\log n_T)$, respectively, where $N$ is a suitable positive integer and $n_T$ is the final number of time steps. Numerical simulations, including the solution of a system of two nonlinear fractional diffusion equations with different fractional orders in two-dimensions, are presented to verify the effectiveness of the semi-implicit methods.

Mitchell T. Griggs, Kevin Burrage, Pamela M. Burrage

This paper develops methods for numerically solving stochastic delay-differential equations (SDDEs) with multiple fixed delays that do not align with a uniform time mesh. We focus on numerical schemes of strong convergence orders $1/2$ and $1$, such as the Euler--Maruyama and Milstein schemes, respectively. Although numerical schemes for SDDEs with delays $τ_1,\ldots,τ_K$ are theoretically established, their implementations require evaluations at both present times such as $t_n$, and also at delayed times such as $t_n-τ_k$ and $t_n-τ_l-τ_k$. As a result, previous simulations of these schemes have been largely restricted to the case of divisible delays. We develop simulation techniques for the general case of indivisible delays where delayed times such as $t_n-τ_k$ are not restricted to a uniform time mesh. To achieve order of convergence (OoC) $1/2$, we implement the schemes with a fixed step size while using linear interpolation to approximate delayed scheme values. To achieve OoC $1$, we construct an augmented time mesh that includes all time points required to evaluate the schemes, which necessitates using a varying step size. We also introduce a technique to simulate delayed iterated stochastic integrals on the augmented time mesh, by extending an established method from the divisible-delays setting. We then confirm that the numerical schemes achieve their theoretical convergence orders with computational examples.

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.

Kevin Burrage, Pamela M. Burrage, Justin N. Kreikemeyer et al.

In this paper, we adapt a two-species agent-based cancer model that describes the interaction between cancer cells and healthy cells on a uniform grid to include the interaction with a third species -- namely immune cells. We run six different scenarios to explore the competition between cancer and immune cells and the initial concentration of the immune cells on cancer dynamics. We then use coupled equation learning to construct a population-based reaction model for each scenario. We show how they can be unified into a single surrogate population-based reaction model, whose underlying three coupled ordinary differential equations are much easier to analyse than the original agent-based model. As an example, by finding the single steady state of the cancer concentration, we are able to find a linear relationship between this concentration and the initial concentration of the immune cells. This then enables us to estimate suitable values for the competition and initial concentration to reduce the cancer substantially without performing additional complex and expensive simulations from an agent-based stochastic model.

Fanhai Zeng, Ian Turner, Kevin Burrage et al.

In this paper, we develop regularized discrete least squares collocation and finite volume methods for solving two-dimensional nonlinear time-dependent partial differential equations on irregular domains. The solution is approximated using tensor product cubic spline basis functions defined on a background rectangular (interpolation) mesh, which leads to high spatial accuracy and straightforward implementation, and establishes a solid base for extending the computational framework to three-dimensional problems. A semi-implicit time-stepping method is employed to transform the nonlinear partial differential equation into a linear boundary value problem. A key finding of our study is that the newly proposed mesh-free finite volume method based on circular control volumes reduces to the collocation method as the radius limits to zero. Both methods produce a large constrained least-squares problem that must be solved at each time step in the advancement of the solution. We have found that regularization yields a relatively well-conditioned system that can be solved accurately using QR factorization. An extensive numerical investigation is performed to illustrate the effectiveness of the present methods, including the application of the new method to a coupled system of time-fractional partial differential equations having different fractional indices in different (irregularly shaped) regions of the solution domain.

Megan E. Farquhar, Timothy J. Moroney, Qianqian Yang et al.

Heart failure is one of the most common causes of death in the western world. Many heart problems are linked to disturbances in cardiac electrical activity, such as wave re-entry caused by ischaemia. In terms of mathematical modelling, the monodomain equation is widely used to model electrical activity in the heart. Recently, Bueno-Orovio et al. [J. R. Soc. Interface 11: 20140352, 2014] pioneered the use of a fractional Laplacian operator in the monodomain equation to account for the complex heterogeneous structures in heart tissue. In this work we consider how to extend this approach to apply to hearts with regions of damaged tissue. This requires the use of a fractional Laplacian operator whose fractional order varies spatially. We develop efficient numerical methods capable of solving this challenging problem on domains ranging from simple one-dimensional intervals with uniform meshes, through to full three-dimensional geometries on unstructured meshes. Results are presented for several test problems in one dimension, demonstrating the effects of different fractional orders in regions of healthy and damaged tissue. Then we showcase some new results for a three-dimensional fractional monodomain equation with a Beeler-Reuter ionic current model on a rabbit heart mesh. These simulation results are found to exhibit wave re-entry behaviour, brought about only by varying the value of the fractional order in a region representing damaged tissue.

Kevin Burrage, Pamela M. Burrage, Ian W. Turner et al.

In this paper we study the class of mixed-index time fractional differential equations in which different components of the problem have different time fractional derivatives on the left hand side. We prove a theorem on the solution of the linear system of equations, which collapses to the well-known Mittag-Leffler solution in the case the indices are the same, and also generalises the solution of the so-called linear sequential class of time fractional problems. We also investigate the asymptotic stability properties of this class of problems using Laplace transforms and show how Laplace transforms can be used to write solutions as linear combinations of generalised Mittag-Leffler functions in some cases. Finally we illustrate our results with some numerical simulations.