Upgrade to Pro — share decks privately, control downloads, hide ads and more …

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

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/

Eec9d25835717f1f1f12a354faf68d87?s=128

PyCon 2016

May 29, 2016
Tweet

More Decks by PyCon 2016

Other Decks in Programming

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()