Shev MacNamara

Arieh Iserles, Shev MacNamara

New directions in research on master equations are showcased by example. Magnus expansions, time-varying rates, and pseudospectra are highlighted. Exact eigenvalues are found and contrasted with the large errors produced by standard numerical methods in some cases. Isomerisation provides a running example and an illustrative application to chemical kinetics. We also give a brief example of the totally asymmetric exclusion process.

Gilbert Strang, Shev MacNamara

When a matrix has a banded inverse there is a remarkable formula that quickly computes that inverse, using only local information in the original matrix. This local inverse formula holds more generally, for matrices with sparsity patterns that are examples of chordal graphs or perfect eliminators. The formula has a long history going back at least as far as the completion problem for covariance matrices with missing data. Maximum entropy estimates, log-determinants, rank conditions, the Nullity Theorem and wavelets are all closely related, and the formula has found wide applications in machine learning and graphical models. We describe that local inverse and explain how it can be understood as a matrix factorization.

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.

Andreas Hellander, Jan Klosa, Per Lötstedt et al.

We analyze the governing partial differential equations of a model of pole-to-pole oscillations of the MinD protein in a bacterial cell. The sensitivity to extrinsic noise in the parameters of the model is explored. Our analysis shows that overall, the oscillations are robust to extrinsic perturbations in the sense that small perturbations in reaction coefficients result in small differences in the frequency and in the amplitude. However, a combination of analysis and simulation also reveals that the oscillations are more sensitive to some extrinsic time-scales than to others.

Katharina Kormann, Shev MacNamara

Error estimates for the numerical solution of the master equation are presented. Estimates are based on adjoint methods. We find that a good estimate can often be computed without spending computational effort on a dual problem. Estimates are applicable to both settings with time-independent, and time-dependent propensity functions. The Finite State Projection algorithm reduces the dimensionality of the problem and time propagation is based on an Arnoldi exponential integrator, which in the time-dependent setting is combined with a Magnus method. Local error estimates are devised for the truncation of both the Magnus expansion and the Krylov subspace in the Arnoldi algorithm. An issue with existing methods is that error estimates for truncation of the state space depend on measuring a loss of probability mass in a way that is not usually compatible with the approximation of the exponential. We suggest an alternative error estimate that is compatible with a Krylov approximation to the matrix exponential. Finally, we apply the new error estimates to develop an adaptive simulation algorithm. Numerical examples demonstrate the benefits of the approach.