Not Another Three-Body Problem Notebook¶

During one of my curiosity-driven deep dives, I was looking at the Wikipedia entry for the three-body problem the other day [1], and the following animation caught my eye:

The image shows a scalene triangle that tracks with the point masses, and shows that the centre of mass stays fixed for three identical bodies with zero initial velocity, acting under the gravitational influence of one another.

That wasn't what I saw, though. What I saw was an equilateral triangle (a 2-simplex) with fixed side length $R$ translating and rotating in 3D space, with one axis stretching through the screen -- that is, the projection or "shadow" of a 2-simplex moving in a higher dimension. As it moves away from us along some hypothetical $z$-axis that stretches through the screen, we see the points come closer together. Similarly, as it moves closer to us, we see those points move further away (this would imply perspective projection, as the size of the object changes with distance). Three-body collisions would occur when the simplex moves away from us towards infinity (a singularity for three point masses), and two body-collisions would occur when the simplex is on its side and rotating along an axis that is orthogonal to the $z$-axis (i.e. a degenerate triangle in 2D space).

I did a bit of reading into the subject, because I was curious, but I couldn't find anything that discussed or visualized what I was imagining (though it's possible that I was looking in the wrong places, or misunderstanding the literature I read). This got me thinking. I'm no physicist or mathematician, but I have enough background knowledge to get myself into trouble. Lagrange famously found a family of solutions where the masses stay oriented on the vertices of an equilateral triangle:

Some simpler examples (see here) can clearly be visualized as a rotating 2-simplex that is aligned with our 2D plane. I wondered if we could view all systems in such a way. That is, even for chaotic systems that don't have stable orbits, can we map them to a 2-simplex, and view its erratic spinning and translation as the chaotic oscillation of this object in 3D space? Similarly, could we also do the inverse and map periodic behaviour in some higher dimensional "simplex space" to a stable solution of some three-body system using a 2D projection?

These are some big questions, and somewhat out of my area of expertise. That said, I think if I had to distill them down to a handful of research questions, they would be:

Can we map any non-collinear configuration of 3 points to a 2-simplex of known size in a higher-dimensional space? What is the behaviour of the 2-simplex (position, velocity, and attitude) when we view known periodic solutions? Related to the above, do (some, or all) periodic solutions correspond to the projection of periodic rotations and translations of the 2-simplex about its centre of mass? Are we able to use this idea to find solutions to the three-body problem (either existing or new)?

That's a fairly ambitious set of questions, and the last one in particular could be quite involved. It's probably not realistic for me to answer them with any kind of rigor or certainty, as I don't really have the toolset to do so. In light of this, I think a realistic goal should be to get the inverse projection of a 2-simplex, look at the behaviour of a few systems, and play around with translations and rotations of the simplex in 3D space to see what we get when we project it onto to a 2D plane.

Full disclaimer: this is very much a learning exercise for myself, to flex a bit of Python, and to practice my experimental design and documentation processes. I can't claim to be an expert in this area, though I've read a few papers on the subject, and understood them as best I can. I'm not an expert on what methods are commonly used to find stable solutions, or where my ideas fit into this landscape (though I roughly understand reduced parameters, the shape sphere, and how they can be used to find solutions). I wouldn't be surprised to learn that I've gone through the process of rediscovering common knowledge within the discipline, but that's okay. I'd geniunely love to hear from people who are experts in this field on this topic, since it's fascinating.

The structure of this notebook is as follows: first, I briefly go over the background of the three-body problem, followed by a quick overview of projections and inverse projections. Then, I'll go through a methodology for obtaining the inverse projection of a given orbital system (using optimization) that I think might work. Then, I'll implement the code (it's all here) and run some experiments to see what we get. I'll conclude with some thoughts on the limitations of the idea, and ways to improve the methodology.

2.1 The Three-Body Problem¶

For the uninitiated, the three-body problem dates back to the time of Sir Isaac Newton, and deals with the dynamics of three bodies acting under the gravitational influence of one another. Whereas the two-body problem is well studied and has known analytic solutions, the three-body problem is highly chaotic, with relatively few stable orbital systems. To date, several families (in fact, 16 [2]) of stable solutions have been discovered using numerical methods, and in 1912, Karl Sundman published a general solution to the problem [1] that takes the form of a convergent series. Unfortunately, Sundman's solution -- while an incredible achievement -- is of little practical value, as it relies on an impossibly large number of terms to be useful.

So what is this problem? From Wikipedia, our model of the three-body system takes the form:

$\frac{dV_1}{dt} = -Gm_2\frac{r_1 - r_2}{|r_1 - r_2|^3}-Gm_3\frac{r_1 - r_3}{|r_1 - r_3|^3}$

$\frac{dV_2}{dt} = -Gm_1\frac{r_2 - r_1}{|r_2 - r_1|^3}-Gm_3\frac{r_2 - r_3}{|r_2 - r_3|^3}$

$\frac{dV_3}{dt} = -Gm_2\frac{r_3 - r_2}{|r_3 - r_2|^3}-Gm_1\frac{r_3 - r_1}{|r_3 - r_1|^3}$

Where $G$ is the gravitational constant, and $m_n$ denotes the mass of the $n_{th}$ body. The presence of the third body makes this problem far more complex than the case for only two bodies. Our hammer for modelling and making predictions about systems with three or more bodies is simulation using numerical integrators -- that is, for each object, we calculate the force acting on it, find its acceleration, and then update its velocity and position using a small time increment. Mathematically:

$V_{t+dt} \approx V_{t} + \frac{dV_{t}}{dt}\Delta t$

$X_{t+dt} \approx X_{t} + \frac{V_{t+dt}}{dt}\Delta t$

Where $\Delta t$ is a small, finite timestep (as opposed to the infinitesimal $dt$). We incur some error with an update over a finitely small time interval, but by making our time step small, this error will be small enough that we can unroll our system over long time horizons and maintain a reasonable degree of accuracy. In general, numerical simulation of orbital systems doesn't conserve energy (due to the aforementioned error); traditionally, this is dealt with by using so-called symplectic integrators, though in more recent works (e.g. [2]), it's common to use adaptive stepsize methods. Since Scipy doesn't include any symplectic integrators (and since the fixed stepsize would smash my poor old macbook) I'll use an adaptive method and lean on the modern miracle that is Scipy.

Finally, I'll make a quick simplification to our system of equations so that I'm not trying to wrangle my way through a second PhD thesis. I'll assume that our three bodies all have the same mass, which -- given that our gravitational constant is (exactly as its name would suggest) constant -- lets us look at the simpler system:

$\frac{dV_1}{dt} = -k\left(\frac{r_1 - r_2}{|r_1 - r_2|^3}+\frac{r_1 - r_3}{|r_1 - r_3|^3}\right)$

$\frac{dV_2}{dt} = -k\left(\frac{r_2 - r_1}{|r_2 - r_1|^3}+\frac{r_2 - r_3}{|r_2 - r_3|^3}\right)$

$\frac{dV_3}{dt} = -k\left(\frac{r_3 - r_2}{|r_3 - r_2|^3}+\frac{r_3 - r_1}{|r_3 - r_1|^3}\right)$

where I'll refer to $k$ as the stiffness. This is a fairly common simplification to make. This way, we're simulating the time evolution of some system without trying to justify what that system is. If we need to recover parameters for a physical system that might exhibit the stable orbits that we find, we can always do that later.

In order to find our hypothetical 2-simplex, we need to first understand what a projection is. We take this for granted, but we see it every day on our phones and laptops -- a projection is an operation that transforms a 3D vector or set of coordinates into a 2D representation that can be shown on a screen. Typically, this is done using the following operation:

$u_w = \begin{bmatrix} x_w \\ y_w \\ z_w \\ w \end{bmatrix} = \mathbf{A} \begin{bmatrix} x_s \\ y_s \\ z_s \\ w_s \end{bmatrix} = \mathbf{A} u_s$

$u_p = \frac{1}{w}\begin{bmatrix} x_w \\ y_w \\ z_w \end{bmatrix}$

Where $\mathbf{A}$ is a $4 \times 4$ projection matrix that encodes our object's rotation and translation (and sometimes shear) relative to the space we're projecting into (though, depending on application, it can also encode other variables, such as camera calibration). When we do a projection, we flatten one of those dimensions (in our case, we'll call it $z$), which is what lets us project it onto our screen. We lose a bit of information in the process (i.e. reversing the operation is more difficult), but that's usually not important in the types of applications where projection is used. Typically, we add a clipping value $w$ into our system, and we project into clip space first. From there, we use this clipping constant to normalize the other coordinates (this is the equivalent of projecting into a coordinate system that is parallel with the camera, and then projecting to the camera itself -- it makes the math a bit more straightforward).

What we're primarily concerned with in this problem is the opposite of a projection (the inverse projection), where we take set of 2D coordinates, and find a matrix transform that maps them to a set of 3D coordinates, which -- as we've already mentioned -- is somewhat more difficult. The reason this is difficult is that for a given set of 2D cordinates for our triangle, there are an infinite number of 2-simplexes of varying dimensions and distances that could be a valid solution to the problem. Consider, for example, that we could have two objects $A$ and $B$, where $A$ is twice the size of $B$. If we position $A$ twice as far from the screen as $B$, it will project the same shadow, making it an equally valid solution. Now consider that there are an infinite combination of sizes and distances that we could choose!

Needless to say, this is a problem, and it barely scratches the surface. For example, we also have to consider that different orientations could also cast the same shadow, and shearing transformations could give us wildly different solutions in 3D space that correspond to the same solution in 2D space. Oof.

Well, let's have a go anyway, since there are ways that we can reduce the space a bit. Since we aren't particularly interested in finding $A$ (what we actually want to find is its inverse $A^{-1}$), we can set this problem up as:

$u_{w^{'}} = \begin{bmatrix} x_{w^{'}} \\ y_{w^{'}} \\ z_{w^{'}} \\ {w^{'}} \end{bmatrix} = \mathbf{A^{-1}} \begin{bmatrix} x_p \\ y_p \\ z_p \\ w_p \end{bmatrix} = \mathbf{A^{-1}} u_p$

$u_s = \frac{1}{w^{'}}\begin{bmatrix} x_{w^{'}} \\ y_{w^{'}} \\ z_{w^{'}} \end{bmatrix}$

Where $w^{'}$ is the inverse clip space (for more info see here). Thankfully, the inverse projection that we find isn't particularly important -- that is, we don't care where in 3D space our object is oriented, because we can always choose a new datum point and translate our object. Rather, what we're interested in is that our inverse projection is consistent over time -- that is, for small changes in our three-body system, the inverse projection shouldn't map to wildly different inverse projections -- so that we can observe the behaviour of our simplex over the same time period.

When we do the opposite, and try to solve the inverse problem of rotating and translating our 2-simplex in 3D space, and and projecting it into 2D, we'll need to do a perspective projection from our 3D space onto a 2D plane. Thankfully, we don't need to find $\mathbf{A}$ for that either. Instead, we will use:

$\begin{bmatrix} x_p \\ y_p \end{bmatrix} = m\begin{bmatrix} x_s \\ y_s \end{bmatrix}$

$m = c/z_s$

for some chosen constant $c$. This will result in an arbitrary projection onto a 2D plane that we can then visualize. As with our inverse projection, the coordinates themselves shouldn't matter, because again, we're mainly interested in the evolution of the system over time. I'm not convinced that if we have what looks to be a valid solution that the coordinates will make a difference. We should be able to find "some" set of parameters (stiffness, initial velocities) to make our system work within this coordinate frame (though maybe someone with a deeper understanding can help expand on this).

Armed with this knowledge, I'll go through the methodology below.

I'm not a 3D graphics expert, so I don't know what techniques are typically used to convert 2D projected points into 3D coordinates for a known object. What I do know, however, is some machine learning. We know the 2D coordinates of our 3 "planets" in the $xy$-plane, and we know that an inverse projection takes the form given previously. What we really need to do is find a set of constraints that will let us find consistent inverse projections over long periods of time, given that we have a previous transform that we think is pretty accurate. The following are probably a good start:

All sides of the triangle have a fixed length $R=1$ (for convenience). The angle between all edges is $\pi/3$. We can pick an arbitrary projection to project our found 3D points back into 2D, which helps constrain the problem because we can map individual vertices to one another in both 2D and 3D. Our solution at each timestep is likely to be similar to the solution for the previous timestep. We can constrain our optimization routine to find answers that are "close" to whatever our previous answer was using the previous surface normal vector.

With this information, I came up with the following methology:

Guess $A^{-1}$, either using a random matrix, or $\mathbf{I}$. If we already have a solution from a previous timestep, we use that instead. If we have a previous solution, assume that the normal vector $\hat{n}$ will be close to the new solution. Add small noise $\epsilon$ to $A^{-1}$ to ensure normal vectors aren't collinear. Calculate $u_f$ using the above equations. Find all edge lengths $e$, and all angles $a$. Project $u_f$ back into 2D using our chosen projection function. Find the cost of our inverse projection matrix using $\sum_i (e_i - e_{g})^2 + (a_i - a_g)^2 + (v_{p,i} - v_{p,i}^{*})^2$, where $e_g$ is our target edge length $R$, $a_g$ is our target angle of $\pi/3$, and $v_{p,i}$ and $v_{p,i}^{*}$ are our 2D projected coordinates from our inverse projected 3D coordinates, and the ground truth, respectively. Use an optimizer to find the inverse projection matrix that minimizes this cost function.

It's not pretty, but it's straightforward, and it should hopefully work. We can also consider adding other constraints on the COM of the simplex to help things out a bit. Since our object is rotating and translating in 3D space, we need to repeat this operation at every timestep in order to get the evolution of its position and attitude. This is actually a good thing; since we're taking small timesteps, our previous solution should be a good guess for our new solution, which will speed things up considerably. Furthermore, it should help ensure that we get solutions that are consistent with the previous timestep. Once we have the 3D coordinates of our simplex's vertices, we can find a surface normal vector by taking the cross product over two of the edge vectors, and normalizing it. This way, we can characterize not only the translation of the object (measured using the position of the centre of mass), but also its attitude.

While we're at it, let's also have a think about the inverse problem: if we have a simplex in some initial state in our 3D space, and we give it some kind of periodic rotation and translation, we can do a straightforward 2D projection to get the behaviour of the vertices in 2D space. Not all of these projections will result in valid systems (I'll expand on a few obvious counterexamples further below), but we should nevertheless get some interesting visuals. Furthermore, some of these projections inarguably do correspond to valid orbital systems, but again, that warrants further discussion below. In theory, we should be able to show this by using optimization to find a set of initialization parameters that results in the system we find, but that's beyond the scope of this study (for now).

The below code imports all of the modules we'll need, and implements our 2D three-body simulation to be called with scipy.integrate.