mjhjmtu

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 (and even if they have, it might not be possible to compute with infinite precision - computers after all represent numbers with finite precision). In general there are trade-offs between accuracy and computational speed.

For the three-body problem we have three point masses in starting locations $mathbfx_i(0)$ with velocities $mathbfv_i(0)$ that we want to calculate for later times $t$. Mathematically we want to find the solution to the system $$mathbfx'_i(t)=mathbfv_i(t),$$ $$mathbfv'_i(t)=mathbff_i(t)/m_i,$$ $$mathbff_i(t)=Gm_i sum_jneq i fracm_j(mathbfx_j-mathbfx_i).$$

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): $$mathbfx_i(t+h)=mathbfx_i(t)+hmathbfv_i(t),$$ $$mathbfv_i(t+h)=mathbfv_i(t)+hmathbff_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. Most serious methods for solving ODEs will also reduce the step-size $h$ when the forces grow big during close encounters to ensure that accuracy remains acceptable. As I said, this is a big topic.

However, for mechanical simulations you may in particular 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): $$mathbfv_i(t+h)=mathbfv_i(t)+hmathbff_i(t),$$ $$mathbfx_i(t+h)=mathbfx_i(t)+hmathbfv_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, which is why it is OK to play around with semi-implicit Euler as long as one is not planning any space mission based on the results.

The "solution" of the three-body problem can be written as the pair of differential equations,
beginalign
vecv&=fracmathrm dvecxmathrm dt\
veca&=fracmathrm dvecvmathrm dt
endalign
where the latter is usually written in terms of the force, $mveca=vecF$. Then using the definition of the derivative,
$$
fracmathrm dfmathrm dx=lim_hto0fracf(x+h)-f(x)h,
$$
the differential equations we started with can be written as the finite difference equations,
beginalign
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
endalign
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).

A reasonably accessible paper on the three body problem, including a discussion on the numerical algorithms, is Musielak & Quarles (2014). For a larger discussion of more-useful higher-order techniques, as suggested by homocomputeris, see Farrés et al (2012)

A few years ago I put together a model of the inner solar system, the sun, Earth's moon, and the planets from Mercury to Mars.

I used a Leapfrog variation of Rk4 and got decent results, orbital periods consistent with initial configurations that stayed stable over time. Even some of the bad things I got I think suggested good numerical practices.

I'd initialize my planets with initial momentum, with the sun having zero initial momentum. The momentum of the planets didn't cancel out, so there was this weird net drift in roughly a straight line direction it took a while to figure out.

In addition to sorting that out, one feature I found important was the order of applying changes to position and speed.

You don't want to iterate through your bodies, calculate the net force in one pass body by body, and then apply your dynamic changes body by body. Otherwise it will, say, update the position of Venus, while keeping the position of Mars unchanged while you are updating the position of Earth. Your time steps get out of sync.

Store somewhere in memory the new displacements and velocities as you iterate through the bodies, and only after finding your dynamics at the present tick, run through and apply your changes. Apologies and congrats if you caught that on your first try. It tripped me up for a bit.

One way to simulate three bodies is using time step increments, in which you pretend that each body accelerates at a constant rate for a brief period of time. You start with the initial position vectors, velocity vectors, and the masses of the bodies, and then use them to calculate the acceleration vector of each body. In Newtonian Physics the gravitational acceleration vector of Body A from Body B is given by the equation $$veca_GA=-fracGleft(vecA-vecBright)M_B$$ with $G$ being the Gravitational Constant, $vecA$ being the position vector of Body A, $vecB$ being the position vector of Body B, $M_B$ being the Mass of Body B, and $veca_GA$ being the Gravitational Acceleration Vector of Body A. In order to find the total gravitational acceleration vector of Body A you add up each of the gravitational acceleration vectors of body A from each of the other bodies.

You can find the position vectors and velocity vectors at the end of time step $n$ using the general formulas $$vecx_n=frac12veca_n-1s^2+vecv_n-1s+vecx_n-1$$ $$vecv_n=veca_n-1s+vecv_n-1$$ with $veca_n-1$ being the acceleration vector for time step $n-1$, $vecv_n-1$ being the velocity vector for time step $n-1$, $vecx_n-1$ being the position vector $n-1$, $s$ being the length of the time step increments, $vecx_n$ being the position vector for time step $n$, and $vecv_n$ being the velocity vector for time step $n$. At each time step you also recalculate the acceleration vectors as they are not constant.

This method of simulating $n$ bodies works for any force law, meaning that you can use it for simulating $n$ bodies for a force law $F=f(r)$, with $f(r)$ being some function of distance, although in the case of the force law $F=f(r)$ the equation I specified previously for calculating the acceleration of Body A from body B becomes $$veca_GA=-fracvecA-vecBGM_Bfleft(||vecA-vecB||right)$$ but the other formulas I mentioned, including the part about adding the acceleration vectors remain the same.

It turns out that knowing the pairwise trajectories one can determine the equation of motion of bodies. The equations of motion are written as
$$fracd^2 vec r^kds^2=-Gsum_n=1,nne k^N fracm_n(vec r^k-vec r^n)tag1$$
We solve the auxiliary problem of pair interaction of bodies with reduced inertial mass
$$fracm_n m_ksum_k=1^N m_k fracd^2 vec R^knds^2=-G fracm_n m_kvec R^kn^3 tag2$$
These equations differ in different initial conditions. Subtract from equation (2) equation (1), we get
$$fracd^2 R_0^kds^2-fracd^2 r^kds^2= Gsum_n=1,nne k^N fracm_n(vec r^k-vec r^n)- Gsum_n=1,nne k^N fracm_nvec R^knvec R^kn=0,vec R^kn=vec r^k-vec r^n$$
Where $vec R_0^k=sum_n=1 n ne k^N fracm_n vec R^knsum_p=1^N m_p$
Get the formula
$$fracd^2 vec R_0^kds^2=fracd^2 vec r^kds^2$$
Integrating this equality, we obtain the equation of motion for each of the N bodies.
$$vec r_0^k=fracdvec r^k(0)dss+vec r^k(0)+vec R_0^k(s)- fracdvec R^k(0)dss-vec R^k(0)= vec R_0^k(s)$$
The movement of each body is determined by the movement of the center of inertia of the pair body system. Initial conditions are determined from the condition $$fracdvec r^k(0)dss+vec r^k(0)= fracdvec R^k(0)dss+vec R^k(0)tag3$$
We select the coordinate system $sum_n=1^N m_n vec r^n(0)=0,sum_n=1^N m_n fracd vec r^n(0)ds=0$ and then condition (3) is satisfied identically. Will get $ vec R^kn(0)=vec r^k(0)- vec r^n(0), fracd vec R^knds(0)=fracd vec r^kds(0)-fracd vec r^nds(0)$ The energy of the pair trajectory of the system are calculated by the formula $$E_kn=fracm_k m_nsum_q=1^N m_q (fracd vec R^knds(0))^2/2-fracG m_k m_n+fracM_kn^2sum_q=1^N m_qR^kn=fracm_k m_nsum_q=1^N m_q (fracd vec R^knds(0))^2/2-fracG m_k m_n+fracm_k m_n^2(fracd vec R^knds(0)times vec R^kn(0))^2$$ The moment of the pulse is determined by the formula vector product $$M_kn=fracm_k m_nsum_q=1^N m_q fracd vec R^knds(0)times vec R^kn(0)$$
In this case, the constructed solution is substituted into the calculated trajectories $vec r_0^k-vec r_0^p $ of the planets and used $vec R^kp=r^k-r^p$, then we get $vec R^kp=vec r_0^k-vec r_0^p$. Indeed
$$ r_0^k-r_0^p =sum_n=1^N fracm_n vec R^knsum_q=1^N m_q-sum_n=1^N fracm_n vec R^pnsum_q=1^N m_q=sum_n=1^N fracm_n (vec r^k-vec r^n)sum_q=1^N m_q-sum_n=1^N fracm_n (vec r^p-vec r^n)sum_q=1^N m_q=vec r^k-vec r^p=R_kp$$ Simplified Formula
$$vec r_0^u=sum_n=1^N fracm_n vec R^un(s)sum_q=1^N m_q fracp_un1+e_uncos(varphi_un-varphi_un0)$$ The angle is calculated by the formula $varphi_un-varphi_un0=int_0^s fracM_un^2du$ The magnitude $$fracR^un(varphi)=fracvec R^un(0)vec R^un(0)cosvarphi+fracvec R^un(0)times vec M_unsinvarphi $$ of the angle $varphi$ common to all pair trajectories corresponds to the direction $fracR^un(varphi)$ with single orths $$fracvec R^un(0)vec R^un(0),fracvec R^un(0)times vec M_un $$. We have the final formula for the trajectory of bodies $$vec r_0^u(varphi)=sum_n=1^N fracm_nsum_q=1^N m_q[fracvec R^un(0)vec R^un(0)cosvarphi+fracvec R^un(0)times vec M_unsinvarphi]fracp_un1+e_uncosvarphi$$ But for the lack of an infinite solution even for a short time, it is necessary to assume $$|vec R^un(varphi)|=fracp_un1+e_uncosvarphi>0,e_un<1$$ This imposes a restriction on the value of the parameters; the energy of pair interaction must be negative. Thus, bodies of small mass are excluded, whose contribution to the energy of the system is not significant. According to the formula of pair interaction, the radius varies according to the law see Landau Lifshits volume 1
$$r=a(ecoshxi-1);t=sqrtfraca^3Gsum_q=1^N m_q(esinhxi-xi),rapprox csqrtfracr_gat$$
After simple transformations, this formula is reduced to $$r=sqrtV^2(0)+frac(vec V(0)times vec R(0))^2R(0)^2-frac2Gsum_q=1^N m_qR(0)t$$ This linear increase in radius shifts the center of gravity of the system. A particle moving to infinity leads to the whole system being carried along and the center of gravity shifts. The second cosmic velocity corresponds to a positive value of the radicand, or positive energy of the system.