test_euler


Modeling and Simulation in Python

Experiments with different ODE solvers

Copyright 2019 Allen Downey

License: Creative Commons Attribution 4.0 International


In [2]:
# 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 [15]:
init = State(y = 2)
system = System(init=init, t_0=1, t_end=3)


Out[15]:
values
init y 2 dtype: int64
t_0 1
t_end 3

In [16]:
def slope_func(state, t, system):
    [y] = state
    dydt = y + t
    return [dydt]

In [17]:
results, details = run_euler(system, slope_func)

In [18]:
get_last_value(results.y)


Out[18]:
24.973714732487576

Glucose minimal model

Read the data.


In [2]:
data = pd.read_csv('data/glucose_insulin.csv', index_col='time');

Interpolate the insulin data.


In [3]:
I = interpolate(data.insulin)


Out[3]:
<scipy.interpolate.interpolate.interp1d at 0x7f0b8f890908>

Initialize the parameters


In [4]:
G0 = 290
k1 = 0.03
k2 = 0.02
k3 = 1e-05


Out[4]:
1e-05

To estimate basal levels, we'll use the concentrations at t=0.


In [5]:
Gb = data.glucose[0]
Ib = data.insulin[0]


Out[5]:
11

Create the initial condtions.


In [6]:
init = State(G=G0, X=0)


Out[6]:
values
G 290
X 0

Make the System object.


In [7]:
t_0 = get_first_label(data)
t_end = get_last_label(data)


Out[7]:
182

In [8]:
system = System(G0=G0, k1=k1, k2=k2, k3=k3,
                init=init, Gb=Gb, Ib=Ib, I=I,
                t_0=t_0, t_end=t_end, dt=2)


Out[8]:
values
G0 290
k1 0.03
k2 0.02
k3 1e-05
init G 290 X 0 dtype: int64
Gb 92
Ib 11
I <scipy.interpolate.interpolate.interp1d object...
t_0 0
t_end 182
dt 2

In [9]:
def update_func(state, t, system):
    """Updates the glucose minimal model.
    
    state: State object
    t: time in min
    system: System object
    
    returns: State object
    """
    G, X = state
    k1, k2, k3 = system.k1, system.k2, system.k3 
    I, Ib, Gb = system.I, system.Ib, system.Gb
    dt = system.dt
    
    dGdt = -k1 * (G - Gb) - X*G
    dXdt = k3 * (I(t) - Ib) - k2 * X
    
    G += dGdt * dt
    X += dXdt * dt

    return State(G=G, X=X)

In [10]:
def run_simulation(system, update_func):
    """Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    """
    t_0, t_end, dt = system.t_0, system.t_end, system.dt
    
    frame = TimeFrame(columns=init.index)
    frame.row[t_0] = init
    ts = linrange(t_0, t_end, dt)
    
    for t in ts:
        frame.row[t+dt] = update_func(frame.row[t], t, system)
    
    return frame

In [11]:
%time results = run_simulation(system, update_func);


CPU times: user 166 ms, sys: 4.08 ms, total: 170 ms
Wall time: 167 ms

In [12]:
results


Out[12]:
G X
0 290 0
2 278.12 0
4 266.953 0.0003
6 256.295 0.002668
8 245.07 0.00404128
10 233.905 0.00467963
12 223.202 0.00525244
14 212.985 0.00572235
16 203.288 0.00609345
18 194.133 0.00632971
20 185.548 0.00648986
22 177.527 0.00661026
24 170.048 0.00672585
26 163.078 0.00681282
28 156.591 0.00687231
30 150.563 0.00692941
32 144.963 0.00700824
34 139.753 0.00710791
36 134.901 0.00717159
38 130.392 0.00720073
40 126.211 0.0071967
42 122.342 0.00716083
44 118.769 0.0070944
46 115.478 0.00700262
48 112.452 0.00688652
50 109.676 0.00674706
52 107.135 0.00658517
54 104.816 0.00640177
56 102.705 0.0062257
58 100.784 0.00605667
... ... ...
124 86.3907 0.00109474
126 86.5381 0.000972951
128 86.6974 0.000858033
130 86.8668 0.000749712
132 87.0445 0.000647723
134 87.2291 0.000551815
136 87.4191 0.000461742
138 87.6132 0.000377272
140 87.8103 0.000298181
142 88.0093 0.000224254
144 88.2093 0.000155284
146 88.4093 8.90726e-05
148 88.609 2.55097e-05
150 88.808 -3.55107e-05
152 89.0058 -9.40903e-05
154 89.2022 -0.000150327
156 89.3969 -0.000204314
158 89.5896 -0.000256141
160 89.7801 -0.000305895
162 89.9682 -0.00035366
164 90.1538 -0.000399513
166 90.3366 -0.000445533
168 90.5169 -0.000491711
170 90.6949 -0.000538043
172 90.8708 -0.000584521
174 91.0448 -0.00063114
176 91.217 -0.000677895
178 91.3877 -0.000724779
180 91.5569 -0.000771788
182 91.7248 -0.000818916

92 rows × 2 columns

Numerical solution

In the previous chapter, we approximated the differential equations with difference equations, and solved them using run_simulation.

In this chapter, we solve the differential equation numerically using run_euler...

Instead of an update function, we provide a slope function that evaluates the right-hand side of the differential equations. We don't have to do the update part; the solver does it for us.


In [13]:
def slope_func(state, t, system):
    """Computes derivatives of the glucose minimal model.
    
    state: State object
    t: time in min
    system: System object
    
    returns: derivatives of G and X
    """
    G, X = state
    k1, k2, k3 = system.k1, system.k2, system.k3 
    I, Ib, Gb = system.I, system.Ib, system.Gb
    
    dGdt = -k1 * (G - Gb) - X*G
    dXdt = k3 * (I(t) - Ib) - k2 * X
    
    return dGdt, dXdt

We can test the slope function with the initial conditions.


In [14]:
slope_func(init, 0, system)


Out[14]:
(-5.9399999999999995, 0.0)

Here's how we run the ODE solver.


In [15]:
system = System(G0=G0, k1=k1, k2=k2, k3=k3,
                init=init, Gb=Gb, Ib=Ib, I=I,
                t_0=t_0, t_end=t_end, dt=1)

%time results2, details = run_euler(system, slope_func)


CPU times: user 276 ms, sys: 368 µs, total: 277 ms
Wall time: 274 ms

results is a TimeFrame with one row for each time step and one column for each state variable:


In [16]:
results2


Out[16]:
G X
0 290 0
1 284.06 0
2 278.298 7.5e-05
3 272.688 0.0002235
4 267.207 0.00088903
5 261.713 0.00206125
6 256.082 0.00298502
7 250.395 0.00366532
8 244.726 0.00416202
9 239.125 0.00447878
10 233.641 0.0047792
11 228.275 0.00506362
12 223.031 0.00532235
13 217.913 0.0055559
14 212.925 0.00576478
15 208.069 0.00594948
16 203.349 0.0061005
17 198.768 0.00621849
18 194.329 0.00631745
19 190.032 0.00639777
20 185.875 0.00645981
21 181.858 0.00652061
22 177.976 0.0065802
23 174.226 0.0066386
24 170.603 0.00668983
25 167.103 0.00673403
26 163.725 0.00677135
27 160.465 0.00680192
28 157.319 0.00682588
29 154.286 0.00685537
... ... ...
153 89.2245 -0.000116549
154 89.3181 -0.000144218
155 89.4115 -0.000171334
156 89.5045 -0.000197907
157 89.597 -0.000223949
158 89.6892 -0.00024947
159 89.7809 -0.000274481
160 89.8721 -0.000298991
161 89.9628 -0.000323011
162 90.053 -0.000346551
163 90.1426 -0.00036962
164 90.2316 -0.000392728
165 90.3201 -0.000415873
166 90.4081 -0.000439056
167 90.4955 -0.000462274
168 90.5825 -0.000485529
169 90.669 -0.000508818
170 90.7551 -0.000532142
171 90.8407 -0.000555499
172 90.926 -0.000578889
173 91.0108 -0.000602311
174 91.0953 -0.000625765
175 91.1795 -0.00064925
176 91.2633 -0.000672765
177 91.3468 -0.00069631
178 91.43 -0.000719883
179 91.5129 -0.000743486
180 91.5955 -0.000767116
181 91.6779 -0.000790774
182 91.7601 -0.000814458

