The Forward Euler Method in the n-body problem

Copyright 2004 S.Edgeworth.
This is version 3, published on 12 Dec 2004. (Version 1 was published 06 Dec 2004).

Abstract:

The Explicit Forward Euler method, applied to a gravitational n-body system with no velocity-dependent forces, is when properly implemented, a 2nd order method and is time-reversible.

Full paper:

We examine the Explicit Forward Euler method as used to calculate trajectories in a gravitational n-body system with no velocity-dependent forces.

Because the names Explicit Euler and Forward Euler are sometimes used inconsistently in the literature, I make it clear that by Explicit Forward Euler I mean the method defined in algorithm A.

Algorithm A: The usual implementation of the Explicit Forward Euler method can be written as

t=0
start loop
a=f(t,p)
v=v+a*h
p=p+v*h
t=t+h
output(t,p)
end loop

where h=timestep,
p=position vector,
v=velocity vector,
a=gravitational acceleration vector calculated by function f at a certain p and t,
output() denotes the place in the loop where output is performed.

Algorithm A is 1st order, and non-reversible.

In the usual implementation (algorithm A), the instantaneous velocity change applied to a body at time t is commonly described in the literature as an approximation of the real-life acceleration over the period from t to (t+h). Hence the name "Forward Euler".

It is more logical to treat the instantaneous velocity change applied to a body at time t as an approximation of the real-life acceleration over the period (t-h/2) to (t+h/2). The name "Forward Euler" is inappropriate. The Euler velocity, before this instantaneous velocity change, is really the average velocity for the period from (t-1) to t. The Euler velocity, after this instantaneous velocity change, is really the average velocity for the period from t to (t+1).

The externally specified initial velocity is the instantaneous velocity at time t=0. But the loop in algorithm A clearly needs to start with the average velocity over the period from -h to 0. So before the start of the loop we must apply the correction v=v-a*h/2. And to obtain the final velocity, after the end of the loop we need to apply the correction v=v+a*h/2.

Algorithm B: The first proposed implementation of the Explicit Forward Euler method is

t=0
a=f(t,p)
v=v-a*h/2
start loop
a=f(t,p)
v=v+a*h
p=p+v*h
t=t+h
output(t,p)
end loop
a=f(t,p)
v=v+a*h/2
output(t,v)

where
h=timestep,
p=position vector,
v=velocity vector,
a=gravitational acceleration vector calculated by function f at a certain p and t,
output() denotes the place in the loop where output is performed.

Algorithm B is 2nd order in p and final v, and time-reversible. Note that the loop is identical to algorithm A. We have simply added an intial and a final velocity correction outside of the loop. When coded efficiently it requires only one evaluation of the gravitational acceleration per step. This implementation is sufficient for applications which require p but not v as output. Those applications can omit the last three lines. The last three lines were included to enable testing for 2nd order in final v, and for time-reversibility. The correct treatment of the initial velocity has been previously stated by other authors.

An empirical test was conducted in which an elliptical orbit was simulated using algorithm B. The test was repeated using various timesteps. Each time the time-step was halved, the errors in p and in final v (versus the analytically known solution) decreased to approximately a quarter, confirming that algorithm B is 2nd-order in p and final v.

Algorithm C: If we also calculate the instantaneous velocity at every evaluated position, we obtain the second proposed implementation of the Explicit Forward Euler method

t=0
p=p0
v=v0
for i = 1 to q+1
{
   a=f(t,p)
   if(i>1)
   {
      v=w+a*h/2
      output(t,p,v)
   }
   if(i<q+1)
   {
      w=v+a*h/2
      p=p+w*h
      t=t+h
   }
}

where t=time
h=timestep,
i=loop counter
q=quantity of steps
p0=initial position vector,
v0=initial velocity vector,
p=position vector,
v=velocity vector at evaluated position,
w=average velocity vector between 2 adjacent evaluated positions,
a=gravitational acceleration vector calculated by function f at a certain p and t,
output() denotes the place in the loop where output is performed.
(Note that the loop is executed q+1 times).

Algorithm C is 2nd order in p and v, and time-reversible. It requires only one evaluation of the gravitational acceleration per step. I believe that it may also be obtained by correct implementation of the Verlet or leapfrog method. The equivalent of this algorithm has been previously stated by other authors.