This notebook is due on Friday, December 9th, 2016 at 11:59 p.m.. Please make sure to get started early, and come by the instructors' office hours if you have any questions. Office hours and locations can be found in the course syllabus. IMPORTANT: While it's fine if you talk to other people in class about this homework - and in fact we encourage it! - you are responsible for creating the solutions for this homework on your own, and each student must submit their own homework assignment.
FOR THIS HOMEWORK: In addition to the correctness of your answers, we will be grading you on:
To that end:
Put your name here!
Emboldened by his success in parachuting from a helicopter that's hovering in place, your professor has decided to start jumping from actual moving airplanes. In this section, we're going to model the behavior of a skydiver that is jumping from a moving plane and opening their parachute at some point above the ground. Your professor has requested your help in trying to figure out where he is going to land, knowing some basic information about how high above the ground he is and how fast he's going when he jumps out of the plane.
An object's position ($\vec{x}$) and velocity ($\vec{v}$) evolve according to a straightforward set of mathematical rules:
$\vec{x} = \vec{x}_0 + \int_0^t \vec{v} dt$
$\vec{v} = \vec{v}_0 + \int_0^t \vec{a} dt$
$\vec{a} = \frac{d\vec{v}}{dt} = \frac{1}{m}\sum \vec{F}$
Note that $\vec{a}$ is the acceleration, m is the mass of the skydiver, and $\sum \vec{F}$ are the sum of all of the forces acting on the sky diver. $\vec{x}_0$ and $\vec{v}_0$ are the position and velocity at $t=0$, respectively. The little arrows on top of letters indicate vector quantities - i.e., quantities that have direction as well as magnitude. With vector quantities, you can separate the components and solve them independently - for example, you only need to know the x-velocity to calculate the x-position (though ultimately you may need some quantities that need information from both - you need the speed, which includes x- and y-velocity, to calculate the acceleration due to air resistance). In the last term of the third equation, the sum is over all forces acting on the object (i.e., the sky diver).
The forces acting on the object are gravity, which points downward with a magnitude of $g=9.81$m s$^{-2}$, and air resistance, which can be written as:
$\vec{F}_{air} = -\frac{1}{2} C \rho A v^2 \hat{v}$
Where C is a constant corresponding to the shape of the object ($C=0.47$ for a sphere, $C=2.0$ for a parachute), $\rho$ is the density of the atmosphere, 1.2 kg/m$^3$, A is the cross-sectional area of the object, v is the magnitude of the velocity, and $\hat{v}$ is the unit vector of the velocity. Note the minus sign - this means that the force of air resistance always opposes the motion of the object, regardless of its direction.
Your task is to create a code that implements the equations presented above to model the position and velocity of the sky diver as a function of time, assuming the sky diver jumps out of the plane at some height above the ground Y and traveling at some horizontal speed V, and that the parachute opens at some time T. You have been provided a function that calculates the acceleration on the sky diver due to both air resistance and gravity; you need to write a function that takes in the height at which sky diver jumps out of the plane, his initial horizontal speed, and the time that the parachute opens, and returns the sky diver's position and velocity over time (so position, velocity, and time). Assume that the sky diver jumps out of the plane 2,000 meters above the ground while traveling at 50 m/s in the horizontal direction, and pulls the rip cord when he is 500 meters above the ground.
After you create the model, plot the sky diver's horizontal and vertical position and velocity as a function of time. Assuming that he jumps out of the plane at $t=0$ and $x=0$, how far does he land from the place he jumped out of the plane (in terms of horizontal distance), and at what time does he reach the ground?
PUT YOUR ANSWSERS TO THE QUESTIONS HERE! (And put your plots below, after your model)
ANSWER: He reaches the ground approximately 101 seconds after he jumps out of the plane, and lands approximatley 244 meters (horizontal distance) from where he jumps.
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
def calc_accel(velocity, parachute_open=False, use_resistance=True):
'''
This function calculates and returns the acceleration vector
for a sky diver, given an input velocity. You can indicate whether the
parachute is open or not (which affects the air resistance). Optionally,
you can turn on and off air resistance for comparison purposes with
resistance=True or False.
Inputs are velocity, a Boolean argument indicating if the parachute is open
(default False, so the parachute is assumed to be closed), and a Boolean argument
indicating if air resistance is to be used.
Output is a 2-element numpy array that indicates velocity.
Note 1: inputs and outputs are in SI units - meters per second for
velocity and meters per second^2 for acceleration.
Note 2: both input velocity and returned acceleration are 2-element
numpy arrays, with the horizontal (x) direction in the first element
and vertical (y) direction in the second element of the array.
'''
C_diver = 0.47 # skydiver, parachute closed
C_chute = 1.0 # skydiver, parachute open
rho_air = 1.225 # kg/m^3
mass_obj = 100.0 # mass of skydiver in kg, including all gear
# speed = magnitude of velocity vector
speed = (velocity[0]**2 + velocity[1]**2)**0.5
accel_vector = np.zeros(2)
# if we have air resistance turned on, compute the force due to air
# resistance and add it to the acceleration.
if use_resistance==True:
'''
Depending on whether the parachute is open or not,
that radically changes our air resistance! Check to
see if parachute's open and set parameters accordingly.
'''
if parachute_open==True:
rad_obj = 3.0
C_obj = C_chute
else:
rad_obj = 0.75
C_obj = C_diver
'''
Calculates acceleration due to air resistance
|F_air| = 0.5 * C * rho * A * v^2
C = constant based on shape
rho = density of air
A = cross-sectional area (depends on whether parachute is open or closed)
v = speed (magnitude of velocity)
Note that we calculate the magnitude in the first step and turn it into a vector in the second step.
'''
accel_air_resist = -0.5 * (C_obj/mass_obj) * rho_air * math.pi * rad_obj**2 * speed**2
# turn this into a vector quantity - velocity[0]/speed gives x-component of velocity unit vector
# and velocity[1]/speed gives y-component of velocity unit vector
accel_vector[0] += accel_air_resist*velocity[0]/speed
accel_vector[1] += accel_air_resist*velocity[1]/speed
# quick check to make sure that the acceleration magnitude isn't too ridiculous
# (this is for numerical stability, and is a slight hack. Forgive me!)
accel_mag = (accel_vector[0]**2 + accel_vector[1]**2)**0.5
max_accel = 40.0
if accel_mag > max_accel:
accel_vector[0] *= max_accel / accel_mag
accel_vector[1] *= max_accel / accel_mag
# we always have gravity, which always points downward.
accel_vector[1] -= 9.81
return accel_vector
In [ ]:
# put your model here
def calc_pos_vel(speed, height, chute_open_height, dt = 0.05, use_resistance=True):
'''
Calculates the flight path of the skydiver, given initial speed (assumed to
be horizontal), initial height, and the height at which the chute opens.
Inputs: initial speed (m/s), initial height (m), height at which the chute
opens (m), and optionally timestep (in seconds) and an option to turn off air
resistance (used for testing).
'''
# create arrays for initial position and velocity
position = np.array([0.0,height])
velocity = np.array([speed, 0.0])
# start at t=0
this_time = 0.0
# lists to keep track of x,y position and x,y velocity and time
pos = [[],[]]
vel = [[],[]]
time = []
# go until we hit the ground! (Er, reach the ground.)
while(position[1]>=0.0):
# need to make sure we keep track of when the parachute opens.
if position[1] > chute_open_height:
chute_open = False
else:
chute_open = True
# calculate acceleration given current velocity and knowledge about whether or not
# the parachute is open.
accel = calc_accel(velocity, parachute_open=chute_open)
# update velocity and position
velocity += accel*dt
position += velocity*dt
# update time
this_time += dt
# keep track of position, velocity, time
pos[0].append(position[0])
pos[1].append(position[1])
vel[0].append(velocity[0])
vel[1].append(velocity[1])
time.append(this_time)
# return position, velocity, time
return pos, vel, time
In [ ]:
# actually run the model!
pos, vel, time = calc_pos_vel(50,2000,500)
In [ ]:
# put the code that generates your plots here
# make some plots about velocity, print out last time.
print("last time is:", time[-1], "seconds")
vy, = plt.plot(time,vel[1],'r-')
vx, = plt.plot(time,vel[0],'b-')
plt.title("velocity components vs. time")
plt.xlabel("time [s]")
plt.ylabel("velocity [m/s]")
plt.legend((vx,vy),('x-velocity','y-velocity'))
In [ ]:
# make plot with position info
print("max x-value is at {:.2f} meters".format(pos[0][-1]))
y, = plt.plot(time,pos[1],'r-')
x, = plt.plot(time,pos[0],'b-')
plt.title("position components vs. time")
plt.xlabel("time [s]")
plt.ylabel("position [m]")
plt.legend((x,y),('x-position','y-position'))
plt.ylim(0,2000)
In this part of the homework, we're going to extend our fire model to make it more realistic. This model, like all models, makes a variety of simplifying assumptions. One that is particularly important is that there are many factors that affect the spread of fire in a forest, such as wind, the type of wood, the amount of dead underbrush, and even how close the branches of adjoining trees are to each other.
Given that there are so many possible complications, we're going to make a small modification of our model to take all of these uncertainties into account. We know that it's not certain that fire will spread, but we don't know exactly what that probability is. So, we're going to make this a new parameter for our system, called prob_spread. By experimenting with values of prob_spread (which is a fraction that should vary between 0 and 1), we can see what effect the probability of the fire spreading from cell to cell has on the overall spread of the fire.
The plan: Modify the fire model you developed in class as follows. In the original model, if a cell has trees and an adjacent cell is on fire, that cell automatically catches on fire in the next turn. Now, instead of it automatically catchin on fire, you generate a random number and compare it to prob_spread. If the number is smaller than prob_spread, the trees catch on fire.
Question: In the original model, you see a rapid transition in the amount of the forest that catches on fire when more than approximately 59% of cells are occupied by trees (this is referred to as the "tipping point" or "critical threshold"). As you try different values of prob_spread (say, values of 0.1 to 1.0 by steps of 0.1), how does the behavior differ from the standard model? At what value of prob_spread is it essentially indistinguishable from the standard model? Answer the question below, and support your answer with plots!
Put your answer here!
ANSWER: When the probability of fire spreading is prob_spread=1.0, the result is indistinguishable from the original model (as expected). As we decrease the probability of fire spreading a bit (to values of prob_spread=0.6-0.9) the critical point moves to the right (i.e., the forest density has to be higher to offset the lower probability the fire will spread from cell to cell). At prob_spread=0.5, the critical point is somewhere around a tree density of 90%, and at lower values of prob_spread (0.1-0.4) the fire basically immediately dies out. So, there is clearly some relationship between the fractional chance that the fire spreads from cell to cell and the forest density. One could imagine figuring this out by doing a much finer grid and figuring out where the critical point is!
In [ ]:
# Put your code here!
'''
NOTE: Virtually all of thise code is copy-pasted directly from the solutions
to the in-class assignment. We've only modified it to allow us to have a probability
of fire transmission in the advance_board() function. So, it looks like a lot of code but
there isn't really a ton of new code put in here. Also, there is now a double loop over both
the probability of fire transmission and the tree density.
'''
# standard includes
import numpy as np
import numpy.random as rand
%matplotlib inline
import matplotlib.pyplot as plt
# Next we are going to import some specific libraries we will
# use to get the animation to work cleanly
from IPython.display import display, clear_output
import time
# ADAPTED FROM THE 2D NUMPY TUTORIAL
# function plotgrid() takes in a 2D array and uses pyplot to make a plot.
# this function returns no values!
def plotgrid(myarray):
# first create two vectors based on the x and y sizes of the grid
x_range = np.linspace(0, myarray.shape[0], myarray.shape[0])
y_range = np.linspace(0, myarray.shape[0], myarray.shape[0])
# use the numpy meshgrid function to create two matrices
# of the same size as myarray with x and y indexes
x_indexes, y_indexes = np.meshgrid(x_range, y_range)
# make a list of all the x and y indexes that are either trees or fire.
tree_x = x_indexes[myarray == 1]; # trees!
tree_y = y_indexes[myarray == 1];
fire_x = x_indexes[myarray == 2]; # fire!
fire_y = y_indexes[myarray == 2];
# plot the trees and fire using appropriate colors.
# make the size of the polygons
# larger than the default so they're easy to see!
plt.plot(tree_x,tree_y, 'gs',markersize=10)
plt.plot(fire_x,fire_y, 'rs',markersize=10)
#Set the x and y limits to include a space of overlap so we don't cut off the shapes
plt.ylim([-1,myarray.shape[0]+1])
plt.xlim([-1,myarray.shape[0]+1])
def set_board(board_size=40,prob_tree=0.5):
'''
Creates the initial game board.
Inputs:
board_size: length of one edge of the board
prob_tree: probability that a given cell is a tree.
Outputs a 2D numpy array with values set to either 0, 1, or 2
(empty, tree, or fire)
'''
# all cells initialized to 'empty' (0) by default
game_board = np.zeros((board_size,board_size),dtype='int64')
# loop over board and roll the dice; if the random number is less
# than prob_tree, make it a tree.
for i in range(board_size):
for j in range(board_size):
if rand.random() <= prob_tree:
game_board[i,j] = 1
# set the whole i=0 edge of the board on fire. We're arsonists!
game_board[0,:] = 2
return game_board
def advance_board(game_board, prob_transmission=1.0):
'''
Advances the game board using the given rules.
Input: the game board, and a probability of fire getting transmitted
from cell to cell.
Output: a new game board
'''
# create a new array that's just like the original one, but initially
# set to all zeros (i.e., totally empty)
new_board = np.zeros_like(game_board)
# loop over each cell in the board and decide what to do.
for i in range(game_board.shape[0]):
for j in range(game_board.shape[1]):
# if the cell was empty last turn, it's still empty.
# if it was on fire last turn, it's now empty.
if game_board[i,j] == 0 or game_board[i,j] == 2:
new_board[i,j] = 0
# now, if it's a tree we have to decide what to do.
if game_board[i,j] == 1:
# initially make it a tree
new_board[i,j] = 1
# If one of the neighboring cells was on fire last turn,
# and a random number is less than the probability of fire transmission,
# this cell is now on fire!
if i > 0:
if game_board[i-1,j] == 2 and rand.random() <= prob_transmission:
new_board[i,j] = 2
if j > 0:
if game_board[i,j-1] == 2 and rand.random() <= prob_transmission:
new_board[i,j] = 2
if i < game_board.shape[0]-1:
if game_board[i+1,j] == 2 and rand.random() <= prob_transmission:
new_board[i,j] = 2
if j < game_board.shape[1]-1:
if game_board[i,j+1] == 2 and rand.random() <= prob_transmission:
new_board[i,j] = 2
# return the new board
return new_board
def calc_stats(game_board):
'''
Calculates the fraction of cells on the game board that are
a tree or are empty.
Input: the game board
Output: fraction that's empty, fraction that's covered in trees.
'''
# use numpy to count up the fractions that are empty
frac_empty = (game_board == 0).sum() / game_board.size
# do the same for trees
frac_tree = (game_board == 1).sum() / game_board.size
# return it!
return frac_empty, frac_tree
In [ ]:
'''
This cell does a quick test run with a user-specified fraction of trees,
and makes an animation. We don't necessarily expect the students to do
something this fancy, but it would certainly be nice. All of the plotting
stuff has been adapted from the Numpy 2D tutorial animation section.
'''
# some initial parameters
prob_tree=0.9
prob_trans=0.9
board_size = 40
# set up a figure
fig = plt.figure(figsize=(6,6))
# set up a game board
game_board = set_board(board_size=board_size, prob_tree=prob_tree)
# plot the initial board
plotgrid(game_board)
# need to decide when to stop - we're just going to ask if something
# is still on fire!
on_fire = True
# while there are fires, keep going
while on_fire == True:
# advance the game board
game_board = advance_board(game_board, prob_trans)
# do a bunch of plotting stuff, using code from the Numpy 2D array tutorial.
plotgrid(game_board)
time.sleep(0.05) # set this to be small because I don't want to wait!
clear_output(wait=True)
display(fig)
fig.clear()
# calculate the fractions that are empty and are trees
frac_empty, frac_trees = calc_stats(game_board)
# if 100% of cells are empty or are trees, there is no fire and we can stop.
if frac_empty + frac_trees == 1.0:
on_fire = False
# clean it up.
plt.close() # Close dynamic display
In [ ]:
'''
This cell is modified somewhat from the in-class solution. It now loops over the
probability of the fire being transmitted from cell to cell in addition to the tree
fraction (i.e., tree density), and stores all of the information in lists-of-lists for
easy plotting later.
'''
board_size = 40
# our lists-of-lists
f_trans_accum = []
f_tree_accum = []
f_burned_accum = []
# loop over probability of fire transmission, from 0.1 to 1.0
for prob_trans in np.arange(0.1,1.1,0.1):
# make empty lists to store various info
f_tree = []
f_burned = []
# print out so we know it's actually doing something.
print("prob_trans:",prob_trans)
# loop over tree fraction and keep track of the fraction of stuff
# burned as a function of the tree fraction
for tree_fraction in np.arange(0.01,1.01,0.01):
game_board = set_board(board_size=board_size, prob_tree=tree_fraction)
on_fire = True
while on_fire == True:
game_board = advance_board(game_board,prob_trans)
frac_empty, frac_trees = calc_stats(game_board)
if frac_empty + frac_trees == 1.0:
on_fire = False
f_tree.append(tree_fraction)
f_burned.append(frac_empty - (1.0-tree_fraction))
# accumulate everything in our lists-of-lists
f_trans_accum.append(prob_trans)
f_tree_accum.append(f_tree)
f_burned_accum.append(f_burned)
In [ ]:
# plot everything out!
# print the last list in our list-of-lists, since that's the baseline model (prob_trans = 1.0)
plt.plot(f_tree_accum[-1],f_burned_accum[-1],'r-',alpha=0.5,linewidth=5)
# then plot the rest of the models so we can see what's going on
for i in range(len(f_trans_accum)-1):
plt.plot(f_tree_accum[i],f_burned_accum[i])
plt.ylabel("burned fraction (normalized)")
plt.ylim(0,1)
In [ ]:
from IPython.display import HTML
HTML(
"""
<iframe
src="https://goo.gl/forms/OU4CVKhXdrdTGVVx1?embedded=true"
width="80%"
height="1200px"
frameborder="0"
marginheight="0"
marginwidth="0">
Loading...
</iframe>
"""
)