Euler-Richardson Method

As can be seen from the proceeding discussion, the algorithm for obtaining a numerical solution of a differential equation is not unique, and there are many algorithms that reduce to the same differential equation in the limit $\Delta t \rightarrow 0$. It might occur that it would be better to compute the velocity at the middle of the interval, rather that at the beginning or at the end of the interval. The Euler-Richardson algorithm is fusion of this idea together with the simple Euler method. This algorithm is particularly useful for velocity-dependent forces, but does as well as other simple algorithms for forces that do not depend on the velocity. The algorithm consists of using the Euler method to find the intermediate position $x_{\rm mid}$ and the velocity $v_{\rm mid}$ at the time $t_{\rm mid}=t+\Delta t/2$. Then we compute the force, $F(x_{\rm mid}, v_{\rm mid}, t_{\rm mid})$, and the acceleration $a_{\rm mid}$ at $t_{\rm mid}$. The new position $x_{n+1}$ and velocity $v_{n+1}$ at time $t_{n+1}$ is found using $v_{\rm mid}$ and $a_{\rm mid}$. We summarize the Euler-Richardson algorithm as follows:
$\displaystyle a_n$ $\textstyle =$ $\displaystyle \frac{F(x_n,v_n,t_n)}{m},$ (13)
$\displaystyle v_{\rm mid}$ $\textstyle =$ $\displaystyle v_n + \frac{1}{2} a_n \Delta t,$ (14)
$\displaystyle a_{\rm mid}$ $\textstyle =$ $\displaystyle x_n + \frac{1}{2} a_n \Delta t,$ (15)
$\displaystyle a_{\rm mid}$ $\textstyle =$ $\displaystyle \frac{F(x_{\rm mid}, v_{\rm mid}, t + \frac{1}{2} \Delta t)}{m},$ (16)

so that the particle velocity and position are obtained from
$\displaystyle v_{n+1}$ $\textstyle =$ $\displaystyle v_n + a_{\rm mid} \Delta t,$ (17)
$\displaystyle x_{n+1}$ $\textstyle =$ $\displaystyle x_n + v_{\rm mid} \Delta t.$ (18)

Even though we need to do twice as many computations per time step, the Euler-Richardson algorithm is much faster because we can make the time step larger and still obtain better accuracy that with either the Euler or Euler-Cromer algorithm. Formal justification of the Euler-Richardson algorithm can be obtained, as usual, from the appropriate Taylor series expansion.