# Modeling and Simulation in Python

Experiments with different ODE solvers

Copyright 2019 Allen Downey



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]:

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

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

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

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]:

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

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

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

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]:

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

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

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

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]:

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

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

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

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]:

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

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

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

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]:

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

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

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

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]:

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

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

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

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]:

.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.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]:

.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.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]:

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

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

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

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]:

.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.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.