Unit 3: Simulations

Lesson 23: Free Fall

Notebook Authors

(fill in your two names here)

Team Roles

Facilitator: (fill in name)
Spokesperson: (fill in name)
Process Analyst: (fill in name)
Quality Control: (fill in name)

If there are only three people in your team, have one person serve as both spokesperson and process analyst for the rest of this activity.

At the end of this Lesson, you will be asked to record how long each Model required for your team. The Facilitator should keep track of time for your team.

Computational Focus: NumPy Arrays

Many science applications require the use of vectors (1D) and matrices (2D, and sometimes more). If you know any linear algebra or have taken physics, some amount of vectors and matrices will be familiar to you.

The NumPy library's main object (remember everything is an object in Python) is the "homogeneous multidimensional array". It is a table of elements (usually numbers), all of the same type, indexed by a tuple of positive integers. In NumPy the dimensions of an array are called axes.

Here are some resources on NumPy:

Model 1: Mathematical Operations (list vs. array)

An array is a variable that stores a set or a sequence of numbers. Numerical arrays are not yet defined in the standard Python language as a "core" data type like integers, floating points and strings. In order to have access to the array type from the NumPy library, we must first import the NumPy library. In this Lesson, we will do this by adding the following line to the start of every notebook, Class, function, or method (as appropriate) in which arrays are used:

from numpy import *

You might recall that in the past we imported NumPy differently (import numpy as np), the difference is that when done as recommended for this lesson (syntax above with *) you do not need to include np. in front of methods from the NumPy library when you call them.

Let's get started exploring the NumPy array and see how it is different from a list.

Type and run the lines of code below each in their own code cell:

list1 = [1,2,3,4,5]  
array1 = array([1,2,3,4,5])  
print(list1 + list1)
print(array1 + array1) 
print(list1 * 3)
print(array1 * 3)
print(list1 + 3)
print(array1 + 3)
print(array1[0])

In [ ]:

Critical Thinking Questions

1. How does defining an array differ from defining a list?

2. Describe how the following mathematical operations differ on an array versus on a list:
2a. +

2b. *

3. Consider adding 10 to all the numbers in a list/array:
3a. Define a new list, write and run the code necessary to add 10 to all numbers in a list.


In [ ]:

3b. Define a new list, write and run the code necessary to add 10 to all numbers in an array.


In [ ]:

3c. What is the advantage of performing this same operation to all the numbers in an array instead of a list?

4. If you only wanted to add 10 to the first number of a sequence, is there any advantage to using a list or an array? Explain.

Model 2: Initialization of Arrays

Another significant distinction between a list and an array is that an array's size is immutable. This means that unlike lists, there is no append method to grow the size of the sequence. The array size and content must be set upon declaration.
Using the zeros() method, we can easily initialize (create) an array and redefine the array elements as desired. You can also use the range method. One final difference between a list and an array is that each element in an array must have the same variable/data type (known as the base type of the array).

Type and run the pairs of lines of code below each pair in their own code cell:

time = zeros(10,int)
print(time)

time = array(range(10),float)
print(time) 

position = time + (9.8/2)*time
print(position) 

mix_list = [1,2.5]
print(mix_list) 

mix_array = array(mix_list)
print(mix_array)

Critical Thinking Questions

5. How does the syntax to create an array of floating point numbers differ from an array of integers?

6. Consider the Python code in the model and identify the line of code that was used to create an array: 6a. using a NumPy array method

6b. from another array

6c. from a list

7. What variable type(s) are the elements inside:
7a. mix_list?

7b. mix_array?

8. How many numbers would be in a sequence of floating-point numbers starting at 0 and ending at 4.0, incrementing in steps of 0.1?

9. Write and run the Python code necessary to create an array of floating-point numbers from 0 to 4.0, incrementing in steps of 0.1. You should not use a loop. (Hint: if you run into errors, call help on the methods that you are calling to better understand their functionalities.)


In [ ]:

Model 3: Physics of 1D free fall

I. Exact Analytical Solution (what you learned in physics)

Imagine an object, such as a ball, rising or falling vertically near the surface of the Earth, i.e., no horizontal motion. Consider the case where the only force acting on the object is the gravitational force by the Earth, i.e., air friction is small compared to the gravitational force and can be ignored. If we consider up to be the positive direction, then the vertical component of the gravitational force on the object by the Earth is given by

$$F^{grav}_y = - m g \hspace{45mm} (1) $$

where $m$ is the mass of the object and $g = -9.8$ N/kg is the magnitude of the local gravitational constant, and the minus sign is because the gravitational force points downward directly toward the Earth. Using Newton’s second law, we can determine the motion of the object,

$$ a_y = \frac{F^{net}_y}{m} \hspace{47mm} (2)$$

where $a_y$ is the vertical component of the acceleration of the object and $F^{net}_y$ is the sum of all the forces acting on the object in the vertical direction. If we set $F^{net}_y = F^{grav}_y$ as we have justified above, and use the definition of acceleration,

$$a_y \equiv \frac{d^2 y(t)}{dt^2} \hspace{45mm} (3)$$

the model for the motion of the object is in the form of the second-order differential equation,

$$ \frac{d^2 y(t)}{dt^2} = \frac{-mg}{m}\hspace{42mm}$$

or more simply (since the $masses$ cancel),

$$ \frac{d^2 y(t)}{dt^2} = -g \hspace{45mm} (4)$$

Here, $y(t)$ is a function representing the vertical position of the object at any particular instant of time, $t$.

Hopefully, you remember from physics the exact analytical solutions to equation (4) for position and velocity as a function of time:

$$y(t)=y_1 + v_1 (t-t_1) - \frac{1}{2}g(t-t_1)^2 \hspace{10mm} (5)$$


$$v(t)=v_1 - g(t-t_1) \hspace{33mm} (6)$$

where $y_1$ is the initial position of the object, $v_1$ is the initial velocity, and $t_1$ is the initial time. For convenience, we can set $t_1 = 0$ s, so the exact analytical solutions can be simplified to:

$$y(t)=y_1 + v_1 t - \frac{1}{2}gt^2 \hspace{26mm} (7)$$


$$v(t)=v_1 - gt \hspace{44mm} (8)$$

II. Approximate Numerical Solutions (useful when analytical solutions are not available)

An example of this idea was used to calculate the trajectory of John Glenn's space capsule in the recent film and book Hidden Figures, which was really good if you didn't see it.

We can also determine an approximate numerical solution for the problem posed above using the finite difference method.
This assumes that the value of the first derivative $f'(x)$ of the function $f(x)$ can be approximated by:

$$f'(x) \approx \frac{f(x+ \Delta x) - f(x)}{\Delta x} \hspace{16mm} (9)$$

Using this definition in conjunction with the appropriate first-order time derivatives of position and velocity:

$$y'(t) = v(t) \hspace{44mm} (10)$$$$v'(t) = a(t) = g \hspace{35mm} (11)$$

The relationships above, allow us to use the Euler method to determine the position and velocity of a free falling object as a function of time. We do this by using an the value for $y$ and $v$ at time $t$ to determine the next $y$ and $v$ at time $t + \Delta t$ using an iterative process. For example if our time step ($\Delta t$) is 1 second, we use our initial values of $y_0$ and $v_0$ at $t_0$ to calculate $y_1$ and $v_1$ at $t_0 + \Delta t$ or $t_1$.

If $\Delta t$ is sufficiently small, we will obtain a numerical answer that is close to the original analytical solution as shown in Equations 7 and 8.

III. Applying the Euler Method to 1D free fall

Now in order to use Euler's Method to approximate free fall, we need to rewrite Equation 9 using Equation 10 and the variables $y$ and $t$ appropriately to calculate an approximate solution for velocity.
since $v(t) = y'(t)$ and we can use the finite difference method (Eq 9):

$$v(t) = \frac{y(t + \Delta t) - y(t)}{\Delta t} \hspace{16mm} (12)$$

But what we really need to know is the position of the object ($y$) at $time = t + \Delta t$ so we need to rearrange the equation to show how to calculate a new position in terms of variables such as position ($y$) and velocity ($v$):

10. Rearrange Equation 12 above so that we can calculate the next position $y(t+ \Delta t)$.

Now we also need the rate of change of the velocity ($v'(t)$), which from Equation 11 is also $g$, so we need to rewrite Equation 9 using the variables $v$ and $t$ appropriately as an approximate solution for the gravity acceleration constant, $g$, since $v'(t)=g$:

$$g = \frac{v(t+ \Delta t)-v(t)}{\Delta t} \hspace{20mm} (13)$$

11. Rearrange Equation 13 above to show how to calculate the next velocity in terms of the variables position ($y$), velocity ($v$), and any constants (to make our calculation more convenient).

Model 3, continued: Iteration of Arrays

When a sequence of numbers is stored in an array, it is important to consider whether or not an operation can be performed on all elements at once, or if the operations should be performed on each element sequentially. For example, an iterative algorithm is a one such that each new solution depends on the previous answer so you would want to perform the calculations sequentially.

Critical Thinking Questions

12. Consider Equation 7 in Part I of Model 3 above. What type of programing construct (sequential, branching, looping) would be required to generate a list of y positions given a list of times?

13. Again considering Equation 7, what type of programing construct (sequential, branching, looping) would be required to generate an array of y positions given an array of times?

14. Considering your answers to Questions 10 and 11 above, what type of programing construct (sequential, branching, looping) would be required to generate an array of y positions using Euler’s method?

Team Programming

Make sure that as you work on these programming tasks, you are using your best pair programming technique (both driver and navigator) and switching drivers frequently.

Your code must have docstrings and comments for you to receive full credit.

Instead of using:

from numpy import *

your code should use:

import numpy as np
def FreeFall:

which does mean that you will need to use the dot notation to call numpy methods, as in:

np.zeros(5)

A. Define a FreeFall class with a constructor that takes four parameters besides self:

  1. the initial height/position of the object (in meters)
  2. the initial velocity (in meters per second)
  3. the total simulation time in seconds and
  4. the time interval/step ($\Delta t$) in seconds

Be careful about storing each value with the proper type.

Run your code several times and demonstrate that it works.


In [ ]:


In [ ]:
## testing

B. Define a free_fall_exact method in the FreeFall class that returns the position (based on Equation 7) at each time interval (starting at t = 0 s) up to the total simulation time. The method does not require any additional parameters besides self. If the object hits the ground before the total simulation time is reached, all the remaining elements in the array (after contact with the ground) should be 0.

Run your code several times and demonstrate that it works.


In [ ]:


In [ ]:
## testing

C. Define a new method, free_fall_approx, in the FreeFall class that uses Euler’s method to iteratively approximate the position of the object (use your answers to Questions 10 and 11). This function should take the same input and produce the same kind of output as free_fall_exact.

Run your code several times and demonstrate that it works.


In [ ]:


In [ ]:
## testing

D. Define a new method, plot_output, in the FreeFall class that plots the output of both the exact and approximate solutions (position with respect to time) on the same axes for comparison. Your plot should only plot the position of the object until it hits the ground, then stop plotting. Run your code several times and demonstrate that it works.


In [ ]:


In [ ]:
## testing

Temporal Analysis Report

How much time did it require for your team to complete each part?

Model 1:

Model 2:

Model 3:

Team Programming: