# Day 9: Pre-class Assignment: Modeling with Ordinary Differential Equations#

✅ Put your name here

### Assignment instructions#

**This assignment is due by 7:59 p.m. the day before class,** and should be uploaded into the appropriate “Pre-class assignments” submission folder. Submission instructions can be found at the end of the notebook.

## Goals for today’s pre-class assignment#

In this assignment, you will:

Create a model for a skydiver falling using differential equations

Setup a

*system*of differential equations (i.e., more than 1).Solve these systems of differential equations using the update method from the previous assignment.

Solve these systems of differential equations using

`solve_ivp`

We will build off the previous assignment where we introduced using update equations to numerically integrate a (mathematical) function. Here we will expand our tool kit by numerically solving *multiple* differential equations simultaneously, and we’ll solve them using a tool built into SciPy called `solve_ivp`

.

But first, let’s talk about the role of differential equations in modeling.

## The Role of Differential Equations in Computational Modeling#

We often use equations in computational modeling, and many of the systems we wish to model contain things that are **changing**:

changing in time (E.g., Temperature or heart rate)

changing in space (E.g., Temperature or elevation)

changing strategies (E.g., deciding whether to invest)

changing numbers (of people, molecules, stars, etc.)

changing prices

and so on…

In mathematics, we describe change with calculus and derivatives in particular. When an equation has a derivative in it, it is called a “differential equation” (as opposed to, say, an algebraic equation). The types of equations can be used to model many different systems, and in this assignment, we’ll see how we can use them to model the free-fall of a skydiver.

✅ **Question**: What is a value or system that you deal with on a daily basis that is changing?

✎ *Put your answer here.*

## Modeling a falling skydiver without air resistance#

In this assignment, we will be modeling the motion of a skydiver. There are two variables that will be changing that we need to consider:

The

*height*of the skydiver (I.e., how far above the ground they are). This value will change with time as the skydiver gets closer and closer to the ground during their fall.The

*speed*(or, really, the*velocity*) of the skydiver * (I.e., how fast they are falling towards the ground). The longer the skydiver falls, the faster they will be falling.

Since we have two values that are changing, we will need two differential equations. We can express these two values using update equations:

The main variables to know are: \(h\) refers to height, and \(v\) refers to velocity.

**Note 1:** \(\frac{dh}{dt} = v\), meaning that velocity is defined as the derivative of height.

**Note 2:** For the first part of this problem, \(\frac{dv}{dt} = g\), where \(g\) is a constant equal to \(-9.81\). (Everything on Earth picks up speed at roughly the same rate, as long as we ignore air resistance.)

### The `derivs`

Function#

In order to solve these problems numerically, we need to specify what \(\frac{dh}{dt}\) and \(\frac{dv}{dt}\) are. We’re going to use a `derivs`

function, just like we did in the previous assignment, except this time it will return the derivatives for *two* equations. Specifically, we want it to return:

Putting it into code, we get the following:

```
# Create a function to compute derivatives of velocity and height
def derivs(v,g):
# derivative of height is velocity
dhdt = v
# derivative of velocity is gravity for this example
dvdt = g
return dhdt, dvdt
```

#### The `derivs`

function is the heart of how we will solve differential equations.#

(It may seem silly to define a function that takes in two values and returns the same two values, but we’ll be making it more complex in a second.)

### Initial Values#

We first will begin by assuming:

\(h_0 = 3000\) meters

\(v_0 = 0\) meters/second

\(g = -9.81\) gravity: meters/second^2

\(t_{max} = 30\) falling duration: seconds

\(dt = 1\) timestep

### Running the Model#

✅ **This code numerically solves for the motion of the skydiver and plots the result. It should be similar to the code you wrote in your previous assignment. Review each line of code an make sure you understand what it’s doing and then execute the code.**

