# Modeling and Simulation in Python

Chapter 22



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



Here's a two dimensional Vector in meters.



In [3]:

A = Vector(3, 4) * m



We can access the elements by name.



In [4]:

A.x




In [5]:

A.y



The magnitude is the length of the vector.



In [6]:

A.mag



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



In [7]:

A.angle



If we make another Vector with the same units,



In [8]:

B = Vector(1, 2) * m



We can add Vector objects like this



In [9]:

A + B



And subtract like this:



In [10]:

A - B



We can compute the Euclidean distance between two Vectors.



In [11]:

A.dist(B)



And the difference in angle



In [12]:

A.diff_angle(B)



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



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)



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



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()



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 goes here



Pint provides units to represent degree and radians.



In [18]:

degree = UNITS.degree



If you have an angle in degrees,



In [19]:

angle = 45 * degree
angle





In [20]:



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



In [21]:



You can also convert from radians to degrees.



In [22]:



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



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 goes here




In [25]:

# Solution goes here



### 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)



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

# 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)



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)



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)



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)



Now we can call run_ode_solver



In [35]:

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



The final label tells us the flight time.



In [36]:

flight_time = get_last_label(results) * s



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



### 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')



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')



### 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)



A ModSimVector is a specialized kind of Pint Quantity.



In [44]:

isinstance(V, Quantity)



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




In [46]:

type(V1)



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



In [47]:

V2 = m * V




In [48]:

type(V2)



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



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



In [50]:

V2[0], V2[1]



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



In [51]:

vector_mag(V2), vector_angle(V2)



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)




In [53]:

# Solution goes here




In [54]:

# Solution goes here




In [55]:

# Solution goes here




In [56]:

# Solution goes here




In [57]:

# Solution goes here



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)




In [59]:

# Solution goes here




In [60]:

# Solution goes here



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

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



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 goes here




In [63]:

# Solution goes here




In [64]:

C_d = drag_interp(43 * m / s)




In [65]:

# Solution goes here




In [66]:

# Solution goes here




In [67]:

# Solution goes here




In [68]:

# Solution goes here




In [69]:

# Solution goes here




In [70]:

# Solution goes here




In [71]:

# Solution goes here




In [72]:

# Solution goes here




In [ ]: