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/

Transcript

1. Computational Physics with Python: Planetary Orbits from Newton to Feynman Pramod Gupta Department of Astronomy, University of Washington
2. 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.
3. Ellipse
4. 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?
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()