Orbital Motion
Introduction: I will simulate an orbit using Newton's equations of
motion to test Kepler's laws of planetary motion. I have used Euler,
Euler-Cromer and Euler-Richardson algorithms to simulate an orbital system that
obeys Newton's universal law of gravitation.
Method: Newton's universal law of gravitation states that the
gravitational force acting on two bodies are directly proportional to their
mass and inversely proportional to their distance from each other:
F = - G M m / r^2. (1a)
(G is the universal gravitational constant (G = 6.67*10-11 m3 / Kg s2).
Since force is a vector, we can rewrite this as:
F= (- G M m / r^3) r (1b)
Combing Newton's law of gravitation and laws of motion, we can produce
various algorithms to simulate an orbit.
Newton's 2nd law of motion states that :
F = ma. (1c)
If we set the coordinate system at mass M, we can write the equation of
motion of the particle of mass m, as in the book (p62), as:
m*d^2r/dt^2 = (-G M m / r^3) r. (1d)
Hence, we see that:
d^2x/dt^2 = (- G M / r^3) x (1e)
and
d^2y/dt^2 = (- G M / r^3) y. (1f)
(d^2r/dt^2 is the second order derivative of distance with respect to time, which
is the acceleration).
The Newton's laws of motion can be written as two coupled first order
differential equations. The numerical solution for the Euler algorithm can be
written in the form:
Vn+1 = vn + an delta_t (2a)
(an is the acceleration from equations 1e and 1f).
yn+1 = yn + vn delta_t. (2b)
The Euler algorithm can be modified slightly to give the Euler-Cromer
algorithm:
Vn+1 = vn + an delta_t (3a)
yn+1 = yn + vn+1 delta_t. (3b)
Or we could use an algorithm that evaluates the acceleration and velocity
at time, tmid = t + delta_t/2. This is called the Euler-Richardson algorithm.
vmid = vn + an delta_t/2 (4a)
ymid = yn + vn delta_t/2 (4b)
amid = F (ymid, vmid, tmid) / m (4c)
therefore,
vn+1 = vn + amid delta_t (5a)
yn+1 = yn + vmid delta_t . (5b)
(This is given in the book in page 42.)
Notice that this is a two-body problem, therefore we will have x and y
components for the acceleration, position and velocity.
Verification: I down-loaded Program planet from the from the Computer
Simulation Laboratory home page. Since this program adopted the Euler-Cromer
algorithm, I modified it to test the Euler and the Euler-Richardson algorithms.
The subroutines are listed below:
SUB Euler_Cromer(pos(),vel(),t,GM,dt) ! Euler-Cromer algorithm
DIM accel(2)
LET r^2 = pos(1)*pos(1) + pos(2)*pos(2)
LET r^3 = r^2*sqr(r^2)
FOR i = 1 to 2
LET accel(i) = -GM*pos(i)/r^3
LET vel(i) = vel(i) + accel(i)*dt
LET pos(i) = pos(i) + vel(i)*dt
NEXT i
LET t = t + dt
END SUB
(notice that you only have to switch lines 8 and 9 around to get the Euler from
the Euler-Cromer algorithm)
SUB Euler(pos(),vel(),t,GM,dt) ! Euler algorithm
DIM accel(2)
LET r^2 = pos(1)*pos(1) + pos(2)*pos(2)
LET r^3 = r^2*sqr(r^2)
FOR i = 1 to 2
LET accel(i) = -GM*pos(i)/r^3
LET pos(i) = pos(i) + vel(i)*dt !for euler put this statement
before
LET vel(i) = vel(i) + accel(i)*dt !defining this
NEXT i
LET t = t + dt
END SUB
SUB Euler_Richardson(pos(),vel(),t,GM,dt) ! Euler_Richardson algorithm
DIM accel(2)
DIM vmid(2)
DIM pmid(2)
DIM amid(2)
LET r^2 = pos(1)*pos(1) + pos(2)*pos(2)
LET r^3 = r^2*sqr(r^2)
FOR i = 1 to 2
LET accel(i) = -GM*pos(i)/r^3
LET vmid(i)=vel(i) + accel(i)*(dt/2) !velocity at middle
LET pmid(i)=pos(i) + vel(i)* (dt/2) !position at middle
NEXT i
LET r^2 = pmid(1)*pmid(1) + pmid(2)*pmid(2)
LET r^3 = r^2*sqr(r^2)
FOR i= 1 to 2
LET amid(i)= -GM*pmid(i)/r^3 !acceleration at middle
LET vel(i) = vel(i) + amid(i)*dt !velocity using euler-richardson
LET pos(i) = pos(i) + vmid(i)*dt !position using euler-richardson
NEXT i
LET t = t + dt
END SUB
To find out which algorithm to adopt, I will compare the outputs of the
algorithm with a known output. We know that the acceleration of a circular
orbit is given by:
a = v2/r (6a)
Since F=ma (by Newton's2nd law), the force is given by:
mv2/r = GMm/r^2 (6b)
therefore,
v=(GM / r)1/2. (6c)
By plugging in an initial value, r=1AU, we can calculate the initial
y-velocity, since we know G M = 4¹2 AU3/yr^2.
vy(initial) = 6.28 AU/yr.
By inputting in 6.28 AU/yr for the initial y velocity and 1AU for the
initial x-position, we should get a circular orbit. We will adopt whichever
algorithm produces the circular orbit that repeats over the most periods.
--Click on Euler, Euler Cromeror Euler Richardson to see the orbit simulated using each
of the algorithms.
Since the Euler-Richardson algorithm produces the most accurate path for a
circular orbit, I will adopt the Euler-Richardson algorithm for the rest of
this report. I will use delta_t = .001 since this gives the least deviation of the
path for a fairly fast speed.
We can also check which algorithm produces the most accurate results by
calculating the energy. Since this is a closed system, the energy should be
conserved due to the Principal of Conservation of Energy. Whichever algorithm
conserves the Energy better, gives the closest output. Energy of the system is
given by:
E = 1/2 m v2 - G M m/ r (7a)
therefore,
E/m = 1/2 (vx2 + vy2) - GM/sqr(x2 + y2) (7b)
It is sufficient to compare the value for E/m since this will also be a
constant.
I will use the following sub-routine to calculate the values for E/m:
SUB Energy(vel(),pos(),GM)
LET r_square=0
LET v_square = 0
FOR i=1 to 2
LET r_square=r_square+pos(i)*pos(i)
LET v_square = v_square +vel(i) *vel(i)
NEXT i
LET r = sqr(r_square)
LET E_over_m = .5*v_square-GM/r
SET CURSOR 1,1
PRINT "Energy/mass =";E_over_m ; " " !clears line to
print over again
END SUB
For the Euler-Cromer algorithm, the value of E/m fluctuated from -19.7592
to -19.7584. The Euler-Richardson gave a constant E/m value of -19.7592. We
can see that the Euler-Richardson conserves energy better than the
Euler-Cromer. This once again proves that Euler-Richardson produces more
accurate results.
Data: My goal is to verify Kepler's 2nd and 3rd laws planetary motion.
Kepler's laws state that:
1) Each planet moves in an elliptical orbit with the Sun located at one of the
foci of the ellipse.
2) The speed of the planet increases as its distance from the sun decreases
such that the line from the sun to the planet sweeps out equal areas in equal
time.
3) The ratio T2/a3 is the same for all the planets that orbit the sun, where T
is the period of the planet and a is the semi-major axis of the ellipse.
To prove Kepler's second law, I have to show that the speed of the planet
is greatest when its distance from the sun is the least and the area swept out
by the planet's path is equal over a time step.
We know that the speed at any given time is:
|v| = sqr (vx2 + vy2) (8a)
I will incorporate this into Planet program to find the positions where the
velocity is at its maximum and minimum. I used the following sub-routine to do
this:
SUB max_min_speed_position(vel(), vmax, vmin, pos(), pos_max(), pos_min())
LET v2 = 0
FOR i = 1 to 2
LET v2 = v2 + vel(i)^2
NEXT i
LET v = sqr (v2)
IF v > vmax THEN
LET vmax = v !determines max speed
FOR i = 1 to 2
LET pos_max(i) = pos(i) !determines position at maximum speed
NEXT i
END IF
IF v < vmin then
LET vmin = v !determines min speed
FOR i = 1 to 2
LET pos_min(i) = pos(i) !determines position at minimum speed
NEXT i
END IF
SET CURSOR 7,1
PRINT "max speed =" ; vmax; " "
SET CURSOR 8,1
PRINT "min speed =" ; vmin; " "
SET CURSOR 9,1
PRINT "position when speed max(x,y) :" ; pos_max(1);","; pos_max(2); " "
SET CURSOR 10,1
PRINT "position when speed min(x,y) :" ; pos_min(1);","; pos_min(2); " "
END SUB
By looking at the output screen (see appendix B), we can see that the
velocity is minimum at the furthest x position in the y-axis (is we take sun to
lie on the y=0 line) from the sun. The velocity is maximum where the planet is
closest to the sun.
We can assume that in one time step, the planet will sweep out a triangle
with area of:
A = 1/2 b h. (9a)
We can rewrite this as,
A = 1/2 delta_t (r * v) (9b)
A= 1/2 delta_t (xvy-yvx). (9c)
We can easily incorporate this into our program:
SUB area (pos(), vel(), dt)
LET A = .5 * dt * (pos(1)*vel(2)-pos(2)*vel(1))
SET CURSOR 11,1
PRINT "area swept out in one time step =" ; A
END SUB
If we run this program (see appendix B), we can see that area is constant over
an equal time step. This verifies Kepler's 2nd Law of Planetary Motion.
To verify Kepler's 3rd Law, we must first calculate the period of the
orbit (T) and its semi-major axis (a). Semi-major axis of an ellipse is just
half the length of the two furthest points of the ellipse. Semi-minor axis (b)
is half the length of the shorter length of an ellipse. We can
calculate them when the sign of the position changes. I have incorporated the following
sub-routine to calculate the semi-major and semi-minor axis of an ellipse:
SUB semimajor (pos(), pos2old, pos1new, pos1newer, a)
IF pos(2)* pos2old <=0 then !calculates semimajor axis when sign changes
IF pos2old>pos(2) then
LET pos1new=pos(1)
ELSE
LET pos1newer = pos(1)
END IF
LET a = abs ((pos1new-pos1newer)/2)
SET CURSOR 4,1
PRINT "semi-major axis ="; a;" "
END IF
END SUB
SUB semiminor (pos(), pos1old, pos2new, pos2newer, b)
IF pos(1) * pos1old <=0 then
IF pos1old>pos(1) then
LET pos2new = pos (2)
ELSE
LET pos2newer = pos(2)
END IF
LET b = abs ((pos2new-pos2newer)/2)
SET CURSOR 5,
PRINT "semi-minor axis =";b;" "
END IF
END SUB
(Note that we must define pos1old and pos2old outside of the if statement in
the do loop.)
The period of the orbit will be twice the time it takes for the position to
change signs. We will use the same trick to find the period as we used to find
the semi-major and semi-minor axis.
SUB PER (pos(), pos2old, tlast, t, period)
IF pos(2) * pos2old <=0 then !to print the period when sign changes
LET period = 2*(t-tlast)
LET tlast = t
SET CURSOR 2,1v
PRINT " "
SET CURSOR 2,1
PRINT "Period =";period !print the period everytime sign changes
END IF
END SUB
Since we have the semi-major and semi-minor axis, we may as well calculate
the eccentricity (e) of the ellipse. Eccentricity is just a property of an
ellipse. It determines the shape of an ellipse (circular ellipse has an
eccentricity of zero). We can calculate eccentricity by:
e = sqr (1 - b2/a2). (10a)
I have collected data of T, a, b and e for many elliptical orbits of convenient
size.
Initial initial E/m Period L/m a
x-position y-velocity (T)
----------------------------------------------------------------------------------
1.0000 4.0000 31.500 0.50000 4.0000 0.63000
1.0000 5.0000 26.900 0.62800 5.0000 0.73000
1.0000 6.0000 21.480 0.88000 6.0000 0.92000
1.0000 7.0000 14.980 1.5000 7.0000 1.3200
2.0000 2.0000 17.710 1.2000 4.0000 1.1100
2.0000 3.0000 15.240 1.4740 6.0000 1.3000
2.0000 4.0000 11.740 2.1800 8.0000 1.6800
2.0000 5.0000 7.2400 4.5000 10.000 2.7300
3.0000 2.0000 11.160 2.5000 6.0000 1.7700
3.0000 3.0000 8.6400 3.4400 9.0000 2.2800
3.0000 4.0000 5.1600 7.5000 12.000 3.8300
b e T2 a3 T2/a3
-----------------------------------------------------------------------------------
0.40000 0.76000 0.25000 0.25000 1.0000
0.63000 0.50000 0.39000 0.40000 0.97500
0.91000 1.2400 0.77000 0.79000 0.97500
1.2400 0.34000 2.2500 2.3000 0.97800
0.41000 0.93000 1.4400 1.3700 1.1000
0.90000 0.71000 2.1700 2.1900 0.99100
1.6200 0.26000 4.7500 4.7400 1.0020
2.5300 0.37000 20.250 20.350 0.99500
0.91000 0.86000 5.5000 5.5500 0.99100
2.1000 0.44000 11.800 11.850 0.99600
3.6500 0.30200 56.250 56.180 1.0010
Table 1. Data of E/m, T, L/m (angular momentum/mass), a, b and e for many
convenient elliptical orbits. Kepler's third law says that T2/a3 is
constant, my values of T2/a3 seems to prove Kepler's 3rd law.
Click here if you want to see the graph of T^2 vs a^3.
Analysis: I used the Euler-Richardson for all the runs since it seemed
to produce the most accurate results. I used Æt = .001, this gave a good orbit
that repeated over most periods. I found that all my data was correct to 3
decimal places and my value for Æt was in the magnitude of 3 decimal places.
Interpretation: All the values of Energy were negative because I was
working in a closed system. This meant that the forces I was dealing with were
attractive since NewtonÕs law of Gravitation defines attractive forces with a
-ve sign. This is obvious because we know that Planets attract one another
causing orbital motion. Energy was slightly less conserved for the elliptical
orbit than the circular orbit because the radius of the elliptical orbit keeps
changing making the mathematics a bit harder causing round off errors.
Critique: I learnt to program in True basic during this lab since there
was so much programming required. I enjoyed simulating motion as large scaled
as planetary movement. I owe some of the credit to Greg Johnson who helped me
find methods to calculate the period and the semi-major axis. I also learnt a
lot about arrays, which was a new concept to me when I started this lab.