We consider the solution of large stiff systems of ordinary differential equations with explicit exponential Runge–Kutta integrators [2] for systems of the form

\begin{equation} \label {eq:ode_system} \frac {\partial \bm {u}(t)}{\partial t} = F(t,\bm {u}(t)) = - \bm {A} \bm {u}(t) + g(t, \bm {u}(t)), \quad \bm {u}(0) = \bm {u}_0, \end{equation}

where we view them as being semi-discretized semi-linear parabolic partial differential equations on continuous domains or on inherently discrete graph domains. This results in \(\bm {A}=\Delta \) (the continuour case) or \(\bm {A}=\bm {L}\) the graph case.

We suggest the use computing linear combinations of \(\varphi \)-functions in exponential integrators [1] to solve the above problem. State-of-the-art computational methods use polynomial Krylov subspaces of adaptive size for this task. These methods often suffer from that the required number of Krylov subspace iterations to obtain a desired tolerance increase drastically with the spectral radius of the system matrix \(\bm {A}\).

We present an approach that leverages rational Krylov subspace methods for this task (cf. [5]). The main workhorse are Runge-Kutta schemes. For these, one can now employ the scheme with \(s\) internal stages \(t+c_1 h_i, \dots , t+c_s h_i\) with \(c_j\in [0,1]\) for \(1 \leq j \leq s\) into the time interval \([t,t+h_i]\) leading to schemes of the form

\begin{align} \bm {u}_{i+1} & = \chi (-h_i\bm {A})\bm {u}_i + h_i \sum _{j=1}^s b_j(-h_i \bm {A})\bm {G}_{ij}.\label {eq:exprk1}\\ \bm {U}_{ij} & = \chi _j(-h_i\bm {A})\bm {u}_i + h_i \sum _{k=1}^s a_{jk}(-h_i \bm {A})\bm {G}_{ik},\label {eq:exprk2}\\ \bm {G}_{ik} & = g(t_i + c_k h_i, \bm {U}_{ik}),\label {eq:exprk3} \end{align} where \(\chi \), \(\chi _j\), \(a_{jk}\), and \(b_j\) are \(\varphi \)-functions (cf. [4]).

Our exponential integrator approach relies on the use of rational Krylov approximations to the matrix functions via

\begin{equation} \label {eq:fAb_rational_krylov} f(\widetilde {\bm {A}})\widetilde {\bm {c}} \approx \|\widetilde {\bm {c}}\|_2 \bm {V}_m f(\bm {H}_m\bm {K}_m^{-1})\bm {e}_1, \end{equation}

where \(\bm {H}_m\) and \(\bm {K}_m\) are small matrices coming from the rational Arnoldi method. We then give a novel a-posteriori error estimate of the rational Krylov approximation to the action of the matrix exponential on vectors for single time points. This bound allows for an adaptive approach similar to existing polynomial Krylov techniques. We then briefly discuss the selection of poles and the need for solving linear systems efficiently. The key to the convincing performance is to construct preconditioners that lead to approximately constant iteration numbers.

Numerical experiments show that our method outperforms the state of the art for sufficiently large spectral radii of the discrete linear differential operators [3]. We focus on well-known nonlinear partial differential equation models allowing the fast simulations of examples including Turing patterns. Additionally, we show that this approach allows the fast simulation of nonlinear network dynamical systems.

References

• [1]  A. H. Al-Mohy and N. J. Higham, Computing the action of the matrix exponential, with an application to exponential integrators, SIAM J. Sci. Comput., 33 (2011), pp. 488–511.

• [2]  R. Alexander, Diagonally implicit Runge–Kutta methods for stiff ODE’s, SIAM J. Numer. Anal., 14 (1977), pp. 1006–1021.

• [3]  S. M. Allen and J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metallurgica, 27 (1979), pp. 1085–1095.

• [4]  K. Bergermann and M. Stoll, Adaptive rational Krylov methods for exponential Runge–Kutta integrators, SIAM J. Matrix Anal. Appl., 45 (2023), pp. 744–770.

• [5]  S. Massei and L. Robol, Rational Krylov for Stieltjes matrix functions: Convergence and pole selection, BIT, 61 (2021), pp. 237–273.