Performance Evaluation of Single-step Explicit Exponential Integration Methods on Stiff Ordinary Differential Equations

TL;DR

The study benchmarks explicit single-step exponential integrators against classical implicit schemes for stiff ODEs, revealing that explicit methods lag in accuracy and scalability. By decomposing dy/dt = Ly + N(t,y) and employing the matrix exponential $e^{Lh}$, the authors assess a range of IF and ETD schemes, as well as EPI-like and SSP variants. Across three stiff models, integrating factor Euler emerges as the only reliably accurate and inexpensive option among exponentials, with higher-order exponential methods failing to surpass first-order performance and often matching or underperforming backward Euler or Radau methods. The results suggest limited practical viability for explicit exponential methods in stiff contexts and motivate the development of more robust solvers or targeted use of IF Euler in neural ODE and Bayesian inference workflows.

Abstract

Stiff systems of ordinary differential equations (ODEs) arise in a wide range of scientific and engineering disciplines and are traditionally solved using implicit integration methods due to their stability and efficiency. However, these methods are computationally expensive, particularly for applications requiring repeated integration, such as parameter estimation, Bayesian inference, neural ODEs, physics-informed neural networks, and MeshGraphNets. Explicit exponential integration methods have been proposed as a potential alternative, leveraging the matrix exponential to address stiffness without requiring nonlinear solvers. This study evaluates several state-of-the-art explicit single-step exponential schemes against classical implicit methods on benchmark stiff ODE problems, analyzing their accuracy, stability, and scalability with step size. Despite their initial appeal, our results reveal that explicit exponential methods significantly lag behind implicit schemes in accuracy and scalability for stiff ODEs. The backward Euler method consistently outperformed higher-order exponential methods in accuracy at small step sizes, with none surpassing the accuracy of the first-order integrating factor Euler method. Exponential methods fail to improve upon first-order accuracy, revealing the integrating factor Euler method as the only reliable choice for repeated, inexpensive integration in applications such as neural ODEs and parameter estimation. This study exposes the limitations of explicit exponential methods and calls for the development of improved algorithms.

Figures (12)

• Figure 1: Comparison of the integration of the stiff van der Pol oscillator with $\mu=1000$ using two different integration schemes: explicit Runge-Kutta-Fehlberg, which is slow with 422,442 time points and 2,956,574 function evaluations, and implicit Radau IIA 5th order, which is faster with only 857 time points and 7,123 function evaluations.
• Figure 2: The stiff Van der Pol model. We have set $\mu=1000$, creating an oscillatory system with high stiffness. The model has sharp jumps in y, followed by stiff areas where the solution evolves more slowly.
• Figure 3: Absolute error for various integration schemes on the stiff Van der Pol model. The step size is also shown for reference. The data shown correspond to $n=1555$ data points across the integration time interval.
• Figure 4: Absolute error for various integration schemes on the stiff Van der Pol model. The step size is also shown for reference. The data shown correspond to $n=24849$ data points across the integration time interval.
• Figure 5: The “High Irradiance RESponse” (HIRES) model