Problem and background: Many problems in science and engineering involve multiple physical processes, where complex interactions between these components can result in dynamics evolving on different time scales (known as stiff systems). They arise from a wide range of fields such as fluid dynamics in meteorology, electrochemistry, combustion, and plasma physics. The evolution of such problems can often be described by parabolic or hyperbolic partial differential equations (PDEs). For instance, parabolic PDEs can be cast in the general form
$\dfrac{\partial u}{\partial t }=- \sum_{i, j= 1}^{d}\frac{\partial }{\partial x_i } \left(a_{ij}(\mathbf{x})\frac{\partial u}{\partial x_j} \right)+ \sum_{i= 1}^{d} b_i(\mathbf{x})\frac{\partial u}{\partial x_i } + c(\mathbf{x}) u+ f(\mathbf{x},t, u, \nabla u), \quad u(\mathbf{x},t_0)=u_0 (\mathbf{x}) $
on some open and bounded domain $\Omega \subset \mathbb{R}^d$ (with smooth boundary $\partial \Omega$ and closure $\overline{\Omega}$), subject to some boundary condition $(B u)(\mathbf{x},t)=g(\mathbf{x},t)$ for all $(\mathbf{x},t) \in \partial\Omega \times (t_0, T)$ ($B$ is some boundary operator e.g. Dirichlet, Neumann, or mixed). Here, $u=u(\mathbf{x},t)$ for $(\mathbf{x},t) \in \Omega \times (t_0, T)$, the coefficient matrix $[a_{ij}(\mathbf{x})]_{i,j=1}^{d}$ is symmetric and strictly positive definite, $a_{ij}(\mathbf{x}) \in C^1 (\overline{\Omega})$, $b_i (\mathbf{x}), c (\mathbf{x}) \in C(\overline{\Omega})$.
Under appropriate regularity conditions on the solution $u$, the function $f$, and suitable initial and boundary conditions, the existence and uniqueness results for the solution have been established both in the classical and weak sense (a well-posed problem), see for instance the book by Taylor (2010).
Typical examples are reaction-diffusion systems describing chemical kinetic mechanisms involved in ignition, which can create thousands of reactions occurring over a wide range of time scales. In magnetohydrodynamic (MHD) equations (describing large scale plasma behavior), stiffness arises from the presence of a wide range of waves encoded in the complex nonlinear terms. Other prominent examples include the heat equation, and the incompressible Navier–Stokes equations, the nonlinear Schrödinger equation, etc.
Numerical approach:
A common way to solve such PDEs numerically is to discretize in space (e.g. using finite difference/element/volume methods or spectral methods), yielding a large system of stiff ODEs of the form
$U'(t) =F(t, U(t)), \quad U(t_0) =U_0$.
Solving this stiff system is challenging since the Jacobian often possesses eigenvalues with large negative real parts or is even an unbounded operator. Classical integrators such as explicit Runge–Kutta methods usually lack stability and are required to use extremely small time steps due to the CFL condition. Thus standard integrators such as implicit Runge–Kutta methods, e.g., Gauss/Radau/BDF/Rosenbrock/W-methods are usually used. Specifically, implicit-explicit (IMEX) or semi-implicit methods applied to the case where $F(t, U)$ can be partitioned into linear and nonlinear parts, i.e.
$U'(t) =LU(t)+N(t,U(t))$
have been widely used with some of the earlier applications coming from fluid dynamics. The downside of these methods, however, is that a nonlinear/linear system of equations has to be solved in each step, which may require a good preconditioner to speed up convergence of an iterative method used in the implicit solver (e.g. Newton-Krylov methods). Depending on the complexity of the differential operator, constructing an efficient preconditioner can be difficult. Moreover, as the stiffness increases such methods may require increased computational effort. Therefore, development of an efficient time integrator that enables simulation in the regimes of interest requires much effort and care.