Paper detail

Wider contours and adaptive contours

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.