# Modeling and Simulation in Python

Chapter 20



In :

# 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 *



### Dropping pennies

I'll start by getting the units we need from Pint.



In :

m = UNITS.meter
s = UNITS.second




Out:

second



And defining the initial state.



In :

init = State(y=381 * m,
v=0 * m/s)




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

y
381 meter

v
0.0 meter / second



Acceleration due to gravity is about 9.8 m / s$^2$.



In :

g = 9.8 * m/s**2




Out:

9.8 meter/second2



I'll start with a duration of 10 seconds and step size 0.1 second.



In :

t_end = 10 * s




Out:

10 second




In :

dt = 0.1 * s




Out:

0.1 second



Now we make a System object.



In :

system = System(init=init, g=g, t_end=t_end, dt=dt)




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

init
y             381 meter
v    0.0 meter / secon...

g
9.8 meter / second ** 2

t_end
10 second

dt
0.1 second



And define the slope function.



In :

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
"""
y, v = state
g = system.g

dydt = v
dvdt = -g

return dydt, dvdt



It's always a good idea to test the slope function with the initial conditions.



In :

dydt, dvdt = slope_func(system.init, 0, system)
print(dydt)
print(dvdt)




0.0 meter / second
-9.8 meter / second ** 2



Now we're ready to call run_ode_solver



In :

results, details = run_ode_solver(system, slope_func)
details




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

success
True

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



Here are the results:



In :




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

y
v

0.0
381 meter
0.0 meter / second

0.1
380.951 meter
-0.9800000000000001 meter / second

0.2
380.80400000000003 meter
-1.9600000000000002 meter / second

0.3
380.559 meter
-2.9400000000000004 meter / second

0.4
380.216 meter
-3.9200000000000004 meter / second




In :

results.tail()




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

y
v

9.6
-70.58399999999979 meter
-94.08000000000003 meter / second

9.7
-80.0409999999998 meter
-95.06000000000003 meter / second

9.8
-89.5959999999998 meter
-96.04000000000003 meter / second

9.9
-99.24899999999981 meter
-97.02000000000004 meter / second

10.0
-108.99999999999982 meter
-98.00000000000004 meter / second



And here's position as a function of time:



In :

def plot_position(results):
plot(results.y, label='y')
decorate(xlabel='Time (s)',
ylabel='Position (m)')

plot_position(results)
savefig('figs/chap20-fig01.pdf')




Saving figure to file figs/chap20-fig01.pdf



### Onto the sidewalk

To figure out when the penny hit the sidewalk, we can use crossings, which finds the times where a Series passes through a given value.



In :

t_crossings = crossings(results.y, 0)




Out:

array([8.81788535])



For this example there should be just one crossing, the time when the penny hits the sidewalk.



In :

t_sidewalk = t_crossings * s




Out:

8.81788534972054 second



We can compare that to the exact result. Without air resistance, we have

$v = -g t$

and

$y = 381 - g t^2 / 2$

Setting $y=0$ and solving for $t$ yields

$t = \sqrt{\frac{2 y_{init}}{g}}$



In :

sqrt(2 * init.y / g)




Out:

8.817885349720552 second



The estimate is accurate to about 9 decimal places.

## Events

Instead of running the simulation until the penny goes through the sidewalk, it would be better to detect the point where the penny hits the sidewalk and stop. run_ode_solver provides exactly the tool we need, event functions.

Here's an event function that returns the height of the penny above the sidewalk:



In :

def event_func(state, t, system):
"""Return the height of the penny above the sidewalk.
"""
y, v = state
return y



And here's how we pass it to run_ode_solver. The solver should run until the event function returns 0, and then terminate.



In :

results, details = run_ode_solver(system, slope_func, events=event_func)
details




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

success
True

message
A termination event occurred.



The message from the solver indicates the solver stopped because the event we wanted to detect happened.

Here are the results:



In :

results.tail()




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

y
v

8.500000
26.97500000000022 meter
-83.29999999999998 meter / second

8.600000
18.596000000000224 meter
-84.27999999999999 meter / second

8.700000
10.119000000000225 meter
-85.25999999999999 meter / second

8.800000
1.5440000000002243 meter
-86.24 meter / second

8.817802
-4.440892098500626e-16 meter
-86.41446327683617 meter / second



With the events option, the solver returns the actual time steps it computed, which are not necessarily equally spaced.

The last time step is when the event occurred:



In :

t_sidewalk = get_last_label(results) * s




Out:

8.81780237518735 second



The result is accurate to about 4 decimal places.

We can also check the velocity of the penny when it hits the sidewalk:



In :

v_sidewalk = get_last_value(results.v)




Out:

-86.41446327683617 meter/second



And convert to kilometers per hour.



In :

km = UNITS.kilometer
h = UNITS.hour
v_sidewalk.to(km / h)




Out:

-311.09206779661025 kilometer/hour



If there were no air resistance, the penny would hit the sidewalk (or someone's head) at more than 300 km/h.

So it's a good thing there is air resistance.

## Under the hood

Here is the source code for crossings so you can see what's happening under the hood:



In :

source_code(crossings)




def crossings(series, value):
"""Find the labels where the series passes through value.

The labels in series must be increasing numerical values.

series: Series
value: number

returns: sequence of labels
"""
values = magnitudes(series - value)
interp = InterpolatedUnivariateSpline(series.index, values)
return interp.roots()



### Exercises

Exercise: Here's a question from the web site Ask an Astronomer:

"If the Earth suddenly stopped orbiting the Sun, I know eventually it would be pulled in by the Sun's gravity and hit it. How long would it take the Earth to hit the Sun? I imagine it would go slowly at first and then pick up speed."

Use run_ode_solver to answer this question.

Here are some suggestions about how to proceed:

1. Look up the Law of Universal Gravitation and any constants you need. I suggest you work entirely in SI units: meters, kilograms, and Newtons.

2. When the distance between the Earth and the Sun gets small, this system behaves badly, so you should use an event function to stop when the surface of Earth reaches the surface of the Sun.

3. Express your answer in days, and plot the results as millions of kilometers versus days.

If you read the reply by Dave Rothstein, you will see other ways to solve the problem, and a good discussion of the modeling decisions behind them.

You might also be interested to know that it's actually not that easy to get to the Sun.



In :

# Solution

N = UNITS.newton
kg = UNITS.kilogram
m = UNITS.meter
AU = UNITS.astronomical_unit




Out:

astronomical_unit




In :

# Solution

r_0 = (1 * AU).to_base_units()
v_0 = 0 * m / s
init = State(r=r_0,
v=v_0)




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

r
149597870691.0 meter

v
0.0 meter / second




In :

# Solution

t_end = 1e7 * s
dt = t_end / 100

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




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

init
r    149597870691.0 meter
v      0.0 meter / s...

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
10000000.0 second

dt
100000.0 second




In :

# Solution

def universal_gravitation(state, system):
"""Computes gravitational force.

state: State object with distance r
system: System object with m1, m2, and G
"""
r, v = state
G, m1, m2 = system.G, system.m1, system.m2

force = G * m1 * m2 / r**2
return force




In :

# Solution

universal_gravitation(init, system)




Out:

3.5423376937972604×1022 newton




In :

# Solution

def slope_func(state, t, system):
"""Compute derivatives of the state.

state: position, velocity
t: time
system: System object containing m2

returns: derivatives of y and v
"""
y, v = state
m2 = system.m2

force = universal_gravitation(state, system)
dydt = v
dvdt = -force / m2

return dydt, dvdt




In :

# Solution

slope_func(system.init, 0, system)




Out:

(0.0 <Unit('meter / second')>,
-0.005931576848287442 <Unit('newton / kilogram')>)




In :

# Solution

def event_func(state, t, system):
r, v = state
return r - system.r_final




In :

# Solution

event_func(init, 0, system)




Out:

148895991691.0 meter




In :

# Solution

results, details = run_ode_solver(system, slope_func, events=event_func)
details




Out:

.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}

.dataframe tbody tr th {
vertical-align: top;
}

text-align: right;
}

values

success
True

message
A termination event occurred.




In :

# Solution

t_event = get_last_label(results)




Out:

5600110.3392641265




In :

# Solution

seconds = t_event * s
days = seconds.to(UNITS.day)




Out:

64.81609188963108 day




In :

# Solution

results.index /= 60 * 60 * 24




In :

# Solution

results.r /= 1e9




In :

# Solution

plot(results.r, label='r')

decorate(xlabel='Time (day)',
ylabel='Distance from sun (million km)')







In [ ]: