PyCon 2016
May 29, 2016
630

# Pramod Gupta - Computational Physics with Python: Planetary Orbits from Newton to Feynman

Newton's explanation of planetary orbits is one of the greatest achievements of science. We will follow Feynman's approach to show how the motion of the planets around the sun can be calculated using computers and without using Newton's advanced mathematics. This talk will convince you that doing physics with Python is way more fun than the way you did physics in high school or university.

https://us.pycon.org/2016/schedule/presentation/1711/

May 29, 2016

## Transcript

1. Computational Physics with
Python: Planetary Orbits
from Newton to Feynman
Pramod Gupta
Department of Astronomy, University of Washington

2. Kepler
Kepler’s ﬁrst law: The planets move
around the sun in an elliptical orbit with
the sun at one of the two foci of the
ellipse.

3. Ellipse

4. Newton
Newton applied his laws of motion and
law of gravity to the orbit of planets and
derived Kepler’s laws.
mathematics.
Is there some other way?

5. Feynman

6. 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.”

7. 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.”

8. Slide Rule

9. Computers
The computer can do arithmetic.
But we must tell it what to compute...

10. 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

11. Newton’s law of Gravity
Direction of the force is from the planet
to the sun

12. 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

13. 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

14. 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))

15. 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)

16. 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)

17. Leapfrog Integration
Feynman uses leapfrog integration
which is more accurate.
It also has some nice properties.
It conserves energy.
It has time reversal symmetry.

18. 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 ...

19. 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)

20. Leapfrog Integration First Iteration
x( ) = x(0) + vx
(
2
)
y( ) = y(0) + vy
(
2
)
r = (x( ))2 + (y( ))2
ax
( ) = −x( )/r3
ay
( ) = −y( )/r3
vx
(
2
+ ) = vx
(
2
) + ax
( )
vy
(
2
+ ) = vy
(
2
) + ay
( )

21. Leapfrog Integration t = n
x(t + ) = x(t) + vx
(t +
2
)
y(t + ) = y(t) + vy
(t +
2
)
r = (x(t + ))2 + (y(t + ))2
ax
(t + ) = −x(t + )/r3
ay
(t + ) = −y(t + )/r3
vx
(t +
2
+ ) = vx
(t +
2
) + ax
(t + )
vy
(t +
2
+ ) = vy
(t +
2
) + ay
(t + )

22. 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

23. 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.

24. Feynman’s solution

25. 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!

26. References
The Feynman Lectures on Physics
Volume 1 Chapter 9
www.feynmanlectures.caltech.edu
Numpy + Scipy + Matplotlib
Anaconda Python/ Canopy Python

27. 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

28. 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

29. Code
#initialize starting values
x[0] = x_initial
y[0] = y_initial
r[0] = math.sqrt(x[0]*x[0]+y[0]*y[0])
r3 = r[0]*r[0]*r[0]
ax[0] = -x[0]/r3
ay[0] = -y[0]/r3
vx[0] = vx_initial + (epsilon/2.0)*ax[0]
vy[0] = vy_initial + (epsilon/2.0)*ay[0]

30. Code
#x[i], y[i], ax[i], ay[i] correspond to same time.
#vx[i] and vy[i] correspond to time
#offset by +epsilon/2 . (Leapfrog integration)
for i in range(1,MAX_STEPS):
x[i] = x[i-1] + epsilon*vx[i-1]
y[i] = y[i-1] + epsilon*vy[i-1]
r[i] = math.sqrt(x[i]*x[i]+y[i]*y[i])
r3 = r[i]*r[i]*r[i]
ax[i] = -x[i]/r3
ay[i] = -y[i]/r3
vx[i] = vx[i-1] + epsilon*ax[i]
vy[i] = vy[i-1] + epsilon*ay[i]

31. Code
#use matplotlib to plot the planet orbit
plt.plot(x,y)
plt.grid()
plt.xlim(-2.0,2.0)
plt.ylim(-2.0,2.0)
plt.show()