Wider contours and adaptive contours
For researchers using contour integral methods for matrix functions in applications with non-normal operators (e.g., finance, biology), this provides practical contour selection guidelines to improve accuracy.
The paper addresses the challenge of computing matrix functions for non-normal matrices using contour integrals, where standard contours can cause accuracy issues due to pseudospectra. It proposes using wider contours and a semi-analytic adaptive contour based on a parabolic bound from the field of values, demonstrating improved accuracy on biological reaction models.
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.