The simulation of a planet orbits around the Sun is a very good practice to understand the solution to partial differential equation.

Starting equation

The balance between gravity and centrifugal force determines the orbit of a planet circling around the Sun. Here the most simple but yet useful Newton’s gravitational law in the 2D vector form is used.

$m\frac{d^2 \vec{r}}{dt^2} = \vec{F} = -G\frac{Mm\vec{r}}{r^3} = \left\{ \begin{array}{l} -G\frac{M m x}{r^3}\vec{e_x}\\ -G\frac{M m y}{r^3}\vec{e_y}\\ \end{array} \right.$

Where $r=\sqrt{x^2 + y^2}$. We further get the momentum equation as:

$\left\{ \begin{array}{l} \frac{dx}{dt} = V_x \\ \frac{dy}{dt} = V_y \\ \frac{dV_x}{dt} = -G\frac{M x}{r^3} \\ \frac{dV_y}{dt} = -G\frac{M y}{r^3} \end{array} \right.$

Further we discrete the above equations with a forward difference form.

$\left\{ \begin{array}{l} x2 - x1 = V_{x1} \Delta t \\ y2 - y1 = V_{y1} \Delta t \\ V_{x2} - V_{x1} = -G\frac{M x1}{r^3}\Delta t \\ V_{y2} - V_{y1} = -G\frac{M y1}{r^3}\Delta t \end{array} \right.$

Using the above discrete form, we should be able to use loop to finish the simulation of planet orbit. Detailed code in python see here:

Orbit simulation in python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

# solve the planet orbit equations
# ------ import necessary modules
import pylab
import numpy
import sys
sys.path.append('../../')

from function_lib import plot_lib


pylab.close('all')
# ------ define difference form of orbit equations
def orbit_diff(x, y, vx, vy, dt):
    GM = 1
    r = numpy.sqrt(x*x+y*y)
    x1 = x + vx*dt
    y1 = y + vy*dt
    vx1 = vx - GM*x/(r**3)*dt
    vy1 = vy - GM*y/(r**3)*dt
    return x1, y1, vx1, vy1


# ------ make a time evolution calculation
# set initial values for these data
t0 = 3
dt = 0.0001
vx0 = 0
vy0 = 1.8*0.7
x0 = 1
y0 = 0

x = x0
y = y0
vx = vx0
vy = vy0
i = 0
N = 1000000
x_array = numpy.linspace(0, 0, N)
y_array = numpy.linspace(0, 0, N)
t = numpy.linspace(0, 0, N)

for i in range(0, N) :
    x_array[i] = x
    y_array[i] = y
    x, y, vx, vy = orbit_diff(x, y, vx, vy, dt)
    t[i] = t0 + i*dt
    i = i + 1
    

t_array = numpy.linspace(0, N - 1, N)*dt


# ------ plot planet orbit out
# pylab.close()
pylab.figure()
ax = pylab.subplot(111)
pylab.plot(x_array, y_array, '-', color = 'green', markersize = 3)
# mark the position of sun
pylab.plot(x0, y0, 'o', markersize = 10)
pylab.plot(0, 0, 'o r', markersize = 30)
pylab.legend(['orbit', 'planet', 'Sun'])
pylab.xlabel('x(A.U)')
pylab.ylabel('y(A.U)')
ymax = numpy.max(numpy.abs(y_array))*1.1
xmax = numpy.max(numpy.abs(x_array))*1.1
pylab.xlim([-xmax, xmax])
pylab.ylim([-ymax, ymax])
plot_lib.single_plot_paras()
ax.set_aspect('equal', 'box')
pylab.show()

Simulation result see this figure

Attention: the choice of dt apparently affect the accuracy of the orbit simulation. Small dt results in good accuracy but increases the calculation time.

Orbit of many planets

In a real star system, there might be mnay planets orbiting around a main star. Thus we should consider the realistic condition of many planets and even the interaction between the planets should be considered.

2 planets orbiting

Here we consider the interaction between two planets but still fix the main star in the system. The forces act on planet 1 and planet 2 are like this:

$\begin{array}{l} m_1\frac{d^2 \vec{r}}{dt^2} = -G\frac{M m_1\vec{r}}{r^3} -G\frac{m_2 m_1\vec{r_{12}}}{r_{12}^3}\\ m_2\frac{d^2 \vec{R}}{dt^2} = -G\frac{M m_2\vec{R}}{R^3} -G\frac{m_1 m_2\vec{r_{21}}}{r_{21}^3} \end{array}$

Then after discretion, we get this form:

$\begin{array}{l} r = \sqrt{x_1^2 + y_1^2}\\ R = \sqrt{X_1^2 + Y_1^2}\\ r12 = \sqrt{(x_1 - X_1)^2 + (y_1 - Y_1)^2}\\ r21 = r12\\ vx2 = vx1 - GM\frac{x1}{r^3} - G m_2\frac{x1 - X1}{r_{12}^3}\Delta t\\ vy2 = vy1 - GM\frac{y1}{r^3} - G m_2\frac{y1 - Y1}{r_{12}^3}\Delta t\\ VX2 = VX2 - GM\frac{X1}{R^3} - G m_1\frac{X1 - x1}{r_{21}^3}\Delta t\\ VY2 = VY2 - GM\frac{Y1}{R^3} - G m_1\frac{Y1 - y1}{r_{21}^3}\Delta t\\ \end{array}$

Then we translate these equations into python codes as:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

def orbit_diff_2(x, y, vx, vy, X, Y, VX, VY, dt) :
    GM = 1
    Gm1 = 10**-2
    Gm2 = 3*10**-2
    r = numpy.sqrt(x*x + y*y)
    R = numpy.sqrt(X*X + Y*Y)
    r12 = numpy.sqrt((x-X)**2 + (y - Y)**2)
    r21 = r12
    x2 = x + vx*dt
    y2 = y + vy*dt
    X2 = X + VX*dt
    Y2 = Y + VY*dt
    vx2 = vx - GM*x/(r**3)*dt - Gm2*(x - X)/(r12**3)*dt
    vy2 = vy - GM*y/(r**3)*dt - Gm2*(y - Y)/(r12**3)*dt
    VX2 = VX - GM*X/(R**3)*dt - Gm1*(X - x)/(r21**3)*dt
    VY2 = VY - GM*Y/(R**3)*dt - Gm1*(Y - y)/(r21**3)*dt
    return x2, y2, vx2, vy2, X2, Y2, VX2, VY2
# set start position and speed of planet 1
x = x0
y = y0
vx = vx0
vy = vy0
i = 0
N = 1000000
x_array = numpy.linspace(0, 0, N)
y_array = numpy.linspace(0, 0, N)
t = numpy.linspace(0, 0, N)
# set start position and speed of planet 2
X0 = 2
Y0 = 0
X = X0
Y = Y0
VX = 0
VY = 0.9
X_array = numpy.linspace(0, 0, N)
Y_array = numpy.linspace(0, 0, N)
# perform loop simulation
for i in range(0, N) :
    x_array[i] = x
    y_array[i] = y
    X_array[i] = X
    Y_array[i] = Y
    x, y, vx, vy, X, Y, VX, VY = orbit_diff_2(x, y, vx, vy, X, Y, VX, VY, dt)
    i = i + 1
t_array = numpy.linspace(0, N - 1, N)*dt
# ------ plot planets orbit out
# pylab.close()
pylab.figure(figsize = (10, 5))
ax = pylab.subplot(111)
pylab.plot(x_array, y_array, '-', color = 'green')
pylab.hold('on')
pylab.plot(X_array, Y_array,'--', color = 'black')
# mark the position of sun
pylab.hold('on')
pylab.plot(x0, y0, 'o', color = 'green', markersize = 10)
pylab.hold('on')
pylab.plot(X0, Y0, 'o', color = 'black', markersize = 15)
pylab.hold('on')
pylab.plot(0, 0, 'o r', markersize = 30)
pylab.legend(['orbit 1', 'orbit 2', 'planet 1', 'planet 2', 'Sun'], frameon = 'false')
pylab.xlabel('x(A.U)')
pylab.ylabel('y(A.U)')
ymax1 = numpy.max(numpy.abs(y_array))*1.1
xmax1 = numpy.max(numpy.abs(x_array))*1.1
ymax2 = numpy.max(numpy.abs(Y_array))*1.1
xmax2 = numpy.max(numpy.abs(X_array))*1.1
Ymax = numpy.max([ymax1, ymax2])
Xmax = numpy.max([xmax1, xmax2])
ylab.xlim([-Xmax, Xmax])
pylab.ylim([-Ymax, Ymax])
plot_lib.single_plot_paras()
ax.set_aspect('equal', 'box')
pylab.tight_layout()
pylab.show()

Simulation results shows that for planet mass much smaller than the central star, the interaction between planets rally change their orbits.

For the case when the planet mass is close to the mass of the central star, the interaction between planets will disturb their orbits.

Simulation of many planets

Under this condition, the interaction between planets will be replaced by a gravity potential.

Newton’s Gravitational Force and Kapler’s Law of Planet Orbit

The above works are just numerical simulations, but the physics behind it is not understand clearly. Eg, why following the Newton’s law, the simulated planet orbit should be an elliptic? The simulation results only tells you the resulted orbit is an elliptic but it does not tell you why the orbit is an elliptic. To solve this puzzle, analytic analysis of orbit equation should be conducted.

Gravitational Force in the polar coordinate

$\frac{d^2\vec{r}}{dt^2} = -GM\frac{\vec{r}}{r^3}$

Then we apply the polar coordinate with $\vec{r} = r\vec{e_r}$, $\frac{d}{dt}\vec{e_r} = \frac{d\theta}{dt}\vec{e_{\theta}} = \omega\vec{e_{\theta}}$ and $\frac{d}{dt}\vec{e_{\theta}} = -\omega\vec{e_r}$. With these rules, the vector form of Gravitation force can be simplified as:

$(\frac{d^2 r}{dt^2} - r\omega^2)\vec{e_r} + (2\frac{dr}{dt}\omega + r\frac{d\omega}{dt})\vec{e_{\theta}} = -GM\frac{1}{r^2}\vec{e_r}$

The above vector equation can be split in two equations. One in the radial direction, one in poloidal direction.

$\left. \begin{array}{l l} \frac{d^2 r}{dt^2} - r\omega^2 = -GM\frac{1}{r^2} && (a)\\ 2\frac{dr}{dt}\omega + r\frac{d\omega}{dt} = 0 && (b) \end{array} \right\} (I)$

Proof of Kapler’s Second Law

From the above equation(b) we can get:
$2\frac{dr}{dt}\omega + r\frac{d\omega}{dt} = 0 \Rightarrow \omega(\omega {dr}{dt} + 2r\frac{d\omega}{dt}) = \frac{d}{dt}(r^2\omega) = 0$, which means $r^2\omega = constant$. This should be the proven of Kapler’s second law, the sweep area by the line linking the planet and the Sun is the same in equal time interval, where $\frac{dS}{dt} = \frac{d(1/2 r^2\omega)}{dt} = 0 \Rightarrow dS/dt = 0$.

{width=200px}

Proof of Kapler’s First Law

Start from equation I(a), $\frac{d^2 r}{dt^2} - r\omega^2 + GM\frac{1}{r^2} = 0$. Using the second law of Kapler $r^2\omega = C$ and replace $\omega$ with $r$, we can get $\frac{d^2 r}{dt^2} = C\frac{1}{r^3} - GM\frac{1}{r^2}$.

If we perform a smart transform as: $\frac{d^2 r}{dt^2} = \frac{dr}{dt}\frac{d}{dr}(\frac{dr}{dt}) = 1/2 \frac{d}{dr}((\frac{dr}{dt})^2)$ , then equation I(a) can be rewritten as: $1/2 \frac{d}{dr}((\frac{dr}{dt})^2) = C\frac{1}{r^3} - GM\frac{1}{r^2}$.

Perform an integration over r we can further get: $(\frac{dr}{dt})^2 = 2GM\frac{1}{r} - C\frac{1}{r^2} + D$.

Then we introduce r as a function of $\theta$ as $r = r(\theta)$ and get: $\frac{dr}{dt} = \frac{dr}{d\theta}\frac{d\theta}{dt} = \omega\frac{dr}{d\theta}$ and consider Kapler’s second law, the previous equation can be rewritten as: $(\frac{dr}{d\theta})^2\frac{C^2}{r^4} = 2GM\frac{1}{r} - C\frac{1}{r^2} + D$.

Then we introduce a transform as $u = 1/r \Rightarrow du = \frac{-1}{r^2}dr$ and the previous equation can be rewrite in this form: $(\frac{du}{d\theta})^2 = A + Bu - Cu^2 \Rightarrow \frac{du}{d\theta} = \pm \sqrt{A + Bu - Cu^2}$. It can be easily found that $\frac{du}{d\theta} < 0$, which means $\frac{du}{d\theta} = -\sqrt{A + Bu - Cu^2}$.

Now we get $\frac{-du}{\sqrt{A + Bu - Cu^2}}=d\theta \Rightarrow \int\frac{-du}{\sqrt{A + Bu - Cu^2}}=\int d\theta$. In the case $\sqrt{A + Bu -Cu^2} = E\sqrt{1 - (\frac{h - u}{E})^2}$. Then the function can be written as:
$\int \frac{d(\frac{h-u}{E})}{\sqrt{1 - (\frac{h - u}{E})^2}} = \int d\theta \Rightarrow arccos(\frac{h-u}{E}) = \theta + H$, considering initial condition we can get $H = 0$. Further perform a cos operation on the whole equation we can get:

$\frac{h - u}{E} = cos(\theta)$, then we replace u with $1/r$ and get:
$r=\frac{1}{h - E cos(\theta)}$ which is just the equation of elliptic in polar coordinate.

Proof of Kapler’s Third Law

Start from equation I(a), $\frac{d^2 r}{dt^2} - r\omega^2 = -GM\frac{1}{r^2}$, we take an extreme condition by take the point near the longer leg of the orbit of a planet, where the change on radius is zero $\frac{dr}{dt} = 0$, which will simplify the above equation as: $- r\omega^2 = -GM\frac{1}{r^2}$, then we can get $r_1^3\omega_1^2 = r_2^3\omega_2^2 = GM \Rightarrow \frac{R_1^3}{T_1^2} = \frac{R_2^3}{T_2^2} = constant$, which is the proof of Kapler’s third law. Be careful that this is only an approximation of the real condition.

{width=200px}