Differential Inversion of the Implicit Euler Method: Symbolic Analysis

The implicit Euler method integrates systems of ordinary differential equations $$\frac{d x}{d t}=G(t,x(t))$$ with differentiable right-hand side $G : R \times R^n \rightarrow R^n$ from an initial state $x=x(0) \in R^n$ to a target time $t \in R$ as $x(t)=E(t,m,x)$ using an equidistant discretization of the time interval $[0,t]$ yielding $m>0$ time steps. We aim to compute the product of its inverse Jacobian $$ (E')^{-1} \equiv \left (\frac{d E}{d x}\right )^{-1} \in R^{n \times n} $$ with a given vector efficiently. We show that the differential inverse $(E')^{-1} \cdot v$ can be evaluated for given $v \in 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).$ The theoretical results are supported by actual run times. A reference implementation is provided.