Computational Physics with
Python: Planetary Orbits
from Newton to Feynman
Pramod Gupta
Department of Astronomy, University of Washington
Slide 2
Slide 2 text
Kepler
Kepler’s first law: The planets move
around the sun in an elliptical orbit with
the sun at one of the two foci of the
ellipse.
Slide 3
Slide 3 text
Ellipse
Slide 4
Slide 4 text
Newton
Newton applied his laws of motion and
law of gravity to the orbit of planets and
derived Kepler’s laws.
The derivation requires advanced
mathematics.
Is there some other way?
Slide 5
Slide 5 text
Feynman
Slide 6
Slide 6 text
Feynman
Feynman: “If force law is known (as
function of positions and velocities)
equations can be solved in close
approximation step by step in time by
arithmetic.”
Slide 7
Slide 7 text
Calculating the orbit
Feynman: “Now how do we make the
calculation......This work can be done
rather easily by using a table of squares,
cubes, and reciprocals: then we need
only multiply x by 1/r3, which we do on
a slide rule.”
Slide 8
Slide 8 text
Slide Rule
Slide 9
Slide 9 text
Computers
The computer can do arithmetic.
But we must tell it what to compute...
Slide 10
Slide 10 text
Newton’s law of Gravity
F =
GMm
r2
G is Newton’s gravitational constant
M is mass of sun
m is mass of the planet
r is the distance of the planet from the
sun
Slide 11
Slide 11 text
Newton’s law of Gravity
Direction of the force is from the planet
to the sun
Slide 12
Slide 12 text
Components of the Force
Fx
= −GMm
r2
x
r
Newton’s second law Fx
= max
Hence, max
= −GMm
r2
x
r
ax
= −GM
r2
x
r
Similarly
ay
= −GM
r2
y
r
Slide 13
Slide 13 text
Components of the Force
Choose units such that GM=1. Hence:
ax
= − 1
r2
x
r
= − x
r3
Similarly
ay
= − 1
r2
y
r
= − y
r3
By Pythagoras theorem:
Distance from sun to the planet
r = x2 + y2
Slide 14
Slide 14 text
Position, Velocity, Accelaration
Position of the planet at time t:
(x(t), y(t))
Velocity of the planet at time t:
(vx
(t), vy
(t))
Acceleration of the planet at time t:
(ax
(t), ay
(t))
Slide 15
Slide 15 text
Position at time t +
is a very small time interval
vx
= ∆x
∆t
vx
(t) = x(t+ )−x(t)
(t+ )−t
= x(t+ )−x(t)
x(t + ) = x(t) + vx
(t)
Similarly:
y(t + ) = y(t) + vy
(t)
Slide 16
Slide 16 text
Velocity at time t +
is a very small time interval
ax
= ∆vx
∆t
ax
(t) = vx(t+ )−vx(t)
(t+ )−t
= vx(t+ )−vx(t)
vx
(t + ) = vx
(t) + ax
(t)
Similarly:
vy
(t + ) = vy
(t) + ay
(t)
Slide 17
Slide 17 text
Leapfrog Integration
Feynman uses leapfrog integration
which is more accurate.
It also has some nice properties.
It conserves energy.
It has time reversal symmetry.
Slide 18
Slide 18 text
Leapfrog Integration
Calculate position (x, y) at
t = 0, , 2 ...
Calculate velocity (vx
, vy
) at
t =
2
, 3
2
, 5
2
...
Calculate acceleration (ax
, ay
) at
t = 0, , 2 ...
Slide 19
Slide 19 text
Leapfrog Integration Initialization
We know the state of the system at
t = 0.
So we know x(0), y(0), vx
(0), vy
(0).
Using x(0), y(0) we can calculate ax
(0),
ay
(0).
vx
(
2
) = vx
(0) +
2
ax
(0)
vy
(
2
) = vy
(0) +
2
ay
(0)
Choose Initial Conditions at t = 0
Feynman’s initial conditions
Just after equation 9.17 in Chapter 9.
x initial = 0.5
y initial = 0.0
vx initial = 0.0
vy initial = 1.63
Slide 23
Slide 23 text
Choose time step size epsilon
Smaller epsilon gives more accurate
results.
Feynman used epsilon=0.1 to minimize
computation.
We will use epsilon=0.001 for his initial
conditions.
For other initial conditions, one may
need an even smaller epsilon.
Slide 24
Slide 24 text
Feynman’s solution
Slide 25
Slide 25 text
Conclusion
Feynman: Now, armed with the tremendous
power of Newton’s laws, we can not only
calculate such simple motions but also,
given only a machine to handle the
arithmetic, even the tremendously complex
motions of the planets, to as high a degree
of precision as we wish!
Slide 26
Slide 26 text
References
The Feynman Lectures on Physics
Volume 1 Chapter 9
www.feynmanlectures.caltech.edu
Numpy + Scipy + Matplotlib
Anaconda Python/ Canopy Python
Slide 27
Slide 27 text
Code
import math
import numpy as np
import matplotlib.pyplot as plt
#Feynman Lectures on Physics Volume 1 Chapter 9
#Smaller epsilon gives more accurate results.
#Feynman used epsilon=0.1 to minimize computation.
#We use epsilon=0.001 for his initial conditions.
#For other initial conditions,
#one may need an even smaller epsilon.
epsilon = 0.001
MAX_STEPS = 10000
Slide 28
Slide 28 text
Code
#Feynman initial conditions just after eq 9.17
x_initial = 0.5
y_initial = 0.0
vx_initial = 0.0
vy_initial = 1.63
#allocate lists (like arrays in other languages)
x = [0.0]*MAX_STEPS
y = [0.0]*MAX_STEPS
vx = [0.0]*MAX_STEPS
vy = [0.0]*MAX_STEPS
ax = [0.0]*MAX_STEPS
ay = [0.0]*MAX_STEPS
r = [0.0]*MAX_STEPS