Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

GitHub Repository: AllenDowney/ModSimPy
Path: blob/master/notebooks/chap22.ipynb
Views: ⁵³¹

Kernel: Python 3

Modeling and Simulation in Python

Chapter 22

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

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

Degrees and radians

Pint provides units to represent degree and radians.

In [18]:

degree = UNITS.degree
radian = UNITS.radian

If you have an angle in degrees,

In [19]:

angle = 45 * degree
angle

You can convert to radians.

In [20]:

angle_rad = angle.to(radian)

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

In [21]:

angle_rad.to(radian)

You can also convert from radians to degrees.

In [22]:

angle_rad.to(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)

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

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

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

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

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

Modeling and Simulation in Python

Vectors

Degrees and radians

Baseball

Visualizing the results

Animation

Under the hood

Exercises

Product

Resources

Company

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more, all in one place.

Modeling and Simulation in Python

Vectors

Degrees and radians

Baseball

Visualizing the results

Animation

Under the hood

Exercises

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.