# 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

## 1 Answer

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.

• Thank you very much for this. Could you please explain a little the benefits of the implicit trapezoidal method compared to the implicit Euler method? – MilkMan May 21 '17 at 8:28
• 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
• Ok. Now we are getting closer... I have added two more external references. There, it is mentioned that implicit multi-stage Runge-Kutta methods (and if I understood correctly, implicit Euler and implicit trapezoidal are belonging to that class of methods) could be used for my task. – MilkMan May 21 '17 at 9:15
• You might want to explore exponential methods (see google/wikipedia) and DIRK, see math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/IRK.pdf – LutzL May 21 '17 at 10:31
• It should not. If you observe that it would be worth opening a new question with a small example and code for both methods that demonstrates that behavior. But it might also be that you observe the usual imprecisions of first order methods, did you try varying the step-size? – LutzL May 22 '17 at 19:56