Modeling and Simulation in Python

In [1]:

    

# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *


    

In [2]:

    

# Here are the units we'll need

s = UNITS.second
N = UNITS.newton
kg = UNITS.kilogram
m = UNITS.meter
year = UNITS.year


    

    
Out[2]:




year

In [3]:

    

# And an inition condition (with everything in SI units)

# distance from the sun to the earth at perihelion
R = Vector(147e9, 0) * m

# initial velocity
V = Vector(0, -30300) * m/s

init = State(R=R, V=V)


    

    
Out[3]:







  

    

      

      
values
    
  
  

    

      
R
      
[147000000000.0 meter, 0.0 meter]
    
    

      
V
      
[0.0 meter / second, -30300.0 meter / second]
    
  

In [4]:

    

# Making a system object

r_earth = 6.371e6 * m
r_sun = 695.508e6 * m

t_end = (1 * year).to_base_units()

system = System(init=init,
                G=6.674e-11 * N / kg**2 * m**2,
                m1=1.989e30 * kg,
                r_final=r_sun + r_earth,
                m2=5.972e24 * kg,
                t_end=t_end)


    

    
Out[4]:







  

    

      

      
values
    
  
  

    

      
init
      
R                [147000000000.0 meter, 0.0 me...
    
    

      
G
      
6.674e-11 meter ** 2 * newton / kilogram ** 2
    
    

      
m1
      
1.989e+30 kilogram
    
    

      
r_final
      
701879000.0 meter
    
    

      
m2
      
5.972e+24 kilogram
    
    

      
t_end
      
31556925.9747 second
    
  

In [5]:

    

# Here's a function that computes the force of gravity

def universal_gravitation(state, system):
    """Computes gravitational force.
    
    state: State object with distance r
    system: System object with m1, m2, and G
    
    returns: Vector
    """
    R, V = state
    G, m1, m2 = system.G, system.m1, system.m2
        
    # make sure the result is a vector, either
    # by forcing it, or by putting v.hat() on
    # the left
    
    force = -G * m1 * m2 / R.mag2 * R.hat()
    return Vector(force)

    force = -R.hat() * G * m1 * m2 / R.mag2
    return force


    

In [6]:

    

universal_gravitation(init, system)


    

    
Out[6]:




\[\begin{pmatrix}-3.6686485997501037×1022 & -0.0\end{pmatrix} newton\]

In [7]:

    

# The slope function

def slope_func(state, t, system):
    """Compute derivatives of the state.
    
    state: position, velocity
    t: time
    system: System object containing `g`
    
    returns: derivatives of y and v
    """
    R, V = state

    F = universal_gravitation(state, system)
    A = F / system.m2
    
    return V, A


    

In [8]:

    

# Always test the slope function!

slope_func(init, 0, system)


    

    
Out[8]:





(array([     0., -30300.]) <Unit('meter / second')>,
 array([-0.00614308, -0.        ]) <Unit('newton / kilogram')>)

In [9]:

    

# Here's an event function that stops the simulation
# before the collision

def event_func(state, t, system):
    R, V = state
    return R.mag - system.r_final


    

In [10]:

    

# Always test the event function!

event_func(init, 0, system)


    

    
Out[10]:




146298121000.0 meter

In [11]:

    

# Finally we can run the simulation

system.set(dt=system.t_end/1000)
results, details = run_euler(system, slope_func)
details


    

    
Out[11]:







  

    

      

      
values
    
  
  

    

      
message
      
Success
    
  

In [12]:

    

def plot_trajectory(R):
    x = R.extract('x') / 1e9
    y = R.extract('y') / 1e9

    plot(x, y)
    plot(0, 0, 'yo')

    decorate(xlabel='x distance (million km)',
             ylabel='y distance (million km)')


    

In [13]:

    

plot_trajectory(results.R)


    

In [14]:

    

error = results.last_row() - system.init


    

    
Out[14]:





R    [1993548154.0907593 meter, 53558254525.37807 m...
V    [9671.881139414403 meter / second, 2805.422229...
dtype: object

In [15]:

    

error.R.mag


    

    
Out[15]:




53595343660.133934 meter

In [16]:

    

error.V.mag


    

    
Out[16]:




10070.535172486145 meter/second

In [17]:

    

# Ralston

system.set(dt=system.t_end/1000)
results, details = run_ralston(system, slope_func)
details


    

    
Out[17]:







  

    

      

      
values
    
  
  

    

      
success
      
True
    
    

      
message
      
The solver successfully reached the end of the...
    
  

In [18]:

    

plot_trajectory(results.R)


    

In [19]:

    

error = results.last_row() - system.init


    

    
Out[19]:





R    [-3723637.7203674316 meter, -1060434466.552573...
V    [-214.86614651603597 meter / second, 0.7785432...
dtype: object

In [20]:

    

error.R.mag


    

    
Out[20]:




1060441004.1725632 meter

In [21]:

    

error.V.mag


    

    
Out[21]:




214.8675569931363 meter/second

In [24]:

    

xs = results.R.extract('x') * 1.1
ys = results.R.extract('y') * 1.1

def draw_func(state, t):
    set_xlim(xs)
    set_ylim(ys)
    x, y = state.R
    plot(x, y, 'b.')
    plot(0, 0, 'yo')
    decorate(xlabel='x distance (million km)',
             ylabel='y distance (million km)')


    

In [25]:

    

system.set(dt=system.t_end/100)
results, details = run_ralston(system, slope_func)
animate(results, draw_func)


    

In [ ]: