# Find a superior alternative to Backward Euler method

Currently, we are using the backward Euler (or implicit Euler) method for the solution of stiff ordinary differential equations during scientific computing.

Assuming a quite performant computer hardware and an identical step size which is smaller than 100us. Are there other stable integration methods that are able to compute y(n+1) in just one time step (real-time) and have lower truncation errors? What are their pros and cons?

I would like to implement the most promising ones and benchmark their results.

External references:

Numerical Solution of Ordinary Diﬀerential Equations

One-Step Methods: Chapter 3.3

John Butcher´s Tutorials

• Runge-Kutta is one of the first traditionally introduced at least when I was a student nl.wikipedia.org/wiki/Runge-Kuttamethode – mathreadler May 20 '17 at 19:08
• The implicit trapezoidal method is only a slight variation of the backward Euler method. Are you sure that you have to exclude it? – LutzL May 20 '17 at 19:54
• Is there a predictor-corrector scheme that meets your requirements? First-order in time is pretty restrictive; there are many second-order methods with superior accuracy and stability using very little added computation. You may also be interested in local-linearization integration. In brief: rather than approximately integrating $\frac{dx}{dt} = f(x)$, we can linearize $f(x) \approx \frac{df}{dx}\Big{|}_{x=x_{i-1}}(x - x_{i-1})$ and then apply one timestep of the analytical solution to $\frac{dx}{dt} = Fx$. – jnez71 May 20 '17 at 19:56
• By "order one in time", I assume you mean a method with first-order accuracy? If so, then by definition every first-order method will have exactly the same behavior of local truncation error, $\tau = O(\Delta t^2)$, so I'm not sure I understand the question. – Rahul May 20 '17 at 20:02
• Gear, William,Numerical Initial Value Problems in Ordunary Differential Equations – Robert Lewis May 21 '17 at 0:15

If you are currently using backwards Euler for a linear ODE $\dot y=Ay+b$, then you are solving in every step $$\Bigl(I-hA(t_n+h)\Bigr)\,y_{n+1} = y_n+hb(t_n+h)$$ If $A$ is constant, that is one matrix factorization at initialization and the corresponding backwards substitutions per step.
Using the implicit trapezoidal formula $$y_{n+1}=y_n+\frac h2(f(t_n,y_n)+f(t_n+h,y_{n+1}))$$ requires to solve for the linear equation the system $$\left(I-\frac h2 A(t_n+h)\right)y_{n+1}=\left(I+\frac h2A(t_n)\right)y_n+\frac h2\left(b(t_n)+b(t_n+h)\right)$$ Compared to the method before, this has one additional matrix-vector multiplication, which is not that much effort, especially if $A$ is sparse. And a bit more organization and storage, which is a concern at coding time and does not materially influence the run-time.
• It is order 2, it is the best you can do with using only values at $t_n$ and $t_{n+1}$, its stability region is still the negative half-plane, it is part of the widely used Crank-Nicholson method,... – LutzL May 21 '17 at 8:31