Differential Inversion of the Implicit Euler Method: Symbolic Analysis

The implicit Euler method integrates systems of ordinary differential equations $$\frac{dx}{dt}=G(t,x(t))$$ with differentiable right-hand side $G:\mathbb{R}\times {\mathbb{R}}^n\to {\mathbb{R}}^n$ from an initial state $x=x(0)\in {\mathbb{R}}^n$ to a target time $t\in \mathbb{R}$ as $x(t)=E(t,m,x)$ using an equidistant discretization of the time interval $[0,t]$ yielding $m>0$ time steps. We present a method for efficiently computing the product of its inverse Jacobian $$(E^{\mathrm{\prime }}{)}^{-1}\equiv {\left(\frac{dE}{dx}\right)}^{-1}\in {\mathbb{R}}^{n\times n}$$ with a given vector $v\in {\mathbb{R}}^n\mathrm{.}$ We show that the differential inverse $(E^{\mathrm{\prime }}{)}^{-1}\cdot v$ can be evaluated for given $v\in {\mathbb{R}}^n$ with a computational cost of $\mathcal{O}(m\cdot n^2)$ as opposed to the standard $\mathcal{O}(m\cdot n^3)$ or, naively, even $\mathcal{O}(m\cdot n^4)\mathrm{.}$ The theoretical results are supported by actual run times. A reference implementation is provided.