Why is this numerical method more accurate and stable?

Question

Suppose I have the current position, velocity and acceleration of a particle. As I see it, the best estimate of its velocity and position at a time $\Delta t$ from now should be $$\vec{v}(t+\Delta t) = \vec{v}(t) + \vec{a}\Delta t$$ $$\vec{s}(t+\Delta t) = \vec{s}(t) + \vec{v}(t)\Delta t + \frac{1}{2} \vec{a} \Delta t^2$$ However, by doing some simulations, I've noticed that an extremely more accurate method would be to consider $$\vec{v}(t+\Delta t) = \vec{v}(t) + \vec{a}\Delta t$$ $$\vec{s}(t+\Delta t) = \vec{s}(t) + \vec{v}(t)\Delta t + \vec{a} \Delta t^2$$ Which is almost the same, except by the missing $1/2$ factor. The position equation is also equivalent to $$\vec{s}(t+\Delta t) = \vec{s}(t) + \vec{v}(t+\Delta t)\Delta t$$ To show how better it is, here is the absolute deviation to the analytical solution with both methods for a simple harmonic oscillator

So the error or deviation of the second method stays almost the same while the first one increases much faster, more specifically it seems the error of the first method increases quadratically and the second linearly. I've also got this result by simulating a gyroscope and a wave propagating through a string.

What explains this counterintuitive result?

Edit

To add more details to the question, here is the implementation for the oscillator example:

• The position and velocity are initially equal to 2 constants x0 and v0 defined in the code
• The acceleration is given by $-\omega^2x$, according to the harmonic oscillator differential equation, where $x$ is the current position and $\omega$ a constant that is related to the mass and elastic constant.
• Each iteration of method 1 is given by

dv = acc(x)*dt
x += dt*(v + 1/2*dv)
v += dv

• Each iteration of method 2 is given by

v += acc(x)*dt
x += v*dt

Answer 1

In order to model the motion and acceleration of a particle you have to use following equations

$$ \begin{align} \text{(i)} \quad \frac{dx}{dt} &= f(t,v),\\ \text{(ii)} \quad \frac{dv}{dt} &= g(t,x) \stackrel{const}{\equiv} a. \end{align} $$

A correct symplectic Euler discretization would result in

$$ \begin{align} \text{(i)} \quad x_{n+1} &= x_n + \Delta t \cdot f(t_n, v_n) ,\\[0.3ex] \text{(ii)} \quad v_{n+1} &= v_n + \Delta t \cdot g(t_n, x_{n+1}) , \end{align} $$

or

$$ \begin{align} \text{(ii)} \quad v_{n+1} &= v_n + \Delta t \cdot g(t_n, x_n) ,\\[0.3em] \text{(i)} \quad x_{n+1} &= x_n + \Delta t \cdot f(t_n, v_{n+1}). \end{align} $$

which is similar to your second method. An alternative would be a Verlet integration method. Your first discretization seems to be wrong or is at least a mix of different discretization methods, as mentioned in the comments.

For more informations read:

• Symplectic integrator
• Symplectic Euler method
• Verlet integration

Edit: (see comments)

Starting form (ii) we can write

$$ \begin{align} \frac{dv}{dt} & = a,\\[0.5em] dv & = a\cdot dt ,\\[0.5em] \int_{v_0}^v dv & =\int_{t_0}^t a \cdot dt , \\[0.5em] \left( v- v_0 \right) & =a \left(t - t_0 \right) , \\[0.5em] (1) \quad \rightarrow \quad v &= a \left(t - t_0 \right) + v_0 , \quad | \quad \text{insert (i)} \\[0.5em] \frac{dx}{dt} &= a \left(t - t_0 \right) + v_0 , \\[0.5em] dx &= \left[ a \left(t - t_0 \right) + v_0\right]dt , \\[0.5em] \int_{x_0}^x dx & = \int_{t_0}^t \left[ a \left(t - t_0 \right)+ v_0\right]dt , \\[0.5em] x -x_0 & = \left[\frac{1}{2} a \left(t - t_0 \right)^2 - a\cdot t_0\left(t - t_0 \right)+ v_0\left(t - t_0 \right) \right] , \\[0.5em] (2) \quad \rightarrow \quad x & = \frac{1}{2} a \left(t - t_0 \right)^2 - a\cdot t_0\left(t - t_0 \right)+ v_0\left(t - t_0 \right) + x_0. \\[0.5em] \end{align} $$

Using your notation

$$ \begin{align} (1) \quad v (t+\Delta t) &= v(t) + a \Delta t \\[0.5em] (2) \quad x (t+\Delta t) & =x(t) + \frac{1}{2} a \Delta t^2 - a \cdot t \Delta t + v(t)\Delta t . \end{align} $$

For $~t_0=0~$ we get

$$ \begin{align} v &= a t + v_0 , \\[0.5em] x & = \frac{1}{2} a t^2 + v_0 t+x_0 .\Box \end{align} $$

This also holds for $t \rightarrow t + \Delta t$.

Answer 2

There is a variant to implement the velocity Verlet method that looks quite similar to method one.

The "half-step" variant of its implementation can look like

x(t+dt/2) = x(t) + v(t)*dt/2
a(t+dt/2) = accel(x(t+dt/2))
v(t+dt/2) = v(t) + a(t+dt/2)*dt/2
v(t+dt) = v(t+dt/2) + a(t+dt/2)dt/2
x(t+dt) = x(t+dt/2) + v(t+dt)dt/2
        = x(t) + v(t)dt/2 + (v(t)+a(t+dt/2)dt)dt/2
        = x(t) + v(t)dt + a(t+dt/2)*dt^2/2

One can now eliminate the middle step as far as possible to get

a = accel(x+0.5*v*dt)
x = x + v*dt + a*dt^2/2
v = v + a*dt

Note that the acceleration is still evaluated at the halfway point of the step. This variant is now second order and should be more precise than the method 2, especially when the energy errors are compared.