High-Order Schemes of Exponential Time Differencing for Stiff Systems with Nondiagonal Linear Part
Evelina V. Permyakova, Denis S. Goldobin
TL;DR
The authors address the challenge of applying exponential time differencing (ETD) to stiff systems with a non-diagonal linear part $\boldsymbol{L}$ by computing ETD coefficients numerically through auxiliary problems, and by rewriting Runge–Kutta ETD schemes in terms of $\mathbf{Q}(\tau)=e^{\boldsymbol{L}\tau}$ and $\mathbf{M}_n(\tau)$. They implement this approach for three test models: the 1D Cahn–Hilliard equation, the Matthews–Cox sixth-order pattern-formation equation with a conservation law, and a Fokker–Planck description of a population of quadratic integrate-and-fire neurons, showing substantial accuracy and, importantly, large speedups (often 2–3 orders) relative to standard predictor–corrector or RK-based methods. The results demonstrate that ETD schemes, particularly ETD4RK, are highly effective for PDEs with high-order spatial derivatives and for long-time simulations in complex, stiff systems, while the auxiliary-problem-based coefficient evaluation remains broadly reusable across problems. The work provides a practical, scalable framework with memory- and accuracy-conscious optimizations, and it sets the stage for extending ETD methods to higher dimensions and more complex heterogeneous models.
Abstract
Exponential time differencing methods is a power tool for high-performance numerical simulation of computationally challenging problems in condensed matter physics, fluid dynamics, chemical and biological physics, where mathematical models often possess fast oscillating or decaying modes -- in other words, are stiff systems. Practical implementation of these methods for the systems with nondiagonal linear part of equations is exacerbated by infeasibility of an analytical calculation of the exponential of a nondiagonal linear operator; in this case, the coefficients of the exponential time differencing scheme cannot be calculated analytically. We suggest an approach, where these coefficients are numerically calculated with auxiliary problems. We rewrite the high-order Runge--Kutta type schemes in terms of the solutions to these auxiliary problems and practically examine the accuracy and computational performance of these methods for a heterogeneous Cahn--Hilliard equation, a sixth-order spatial derivative equation governing pattern formation in the presence of an additional conservation law, and a Fokker--Planck equation governing macroscopic dynamics of a network of neurons.