Introduction

Differential Equations are used as modeling tools in several scientific disciplines. Since the advent of affordable computing technologies, a prevalent theme in these fields is to numerically approximate the solutions to these differential equations using numerical integration.

Consider a first order differential equation, formulated as an initial value problem.

\[\begin{align} \label{eq:IVP} \frac{d}{dx}y(x) = f(y(x),x) \\ y(x=0) = y_0 \in \mathbb{R}^n \end{align}\]

The central theme of numerical integration is to transform the solution of the IVP in (\ref{eq:IVP}) to the solution of the following quadrature problem.

\[\begin{equation} \label{eq:quad} y(x) = y_0 + \int_{x_0}^{x} dz\:f(z) = y_0 + \frac{1}{x-x_0} \int_{0}^{1} dq\:f(q) \end{equation}\]

Admittedly, there are several techniques to approximate the integral in (\ref{eq:quad}). Here we are interested in a class of integrators that are suitable for stiff differential equations.

Stiff Differential Equations

Stiffness of differential equations is an ambiguous term. Hairer and Wanner famouosly write about stiff differential equations, as equations where evaluation of solutions using explicit integration methods is too inefficient. There have been other attempts to formalize this “definition” using the eigenfunctions of the differential operator. Intuitively, stiff differential equations are those with \(f(y(x),x)\) evolving in disparate scales along its different bases directions.

Exponential Integrators

The question of resolving the disparate timescales in \(f(y(x),x)\) is addressed in exponential integrators as follows.

\[\begin{align} \label{eq:exp-int} \frac{d}{dx}y(x) = Ay(x) + g(y(x),x) \textrm{ [Linearization] } \end{align}\]

The solution to the homogenous problem in (\ref{eq:exp-int}), is

\[\begin{align} \label{eq:exact-sol} y(x) = \exp^{xA} y(0) \end{align}\] \[\begin{equation} \label{eq:exp} \exp^{zA} = \mathcal{I} + \frac{zA}{1!} + \frac{(zA)^2}{2!} + \frac{(zA)^3}{3!} + ... \end{equation}\]

For a generic heterogenity, \(g(x)\), the solution is expressed in terms of its propagator.

\[\begin{equation} \label{eq:hetero} y(x) = \exp^{xA} y(0) + \int dz \: \exp^{(x-z)A} g(y(z),z) \end{equation}\]