Category: programming

You are sitting in your chair and reading this because of:

$\mathbf{F}_{ij} = \frac{Gm_im_j(q_j - q_i)}{||q_j - q_i||^3}$

I am writing this blog entry for my astronomy term project (not saying which one because the professor my be reading this, but after I finish the course I will remove this note and say a few good words).

Generations of mathematicians have been frustrated because they failed to find an analytic solution to the $n$-body ($n \ge 3$) problem represented by the above equation.

The $n$-body problem is easy to solve by using numerical methods. In order to use numerical integration methods one has to convert the original problem to an Ordinary Differential Equation (ODE). Here is how Newton's law of gravity looks for $n$ bodies as an ODE:

$\frac{d^2x_i(t)}{dt^2} = G\sum_{k = 1, k \ne i}^{n}{\frac{m_k(x_k(t) - x_i(t))}{||x_k(t) - x_i(t)||^3}}$

This is a second order ODE. To solve it with Python we first have to convert the second order ODE to a system of first-order ODEs. Let's substitute $\frac{dx}{dt}$ with $p$ and $x$ with $q$. We now get:

$\frac{dq}{dt} = p$

$\frac{dp}{dt} = G\sum_{k = 1, k \ne i}^{n}{\frac{m_k(q_k(t) - q_i(t))}{||q_k(t) - q_i(t)||^3}}$

Let us first solve the two-body problem in two dimensions. This, of course, can be done analytically but we want to go through the steps of writing and evolving a Python program. So, here are the four equations:

$\frac{dq_1}{dt} = p_1$

$\frac{dq_2}{dt} = p_2$

$\frac{dp_1}{dt} = \frac{Gm_1(q_2 - q_1)}{||q_2 - q_1||}$

$\frac{dp_2}{dt} = \frac{Gm_2(q_1 - q_2)}{||q_1 - q_2||}$

And this is the Python code that can solve the equations:

#!/usr/bin/env python

import numpy
import scipy
import sys
import matplotlib.pyplot

import scipy.integrate

G = 6.674e-11
m1 = 1e2
m2 = 1e2

def func(Y, t):    
    p1 = numpy.array([ Y[0], Y[1] ])
    p2 = numpy.array([ Y[2], Y[3] ])
    q1 = numpy.array([ Y[4], Y[5] ])
    q2 = numpy.array([ Y[6], Y[7] ])

    p1_prime = G * m2 * (q2 - q1) / \
               numpy.sum((q2 - q1)**2)
    p2_prime = G * m1 * (q1 - q2) / \
               numpy.sum((q1 - q2)**2)
    q1_prime = p1
    q2_prime = p2

    return [ p1_prime[0], p1_prime[1],
             p2_prime[0], p2_prime[1],
             q1_prime[0], q1_prime[1],
             q2_prime[0], q2_prime[1] ]

def main():
    t = numpy.arange(0, 20000.0, 1)
    init_p1 = [ 0, 1 ]
    init_p2 = [ 0, -1 ]
    init_q1 = [ 0, 0 ]
    init_q2 = [ 10, 0 ]
    iv = init_p1 + init_p2 + init_q1 + init_q2
    result = scipy.integrate.odeint(func, iv, t)

    matplotlib.pyplot.scatter(result[:, 4],
                              result[:, 5],
                              s = 1,
                              lw = 0,
                              c = 'blue')
    matplotlib.pyplot.scatter(result[:, 6],
                              result[:, 7],
                              s = 1,
                              lw = 0,
                              c = 'red')
    matplotlib.pyplot.show()

if __name__ == '__main__':
    main()

So what do we set for the masses of the objects ($m_1$ and $m_2$) and for the initial position and velocity. Let's try with very small masses ($m_1 = m_2 = 10$) and the two bodies going in different directions with the same speed. This is the result:

So, given the initial vector, the above plot looks correct. The two bodies have small mass and as they keep flying in the direction in which they were flying initially. The next check is to see if when we increase the mass of the two bodies they will start attracting each other:

Now we are going to change the initial velocity of one of the two bodies and this is the result:

Finally, we can play some more and get this really complicated looking motion:

Why is all this interesting? Simulation is a very important tool in astronomy and the $n$-body problem served as a fundamental building block of today's understanding of the formation of the Solar system. Of course, Newton's formulas work only for small speeds and we have to add relativity into the picture if we want a more accurate prediction. But before updating our equations we want to generalize this introductory example first into 2-bodies in three dimensions and then to $n$ bodies in three dimensions. This is going to happen in our next post.