Exponential Integrators

This module provides exponential integrators for ODEs.

Exponential integrators are a large class of methods from numerical analysis is based on the exact integration of the linear part of the initial value problem. Because the linear part is integrated exactly, this can help to mitigate the stiffness of a differential equation.

We consider initial value problems of the form,

\[u'(t)=f(u(t)),\qquad u(t_{0})=u_{0},\]

which can be decomposed of

\[u'(t)=Lu(t)+N(u(t)),\qquad u(t_{0})=u_{0},\]

where \(L={\frac {\partial f}{\partial u}}\) (the Jacobian of f) is composed of linear terms, and \(N=f(u)-Lu\) is composed of the non-linear terms.

This procedure enjoys the advantage, in each step, that \({\frac {\partial N_{n}}{\partial u}}(u_{n})=0\). This considerably simplifies the derivation of the order conditions and improves the stability when integrating the nonlinearity \(N(u(t))\).

Exact integration of this problem from time 0 to a later time \(t\) can be performed using matrix exponentials to define an integral equation for the exact solution:

\[u(t)=e^{Lt}u_{0}+\int _{0}^{t}e^{L(t-\tau )}N\left(t+\tau, u\left(\tau \right)\right)\,d\tau .\]

This representation of the exact solution is also called as variation-of-constant formula. In the case of \(N\equiv 0\), this formulation is the exact solution to the linear differential equation.

Exponential Rosenbrock methods

Exponential Rosenbrock methods were shown to be very efficient in solving large systems of stiff ODEs. Applying the variation-of-constants formula gives the exact solution at time \(t_{n+1}\) with the numerical solution \(u_n\) as

(1)\[u(t_{n+1})=e^{h_{n}L}u(t_{n})+\int _{0}^{h_{n}}e^{(h_{n}-\tau )L}N(t_n+\tau, u(t_{n}+\tau ))d\tau .\]

where \(h_n=t_{n+1}-t_n\).

The idea now is to approximate the integral in (1) by some quadrature rule with nodes \(c_{i}\) and weights \(b_{i}(h_{n}L)\) (\(1\leq i\leq s\)). This yields the following class of s-stage explicit exponential Rosenbrock methods:

\[\begin{split}\begin{align} U_{ni}=&e^{c_{i}h_{n}L}u_n+h_{n}\sum_{j=1}^{i-1}a_{ij}(h_{n}L)N(U_{nj}), \\ u_{n+1}=&e^{h_{n}L}u_n+h_{n}\sum_{i=1}^{s}b_{i}(h_{n}L)N(U_{ni}) \end{align}\end{split}\]

where \(U_{ni}\approx u(t_{n}+c_{i}h_{n})\).

The coefficients \(a_{ij}(z),b_{i}(z)\) are usually chosen as linear combinations of the entire functions \(\varphi _{k}(c_{i}z),\varphi _{k}(z)\), respectively, where

\[\begin{split}\begin{align} \varphi _{k}(z)=&\int _{0}^{1}e^{(1-\theta )z}{\frac {\theta ^{k-1}}{(k-1)!}}d\theta ,\quad k\geq 1, \\ \varphi _{0}(z)=&e^{z},\\ \varphi _{k+1}(z)=&{\frac {\varphi_{k}(z)-\varphi _{k}(0)}{z}},\ k\geq 0. \end{align}\end{split}\]

By introducing the difference \(D_{ni}=N(U_{ni})-N(u_{n})\), they can be reformulated in a more efficient way for implementation as

\[\begin{split}\begin{align} U_{ni}=&u_{n}+c_{i}h_{n}\varphi _{1}(c_{i}h_{n}L)f(u_{n})+h_{n}\sum _{j=2}^{i-1}a_{ij}(h_{n}L)D_{nj}, \\ u_{n+1}=&u_{n}+h_{n}\varphi _{1}(h_{n}L)f(u_{n})+h_{n}\sum _{i=2}^{s}b_{i}(h_{n}L)D_{ni}. \end{align}\end{split}\]

where \(\varphi_{1}(z)=\frac{e^z-1}{z}\).

In order to implement this scheme with adaptive step size, one can consider, for the purpose of local error estimation, the following embedded methods

\[{\bar {u}}_{n+1}=u_{n}+h_{n}\varphi _{1}(h_{n}L)f(u_{n})+h_{n}\sum _{i=2}^{s}{\bar {b}}_{i}(h_{n}L)D_{ni},\]

which use the same stages \(U_{ni}\) but with weights \({\bar {b}}_{i}\).

For convenience, the coefficients of the explicit exponential Rosenbrock methods together with their embedded methods can be represented by using the so-called reduced Butcher tableau as follows:

\[\begin{split}\begin{array}{c|ccccc} c_{2} & & & & & \\ c_{3} & a_{32} & & & & \\ \vdots & \vdots & & \ddots & & \\ c_{s} & a_{s 2} & a_{s 3} & \cdots & a_{s, s-1} \\ \hline & b_{2} & b_{3} & \cdots & b_{s-1} & b_{s} \\ & \bar{b}_{2} & \bar{b}_{3} & \cdots & \bar{b}_{s-1} & \bar{b}_{s} \end{array}\end{split}\]



Hochbruck, M., & Ostermann, A. (2010). Exponential integrators. Acta Numerica, 19, 209-286.

ExponentialEuler(f[, var_type, dt, name, ...])

The exponential Euler method for ODEs.

class brainpy.integrators.ode.exponential.ExponentialEuler(f, var_type=None, dt=None, name=None, show_code=False)[source]

The exponential Euler method for ODEs.

The simplest exponential Rosenbrock method is the exponential Rosenbrock–Euler scheme, which has order 2.

For an ODE equation of the form

\[u^{\prime}=f(u), \quad u(0)=u_{0}\]

its schema is given by

\[u_{n+1}= u_{n}+h \varphi(hL) f (u_{n})\]

where \(L=f^{\prime}(u_{n})\) and \(\varphi(z)=\frac{e^{z}-1}{z}\).

For a linear ODE system: \(u^{\prime} = Ay + B\), the above equation is equal to \(u_{n+1}= u_{n}e^{hA}-B/A(1-e^{hA})\), which is the exact solution for this ODE system.

  • f (function) – The derivative function.

  • dt (optional, float) – The numerical precision.

  • var_type (optional, str) – The variable type.

  • show_code (bool) – Whether show the code.