CoCalc -- Week-06-projectile.ipynb

Jupyter notebook Week-06-Projectile/Week-06-projectile.ipynb

Path: Week-06-Projectile/Week-06-projectile.ipynb

Views: ¹⁶⁹⁷⁶

Kernel: Python 2 (SageMath)

Python exercise for projectile motion

In this exercise, we will use a simple script to plot the trajectory of a projectile with inclusion of air resistance. We will also calculate the range of the projectile.

The air resistance drag force $F_d$ in $N$ is approximated by the following relation.

F_d = C_d v^2

Here, $v$ is the magnitude of velocity in $m/s$ and $C_d$ is the coefficient of drag.

The following is a step-by-step guide to writing up the script. Some parts of the script have been left out for you to fill in.

Define the parameters in the problem

First we import the modules we need: numpy to provide some mathematics functions, and matplotlib for plotting.

To start describing the problem, we need to define the parameters in the problem.

the mass of the projectile
acceleration due to gravity
initial velocity (magnitude and direction)
drag coefficient

We also set the time step to track the trajectory.

In [2]:

import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline

# Model parameters
M = 1.0          # Mass of projectile in kg
g = 9.8          # Acceleration due to gravity (m/s^2)
V = 80           # Initial velocity in m/s
ang = 60.0       # Angle of initial velocity in degrees
Cd = 0.005       # Drag coefficient
dt = 0.5         # time step in s

# You can check the variables by printing them out
print V, ang

80 60.0

Set up the variables at time zero

Next, we generate velocity of the projectile as a function of time. To do that, we create the following lists.

list to store the values of time $t$
list to store x-component of velocity $v_x$
list to store y-component of velocity $v_y$
list to store x-component of acceleration $a_x$
list to store y-component of acceleration $a_y$

We start by putting in the initial velocity components $v_x = V \cos\theta$ and $v_y = V \sin\theta$ at $t=0$ .

To get the acceleration, we need to find the resultant force, which consists of the air drag $F_d$ and weight $W$ .

W = M g

and

F_d = C_d v^2

The x and y components of the resultant force are given by

F_x = - F_d \cos\theta \qquad \text{and} \qquad F_y = -Mg - F_d \sin\theta

where $\theta$ is the angle that the velocity forms with the positive x-axis. The force components follow the convention that they are positive when pointing to the right and upwards.

By Newton's second law, the acceleration of the projectile is thus given by

a_x = - (F_d \cos\theta)/M \qquad \text{and} \qquad a_y = -g - (F_d \sin\theta)/M

Let's implement these steps to get the velocity, drag force, and acceleration at $t=0$ .

In [3]:

# Set up the lists to store variables
# Start by putting the initial velocities at t=0
t = [0]                         # list to keep track of time
vx = [V*np.cos(ang/180*np.pi)]  # list for velocity x and y components
vy = [V*np.sin(ang/180*np.pi)]

# Drag force
drag = Cd*V**2                      # drag force 

# Create the lists for acceleration components
ax = [-(drag*np.cos(ang/180*np.pi))/M]        
ay = [-g-(drag*np.sin(ang/180*np.pi)/M)]

# Print out some values to check
print ax[0]
print ay[0]
print vx[0]
print vy[0]

-16.0
-37.5128129211
40.0
69.2820323028

Update the velocity for every time-step

To get the velocity at the next time step, we make use of the following approximation.

a(t_{n}) = \frac{dv(t_{n})}{dt} \approx \frac{v(t_{n+1}) - v(t_n)}{\Delta t}

v(t_{n+1}) \approx v(t_n) + a(t_n) \Delta t

This is just a very primitive way of doing integration of acceleration to get velocity. In the literature, this is known as Euler's method.

Let us calculate 10 sets of velocities for 10 time-steps. We will use a while-loop and keep a counter to count to 10.

In [4]:

# Use Euler method to update variables
counter = 0
while (counter < 10):
    t.append(t[counter]+dt)                # increment by dt and add to the list of time 
    vx.append(vx[counter]+dt*ax[counter])  # Update the velocity
    vy.append(vy[counter]+dt*ay[counter])  

    # With the new velocity calculate the drag force
    vel = np.sqrt(vx[counter+1]**2 + vy[counter+1]**2)   # magnitude of velocity
    drag = Cd*vel**2                                   # drag force 
    ax.append(-(drag*np.cos(ang/180*np.pi))/M)     
    ay.append(-g-(drag*np.sin(ang/180*np.pi)/M))
    
    # Increment the counter by 1
    counter = counter +1

# Print the values to check
print "t=", t
print "vx=", vx

# Let's plot the velocity against time
plt.plot(t,vx)
plt.ylabel("x-velocity (m/s)")
plt.xlabel("time (s)")

t= [0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]
vx= [40.000000000000007, 32.000000000000007, 27.528951416567008, 24.787883378136915, 23.023394914592512, 21.846818145611152, 21.027776019653995, 20.413445666052425, 19.891162819384121, 19.368416481078107, 18.760291921319446]

<matplotlib.text.Text at 0x7f351488ab50>

Update the positions

After we have obtained the velocities, we can get the $x$ and $y$ positions of the projectile using a similar concept.

v_x(t_{n}) = \frac{dx(t_{n})}{dt} \approx \frac{x(t_{n+1}) - x(t_n)}{\Delta t} \qquad \text{and} \qquad v_y(t_{n}) = \frac{dy(t_{n})}{dt} \approx \frac{y(t_{n+1}) - y(t_n)}{\Delta t}

x(t_{n+1}) \approx x(t_n) + v_x(t_n) \Delta t \qquad \text{and} \qquad y(t_{n+1}) \approx y(t_n) + v_y(t_n) \Delta t

Starting with $x=0$ and $y=0$ at $t=0$ , we can use Euler's method to get the positions for the 10 time-steps using the velocities found previously.

In [6]:

# Initialise the lists for x and y
x = [0]
y = [0]

# Use Euler method to update variables
counter = 0
while (counter < 10):
    # Update the positions x and y
    x.append(x[counter]+dt*vx[counter])    
    y.append(y[counter]+dt*vy[counter])    
    # Increment the counter by 1
    counter = counter +1


# Let's plot the trajectory
plt.plot(x,y,'ro')
plt.ylabel("y (m)")
plt.xlabel("x (m)")
print "Range of projectile is {:3.1f} m".format(x[counter])

Range of projectile is 124.9 m

Let's put everything together

How about calculating the velocity and position all together in the same while-loop?

Also, let's stop the while-loop after the projectile drops back to the ground again. One simple way to detect this is to check whether the value of $y$ position falls below zero. So, we keep running the while-loop while $y \ge 0$ .

In [2]:

import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline

# Model parameters
M = 1.0          # Mass of projectile in kg
g = 9.8          # Acceleration due to gravity (m/s^2)
V = 80           # Initial velocity in m/s
ang = 60.0       # Angle of initial velocity in degrees
Cd = 0.005       # Drag coefficient
dt = 0.5         # time step in s

# You can check the variables by printing them out
print V, ang

# Set up the lists to store variables
# Initialize the velocity and position at t=0
t = [0]                         # list to keep track of time
vx = [V*np.cos(ang/180*np.pi)]  # list for velocity x and y components
vy = [V*np.sin(ang/180*np.pi)]
x = [0]                         # list for x and y position
y = [0]

# Drag force
drag=Cd*V**2                      # drag force 

# Acceleration components
ax = [-(drag*np.cos(ang/180*np.pi))/M ]          
ay = [-g-(drag*np.sin(ang/180*np.pi)/M) ]

## Leave this out for students to try
# We can choose to have better control of the time-step here
dt = 0.2

# Use Euler method to update variables
counter = 0
while (y[counter] >= 0):                   # Check that the last value of y is >= 0
    t.append(t[counter]+dt)                # increment by dt and add to the list of time 
    
    
    
    # Update velocity
    vx.append(vx[counter]+dt*ax[counter])  # Update the velocity
    vy.append(vy[counter]+dt*ay[counter])

    # Update position
    x.append(x[counter]+dt*vx[counter])    
    y.append(y[counter]+dt*vy[counter])    

    # With the new velocity calculate the drag force and update acceleration
    vel = np.sqrt(vx[counter+1]**2 + vy[counter+1]**2)   # magnitude of velocity
    drag = Cd*vel**2                                   # drag force 
    ax.append(-(drag*np.cos(ang/180*np.pi))/M)     
    ay.append(-g-(drag*np.sin(ang/180*np.pi)/M))
    
    # Increment the counter by 1
    counter = counter +1

# Let's plot the trajectory
plt.plot(x,y,'ro')
plt.ylabel("y (m)")
plt.xlabel("x (m)")
   
# The last value of x should give the range of the projectile approximately.

print "Range of projectile is {:3.1f} m".format(x[counter])

80 60.0
Range of projectile is 174.4 m

In [ ]:

In [ ]:

In [ ]:

How far can the projectile go?

Try changing the initial velocity of the projectile and the drag coefficient to see how the trajectory changes.

Keeping the initial velocity magnitude fixed, find the initial velocity angle that will give the maximum range for a specific drag coefficient!

In [20]:

Trivial

If you set the coefficient of drag to zero, can you get back the same trajectory you learned in Physics? Try checking the range of the projectile against the closed-form solution.