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 *
In [2]:
m = UNITS.meter
s = UNITS.second
Out[2]:
And defining the initial state.
In [3]:
init = State(y=381 * m,
v=0 * m/s)
Out[3]:
Acceleration due to gravity is about 9.8 m / s$^2$.
In [4]:
g = 9.8 * m/s**2
Out[4]:
I'll start with a duration of 10 seconds and step size 0.1 second.
In [5]:
t_end = 10 * s
Out[5]:
In [6]:
dt = 0.1 * s
Out[6]:
Now we make a System object.
In [7]:
system = System(init=init, g=g, t_end=t_end, dt=dt)
Out[7]:
And define the slope function.
In [8]:
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 [9]:
dydt, dvdt = slope_func(system.init, 0, system)
print(dydt)
print(dvdt)
Now we're ready to call run_ode_solver
In [10]:
results, details = run_ode_solver(system, slope_func)
details
Out[10]:
Here are the results:
In [11]:
results.head()
Out[11]:
In [12]:
results.tail()
Out[12]:
And here's position as a function of time:
In [13]:
def plot_position(results):
plot(results.y, label='y')
decorate(xlabel='Time (s)',
ylabel='Position (m)')
plot_position(results)
savefig('figs/chap20-fig01.pdf')
In [14]:
t_crossings = crossings(results.y, 0)
Out[14]:
For this example there should be just one crossing, the time when the penny hits the sidewalk.
In [15]:
t_sidewalk = t_crossings[0] * s
Out[15]:
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 [16]:
sqrt(2 * init.y / g)
Out[16]:
The estimate is accurate to about 9 decimal places.
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 [17]:
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 [18]:
results, details = run_ode_solver(system, slope_func, events=event_func)
details
Out[18]:
The message from the solver indicates the solver stopped because the event we wanted to detect happened.
Here are the results:
In [19]:
results.tail()
Out[19]:
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 [20]:
t_sidewalk = get_last_label(results) * s
Out[20]:
The result is accurate to about 4 decimal places.
We can also check the velocity of the penny when it hits the sidewalk:
In [21]:
v_sidewalk = get_last_value(results.v)
Out[21]:
And convert to kilometers per hour.
In [22]:
km = UNITS.kilometer
h = UNITS.hour
v_sidewalk.to(km / h)
Out[22]:
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.
In [23]:
source_code(crossings)
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:
Look up the Law of Universal Gravitation and any constants you need. I suggest you work entirely in SI units: meters, kilograms, and Newtons.
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.
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 [24]:
# Solution
N = UNITS.newton
kg = UNITS.kilogram
m = UNITS.meter
AU = UNITS.astronomical_unit
Out[24]:
In [25]:
# Solution
r_0 = (1 * AU).to_base_units()
v_0 = 0 * m / s
init = State(r=r_0,
v=v_0)
Out[25]:
In [26]:
# Solution
radius_earth = 6.371e6 * m
radius_sun = 695.508e6 * m
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,
r_final=radius_sun + radius_earth,
m2=5.972e24 * kg,
t_end=t_end,
dt=dt)
Out[26]:
In [27]:
# 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 [28]:
# Solution
universal_gravitation(init, system)
Out[28]:
In [29]:
# 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 [30]:
# Solution
slope_func(system.init, 0, system)
Out[30]:
In [31]:
# Solution
def event_func(state, t, system):
r, v = state
return r - system.r_final
In [32]:
# Solution
event_func(init, 0, system)
Out[32]:
In [33]:
# Solution
results, details = run_ode_solver(system, slope_func, events=event_func)
details
Out[33]:
In [34]:
# Solution
t_event = get_last_label(results)
Out[34]:
In [35]:
# Solution
seconds = t_event * s
days = seconds.to(UNITS.day)
Out[35]:
In [36]:
# Solution
results.index /= 60 * 60 * 24
In [37]:
# Solution
results.r /= 1e9
In [38]:
# Solution
plot(results.r, label='r')
decorate(xlabel='Time (day)',
ylabel='Distance from sun (million km)')
In [ ]: