I like to classify solutions of the problem of the time evolution of the
complete initial state of a set of objects at some epoch time, where the
objects are subject to Newtonian gravitation into two main groups.
One approach is to use orbital elements of some sort.
The other is to use a numerical initial value problem solver,
aka a numerical integrator. The latter is the primary subject of this answer.
Note well: This classification isn't quite perfect as hybrid approaches are
also possible, wherein one uses a numerical integrator to integrate
time-varying orbital elements.
Keplerian orbital elements work quite nicely in the case of two point
masses, or more generally, two objects with a spherical mass distribution.
The anomaly is the only Keplerian orbital element that changes over time.
Keplerian elements can be used in situations where the underlying
assumptions are approximately correct by developing a model of how
those supposedly unchanging elements vary with time. One way of doing
this is to use Lagrange's planetary equations. (There are other related
approaches such as Gauss' planetary equations, Delaunay's planetary
equations, etc.)
Lagrange's planetary equations yield expressions for how Keplerian orbital
elements vary over time given a set of perturbing forces.
Another approach is to use something akin to those Keplerian elements
such as Delaunay elements), coupled with planetary equations for those
alternative elements. Yet another approach is to use orbital elements
(e.g., Brouwer-Lyddane elements, SGP4 elements) in which the planetary
equations are embedded in the orbital element to Cartesian state
transformation algorithm. This final approach is used to this day to
describe vehicles in Earth orbit.
The other approach is to use numerical integration.
I'll start with a discussion of how to solve for the value of a scalar
function $x(t)$ at some time $t_1$ given an initial value $x(t_0) = x_0$
and some well-behaved (continuous and bounded) derivative function
$f(x(t),t) = dx(t)/dt$ that describes the time evolution of $x$.
This falls in the very broad category of initial value problems.
Suppose the ordinary differential equation cannot be solved analytically
and cannot be expressed terms of a useful power series. This doesn't mean
nothing can be done. There are a number of techniques for solving this
problem numerically.
Note the dependence of the derivative function on the dependent variable $x$.
This becomes the much simpler problem of numerical quadrature if the
derivative function can be expressed independent of $x$.
The discussion that follows assumes that the derivative function $f$ does
indeed depend on $x$.
Note well: Newtonian gravitation falls in this category.
It also assumes the derivative function is well-behaved.
Numerically integrating across a discontinuity is a bad idea.
The foundation of the integration-based techniques for a scalar function
is the mean value theorem, which says that at some time $t_c$ between
$t_0$ and $t_1$, the value at $t_1$ is exactly
$x(t_1) = x(t_0) + (t_1-t_0)\,f(x(t_c),t_c)$. If only we could find that
magical $t_c$ and the derivative at that point. There's a chicken and
egg problem here: that magical point in time is not known. Even if it was,
the derivative function depends on state, and that too isn't known.
A very simple approach around this problem is to assume that this magical
point is the initial point:
$$x(t_1) = x(t_0) + (t_1-t_0)f(x(t_0),t0)$$
This works quite nicely for values of $t_1$ that are very close to $t_0$.
It doesn't work very well at all where $|t_1 - t_0|$ is not small.
This suggests splitting the interval $(t_0, t_1)$ into a number of smaller
intervals. This results in Euler's method: Apply the above to advance
state to time $t_0+\Delta t$, then to $t_0+2\Delta t$, and so on, eventually
reaching the desired time. Euler's method is rather lousy, even for a simple
first order scalar ODE. We can do much better than this. The key reason for
discussing Euler's method is that it is the basis for many other integration
techniques. Learn how it works, then toss it.
One approach to approving on Euler's method is to somehow correct the
result from Euler's method. For example, take an Euler step and compute
the derivative at the end point. Then use the average of those two
derivative values (the original value used to make the Euler step, and
the other from the end of the Euler step) to recompute the step from
$t$ to $t+\Delta t$. THis is Heun's method.
Another approach is to guess that the magical point $t_c$ lies
somewhere between $t$ and $t+\Delta t$. Perhaps the middle? We can use
Euler's method to advance state to the midpoint, and then use the
derivative at that point to advance state from $t$ to $t+\Delta t$.
This is the midpoint method.
Both Heun's method and the midpoint method appear to be steps backwards,
computationally. While Euler's method requires but one evaluation of the
derivative function per time step, these improved methods require two.
However, the error growth is in general so much smaller with either
Heun's method or the midpoint method compared to that from Euler's method.
This means that those "improvements" most definitely are improvements.
The expense of calling the derivative function twice per step is more than
offset by the fact that imprpovements enables take steps that are orders
of magnitude larger than one can make with Euler's method.
Both Heun's method and the midpoint method are simple improvements.
This problem has been studied in many guises. There are many more advanced
techniques. One is the class of Runge-Kutta integration techniques.
Both Heun's method and the midpoint method fall into this class.
The most popular of these, classical Runge-Kutta 4, is a significant
improvement on those two methods. There are even higher order Runge-Kutta
integrators than RK4.
Heun's method also falls into the broad class of predictor-correctors,
wherein one method (the predictor) advances state to the end of the interval
and another method (the corrector) uses the derivative at this approximate
endpoint to correct the guess made by the predictor.
The above focused solely on first order ODEs involving a scalar function.
What if the problem is multidimensional or involves higher order derivatives?
The mathematics of the techniques described above can easily handle
multidimensional data: Simply use the vector-valued time derivative.
Since a higher order ODE can be converted to a first order ODE via an
augmented vector-valued state, the same approaches used to address
multidimensional data can also be employed to address higher order ODEs.
There's a problem with doing this: It throws out geometry. For example,
consider the rather simple first order ODE $\dot x = -y, \dot y = x$.
The solution to this multivariate ODE is uniform circular motion.
Applying Euler's method to this results in
$$\begin{aligned}
x(t+\Delta t) &= x(t) - \Delta t\,y(t) \\
y(t+\Delta t) &= y(t) + \Delta t\,x(t)
\end{aligned}$$
The square magnitude of this new vector is
$(x(t)^2+y(t)^2)(1+\Delta t^2)$, which is always greater than magnitude
of the vector at the start of the step. This is not uniform circular motion.
The solution obtained via Euler's method instead spirals out. Other techniques
spiral inward. A geometric integrator on the other hand will somehow
maintain the constraint that $x^2+y^2$ is a constant of motion.
The above example showed why we don't want to toss geometry in a very
simple problem. The geometry of Newtonian gravitation, along with much
classical mechanics in general, is symplectic geometry.
This is why symplectic integrators are of great concern. A simple example
again, with Euler's method: Suppose the second derivative of position
is given by some function $\ddot {\vec x}(t) = \vec f(x(t),t)$.
Applying the basic Euler method against 3+3 dimensional phase space
dictates that
$$\begin{aligned}
\vec x(t+\Delta t) &= \vec x(t) + \Delta t \vec v(t) \\
\vec v(t+\Delta t) &= \vec v(t) + \Delta t \vec f(\vec x(t),t)
\end{aligned}$$
A simple change makes this symplectic:
$$\begin{aligned}
\vec v(t+\Delta t) &= \vec v(t) + \Delta t \vec f(\vec x(t),t) \\
\vec x(t+\Delta t) &= \vec x(t) + \Delta t \vec v(t+\Delta t)
\end{aligned}$$
As is the case with the scalar techniques discussed at the start,
Euler's method is a starting point rather than the end with
regard to symplectic integration techniques. Symplectic Euler's method
is rather lousy. But at least orbits don't spiral outward.
In making an N-body gravitation simulation, the size of N (the number of
bodies) is a key concern. Simulating a galaxy is a very different concern
from simulating a star system. The techniques used in simulating the
formation of a galaxy are very different from those used to develop a
solar system ephemeris. Galactic scale simulations cannot afford to
calculate all of the N2 gravitational interactions amongst all
the particles, and because N is so large, it cannot afford anything
more complex than very simple integrators. A star system model that does
not calculate all N2 of the gravitational interactions
or that uses very a simple integrator will by viewed in disdain.
