The three-body problem

The publication of the novel “The three body problem” written by Chinese author Cixin Liu rekindled my interest about this long-standing issue in classical dynamics. As also pointed out in the novel, an empty space is stationary (and boring). An empty space plus a spherical body is also stationary (still boring). An empty space plus two bodies is a little more interesting, as the two will orbit around their center of mass. But they will do so forever. It’s only when you add a third object that things start to be very interesting. A complex dynamics arises from Newtons’ three laws of motions and the law of universal gravitational.

It was already Newton who attacked for the first time this problem, as it is relevant in the study of the motion of the Sun-Earth-Moon system. However, the problem first gain its “official” name in 1747 thanks to the work of the French mathematician Jean d’Alembert, although he’s mostly known for having found the solution to the wave equation.

If you are not familiar with differential equations, finding a solution means to discover an equation that satisfies all the constraints and gives a valid result  for every possible value of the independent variable (usually time or spatial coordinates).

The most famous result about the three-body problem is probably due to Henri Poincaré. He showed that there is no general analytical solution for the three-body problem given by algebraic expressions and integrals. The motion of three bodies is generally non-repeating, except in special cases. The most relevant is probably the case in which the three masses lay on the same plane, as may happen in a solar system. More complex solutions have been found for 3-dimensional arrangements, but the values of the masses have to appear in peculiar ratios, or the initial positions have to coincide with the vertex of a 3:4:5 triangle. As a matter of fact, the 3-body problem can be considered as a special case of the N-body problem, which has a solution found by Karl Sundman, but the convergence of the series is so slow that has no practical application.

Another possibility is actually to give up finding an analytical solution, and instead integrate numerically the differential equation. This usually requires large computational power if the equation is very complex or if the required level of accuracy must be very high. Especially when dealing with planetary orbits, it is not uncommon to run simulations that account for a time span equivalent to a few billion of years (our solar system is in fact about 5 billion years old). Small uncertainties accrue over the iterations and can potentially lead to highly divergent solutions when the calculation is repeated for a very large number of times (see also The Copernicus complex by Caleb Scarf). The basic method is called Euler Forward Method. Improvements to this calculation lead to the so-called Runge-Kutta family of methods.

I’m showing here a program I created using Processing that you can use to play around with the three-body problem. It is based on the Solar System model shown in previous posts in this blog, and performs a numerical integration based on the Runge-Kutta method. In this particular case, I defined three “suns” of mass 0.5, 1.0 and 1.5 solar masses, plus a planet with mass equal to one Earth mass. Every time you run the simulation, the initial state of the suns and the planet will be different. As you can see, in some cases the planet ends up having almost stable orbits, but most of the times an instability builds up and the planet is expelled from the system. You can see chaos theory acting in front you!

If you think this seems an uncommon scenario, think twice: observations suggest that half of the exoplanets discovered so far actually orbit around binary-systems! The scenario depicted in the novel  The three body problem may not be so far fetched after all.

import java.io.*;
import traer.physics.*;

Boolean recording = false;
String outdir = "/Users/disipio/Documents/Processing/three_body/";

float AU = 1.496e11; // m

float G_N = 6.67e-11;
float m_sun = 1.989e30;
float m_earth = 5.9e24;
float f1 = 0.5;
float f2 = 1.;
float f3 = 2;

// standard grav parameter mu=GM
float mu_solar_system = G_N * m_sun;
float mu_solar_system_reduced = mu_solar_system / m_earth;

float r_earth = 1.; //AU

ParticleSystem physics;
Particle sun1;
Particle sun2;
Particle sun3;
Particle earth; //earth

void setup()
{
 size( 1000, 1000 );
 smooth();
 
 float r_0 = ( width/20 );
 fill( 0 );
 ellipseMode( CENTER );
 
 println( "mu of the solar system " + mu_solar_system );
 println( "speed of earth at perihelion = " + sqrt(mu_solar_system/AU) + " m/s");

println(" width = " + width + " r_0 = " + r_0 );
 
 physics = new ParticleSystem( 0., 0. );
 physics.setIntegrator( ParticleSystem.RUNGE_KUTTA );
 
 float amax = 0.2;
 float amin = -0.2;
 float x01 = width/2 + width * random( amin, amax );
 float y01 = height/2 + height * random( amin, amax );
 float z01 = 0.;
 float x02 = width/2; //+ width * random( amin, amax );
 float y02 = height/2; //+ height * random( amin, amax );
 float z02 = 0.;
 float x03 = width/2 + width * random( amin, amax );
 float y03 = height/2 + height * random( amin, amax );
 float z03 = 0.;

sun1 = physics.makeParticle( f1*m_sun/m_earth, x01, y01, z01 );
 sun2 = physics.makeParticle( f2*m_sun/m_earth, x02, y02, z02 );
 sun3 = physics.makeParticle( f3*m_sun/m_earth, x03, y03, z03 );

sun1.velocity().set( 0.3*sqrt(mu_solar_system_reduced/r_0), 
 0.3*sqrt(mu_solar_system_reduced/r_0), 
 -0.3*sqrt(mu_solar_system_reduced/r_0) );
 //sun2.velocity().set( 0.5*sqrt(mu_solar_system_reduced/r_0), 
 // 0.5*sqrt(mu_solar_system_reduced/r_0), 
 // 0.);
 sun2.velocity().set( 0., 0., 0. );
 sun3.velocity().set( -0.5*sqrt(mu_solar_system_reduced/r_0), 
 -0.5*sqrt(mu_solar_system_reduced/r_0), 
 0.5*sqrt(mu_solar_system_reduced/r_0) );
 
 earth = physics.makeParticle( 1, width/2 + r_0*r_earth, height/2, 0 );
 earth.velocity().set( 0, sqrt(mu_solar_system_reduced/r_0), 0. );
 
 println( "Start:" );
 
 physics.makeAttraction( sun1, earth, G_N, r_0/200 );
 physics.makeAttraction( sun2, earth, G_N, r_0/200 );
 physics.makeAttraction( sun3, earth, G_N, r_0/200 );
 
 physics.makeAttraction( sun1, sun2, G_N, r_0/200 );
 physics.makeAttraction( sun1, sun3, G_N, r_0/200 );
 physics.makeAttraction( sun2, sun3, G_N, r_0/200 );

}

void mousePressed()
{
 
}

void mouseDragged()
{
 
}

void mouseReleased()
{
 
}

void keyPressed() { 
 // If we press r, start or stop recording!
 if (key == 'r' || key == 'R') {
 recording = !recording;
 }
}

void draw()
{
 physics.tick(1000);
 
 background( 0 );
 
 color y = color( 255, 204, 0);
 fill(y);
 noStroke();
 ellipse( sun1.position().x(), sun1.position().y(), f1*12, f1*12 );
 ellipse( sun2.position().x(), sun2.position().y(), f2*12, f2*12 );
 ellipse( sun3.position().x(), sun3.position().y(), f3*12, f3*12 );
 
 color b = color( 0, 0, 255);
 fill(b);
 ellipse( earth.position().x(), earth.position().y(), 5, 5 );
 
 if (recording) {
 saveFrame( outdir + "frames/frames####.png");
 }
}

Advertisement