The motion of the sun and the earth is an example of a “two-body problem”. This is a relatively simple problem that can be solved analytically (it interesting to notice here, however, that adding a new object to the problem, the moon for instance, make it completely intractable. This is the famous “three body problem”). We can assume that, to a good approximation, the sun is stationary and is a convenient origin of our coordinate system. This is equivalent to changing to a “center of mass” coordinate system, where most of the mass is concentrated in the sun. The problem can be reduced to an equivalent one body problem involving an object of reduced mass $\mu$ given by $$\mu=\frac{mM}{m+M}$$ Since the mass of the earth is $m=5.99\times 10^{24}$ kg and the mass of the sun is $M=1.99\times 10^{30}$ kg we find that for most practical purposes, the reduced mass of the earth-sun system is that of the earth. Hence, in the following we are going to consider the problem of a single particle of mass $m$ moving about a fixed center of force, which we take as the origin of the coordinate system. The gravitational force on the particle m is given by $${\mathbf F}=-\frac{GMm}{r^3}{\mathbf r},$$ where the vector ${\mathbf r}$ is directed from $M$ to $m$, and $G$ is the gravitation constant $$G=6.67\times 10^{-11} \frac{m^3}{kg.s^2}$$ The negative sign implies that the gravitational force is attractive, and decreases with the separation $r$. The gravitational force is a “central force”: its magnitude depends on the separation between the particles and its direction is along the line that connects them. The assumption is that the motion is confined to the $xy$ plane. The angular momentum ${\mathbf L}$ lies on the third direction $z$ and is a constant of motion, *i.e.* it is conserved: $$L_z=({\mathbf r}\times m{\mathbf v})_z=m(xv_y-yv_x)=\mathrm{const.}$$ An additional constant of motion ins the total energy $E$ given by $$E=\frac{1}{2}mv^2-\frac{GmM}{r}$$ If we fix the coordinate system in the sun, the equation of motion is $$m\frac{d^2{\mathbf r}}{dt^2}=-\frac{mMG}{r^3}{\mathbf r}$$ For computational purposes it is convenient to write it down in cartesian components: $$\begin{aligned} && F_x=-\frac{GMm}{r^2}\cos{\theta}=-\frac{GMm}{r^3}x, \\ && F_y=-\frac{GMm}{r^2}\sin{\theta}=-\frac{GMm}{r^3}y.\end{aligned}$$ Hence, the equations of motions in cartesian coordinates are: $$\begin{aligned} && \frac{d^2x}{dt}=-\frac{GM}{r^3}x, \\ && \frac{d^2y}{dt}=-\frac{GM}{r^3}y, \end{aligned}$$ where $r^2=x^2+y^2$. These are coupled differential equations, since each differential equation contains both $x$ and $y$.
Since many planetary orbits are nearly circular, it is useful to obtain the condition for a circular orbit. In this case, the magnitude of the acceleration ${\mathbf a}$ is related to the radius by $$a=\frac{v^2}{r}$$ where $v$ is the speed of the object. The acceleration is always directed toward the center. Hence $$\frac{mv^2}{r}=\frac{mMG}{r^2}$$ or $$v=(\frac{MG}{r})^{1/2}. $$ This is a general condition for the circular orbit. We can also find the dependence of the period $T$ on the radius of a circular orbit. Using the relation $$T=\frac{2\pi r}{v},$$ we obtain $$T^2=\frac{4\pi^2r^3}{GM}$$
An ellipse has two foci $F_1$ and $F_2$, and has the property that for any point the distance $F_1P+F_2P$ is a constant. It also has a horizontal semi-axis $a$ and a vertical $b$. It is common in astronomy to characterize an orbit by its “eccentricity” $e$, given by the ratio of the distance between the foci, and the length of the major axis $2a$. Since $F_1P+F_2P=2a$, it is easy to show that (consider a point $P$ at $x=0$,$y=b$) $$e=\sqrt{1-\frac{b^2}{a^2}},$$ with $0<e<1$. A special case is $a=b$ for which the ellipse reduces to a circle and $e=0$. The earth orbit has eccentricity $e=0.0167$.
It is useful to choose a system of units where the product $GM$ is of the order of unity. To describe the earth’s motion, the convention is to choose the earth’s semi-major axis as the unit of length, called “astronomical unit” (AU) and is $$1AU=1.496 \times 10^{11}m.$$ The unit of time is taken to be “one year”, or $3.15 \times 10^7$s. In these units, $T=1$yr, $a=1AU$, and we can write $$GM=\frac{4\pi ^2a^3}{T^2}=4\pi ^2 AU^3/yr^2.$$
Write a program to simulate motion in a central force field. Verify the case of circular orbit using (in astronomical units) ($x_0=1$, $y_0=0$ and $v_x(t=0)=0$. Use the condition ([circular] to calculate $v_y(t=0)$ for a circular orbit. Choose a value of $\Delta t$ such that to a good approximation the total energy $E$ is conserved. Is your value of $\Delta t$ small enough to reproduce the orbit over several periods?
Run the program for different sets of initial conditions $x_0$ and $v_y(t=0)$ consistent with the condition for a circular orbit. Set $y_0=0$ and $v_x(t=0)=0$. For each orbit, measure the radius and the period to verify Kepler’s third law ($T^2/a^3=\mathrm{const.}$). Think of a simple condition which allows you to find the numerical value of the period.
Show that Euler’s method does not shield stable orbits for the same choice of $\Delta t$ used in the previous items. Is it sufficient to simply choose a smaller $\Delta t$ or Euler’s method is no stable for this dynamical system? Use the average velocity $1/2(v_n+v_{n+1})$ to obtain $x_{n+1}$. Are the results any better?
Set $y_0=0$ and $v_x(t=0)=0$. By trial and error find several choices of $x_0$ and $v_y(t=0)$ which yield coonvinient elliptical orbits. Determine total energy, angular momentum, semi-major and semi-nimor axes, eccentricity, and period for each orbit.
You probably noticed that Euler’s algorithm with a fixed $\Delta t$ breaks down if you get to close to the sun. How are you able to visually confirm this? What is the cause of the failure of the method? Think of a simple modification of your program that can improve your results.
In [17]:
class particle2(object):
def __init__(self, mass=1., x=0., y=0., vx=0., vy=0.):
self.mass = mass
self.x = x
self.y = y
self.vx = vx
self.vy = vy
def euler(self, fx, fy, dt):
self.vx = self.vx + fx/self.mass*dt
self.vy = self.vy + fy/self.mass*dt
self.x = self.x + self.vx*dt
self.y = self.y + self.vy*dt
In [33]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot
from matplotlib.colors import ColorConverter as cc
import math
GM = 4*math.pi**2
r = 1. # radius of the orbit
v0 = math.sqrt(GM/r) # This is the condition for circular orbits
x0 = r # initial position
y0 = 0. # we asume we start from the x axis
v0x = 0. # and the initial velocity
v0y = v0 # is perpendicular to the vector position.
dt = 0.001 # time step
tmax = 4.
nsteps = int(tmax/dt)
x = np.zeros(nsteps)
y = np.zeros(nsteps)
vx = np.zeros(nsteps)
vy = np.zeros(nsteps)
energy = np.zeros(nsteps)
x[0] = x0
y[0] = y0
vx[0] = v0x
vy[0] = v0y
energy[0] = 0.5*(v0y**2+v0x**2) - GM/r
p = particle2(1., x0, y0, v0x, v0y)
for i in range(1,nsteps):
r = math.sqrt(p.x*p.x+p.y*p.y)
r3 = r * r * r
fx = -GM*p.x/r3
fy = -GM*p.y/r3
p.euler(fx, fy, dt)
x[i] = p.x
y[i] = p.y
vx[i] = p.vx
vy[i] = p.vy
energy[i] = 0.5*(p.vx**2+p.vy**2) - GM/r
t = np.linspace(0.,tmax,nsteps)
pyplot.plot(t, energy, color='green', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('Energy');
In [34]:
pyplot.plot(t, x, color='red', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('position x');
In [35]:
pyplot.plot(x, y, color='blue', ls='-', lw=3)
pyplot.xlabel('position x')
pyplot.ylabel('position y');
The presence of other planets implies that the total force on a planet is no longer a central force. Furthermore, since the orbits are not exactly on the same plane, the analysis must be extended to 3D. However, for simplicity, we are going to consider a two-dimensional solar system, with two planets in orbit around the sun.
The equations of motion of the two planets of mass $m_1$ and $m_2$ can be written in vector form as $$\begin{aligned} && m_1\frac{d^2 {\mathbf r}_1}{dt^2}=-\frac{m_1MG}{r_1^3}{\mathbf r}_1+\frac{m_1m_2G}{r_{21}^3}{\mathbf r}_{21}, \\ && m_2\frac{d^2 {\mathbf r}_2}{dt^2}=-\frac{m_2MG}{r_2^3}{\mathbf r}_2+\frac{m_1m_2G}{r_{21}^3}{\mathbf r}_{21},\end{aligned}$$ where ${\mathbf r}_1$ and ${\mathbf r}_2$ are directed form the sun to the planets, and ${\mathbf r}_{21}={\mathbf r}_2-{\mathbf r}_1$ is the vector from planet 1 to planet 2. This is a problem with no analytical solution, but its numerical solution can be obtained extending our previous analysis for the two-body problem.
Let us consider astronomical units, and values for the masses $m_1/M=0.001$ and $m_2/M=0.01$. Consider initial positions $r_1=1$ and $r_2=4/3$ and velocities ${\mathbf v}_{1,2}=(0,\sqrt{GM/r_{1,2}})$.
Write a program to calculate the trajectories of the two planets, and plot them.
What would the shape and the periods of the orbits be if the don’t interact? What is the qualitative effect of the interaction. Why is one planet affected more by the interaction that the other? Are the angular momentum and energy of planet 1 conserved? Is the total momentum an energy of the two planets conserved?
In [38]:
NPLANETS = 2
m1 = 0.1
r1 = 1. # radius of the orbit
v1 = math.sqrt(GM/r1) # This is the condition for circular orbits
m2 = 0.01
r2 = 4./3. # radius of the orbit
v2 = math.sqrt(GM/r2) # This is the condition for circular orbits
dt = 0.001 # time step
tmax = 2.
nsteps = int(tmax/dt)
x = np.zeros(shape=(nsteps,NPLANETS))
y = np.zeros(shape=(nsteps,NPLANETS))
vx = np.zeros(shape=(nsteps,NPLANETS))
vy = np.zeros(shape=(nsteps,NPLANETS))
energy = np.zeros(shape=(nsteps,NPLANETS))
x[0,0] = r1
y[0,0] = 0.
vx[0,0] = 0.
vy[0,0] = v1
energy[0][0] = 0.5*v1**2 - GM/r1
x[0,1] = r2
y[0,1] = 0.
vx[0,1] = 0.
vy[0,1] = v2
energy[0][1] = 0.5*v2**2 - GM/r2
planets = [particle2]
for i in range(1,NPLANETS): # create a list of NPLANETS particle2's
planets.append([particle2])
m = [m1, m2]
for i in range(NPLANETS):
planets[i] = particle2(m[i], x[0,i], y[0,i], vx[0,i], vy[0,i])
for i in range(1,nsteps):
for n in range(0,NPLANETS):
r = math.sqrt(planets[n].x*planets[n].x+planets[n].y*planets[n].y);
r3 = r * r * r;
fx = -GM*planets[n].mass*planets[n].x/r3;
fy = -GM*planets[n].mass*planets[n].y/r3;
planets[n].euler(fx, fy, dt)
x[i,n] = planets[n].x
y[i,n] = planets[n].y
vx[i,n] = planets[n].vx
vy[i,n] = planets[n].vy
energy[i,n] = 0.5*(planets[n].vx**2+planets[n].vy**2) - GM/r
t = np.linspace(0.,tmax,nsteps)
In [39]:
pyplot.plot(x[:,0], y[:,0], color='blue', ls='-', lw=3)
pyplot.plot(x[:,1], y[:,1], color='red', ls='-', lw=3)
pyplot.xlabel('position x')
pyplot.ylabel('position y');
In [40]:
pyplot.plot(t, energy[:,0], color='blue', ls='-', lw=3, label='Planet 1')
pyplot.plot(t, energy[:,1], color='red', ls='-', lw=3, label = 'Planet 2')
pyplot.plot(t, energy[:,0]+energy[:,1], color='green', ls='-', lw=3, label = 'Both planets')
pyplot.xlabel('Energy')
pyplot.ylabel('time');
In this particular case, since both planets do not interact, the individual energies are conserved.