183 rows × 2 columns

Plotting the results from run_simulation and run_euler, we can see that they are not very different.


In [17]:
plot(results.G, '-')
plot(results2.G, '-')
plot(data.glucose, 'bo')


Out[17]:
[<matplotlib.lines.Line2D at 0x7f0b8f7954e0>]

The differences in G are less than 1%.


In [18]:
diff = results.G - results2.G
percent_diff = diff / results2.G * 100

max(abs(percent_diff.dropna()))


Out[18]:
0.970403647921524

In [ ]:


In [ ]:

Dropping pennies

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


In [19]:
m = UNITS.meter
s = UNITS.second


Out[19]:
second

And defining the initial state.


In [20]:
init = State(y=381 * m, 
             v=0 * m/s)


Out[20]:
values
y 381 meter
v 0.0 meter / second

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


In [21]:
g = 9.8 * m/s**2


Out[21]:
9.8 meter/second2

When we call odeint, we need an array of timestamps where we want to compute the solution.

I'll start with a duration of 10 seconds.


In [22]:
t_end = 10 * s


Out[22]:
10 second

Now we make a System object.


In [23]:
system = System(init=init, g=g, t_end=t_end)


Out[23]:
values
init y 381 meter v 0.0 meter / secon...
g 9.8 meter / second ** 2
t_end 10 second

And define the slope function.


In [24]:
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 [25]:
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_euler


In [26]:
system.set(dt=0.1*s)
results, details = run_euler(system, slope_func, max_step=0.5)
details.message


Out[26]:
'Success'

In [27]:
results


Out[27]:
y v
0.0 381 meter 0.0 meter / second
0.1 381.0 meter -0.9800000000000001 meter / second
0.2 380.902 meter -1.9600000000000002 meter / second
0.3 380.70599999999996 meter -2.9400000000000004 meter / second
0.4 380.412 meter -3.9200000000000004 meter / second
0.5 380.02 meter -4.9 meter / second
0.6 379.53 meter -5.880000000000001 meter / second
0.7 378.94199999999995 meter -6.860000000000001 meter / second
0.8 378.256 meter -7.840000000000002 meter / second
0.9 377.472 meter -8.820000000000002 meter / second
1.0 376.59 meter -9.800000000000002 meter / second
1.1 375.60999999999996 meter -10.780000000000003 meter / second
1.2 374.532 meter -11.760000000000003 meter / second
1.3 373.356 meter -12.740000000000004 meter / second
1.4 372.082 meter -13.720000000000004 meter / second
1.5 370.71 meter -14.700000000000005 meter / second
1.6 369.23999999999995 meter -15.680000000000005 meter / second
1.7 367.67199999999997 meter -16.660000000000004 meter / second
1.8 366.006 meter -17.640000000000004 meter / second
1.9 364.24199999999996 meter -18.620000000000005 meter / second
2.0 362.37999999999994 meter -19.600000000000005 meter / second
2.1 360.41999999999996 meter -20.580000000000005 meter / second
2.2 358.36199999999997 meter -21.560000000000006 meter / second
2.3 356.20599999999996 meter -22.540000000000006 meter / second
2.4 353.95199999999994 meter -23.520000000000007 meter / second
2.5 351.59999999999997 meter -24.500000000000007 meter / second
2.6 349.15 meter -25.480000000000008 meter / second
2.7 346.602 meter -26.460000000000008 meter / second
2.8 343.95599999999996 meter -27.44000000000001 meter / second
2.9 341.21199999999993 meter -28.42000000000001 meter / second
... ... ...
7.1 137.4700000000001 meter -69.57999999999993 meter / second
7.2 130.5120000000001 meter -70.55999999999993 meter / second
7.3 123.45600000000012 meter -71.53999999999994 meter / second
7.4 116.30200000000012 meter -72.51999999999994 meter / second
7.5 109.05000000000013 meter -73.49999999999994 meter / second
7.6 101.70000000000013 meter -74.47999999999995 meter / second
7.7 94.25200000000014 meter -75.45999999999995 meter / second
7.8 86.70600000000015 meter -76.43999999999996 meter / second
7.9 79.06200000000015 meter -77.41999999999996 meter / second
8.0 71.32000000000016 meter -78.39999999999996 meter / second
8.1 63.48000000000017 meter -79.37999999999997 meter / second
8.2 55.54200000000017 meter -80.35999999999997 meter / second
8.3 47.50600000000017 meter -81.33999999999997 meter / second
8.4 39.37200000000017 meter -82.31999999999998 meter / second
8.5 31.14000000000017 meter -83.29999999999998 meter / second
8.6 22.810000000000173 meter -84.27999999999999 meter / second
8.7 14.382000000000174 meter -85.25999999999999 meter / second
8.8 5.856000000000174 meter -86.24 meter / second
8.9 -2.7679999999998266 meter -87.22 meter / second
9.0 -11.489999999999826 meter -88.2 meter / second
9.1 -20.309999999999825 meter -89.18 meter / second
9.2 -29.227999999999824 meter -90.16000000000001 meter / second
9.3 -38.24399999999983 meter -91.14000000000001 meter / second
9.4 -47.357999999999834 meter -92.12000000000002 meter / second
9.5 -56.56999999999984 meter -93.10000000000002 meter / second
9.6 -65.87999999999984 meter -94.08000000000003 meter / second
9.7 -75.28799999999984 meter -95.06000000000003 meter / second
9.8 -84.79399999999984 meter -96.04000000000003 meter / second
9.9 -94.39799999999984 meter -97.02000000000004 meter / second
10.0 -104.09999999999985 meter -98.00000000000004 meter / second

101 rows × 2 columns


In [28]:
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
    """
    units = get_units(series.values[0])
    values = magnitudes(series - value)
    interp = InterpolatedUnivariateSpline(series.index, values)
    return interp.roots()

In [29]:
t_crossings = crossings(results.y, 0)


Out[29]:
array([8.86802711])

In [30]:
system.set(dt=0.1*s)
results, details = run_ralston(system, slope_func, max_step=0.5)
details.message


Out[30]:
'Success'

In [31]:
t_crossings = crossings(results.y, 0)


Out[31]:
array([8.81788535])

Here are the results:


In [32]:
results


Out[32]:
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
0.5 379.77500000000003 meter -4.9 meter / second
0.6 379.23600000000005 meter -5.880000000000001 meter / second
0.7 378.59900000000005 meter -6.860000000000001 meter / second
0.8 377.86400000000003 meter -7.840000000000002 meter / second
0.9 377.031 meter -8.820000000000002 meter / second
1.0 376.1 meter -9.800000000000002 meter / second
1.1 375.071 meter -10.780000000000003 meter / second
1.2 373.944 meter -11.760000000000003 meter / second
1.3 372.719 meter -12.740000000000004 meter / second
1.4 371.396 meter -13.720000000000004 meter / second
1.5 369.975 meter -14.700000000000005 meter / second
1.6 368.456 meter -15.680000000000005 meter / second
1.7 366.839 meter -16.660000000000004 meter / second
1.8 365.124 meter -17.640000000000004 meter / second
1.9 363.31100000000004 meter -18.620000000000005 meter / second
2.0 361.40000000000003 meter -19.600000000000005 meter / second
2.1 359.391 meter -20.580000000000005 meter / second
2.2 357.284 meter -21.560000000000006 meter / second
2.3 355.079 meter -22.540000000000006 meter / second
2.4 352.776 meter -23.520000000000007 meter / second
2.5 350.375 meter -24.500000000000007 meter / second
2.6 347.876 meter -25.480000000000008 meter / second
2.7 345.279 meter -26.460000000000008 meter / second
2.8 342.584 meter -27.44000000000001 meter / second
2.9 339.791 meter -28.42000000000001 meter / second
... ... ...
7.1 133.99100000000016 meter -69.57999999999993 meter / second
7.2 126.98400000000017 meter -70.55999999999993 meter / second
7.3 119.87900000000018 meter -71.53999999999994 meter / second
7.4 112.67600000000019 meter -72.51999999999994 meter / second
7.5 105.3750000000002 meter -73.49999999999994 meter / second
7.6 97.9760000000002 meter -74.47999999999995 meter / second
7.7 90.4790000000002 meter -75.45999999999995 meter / second
7.8 82.8840000000002 meter -76.43999999999996 meter / second
7.9 75.1910000000002 meter -77.41999999999996 meter / second
8.0 67.4000000000002 meter -78.39999999999996 meter / second
8.1 59.51100000000021 meter -79.37999999999997 meter / second
8.2 51.524000000000214 meter -80.35999999999997 meter / second
8.3 43.43900000000022 meter -81.33999999999997 meter / second
8.4 35.25600000000022 meter -82.31999999999998 meter / second
8.5 26.97500000000022 meter -83.29999999999998 meter / second
8.6 18.596000000000224 meter -84.27999999999999 meter / second
8.7 10.119000000000225 meter -85.25999999999999 meter / second
8.8 1.5440000000002243 meter -86.24 meter / second
8.9 -7.128999999999776 meter -87.22 meter / second
9.0 -15.899999999999777 meter -88.2 meter / second
9.1 -24.768999999999778 meter -89.18 meter / second
9.2 -33.73599999999978 meter -90.16000000000001 meter / second
9.3 -42.80099999999978 meter -91.14000000000001 meter / second
9.4 -51.963999999999785 meter -92.12000000000002 meter / second
9.5 -61.22499999999979 meter -93.10000000000002 meter / second
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

101 rows × 2 columns

And here's position as a function of time:


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

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


Saving figure to file figs/chap09-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 [34]:
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
    """
    units = get_units(series.values[0])
    values = magnitudes(series - value)
    interp = InterpolatedUnivariateSpline(series.index, values)
    return interp.roots()

In [35]:
t_crossings = crossings(results.y, 0)


Out[35]:
array([8.81788535])

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


In [36]:
t_sidewalk = t_crossings[0] * s


Out[36]:
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 [37]:
sqrt(2 * init.y / g)


Out[37]:
8.817885349720552 second

The estimate is accurate to about 10 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_ralston provides exactly the tool we need, event functions.

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


In [38]:
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_ralston. The solver should run until the event function returns 0, and then terminate.


In [39]:
results, details = run_ralston(system, slope_func, events=event_func)
details


Out[39]:
values
message Success

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

Here are the results:


In [40]:
results


