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
   

   View Slide

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.
   

   View Slide

3. Ellipse
   

   View Slide

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?
   

   View Slide

5. Feynman
   

   View Slide

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

   View Slide

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

   View Slide

8. Slide Rule
   

   View Slide

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

   View Slide

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
    

    View Slide

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

    View Slide

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
    

    View Slide

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
    

    View Slide

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

    View Slide

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)
    

    View Slide

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)
    

    View Slide

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.
    

    View Slide

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

    View Slide

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)
    

    View Slide

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

    View Slide

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

    View Slide

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
    

    View Slide

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.
    

    View Slide

24. Feynman’s solution
    

    View Slide

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!
    

    View Slide

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

    View Slide

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
    

    View Slide

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
    

    View Slide

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]
    

    View Slide

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]
    

    View Slide

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

    View Slide