Stochastic Picard-Runge-Kutta Solvers for Large Systems of Autonomous Ordinary Differential Equations

Abstract

We present an explicit stochastic scheme for approximating solutions of autonomous systems ordinary differential equations based on the direct simulation of paths of appropriate Markov jump processes. Its features make it efficient especially for large systems with a sparse incidence matrix, i.e. typically for spatially discretized partial differential equations. The full path of the Markov jump process is simulated by a very simple method and serves as a predictor for improved approximations. One possibility to obtain more precise values is to employ a Picard–iteration over small time intervals with length h. This amounts to the computation of the integral of a step–function, which can be done explicitly. Alternatively, or as a further step, one can apply a Runge–Kutta principle. For this one needs approximations of the integrand at some intermediate equidistant points in the time discretization interval. As suitable values can be taken either those given by the above described predictor, or their improvement by Picard iterations. The corresponding integral is then approximated by a quadrature formula from the Newton–Cotes family. For example we can use Simpson quadrature formulae (1/3 rule for the step h/2 and 3/8 rule for the step h/3). These approaches mimic respectively the well–known Runge–Kutta schemes of order 3, and the 3/8 scheme of order 4. The main difference to the classical Runge–Kutta method is that the intermediate values are computed here by performing a stochastic simulation, possibly followed by a Picard iteration. The described approach has as consequence an increased order of convergence. We also introduce time-adaptive versions of the stochastic Runge–Kutta method and illustrate the features of all considered schemes at a standard benchmark problem, a reaction–diffusion equation modeling a combustion process in one space dimension.