Out[40]:
y v
0.000000 381 meter 0.0 meter / second
0.100000 380.951 meter -0.9800000000000001 meter / second
0.200000 380.80400000000003 meter -1.9600000000000002 meter / second
0.300000 380.559 meter -2.9400000000000004 meter / second
0.400000 380.216 meter -3.9200000000000004 meter / second
0.500000 379.77500000000003 meter -4.9 meter / second
0.600000 379.23600000000005 meter -5.880000000000001 meter / second
0.700000 378.59900000000005 meter -6.860000000000001 meter / second
0.800000 377.86400000000003 meter -7.840000000000002 meter / second
0.900000 377.031 meter -8.820000000000002 meter / second
1.000000 376.1 meter -9.800000000000002 meter / second
1.100000 375.071 meter -10.780000000000003 meter / second
1.200000 373.944 meter -11.760000000000003 meter / second
1.300000 372.719 meter -12.740000000000004 meter / second
1.400000 371.396 meter -13.720000000000004 meter / second
1.500000 369.975 meter -14.700000000000005 meter / second
1.600000 368.456 meter -15.680000000000005 meter / second
1.700000 366.839 meter -16.660000000000004 meter / second
1.800000 365.124 meter -17.640000000000004 meter / second
1.900000 363.31100000000004 meter -18.620000000000005 meter / second
2.000000 361.40000000000003 meter -19.600000000000005 meter / second
2.100000 359.391 meter -20.580000000000005 meter / second
2.200000 357.284 meter -21.560000000000006 meter / second
2.300000 355.079 meter -22.540000000000006 meter / second
2.400000 352.776 meter -23.520000000000007 meter / second
2.500000 350.375 meter -24.500000000000007 meter / second
2.600000 347.876 meter -25.480000000000008 meter / second
2.700000 345.279 meter -26.460000000000008 meter / second
2.800000 342.584 meter -27.44000000000001 meter / second
2.900000 339.791 meter -28.42000000000001 meter / second
... ... ...
6.000000 204.60000000000005 meter -58.799999999999926 meter / second
6.100000 198.67100000000005 meter -59.77999999999992 meter / second
6.200000 192.64400000000006 meter -60.75999999999992 meter / second
6.300000 186.51900000000006 meter -61.73999999999992 meter / second
6.400000 180.29600000000008 meter -62.719999999999914 meter / second
6.500000 173.97500000000008 meter -63.69999999999991 meter / second
6.600000 167.5560000000001 meter -64.67999999999991 meter / second
6.700000 161.0390000000001 meter -65.65999999999991 meter / second
6.800000 154.42400000000012 meter -66.63999999999992 meter / second
6.900000 147.71100000000013 meter -67.61999999999992 meter / second
7.000000 140.90000000000015 meter -68.59999999999992 meter / second
7.100000 133.99100000000016 meter -69.57999999999993 meter / second
7.200000 126.98400000000017 meter -70.55999999999993 meter / second
7.300000 119.87900000000018 meter -71.53999999999994 meter / second
7.400000 112.67600000000019 meter -72.51999999999994 meter / second
7.500000 105.3750000000002 meter -73.49999999999994 meter / second
7.600000 97.9760000000002 meter -74.47999999999995 meter / second
7.700000 90.4790000000002 meter -75.45999999999995 meter / second
7.800000 82.8840000000002 meter -76.43999999999996 meter / second
7.900000 75.1910000000002 meter -77.41999999999996 meter / second
8.000000 67.4000000000002 meter -78.39999999999996 meter / second
8.100000 59.51100000000021 meter -79.37999999999997 meter / second
8.200000 51.524000000000214 meter -80.35999999999997 meter / second
8.300000 43.43900000000022 meter -81.33999999999997 meter / second
8.400000 35.25600000000022 meter -82.31999999999998 meter / second
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

90 rows × 2 columns

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 [41]:
t_sidewalk = get_last_label(results) * s


Out[41]:
8.81780237518735 second

The result is accurate to about 15 decimal places.

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


In [42]:
v_sidewalk = get_last_value(results.v)


Out[42]:
-86.41446327683617 meter/second

And convert to kilometers per hour.


In [43]:
km = UNITS.kilometer
h = UNITS.hour
v_sidewalk.to(km / h)


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