MARS : a Method for the Adaptive Removal of Stiffness in PDEs

The E(xplicit)I(implicit)N(null) method was developed recently to remove numerical instability from PDEs, adding and subtracting an operator $\mathcal{D}$ of arbitrary structure, treating the operator implicitly in one case, and explicitly in the other. Here we extend this idea by devising an adaptive procedure to find an optimal approximation for $\mathcal{D}$. We propose a measure of the numerical error which detects numerical instabilities across all wavelengths, and adjust each Fourier component of $\mathcal{D}$ to the smallest value such that numerical instability is suppressed. We show that for a number of nonlinear and non-local PDEs, in one and two dimensions, the spectrum of $\mathcal{D}$ adapts automatically and dynamically to the theoretical result for marginal stability. Our method thus has the same stability properties as a fully implicit method, while only requiring the computational cost comparable to an explicit solver. The adaptive implicit part is diagonal in Fourier space, and thus leads to minimal overhead compared to the explicit method.