# Modeling and Simulation in Python

Chapter 22

Copyright 2017 Allen Downey



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