Universal Gravitation
Two bodies are mass \( m \) and \( M \), where the first body is a circular orbit with respect to the second.
The centripetal force that acts on \(m\) is:
\[
F\_{m} = m\frac{v^{2}}{r}
\]
We know that for a circular orbit \(v = \frac{2\pi r}{T}\), using this expression in the third law of Kepler \(T^{2}=Kr^{3}\), we have:
\[
v^{2} = \frac{4\pi^{2} r^{2}}{Kr^{3}} = \frac{4\pi^{2}}{Kr} \Rightarrow v^{2} \propto \frac{1}{r}
\]
Substituting (2) into (1):
\[
F\_{m} = \frac{m4\pi^{2}}{Kr^{2}} \Rightarrow F\_{m} \propto \frac{m}{r^{2}}
\]
By Newton’s third law, the body of mass \(m\) exerts an equal and opposite force on the mass \(M\). The centripetal force acting on the body mass \(M\) is given by:
\[
F\_{M} \propto \frac{M}{r^{2}}
\]
As \(m|F\_{M}| = M|F\_{m}|\), Newton deduced that:
\[
F\_{grav} = \frac{GMm}{r^{2}}
\]
Using Newton’s second law (3):
\[
m\frac{d^{2}r}{dt^{2}} = \frac{GMm}{r^{2}}
\]
Despite the gravitational force have been obtained for a system with two bodies and circular orbit, it is more general, and is valid for systems with multiple bodies gravitationally interacting.
Its final form is:
\[
m_{i}\frac{d^{2}\vec{r}_{i}}{dt^{2}} = -\sum\_{j=1,j\neq i}^{N} Gm_{i}m_{j}\frac{\vec{r_{i}}-\vec{r_{j}}}{|\vec{r_{i}}-\vec{r_{j}}|^{3}}
\]
For i = 1, …, n.
Which can be separated for each coordinate x, y and z:
\[
m_{i}\frac{d^{2}x_{i}}{dt^{2}} = -\sum_{j=1,j\neq i}^{N} \frac{Gm_{i}m_{j}(x_{i}-x_{j})}{r_{ij}^{3}}
\]
\[
m_{i}\frac{d^{2}y_{i}}{dt^{2}} = -\sum_{j=1,j\neq i}^{N} \frac{Gm_{i}m_{j}(y_{i}-y_{j})}{r_{ij}^{3}}
\]
\[
m_{i}\frac{d^{2}z_{i}}{dt^{2}} = -\sum_{j=1,j\neq i}^{N} \frac{Gm_{i}m_{j}(z_{i}-z_{j})}{r_{ij}^{3}}
\]
Where \(r_{ij} = \sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}+(z_{i}-z_{j})^{2}}\).
Numerical solutions of ODE
There are several types of numerical methods for solving ordinary differential equations, we see two of them, the improved Euler method and the Runge-Kutta fourth order, they were written to solve ODEs of the first order, so to solve the 2nd order we have to separate them in a system of ODEs of the first order.
Euler method Improved
This method is an extension of Euler’s method, it is relatively simple, and very quickly. Its general form is:
\[
y_{i+1}^{1}=y_{i}+hf(x_{i},y_{i})
\]
\[
y_{i+1}=y_{i}+h\frac{(f(x_{i},y_{i})+f(x_{i+1},y_{i+1}^{1}))}{2}
\]
The algorithm for EDO second order:
a1 = f(t,x0,v0) x1 = x0 + h*v0 v1 = v0 + h*a1 a = f(t,x1,v1) x = x0 + h/2*(v0+v1) v = v0 + h/2*(a1+a) t = t+h
Runge-Kutta method of 4th order
The Runge-Kutta 4th order is much more accurate than the improved Euler, and also slower. Its general form is:
\[
k1=hf(x_{i},y_{i})
\]
\[
k2=hf(x_{i}+\frac{h}{2},y_{i}+\frac{k1}{2})
\]
\[
k3=hf(x_{i}+\frac{h}{2},y_{i}+\frac{k2}{2})
\]
\[
k4=hf(x_{i}+h,y_{i}+k3)
\]
\[
y_{i+1}=y_{i}+\frac{1}{6}(k1+2k2+2k3+k4)
\]
\[
x_{i+1}=x_{i}+h
\]
The algorithm for EDO second order:
v1 = v a1 = f(t,x,v) v2 = v+a1*h/2 a2 = f(t+h/2,x+v1*h/2,v2) v3 = v+a2*h/2 a3 = f(t+h/2,x+v2*h/2,v3) v4 = v+h*a3 a4 = f(t+h,x+v3*h,v4) x = x + h*(v1+2*v2+2*v3+v4)/6 v = v + h*(a1+2*a2+2*a3+a4)/6 t = t+h
Python
The algorithm was developed using version 3.7 of Python, and the only dependencies are the random package that generates random numbers, and the NumPy that supports arrays and multidimensional arrays. To plot the data obtained was used Matplotlib package.
The file name is nbody.py, it can be imported with the following command:
from nbody import *
The class methods are:
def __init__(_, G=6.674287e-11): def f(_, i, m, x): def eulerm(_, f, m, x0, v0, dt): def rk4(_, f, m, x0, v0, dt): def calculate(_, m, x0, v0, dt, T, method):
init(..) is the class constructor, if not parameter passed it assumes G = 6.674287e-11, f(..) calculate the acceleration, eulerm(..) is the improved Euler method used to calculate speed and position RK4(..) is the Runge-Kutta method of 4th order and have the same purpose as the improved Euler method, and calculate(..) is used to perform all calculations, and returns as output, the list (t, x, v, a), where \(x \Rightarrow x(i,n,c)\), i is the vector index of all positions for the same point in time, n identifies the body, c identifies the component (x = 0, y = 1, z = 2).
Files download: https://github.com/clnrp/nbody
Two bodies gravitationally interacting
This is a problem that can usually be solved analytically, and consider two typical cases.
Law of elliptical orbits:
\[
r=\frac{a(1-\epsilon ^{2})}{1\pm \epsilon cos(\theta)}
\]
The every body equations of motion are:
\[
m_{1}\frac{d^{2}\vec{r}_{1}}{dt^{2}} = -Gm_{1}m_{2}\frac{\vec{r_{1}}-\vec{r_{2}}}{|\vec{r_{1}}-\vec{r_{2}}|^{3}}
\]
\[
m_{2}\frac{d^{2}\vec{r}_{2}}{dt^{2}} = -Gm_{2}m_{1}\frac{\vec{r_{2}}-\vec{r_{1}}}{|\vec{r_{2}}-\vec{r_{1}}|^{3}}
\]
Body in free fall on the surface of a planet
To solve this problem we must make some simplifications first if we define that \(m_{1}\) is the mass of the planet, and \(m_{2}\) is the mass of a body on its surface, so we will have \(m_{1} >> m_{2}\), and we can consider \(\frac{d^{2}\vec{r}_{1}}{dt^{2}} \simeq 0\), ie we need only solve the equation (6).
Second, we can assume that \(\vec{r}_{1}\) this at the origin (0,0,0), and \(\vec{r}_{2}\) this near the surface of the planet we have \(|\vec{r}_{2}| \simeq R\).
And if we assume that local displacements on the surface of the planet can be approached to in a Euclidean plane “xy” and that gravity acts in the direction perpendicular to this plane “z coordinate,” we will have the left of (5), (6) and (7):
\[
\frac{d^{2}x_{2}}{dt^{2}} = \frac{d^{2}y_{2}}{dt^{2}} = 0
\]
\[
\frac{d^{2}z_{2}}{dt^{2}} = \frac{-GM}{R^{2}}
\]
And in the case considered \(\frac{GM}{R^{2}}=g\), is the planet’s gravity on the surface.
\[
\frac{d^{2}z_{2}}{dt^{2}} = -g
\]
Integrating (11):
\[
x_{2}(t) = x_{0} + v_{0x}t
\]
\[
y_{2}(t) = y_{0} + v_{0y}t
\]
Integrating (12):
\[
\frac{d}{dt}(\frac{dz_{2}}{dt}) = -g
\]
\[
z_{2}(t) = z_{0} + v_{0z}t -\frac{gt^{2}}{2}
\]
Example: A body is launched from the earth’s surface with an initial velocity of 50 m / s at an angle relative to the ground 60 degrees.
The time for an object thrown reach the maximum height “\(h_{max}\)” occurs when \(\frac{dz_{2}(t)}{dt}=0\), the \(t_{h_{max}}=\frac{v_{0z}}{g}\), and \(h_{max}=\frac{v_{0z}^{2}}{2g}\). The total time is \(T=2t_{h_{max}}\), and admitted that the movement takes place in the plane “yz”, the horizontal offset “\(D\)” is obtained by replacing \( T \) in (14), \(D=v_{0y}T\).
The total time is used to set the number of numerical calculation steps, and the maximum height and horizontal displacement to set the axes of the graph generated.
Code of numerical solution:
from nbody import * # international system of units G=6.674287e-11 # m^3 * kg^-1 * s^-2 M=5.9722e24 # kg R=6400e3 V=50 # m/s phi=pi/3 # ang rad g=G*M/R**2 T=2*V*sin(phi)/g # time (s) dt=0.01 Vy=V*cos(phi) Vz=V*sin(phi) m=[] x0=[] v0=[] m.append(M) # m1 m.append(5e2) # m2 x0.append(array([0.0,0.0,0.0])) # x1 x0.append(array([0,0,R])) # x2 v0.append(array([0.0,0.0,0.0])) # v1 v0.append(array([0.0,Vy,Vz])) # v2 nbody = nbody(G) [t1,x1,v1,a1]=nbody.calculate(m,x0,v0,dt,T,'euler') [t2,x2,v2,a2]=nbody.calculate(m,x0,v0,dt,T,'rk4')
To measure the error, we will use the mean absolute percentage error “MAPE”.
EMPA data obtained through the Euler improved methods and Runge-Kutta in relation to equations (13), (14) and (15), is 0.000333592%, 0.000333591%, respectively.
Satellite in an elliptical orbit
The total energy of the system in polar coordinates is:
\[
E=\frac{m}{2}(\dot{r}^{2}+r^{2}\dot{\theta}^{2})-\frac{GMm}{r}
\]
The Lagrangian that \(L=\frac{m}{2}(\dot{r}^{2}+r^{2}\dot{\theta}^{2})+\frac{GMm}{r}\), we obtain the angular momentum is \(l=\frac{\partial L}{\partial \dot{\theta} }=m r^{2}\dot{\theta }\) and to \(\frac{d}{dt}(\frac{\partial L}{\partial \dot{\theta}})=0\), we have that \(l\) is a constant movement.
Replacing \(l\) in (16):
\[
E=\frac{m\dot{r}^{2}}{2} + \frac{l^{2}}{2mr^{2}} -\frac{GMm}{r}
\]
isolating \(\dot{r}\):
\[
\frac{dr}{dt }=\sqrt{\frac{2}{m}[E – \frac{l^{2}}{2m r^{2}} +\frac{GMm}{r}]}
\]
\[
\frac{dr}{d\theta }=\frac{dr}{dt}\frac{dt}{d\theta }=\frac{mr^{2}}{l}\frac{dr}{dt}
\]
\[
\frac{dr}{d\theta }=\frac{m r^{2}}{l}\sqrt{\frac{2}{m}[E – \frac{l^{2}}{2m r^{2}} +\frac{GMm}{r}]}
\]
\[
\frac{dr}{d\theta }=r^{2}\sqrt{\frac{2mE}{l^{2}} – \frac{1}{r^{2}} +\frac{2GMm^{2}}{l^{2}r}}
\]
Making \(u=1/r\), we obtain \(dr=r^{2}du\) replacing:
\[
\frac{du}{d\theta }=-\sqrt{\frac{2mE}{l^{2}} +\frac{2GMm^{2}}{l^{2}}u – u^{2}}
\]
\[
\theta – \theta _{0} = -\int _{u_{0}}^{u}\frac{du}{\sqrt{\frac{2mE}{l^{2}} +\frac{2GMm^{2}}{l^{2}}u – u^{2}}}
\]
Which is an integral type \(\int \frac{dx}{\sqrt{a + bx + cx^{2}}}=\frac{1}{\sqrt{-c}}arc cos(-\frac{b+2cx}{\sqrt{b^{2}-4ac}})+k\).
integrating:
\[
\frac{1}{r} = \frac{GMm^{2}}{l^{2}}[1+\sqrt{1+\frac{2El^{2}}{G^{2}M^{2}m^{3}}}cos(\theta – \theta _{0})]
\]
And it can be rewritten as:
\[
\frac{1}{r} = C[1+\epsilon cos(\theta – \theta _{0})]
\]
What it is a conic with one focus at the origin, and its kind of eccentricity of overhangs \(\epsilon = \sqrt{1-\frac{2El^{2}}{G^{2}M^{2}m^{3}}}\).
EccentricityEnergyType \(\epsilon = 0\)$E = -\frac{m k^{2}}{2l^{2}}$Circular$0 < \epsilon < 1$E < 0Elliptical$\epsilon = 1$E = 0Satellite dish$\epsilon > 1$E > 1Hyperbolic
For an ellipse, the equation (18) is equivalent to (8), but in Cartesian coordinates as follows:
\[
x=a(cos(u)-\epsilon)
\]
\[
y=asin(u)\sqrt{1-\epsilon^{2}}
\]
Kepler’s equation provides a relationship between the eccentric anomaly “\(u\)” and time:
\[
u-\epsilon sen(u)=M=2\pi \frac{t}{T}
\]
Where M is known as the mean anomaly.
From the orbital energy, we can obtain an equation for calculating the orbital velocity, and it is known as vis-viva:
\[
v^{2}=GM\left(\frac{2}{r}-\frac{1}{a}\right)
\]
Example: A satellite in Earth orbit with eccentricity 0.5 and 2 hours.
Through the third Kepler law, we calculate the value of semi-major axis “\(a\)”:
\[
a=\left(\frac{T^{2}}{K}\right)^{1/3}
\]
From equation (8), we get \(r\) at perihelion \(\theta = 0\):
\[
r=\frac{a(1-\epsilon^{2})}{1+\epsilon cos(0)}
\]
For the speed at perihelion, replace the value of \(a\) and \(r\) the equation vis-viva.
Code of numerical solution:
from nbody import * # system units # u.t. hours # u.c. kilometers # u.m. kg d_km=1000 # 1km = 1000m t_hour=3600 # 1h = 3600s G=6.674287e-11 # m^3 * kg^-1 * s^-2 G=G*(t_hour**2)/(d_km**3) # m^3 * kg^-1 * hour^-2 M=5.9722e24 # earth mass e=0.5 # eccentricity P=2 # period in hours K=(4*pi**2)/(G*M) a=(P**2/K)**(1/3.0) #eq ellipse, 0 perihelium r = a*(1-e**2)/(1+e*cos(0)) #vis-viva v = sqrt(G*M*(2/r-1/a)) Dp=r Vp=v T=P dt=0.002 D=Dp m=[] x0=[] v0=[] m.append(M) # m1 m.append(5e2) # m2 x0.append(array([0.0,0.0,0.0])) # x1 x0.append(array([D,0,0])) # x2 v0.append(array([Vp*0.00,Vp*0.00,Vp*0.00])) # v1 v0.append(array([0.0,Vp,0.0])) # v2 nbody = nbody(G) [t1,x1,v1,a1]=nbody.calculate(m,x0,v0,dt,T,'euler') [t2,x2,v2,a2]=nbody.calculate(m,x0,v0,dt,T,'rk4')
MAPE numerical data obtained in relation to equations (19) and (20) was 1.1% for improved Euler method and 5.1e-06% for the Runge-Kutta fourth order.
More than two bodies interacting gravitationally
This is a type of problem that usually can not be solved analytically, and consider two hypothetical situations where the masses, positions and bodies of the speeds will be generated to form systems with a certain stability, so we use for the initial time a model in which the orbit it is given by a closed conical as seen in the case of a satellite in an elliptical orbit, but in this case and the eccentricity is zero, i.e. is initially circular orbit. As the mass of bodies are not negligible in relation to the greater body mass, the system dynamics over time is chaotic.
Recalling that:
Orbital velocity, known as vis-viva equation:
\[
v^{2}=GM\left(\frac{2}{r}-\frac{1}{a}\right)
\]
Her left we get the speed for each body:
\[
v=\sqrt{\frac{GM_{p}}{R}}
\]
\(M_{p}\) is greater value the mass, and R is the distance of the bodies in relation to body mass \(M_ {p}\), being orthogonal to the velocity vector \(\vec{R}\). The initial conditions will be obtained through the code below:
def val(value):
return (-1)**int(random()*100) *random()*value
def mult_star(Pp,Vp,Mp,Dp,N,sn):
[m,x0,v0]=sn
# body M0
m.append(Mp*.8)
x0.append(Pp)
v0.append(Vp)
# other bodies
for i in range(1,N):
m.append(0.15*Mp*abs(val(1))) # m
x0_a=array([val(Dp),val(Dp),val(Dp)])+Pp
x0.append(x0_a) # position
rp=x0_a-Pp # position relative to M0
R=linalg.norm(rp)
Vn=sqrt(G*Mp/R)
# the vector velocity and orthogonal to rp
vt=array([rp[1], -rp[0], 0])
vtn=linalg.norm(vt)
vel=Vn*vt/vtn+0.2*Vp
v0.append(vel) # velocity
return [m,x0,v0]
Multiple star systems
A multiple star system consists of a small number of stars, the most common are composed of three stars, and generally have a close binary and a more distant orbital companion. These systems may be unstable, possibly resulting in the ejection of one or more stars. Usually they have more difficulty in modeling such systems than the binary because of its chaotic nature.
Example: A 5-star system with an average distance of 12 ua, the largest mass is approximately twice that of the sun, and the evolution period is 60 years.
Example code:
from nbody import * from random import * # system units # u.t. years # u.l. astronomical unit (au) # u.m. sun mass (msun) d_au=149597870700 t_year=31536000 m_sun=1.989e30 G=6.674287e-11 # m^3 * kg^-1 * s^-2 # au^3 * msun^-1 * years^-2 G=G*(t_year**2)*(m_sun)/(d_au**3) Mp=3 # 3 sun mass N=5 T=60 # 60 years dt=0.005 Dp=12 # 12 au D=Dp*3 m=[] x0=[] v0=[] sn=[m,x0,v0] Pp=array([0,0,0]) Vp=array([0.8,0,0]) [m,x0,v0] = mult_star(Pp,Vp,Mp,Dp,N,sn) nbody = nbody(G) [t1,x1,v1,a1]=nbody.calculate(m,x0,v0,dt,T,'euler') [t2,x2,v2,a2]=nbody.calculate(m,x0,v0,dt,T,'rk4') Initial data: m=[2.4000000000000004, 0.15810696549158487, 0.16559642134905525, 0.25212587182456586, 0.0650608626844281] x0=[array([0, 0, 0]), array([-5.9679843, 1.71668818, -6.57311121]), array([0.01378749, -4.09105872, 6.19174474]), array([-1.64727009, -8.87201438, 9.03102758]), array([-7.88436924, -8.77524908, 8.25068406])] v0=[array([0.8, 0., 0.]), array([0.99988975, 3.47606884, 0.]), array([-3.99262512, -0.01345575, 0.]), array([-2.99296721, 0.55570529, 0.]), array([-2.13240063, 1.91591529, 0.])]
The EMPA the improved Euler method over the Runge-Kutta 4th order was 1.3%.