```
# Import modules
import matplotlib.pyplot as plt
# Initialize variables
h = 3000 # initial height: m
v = 0 # initial velocoty: m/s
g = -9.81 # gravity: m/s^2
t = 0 # initial time
tmax = 30 # Falling time
dt = 1 # timestep
# Initialize lists for storing data
height = []
velocity = []
time = []
# Append initial values to lists
height.append(h)
velocity.append(v)
time.append(t)
# Create a time loop that will update the skydiver over time
# Use a while loop that will loop until t > tmax
while t <= tmax:
# Compute derivatives for use in update equations
dhdt, dvdt = derivs(v,g)
# Update Equations
h_new = h + dhdt*dt # new height
v_new = v + dvdt*dt # new velocity
# Append values to height and velocity lists
height.append(h_new)
velocity.append(v_new)
# Update old height/velocity with new height
h = h_new
v = v_new
# Increase time
t += dt # t = t + dt
# Update time list
time.append(t)
# Plotting height/velocity vs time
plt.figure(1)
plt.plot(time, height, color = 'green')
plt.xlabel('Time [seconds]')
plt.ylabel('Height [meters]')
plt.grid()
plt.figure(2)
plt.plot(time, velocity, color = 'purple')
plt.xlabel('Time [seconds]')
plt.ylabel('Veloctiy [meters/second]')
plt.grid()
```

✅ **Question:** If you want the model to stop when the skydiver hits the ground (going far too fast at the moment!), how might you build that into the model?

✎ *Put your answer here.*

## Adding air resistance to the model#

As we can see from our plots above, without air resistance the skydiver continously decreases velocity (or, “accelerates”) under the force of gravity. If they start high enough, they could reach *very big* velocities! In reality, this isn’t how skydiving works, eventually the skydiver reaches a constant **“terminal velocity”**. Check out this video to learn more about air resistance and terminal velocity if you’re interested.

In this case, the derivative of velocity (also known as acceleration) becomes dependent on the skydiver’s velocity. It is defined to be:

In this equation, \(A\) is the projected area of the skydiver, \(\left|v\right|\) is the absolute value of the velocity, and the \(m\) is the mass of skydiver. These make up the left term, or the acceleration due to air resistance. The right term, \(g\), is the acceleration from gravity, which was used in the code above.

As before, you are going to **model the free-fall of the sky-diver for a total of 30 seconds**, but now you need to take into account air resistance, using the formula above.

You can assume that:

The projected area (A) of the skydiver is \(0.4\, \mathrm{m}^2\)

Their mass (m) is 80 kg

Their initial height is 2000 m

Finally,

**Use a timestep (dt) of 0.1 seconds**.

✅ **Copy an paste the code from above and modify it to use the new timestep and initial height and dvdt in the derivs function to account for air resistance.** Think about how a skydiver might move in the air and see if your plots match those ideas. You should see that the velocity flattens out to a constant value, the terminal velocity.

```
#Put your code here
```

✅ **Question**: Based on your results, what does the terminal velocity of the skydiver appear to be?

✎ *Put your answer here.*

## Modeling the falling skydiver using `solve_ivp`

#

As with many computational models, there is often more than one way to find a solution. In fact, Python has a software library called SciPy that includes a tool for solving sets of differential questions called `solve_ivp`

.

Watch the following video to learn how we can use this tool to model the motion of the skydiver. The **SciPy documentation highlighted in the video can be found here.**

If the YouTube video doesn’t work, try this MediaSpace link.

```
from IPython.display import YouTubeVideo
# You may want to watch the video in full screen to make sure you can see everything!
YouTubeVideo("QGmJrEUNXA8",width=640,height=360)
```

## Reviewing the `solve_ivp`

solution for the skydiver problem.#

The code below is the same code used in the video to model the skydiver using the SciPy `solve_ivp`

tool.

✅ **Review the code and make sure you understand each line then execute the code.**

```
# Derivative function
def derivs(time,curr_vals):
# Declare parameters
g = -9.81 # m/s^2
# Unpack the current values of the variables we wish to "update" from the curr_vals list
h, v = curr_vals
# Right-hand side of odes, which are used to computer the derivative
dhdt = v
dvdt = g
return dhdt, dvdt
```

```
# Import commands
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp # This one is new to you!
# Declare Variables for initial conditions
h0 = 3000 # meters
v0 = 0 # m/s
tmax = 10 # seconds
dt = 1 # seconds
# Define the time array
time = np.arange(0, tmax + dt, dt)
# Store the initial values in a list
init = [h0, v0]
# Solve the odes with solve_ivp
sol = solve_ivp(derivs,(0,tmax),init,t_eval = time)
# Plot the results using the values stored in the solution variable, "sol"
# Plot the height using the "0" element from the solution
plt.figure(1)
plt.plot(time, sol.y[0,:],color="green")
plt.xlabel('Time [s]')
plt.ylabel('Height [m]')
plt.grid()
# Plot the velocity using the \"1\" element from the solution
plt.figure(2)
plt.plot(time, sol.y[1,:],color="purple")
plt.xlabel('Time [s]')
plt.ylabel('Velocity [m/s]')
plt.grid()
```

✅ **In your own words, describe what the variable curr_vals holds? You may need use the solve_ivp documentation mentioned above.**

✎ *Put your answer here.*

✅ **Try changing the order of the arguments in the function ( i.e., curr_vals and time). Does it still work?**

✎ *Put your answer here.*

## Adding air resistance to the `solve_ivp`

solution#

✅ Now you’re going to try and **modify the provided code that uses solve_ivp** to include air resistance

**using the same parameters you used when you added air resistance previously**. Make sure that the solution is computed for

**30 seconds**of skydiver freefall time and check if you get plots that look similar to your previous solution.

```
#Put your code here
# Import commands
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp # This one is new to you!
# Derivative function
def derivs(time,curr_vals):
# Declare parameters
g = -9.81 # m/s^2
A = ?? # m^2 ### MODIFY ME PLEASE!
m = ?? # kg ### MODIFY ME PLEASE!
# Unpack the current values of the variables we wish to "update" from the curr_vals list
h, v = curr_vals
# Right-hand side of odes, which are used to compute the derivative
dhdt = v
dvdt = g + ?? # MODIFY ME PLEASE! HINT: Look at your dvdt with air resistance from the update equations solution
return dhdt, dvdt
# Declare Variables for initial conditions
h0 = ?? # meters ### MODIFY ME PLEASE!
v0 = 0 # m/s
g = -9.81 # m/s^2
tmax = ?? # seconds ### MODIFY ME PLEASE!
dt = ?? # seconds ### MODIFY ME PLEASE!
# Define the time array
time = np.arange(0, tmax + dt, dt)
# Store the initial values in a list
init = [h0, v0]
# Solve the odes with solve_ivp
sol = solve_ivp(derivs, (0,tmax),init,t_eval = time)
# Unpack the results stored in the solution variable, "sol"
h = sol.y[0,:]
v = sol.y[1,:]
t = sol.t
# Plot the results we unpacked from "sol"
plt.figure(1)
plt.plot(t,h,color = 'green')
plt.xlabel('Time [s]')
plt.ylabel('Height [m]')
plt.grid()
plt.figure(2)
plt.plot(t,v, color = 'purple')
plt.xlabel('Time [s]')
plt.ylabel('Velocity [m/s]')
plt.grid()
```

✅ **Question**: Do your plots of the height and velocity of the skydiver from `solve_ivp`

agree with your solution when you used the update equations?

✎ *Put your answer here.*

### Assignment wrapup#

Please fill out the form that appears when you run the code below. **You must completely fill this out in order to receive credit for the assignment!**

```
from IPython.display import HTML
HTML(
"""
<iframe
src="https://cmse.msu.edu/cmse801-pc-survey"
width="800px"
height="600px"
frameborder="0"
marginheight="0"
marginwidth="0">
Loading...
</iframe>
"""
)
```

### Congratulations, you’re done!#

Submit this assignment by uploading it to the course D2L web page. Go to the “Pre-class assignments” folder, find the appropriate submission link, and upload it there.

See you in class!

© Copyright 2023, Michigan State University Board of Trustees