Fdzfgy

I've done a bit of research, and have learned that computers "solve" the three-body-problem by using "Numerical methods for ordinary differential equations", but I can't really find anything about it other then Wikipedia. Does anyone have any good sources about this topic that isn't any kind of Wikipedia?

My thoughts:

Currently I'm using simulations of three bodies flying around each other, using Newton's gravitational law, and at a random time in the simulation, everything goes chaotic. I though that this was the only way to kind of "solve" it, but how does this "Numerical methods for ordinary differential equations" method work? And what does the computer actually do?

Numerical analysis is used to calculate approximations to things: the value of a function at a certain point, where a root of an equation is, or the solutions to a set of differential equations. It is a huge and important topic since in practice most real problems in mathematics, science and technology will not have an explicit closed-form solution. In general there are trade-offs between accuracy and computational speed.

For the three-body problem we have three point masses in starting locations $mathbf{x}_i(0)$ with velocities $mathbf{v}_i(0)$ that we want to calculate for later times $t$. Mathematically we want to find the solution to the system $$mathbf{x}'_i(t)=mathbf{v}_i(t),$$ $$mathbf{v}'_i(t)=mathbf{f}_i(t)/m_i,$$ $$mathbf{f}_i(t)=sum_{jneq i} Gm_i m_j/||mathbf{x}_i-mathbf{x}_j||.$$

The obvious method is to think "if we move forward a tiny step $h$ in time, we can approximate everything to be linear", so we make a formula where we calculate the state at time $t+h$ from the state at time $t$ (and so on for $t+2h$ and onwards): $$mathbf{x}_i(t+h)=mathbf{x}_i(t)+hmathbf{v}_i(t),$$ $$mathbf{v}_i(t+h)=mathbf{v}_i(t)+hmathbf{f}_i(t).$$ This is called Euler's method. It is simple but tends to be inaccurate; the error per step is $approx O(h^2)$ and they tend to build up. If you try it for a two body problem it will make the orbiting masses perform a precessing rosette orbit because of the error build-up, especially when they get close to each other.

There is a menagerie of methods for solving ODEs numerically. One can use higher order methods that sample the functions in more points and hence approximate them better. There are implicit methods that instead of trying to find a state at a later time only based on the current state look for a self-consistent late and intermediate state. As I said, this is a big topic.

However, for mechanical simulations you may want to look at methods designed to preserve energy and other conserved quantities (symplectic methods - these are the ones used by professionals for long-run orbit calculations). Perhaps the simplest is the semi-implicit Euler method. There is also the Verlet method and leapfrog integration. I like the semi-implicit Euler method because it is super-simple (but being a first order-method it is still not terribly accurate): $$mathbf{v}_i(t+h)=mathbf{v}_i(t)+hmathbf{f}_i(t),$$ $$mathbf{x}_i(t+h)=mathbf{x}_i(t)+hmathbf{v}_i(t+h).$$ Do you see the difference? You calculate the updated velocity first, and then use it to update the positions - a tiny trick, but suddenly 2-body orbits are well behaved.

The three body problem is chaotic in a true mathematical sense. We know there are situations where tiny differences in initial conditions get scaled up to arbitrarily large differences in later positions (even if we rule out super-close passes between masses). So even with an arbitrarily fine numerical precision there will be a point in time where our calculated orbits will be totally wrong. The overall "style" of trajectory may still be correct.

The "solution" of the three-body problem can be written as the pair of differential equations,
begin{align}
vec{v}&=frac{mathrm dvec{x}}{mathrm dt}\
vec{a}&=frac{mathrm dvec{v}}{mathrm dt}
end{align}
where the latter is usually written in terms of the force, $mvec{a}=vec{F}$. Then using the definition of the derivative,
$$
frac{mathrm df}{mathrm dx}=lim_{hto0}frac{f(x+h)-f(x)}{h},
$$
the differential equations we started with can be written as the finite difference equations,
begin{align}
vec x(t+mathrm dt)&simeq vec x(t)+vec v,mathrm dt \
vec v(t+mathrm dt)&simeq vec v(t)+vec a,mathrm dt
end{align}
which is done assuming that $mathrm dt$ is "just small" rather than infinitesimal.

It is these equations that a computer uses to "solve" the three body problem: given an initial condition for each of the bodies, the computer iteratively steps forward in time using the above finite difference equations with the force being the $n$-body force.

As mentioned in the comments, the simplistic model above is called Euler integration which is not at all well-suited for this problem because it does not conserve energy. A better option is called velocity Verlet, which is disucssed in other posts on physics.SE (I would recommend this post of mine as it gives some decent, but brief, details of the implementation).