Modeling and Simulation in Python

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

Copyright 2018 Allen Downey

License: Creative Commons Attribution 4.0 International


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]:
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]:
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]:
values
success True
message A termination event occurred.

Here are the results.


In [10]:
results


Out[10]:
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 [ ]: