chap22soln


Modeling and Simulation in Python

Chapter 22

Copyright 2017 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 *

Vectors

A Vector object represents a vector quantity. In the context of mechanics, vector quantities include position, velocity, acceleration, and force, all of which might be in 2D or 3D.

You can define a Vector object without units, but if it represents a physical quantity, you will often want to attach units to it.

I'll start by grabbing the units we'll need.


In [2]:
m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram


Out[2]:
kilogram

Here's a two dimensional Vector in meters.


In [3]:
A = Vector(3, 4) * m


Out[3]:
\[\begin{pmatrix}3.0 & 4.0\end{pmatrix} meter\]

We can access the elements by name.


In [4]:
A.x


Out[4]:
3.0 meter

In [5]:
A.y


Out[5]:
4.0 meter

The magnitude is the length of the vector.


In [6]:
A.mag


Out[6]:
5.0 meter

The angle is the number of radians between the vector and the positive x axis.


In [7]:
A.angle


Out[7]:
0.9272952180016122 radian

If we make another Vector with the same units,


In [8]:
B = Vector(1, 2) * m


Out[8]:
\[\begin{pmatrix}1.0 & 2.0\end{pmatrix} meter\]

We can add Vector objects like this


In [9]:
A + B


Out[9]:
\[\begin{pmatrix}4.0 & 6.0\end{pmatrix} meter\]

And subtract like this:


In [10]:
A - B


Out[10]:
\[\begin{pmatrix}2.0 & 2.0\end{pmatrix} meter\]

We can compute the Euclidean distance between two Vectors.


In [11]:
A.dist(B)


Out[11]:
2.8284271247461903 meter

And the difference in angle


In [12]:
A.diff_angle(B)


Out[12]:
-0.17985349979247822 radian

If we are given the magnitude and angle of a vector, what we have is the representation of the vector in polar coordinates.


In [13]:
mag = A.mag
angle = A.angle


Out[13]:
0.9272952180016122 radian

We can use pol2cart to convert from polar to Cartesian coordinates, and then use the Cartesian coordinates to make a Vector object.

In this example, the Vector we get should have the same components as A.


In [14]:
x, y = pol2cart(angle, mag)
Vector(x, y)


Out[14]:
\[\begin{pmatrix}3.0000000000000004 & 3.9999999999999996\end{pmatrix} meter\]

Another way to represent the direction of A is a unit vector, which is a vector with magnitude 1 that points in the same direction as A. You can compute a unit vector by dividing a vector by its magnitude:


In [15]:
A / A.mag


Out[15]:
\[\begin{pmatrix}0.6 & 0.8\end{pmatrix} dimensionless\]

Or by using the hat function, so named because unit vectors are conventionally decorated with a hat, like this: $\hat{A}$:


In [16]:
A.hat()


Out[16]:
\[\begin{pmatrix}0.6 & 0.8\end{pmatrix} dimensionless\]

Exercise: Create a Vector named a_grav that represents acceleration due to gravity, with x component 0 and y component $-9.8$ meters / second$^2$.


In [17]:
# Solution

a_grav = Vector(0, -9.8) * m / s**2


Out[17]:
\[\begin{pmatrix}0.0 & -9.8\end{pmatrix} meter/second2\]

Degrees and radians

Pint provides units to represent degree and radians.


In [18]:
degree = UNITS.degree
radian = UNITS.radian


Out[18]:
radian

If you have an angle in degrees,


In [19]:
angle = 45 * degree
angle


Out[19]:
45 degree

You can convert to radians.


In [20]:
angle_rad = angle.to(radian)


Out[20]:
0.7853981633974483 radian

If it's already in radians, to does the right thing.


In [21]:
angle_rad.to(radian)


Out[21]:
0.7853981633974483 radian

You can also convert from radians to degrees.


In [22]:
angle_rad.to(degree)


Out[22]:
45.0 degree

As an alterative, you can use np.deg2rad, which works with Pint quantities, but it also works with simple numbers and NumPy arrays:


In [23]:
np.deg2rad(angle)


Out[23]:
0.7853981633974483 radian

Exercise: Create a Vector named a_force that represents acceleration due to a force of 0.5 Newton applied to an object with mass 0.3 kilograms, in a direction 45 degrees up from the positive x-axis.

Add a_force to a_grav from the previous exercise. If that addition succeeds, that means that the units are compatible. Confirm that the total acceleration seems to make sense.


In [24]:
# Solution

N = UNITS.newton
mag = 0.5 * N
angle = 45 * degree
theta = angle.to(radian)
x, y = pol2cart(theta, mag)
force = Vector(x, y)

mass = 0.3 * kg
a_force = force / mass
a_force


Out[24]:
\[\begin{pmatrix}1.1785113019775793 & 1.1785113019775793\end{pmatrix} newton/kilogram\]

In [25]:
# Solution

a_force + a_grav


Out[25]:
\[\begin{pmatrix}1.1785113019775793 & -8.621488698022421\end{pmatrix} newton/kilogram\]

Baseball

Here's a Params object that contains parameters for the flight of a baseball.


In [26]:
t_end = 10 * s
dt = t_end / 100

params = Params(x = 0 * m, 
                y = 1 * m,
                g = 9.8 * m/s**2,
                mass = 145e-3 * kg,
                diameter = 73e-3 * m,
                rho = 1.2 * kg/m**3,
                C_d = 0.33,
                angle = 45 * degree,
                velocity = 40 * m / s,
                t_end=t_end, dt=dt)


Out[26]:
values
x 0 meter
y 1 meter
g 9.8 meter / second ** 2
mass 0.145 kilogram
diameter 0.073 meter
rho 1.2 kilogram / meter ** 3
C_d 0.33
angle 45 degree
velocity 40.0 meter / second
t_end 10 second
dt 0.1 second

And here's the function that uses the Params object to make a System object.


In [27]:
def make_system(params):
    """Make a system object.
    
    params: Params object with angle, velocity, x, y,
               diameter, duration, g, mass, rho, and C_d
               
    returns: System object
    """
    angle, velocity = params.angle, params.velocity
    
    # convert angle to degrees
    theta = np.deg2rad(angle)
    
    # compute x and y components of velocity
    vx, vy = pol2cart(theta, velocity)
    
    # make the initial state
    R = Vector(params.x, params.y)
    V = Vector(vx, vy)
    init = State(R=R, V=V)
    
    # compute area from diameter
    diameter = params.diameter
    area = np.pi * (diameter/2)**2
    
    return System(params, init=init, area=area)

Here's how we use it:


In [28]:
system = make_system(params)


Out[28]:
values
x 0 meter
y 1 meter
g 9.8 meter / second ** 2
mass 0.145 kilogram
diameter 0.073 meter
rho 1.2 kilogram / meter ** 3
C_d 0.33
angle 45 degree
velocity 40.0 meter / second
t_end 10 second
dt 0.1 second
init R [0 meter, ...
area 0.004185386812745002 meter ** 2

Here's a function that computes drag force using vectors:


In [29]:
def drag_force(V, system):
    """Computes drag force in the opposite direction of `v`.
    
    V: velocity Vector
    system: System object with rho, C_d, area
    
    returns: Vector drag force
    """
    rho, C_d, area = system.rho, system.C_d, system.area
    
    mag = rho * V.mag**2 * C_d * area / 2
    direction = -V.hat()
    f_drag = direction * mag
    return f_drag

We can test it like this.


In [30]:
V_test = Vector(10, 10) * m/s
drag_force(V_test, system)


Out[30]:
\[\begin{pmatrix}-0.11719680972835739 & -0.11719680972835739\end{pmatrix} kilogram meter/second2\]

Here's the slope function that computes acceleration due to gravity and drag.


In [31]:
def slope_func(state, t, system):
    """Computes derivatives of the state variables.
    
    state: State (x, y, x velocity, y velocity)
    t: time
    system: System object with g, rho, C_d, area, mass
    
    returns: sequence (vx, vy, ax, ay)
    """
    R, V = state
    mass, g = system.mass, system.g
    
    a_drag = drag_force(V, system) / mass
    a_grav = Vector(0, -g)
    
    A = a_grav + a_drag
    
    return V, A

Always test the slope function with the initial conditions.


In [32]:
slope_func(system.init, 0, system)


Out[32]:
(array([28.28427125, 28.28427125]) <Unit('meter / second')>,
 array([ -6.46603088, -16.26603088]) <Unit('meter / second ** 2')>)

We can use an event function to stop the simulation when the ball hits the ground:


In [33]:
def event_func(state, t, system):
    """Stop when the y coordinate is 0.
    
    state: State object
    t: time
    system: System object
    
    returns: y coordinate
    """
    R, V = state
    return R.y

In [34]:
event_func(system.init, 0, system)


Out[34]:
1 meter

Now we can call run_ode_solver


In [35]:
results, details = run_ode_solver(system, slope_func, events=event_func)
details


Out[35]:
values
success True
message A termination event occurred.

The final label tells us the flight time.


In [36]:
flight_time = get_last_label(results) * s


Out[36]:
5.004505488017051 second

The final value of x tells us the how far the ball landed from home plate:


In [37]:
R_final = get_last_value(results.R)
x_dist = R_final.x


Out[37]:
99.30497406350605 meter

Visualizing the results

The simplest way to visualize the results is to plot x and y as functions of time.


In [38]:
xs = results.R.extract('x')
ys = results.R.extract('y')

xs.plot()
ys.plot()

decorate(xlabel='Time (s)',
         ylabel='Position (m)')

savefig('figs/chap22-fig01.pdf')


Saving figure to file figs/chap22-fig01.pdf

We can plot the velocities the same way.


In [39]:
vx = results.V.extract('x')
vy = results.V.extract('y')

vx.plot(label='vx')
vy.plot(label='vy')

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


The x velocity slows down due to drag.

The y velocity drops quickly while drag and gravity are in the same direction, then more slowly after the ball starts to fall.

Another way to visualize the results is to plot y versus x. The result is the trajectory of the ball through its plane of motion.


In [40]:
def plot_trajectory(results):
    xs = results.R.extract('x')
    ys = results.R.extract('y')
    plot(xs, ys, color='C2', label='trajectory')

    decorate(xlabel='x position (m)',
             ylabel='y position (m)')

plot_trajectory(results)
savefig('figs/chap22-fig02.pdf')


Saving figure to file figs/chap22-fig02.pdf

Animation

One of the best ways to visualize the results of a physical model is animation. If there are problems with the model, animation can make them apparent.

The ModSimPy library provides animate, which takes as parameters a TimeSeries and a draw function.

The draw function should take as parameters a State object and the time. It should draw a single frame of the animation.

Inside the draw function, you almost always have to call set_xlim and set_ylim. Otherwise matplotlib auto-scales the axes, which is usually not what you want.


In [41]:
xs = results.R.extract('x')
ys = results.R.extract('y')

def draw_func(state, t):
    set_xlim(xs)
    set_ylim(ys)
    x, y = state.R
    plot(x, y, 'bo')
    decorate(xlabel='x position (m)',
             ylabel='y position (m)')

In [42]:
animate(results, draw_func)


Exercise: Delete the lines that set the x and y axes (or comment them out) and see what the animation does.

Under the hood

Vector is a function that returns a ModSimVector object.


In [43]:
V = Vector(3, 4)
type(V)


Out[43]:
modsim.modsim.ModSimVector

A ModSimVector is a specialized kind of Pint Quantity.


In [44]:
isinstance(V, Quantity)


Out[44]:
True

There's one gotcha you might run into with Vectors and Quantities. If you multiply a ModSimVector and a Quantity, you get a ModSimVector:


In [45]:
V1 = V * m


Out[45]:
\[\begin{pmatrix}3.0 & 4.0\end{pmatrix} meter\]

In [46]:
type(V1)


Out[46]:
modsim.modsim.ModSimVector

But if you multiply a Quantity and a Vector, you get a Quantity:


In [47]:
V2 = m * V


Out[47]:
\[\begin{pmatrix}3.0 & 4.0\end{pmatrix} meter\]

In [48]:
type(V2)


Out[48]:
pint.quantity.build_quantity_class.<locals>.Quantity

With a ModSimVector you can get the coordinates using dot notation, as well as mag, mag2, and angle:


In [49]:
V1.x, V1.y, V1.mag, V1.angle


Out[49]:
(3.0 <Unit('meter')>,
 4.0 <Unit('meter')>,
 5.0 <Unit('meter')>,
 0.9272952180016122 <Unit('radian')>)

With a Quantity, you can't. But you can use indexing to get the coordinates:


In [50]:
V2[0], V2[1]


Out[50]:
(3.0 <Unit('meter')>, 4.0 <Unit('meter')>)

And you can use vector functions to get the magnitude and angle.


In [51]:
vector_mag(V2), vector_angle(V2)


Out[51]:
(5.0 <Unit('meter')>, 0.9272952180016122 <Unit('radian')>)

And often you can avoid the whole issue by doing the multiplication with the ModSimVector on the left.

Exercises

Exercise: Run the simulation with and without air resistance. How wrong would we be if we ignored drag?


In [52]:
# Hint

system_no_drag = System(system, C_d=0)


Out[52]:
values
x 0 meter
y 1 meter
g 9.8 meter / second ** 2
mass 0.145 kilogram
diameter 0.073 meter
rho 1.2 kilogram / meter ** 3
C_d 0
angle 45 degree
velocity 40.0 meter / second
t_end 10 second
dt 0.1 second
init R [0 meter, ...
area 0.004185386812745002 meter ** 2

In [53]:
# Solution

results_no_drag, details = run_ode_solver(system_no_drag, slope_func, events=event_func)
details


Out[53]:
values
success True
message A termination event occurred.

In [54]:
# Solution

plot_trajectory(results)
plot_trajectory(results_no_drag)



In [55]:
# Solution

x_dist = get_last_value(results.R).x


Out[55]:
99.30497406350605 meter

In [56]:
# Solution

xdist_no_drag = get_last_value(results_no_drag.R).x


Out[56]:
164.25596844639244 meter

In [57]:
# Solution

xdist_no_drag - x_dist


Out[57]:
64.95099438288639 meter

Exercise: The baseball stadium in Denver, Colorado is 1,580 meters above sea level, where the density of air is about 1.0 kg / meter$^3$. How much farther would a ball hit with the same velocity and launch angle travel?


In [58]:
# Hint

system2 = System(system, rho=1.0*kg/m**3)


Out[58]:
values
x 0 meter
y 1 meter
g 9.8 meter / second ** 2
mass 0.145 kilogram
diameter 0.073 meter
rho 1.0 kilogram / meter ** 3
C_d 0.33
angle 45 degree
velocity 40.0 meter / second
t_end 10 second
dt 0.1 second
init R [0 meter, ...
area 0.004185386812745002 meter ** 2

In [59]:
# Solution

results2, details2 = run_ode_solver(system2, slope_func, events=event_func)
x = results2.R.extract('x')
x_dist2 = get_last_value(x)


Out[59]:
105.77787365390012 meter

In [60]:
# Solution

x_dist2 - x_dist


Out[60]:
6.472899590394064 meter

Exercise: The model so far is based on the assumption that coefficient of drag does not depend on velocity, but in reality it does. The following figure, from Adair, The Physics of Baseball, shows coefficient of drag as a function of velocity.

I used an online graph digitizer to extract the data and save it in a CSV file. Here's how we can read it:


In [61]:
baseball_drag = pd.read_csv('data/baseball_drag.csv')
mph = Quantity(baseball_drag['Velocity in mph'], UNITS.mph)
mps = mph.to(m/s)
baseball_drag.index = magnitude(mps)
baseball_drag.index.name = 'Velocity in meters per second'
baseball_drag


Out[61]:
Velocity in mph Drag coefficient
Velocity in meters per second
0.026146 0.058486 0.49965
8.871509 19.845000 0.49878
17.647351 39.476000 0.49704
22.432914 50.181000 0.48225
26.882303 60.134000 0.45004
30.636992 68.533000 0.40914
32.977694 73.769000 0.38042
34.604472 77.408000 0.36562
37.497268 83.879000 0.34822
40.460249 90.507000 0.33081
43.492522 97.290000 0.31427
46.733562 104.540000 0.30035
50.886563 113.830000 0.28816
54.047136 120.900000 0.28381
56.926074 127.340000 0.28033
60.086646 134.410000 0.28207

Modify the model to include the dependence of C_d on velocity, and see how much it affects the results. Hint: use interpolate.


In [62]:
# Solution

drag_interp = interpolate(baseball_drag['Drag coefficient'])
vs = linspace(0, 60)
cds = drag_interp(vs)
plot(vs, cds)
decorate(xlabel='Velocity (m/s)', ylabel='C_d')



In [63]:
# Solution

def drag_force(V, system):
    """Computes drag force in the opposite direction of `v`.
    
    v: velocity
    system: System object with rho, C_d, area
    
    returns: Vector drag force
    """
    rho, C_d, area = system.rho, system.C_d, system.area
    
    C_d = drag_interp(V.mag)
    mag = -rho * V.mag**2 * C_d * area / 2
    direction = V.hat()
    f_drag = direction * mag
    return f_drag

In [64]:
C_d = drag_interp(43 * m / s)


Out[64]:
0.31695653551010883

In [65]:
# Solution

system = System(system, drag_interp=drag_interp)
V = Vector(30, 30) * m/s
f_drag = drag_force(V, system)


Out[65]:
\[\begin{pmatrix}-1.023081126128183 & -1.023081126128183\end{pmatrix} kilogram meter/second2\]

In [66]:
# Solution

slope_func(system.init, 0, system)


Out[66]:
(array([28.28427125, 28.28427125]) <Unit('meter / second')>,
 array([ -6.53489118, -16.33489118]) <Unit('meter / second ** 2')>)

In [67]:
# Solution

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


Out[67]:
values
success True
message A termination event occurred.

In [68]:
# Solution

results.tail()


Out[68]:
R V
4.600000 [86.37265447523836 meter, 6.352840864406531 me... [12.388380175637913 meter / second, -18.940140...
4.700000 [87.59981910596807 meter, 4.4276738010043495 m... [12.155784727055282 meter / second, -19.555228...
4.800000 [88.8038132288331 meter, 2.4417868717825417 me... [11.925028455712733 meter / second, -20.154598...
4.900000 [89.98482652030708 meter, 0.396745649345958 me... [11.696226535097686 meter / second, -20.738387...
4.918869 [90.20337227180956 meter, 0.0 meter] [11.653444424369289 meter / second, -20.845633...

In [69]:
# Solution

x = results.R.extract('x')
x_dist3 = get_last_value(x)


Out[69]:
90.20337227180956 meter

In [70]:
# Solution

x_dist - x_dist3


Out[70]:
9.101601791696496 meter

In [71]:
# Solution

# Here are the highest and lowest speeds

vs = results.V.extract('mag')
interval = min(vs), max(vs)


Out[71]:
(17.461046752327743 <Unit('meter / second')>, 40.0 <Unit('meter / second')>)

In [72]:
# Solution

# And here are the drag coefficients at the highest and lowest speed.
# They are substantially different.

drag_interp(interval)


Out[72]:
array([0.49707694, 0.33351435])

In [ ]: