# Modeling and Simulation in Python

Comparison of the penny models from Chapters 1, 20, and 21

Copyright 2018 Allen Downey

``````

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 *

``````

### With air resistance

Next we'll add air resistance using the drag equation

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

``````

In [2]:

m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram

``````
``````

Out[2]:

kilogram

``````

Now I'll create a `Params` object to contain the quantities we need. Using a Params object is convenient for grouping the system parameters in a way that's easy to read (and double-check).

``````

In [3]:

params = Params(height = 381 * m,
v_init = 0 * m / s,
g = 9.8 * m/s**2,
mass = 2.5e-3 * kg,
diameter = 19e-3 * m,
rho = 1.2 * kg/m**3,
v_term = 18 * m / s)

``````
``````

Out[3]:

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

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

.dataframe thead th {
text-align: right;
}

values

height
381 meter

v_init
0.0 meter / second

g
9.8 meter / second ** 2

mass
0.0025 kilogram

diameter
0.019 meter

rho
1.2 kilogram / meter ** 3

v_term
18.0 meter / second

``````

Now we can pass the `Params` object `make_system` which computes some additional parameters and defines `init`.

`make_system` uses the given radius to compute `area` and the given `v_term` to compute the drag coefficient `C_d`.

``````

In [4]:

def make_system(params):
"""Makes a System object for the given conditions.

params: Params object

returns: System object
"""
unpack(params)

area = np.pi * (diameter/2)**2
C_d = 2 * mass * g / (rho * area * v_term**2)
init = State(y=height, v=v_init)
t_end = 30 * s

return System(params, area=area, C_d=C_d,
init=init, t_end=t_end)

``````

Let's make a `System`

``````

In [5]:

system = make_system(params)

``````
``````

Out[5]:

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

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

.dataframe thead th {
text-align: right;
}

values

height
381 meter

v_init
0.0 meter / second

g
9.8 meter / second ** 2

mass
0.0025 kilogram

diameter
0.019 meter

rho
1.2 kilogram / meter ** 3

v_term
18.0 meter / second

area
0.0002835287369864788 meter ** 2

C_d
0.4445009981135434 dimensionless

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

t_end
30 second

``````

Here's the slope function, including acceleration due to gravity and drag.

``````

In [6]:

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

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

returns: derivatives of y and v
"""
y, v = state
rho, C_d, area = system.rho, system.C_d, system.area
mass = system.mass
g = system.g

f_drag = rho * v**2 * C_d * area / 2
a_drag = f_drag / mass

dydt = v
dvdt = -g + a_drag

return dydt, dvdt

``````

As always, let's test the slope function with the initial conditions.

``````

In [7]:

slope_func(system.init, 0, system)

``````
``````

Out[7]:

(0.0 <Unit('meter / second')>, -9.8 <Unit('meter / second ** 2')>)

``````

We can use the same event function as in the previous chapter.

``````

In [8]:

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

``````

And then run the simulation.

``````

In [9]:

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

``````
``````

Out[9]:

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

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

.dataframe thead th {
text-align: right;
}

values

success
True

message
A termination event occurred.

``````

Here are the results.

``````

In [10]:

results

``````
``````

Out[10]:

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

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

.dataframe thead th {
text-align: right;
}

y
v

0.000000
381 meter
0.0 meter / second

0.300000
380.559 meter
-2.913855777777778 meter / second

0.600000
379.2553998560887 meter
-5.676321799192868 meter / second

0.900000
377.15535917269835 meter
-8.166374221525693 meter / second

1.200000
374.35521895425103 meter
-10.311720070623485 meter / second

1.500000
370.96543201556204 meter
-12.090247228814045 meter / second

1.800000
367.09631700871336 meter
-13.518952241005275 meter / second

2.100000
362.8483908201627 meter
-14.638370780659383 meter / second

2.400000
358.3075410597563 meter
-15.498718976734299 meter / second

2.700000
353.54387826134905 meter
-16.150356583635052 meter / second

3.000000
348.6127953660127 meter
-16.638537486499228 meter / second

3.300000
343.5570453854738 meter
-17.001300692432412 meter / second

3.600000
338.4090764843133 meter
-17.269254229780273 meter / second

3.900000
333.1932204915177 meter
-17.466306425809933 meter / second

4.200000
327.92756526232625 meter
-17.610750940510915 meter / second

4.500000
322.62547300477706 meter
-17.716383148064374 meter / second

4.800000
317.2967703203745 meter
-17.793499099432996 meter / second

5.100000
311.94866008776336 meter
-17.84972641814655 meter / second

5.400000
306.586409493623 meter
-17.890685845112163 meter / second

5.700000
301.21386361121205 meter
-17.920503377109643 meter / second

6.000000
295.8338258653892 meter
-17.94219936904363 meter / second

6.300000
290.448338371113 meter
-17.957980371929946 meter / second

6.600000
285.05888770100324 meter
-17.96945605745975 meter / second

6.900000
279.66655550040554 meter
-17.977799435972766 meter / second

7.200000
274.2721285128205 meter
-17.9838646561385 meter / second

7.500000
268.87617883849407 meter
-17.988273335755586 meter / second

7.800000
263.47912241839214 meter
-17.99147768251337 meter / second

8.100000
258.08126161893864 meter
-17.99380656697738 meter / second

8.400000
252.68281622283766 meter
-17.99549911112701 meter / second

8.700000
247.28394597351817 meter
-17.996729153771838 meter / second

...
...
...

13.800000
155.48693674391328 meter
-17.99998562459072 meter / second

14.100000
150.08694035214128 meter
-17.99998955354878 meter / second

14.400000
144.6869429742007 meter
-17.99999240867943 meter / second

14.700000
139.28694487962224 meter
-17.999994483471404 meter / second

15.000000
133.88694626427096 meter
-17.999995991199814 meter / second

15.300000
128.48694727047982 meter
-17.99999708684937 meter / second

15.600000
123.08694800168064 meter
-17.99999788304576 meter / second

15.900000
117.68694853303616 meter
-17.999998461632856 meter / second

16.200000
112.28694891916632 meter
-17.999998882085677 meter / second

16.500000
106.88694919976282 meter
-17.999999187624077 meter / second

16.800000
101.48694940366919 meter
-17.99999940965544 meter / second

17.100000
96.08694955184568 meter
-17.999999571003165 meter / second

17.400000
90.68694965952389 meter
-17.999999688252768 meter / second

17.700000
85.28694973777245 meter
-17.99999977345675 meter / second

18.000000
79.8869497946348 meter
-17.999999835373536 meter / second

18.300000
74.48694983595604 meter
-17.999999880367778 meter / second

18.600000
69.08694986598373 meter
-17.999999913064592 meter / second

18.900000
63.686949887804516 meter
-17.999999936825006 meter / second

19.200000
58.28694990366144 meter
-17.999999954091432 meter / second

19.500000
52.88694991518449 meter
-17.999999966638754 meter / second

19.800000
47.48694992355816 meter
-17.999999975756754 meter / second

20.100000
42.08694992964321 meter
-17.999999982382704 meter / second

20.400000
36.686949934065154 meter
-17.999999987197707 meter / second

20.700000
31.28694993727853 meter
-17.999999990696715 meter / second

21.000000
25.886949939613654 meter
-17.999999993239406 meter / second

21.300000
20.486949941310563 meter
-17.99999999508715 meter / second

21.600000
15.086949942543686 meter
-17.99999999642989 meter / second

21.900000
9.686949943439785 meter
-17.99999999740564 meter / second

22.200000
4.286949944090969 meter
-17.999999998114706 meter / second

22.438164
0.0 meter
-17.99999999852377 meter / second

76 rows × 2 columns

``````

The final height is close to 0, as expected.

Interestingly, the final velocity is not exactly terminal velocity, which suggests that there are some numerical errors.

We can get the flight time from `results`.

``````

In [11]:

t_sidewalk = get_last_label(results)

``````
``````

Out[11]:

22.438163885803732

``````

Here's the plot of position as a function of time.

``````

In [13]:

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

plot_position(results)

``````
``````

``````

And velocity as a function of time:

``````

In [14]:

def plot_velocity(results):
plot(results.v, color='C1', label='v')

decorate(xlabel='Time (s)',
ylabel='Velocity (m/s)')

plot_velocity(results)

``````
``````

``````

From an initial velocity of 0, the penny accelerates downward until it reaches terminal velocity; after that, velocity is constant.

### Back to Chapter 1

We have now considered three models of a falling penny:

1. In Chapter 1, we started with the simplest model, which includes gravity and ignores drag.

2. As an exercise in Chapter 1, we did a "back of the envelope" calculation assuming constant acceleration until the penny reaches terminal velocity, and then constant velocity until it reaches the sidewalk.

3. In this chapter, we model the interaction of gravity and drag during the acceleration phase.

We can compare the models by plotting velocity as a function of time.

``````

In [15]:

g = 9.8
v_term = 18
t_end = 22.4

``````
``````

Out[15]:

22.4

``````
``````

In [16]:

ts = linspace(0, t_end, 201)
model1 = -g * ts;

``````
``````

In [17]:

model2 = TimeSeries()
for t in ts:
v = -g * t
if v < -v_term:
model2[t] = -v_term
else:
model2[t] = v

``````
``````

In [18]:

results, details = run_ode_solver(system, slope_func, events=event_func, max_step=0.5)
model3 = results.v;

``````
``````

In [19]:

plot(ts, model1, label='model1', color='gray')
plot(model2, label='model2', color='C0')
plot(model3, label='model3', color='C1')

decorate(xlabel='Time (s)',
ylabel='Velocity (m/s)')

``````
``````

``````
``````

In [20]:

plot(model2, label='model2', color='C0')
plot(results.v, label='model3', color='C1')

decorate(xlabel='Time (s)',
ylabel='Velocity (m/s)')

``````
``````

``````

Clearly Model 1 is very different from the other two, which are almost identical except during the transition from constant acceleration to constant velocity.

We can also compare the predictions:

• According to Model 1, the penny takes 8.8 seconds to reach the sidewalk, and lands at 86 meters per second.

• According to Model 2, the penny takes 22.1 seconds and lands at terminal velocity, 18 m/s.

• According to Model 3, the penny takes 22.4 seconds and lands at terminal velocity.

So what can we conclude? The results from Model 1 are clearly unrealistic; it is probably not a useful model of this system. The results from Model 2 are off by about 1%, which is probably good enough for most purposes.

In fact, our estimate of the terminal velocity could by off by 10% or more, and the figure we are using for the height of the Empire State Building is not precise either.

So the difference between Models 2 and 3 is swamped by other uncertainties.

``````

In [ ]:

``````