Workshop: Projectile Motion

This tutorial is adapted from a presentation held by Jonas van den Brink at iCSE4schools’ workshop in Katowice January 29, 2015.

This workshop focuses on introducing computations to introductory physics. Introducing computations should lead to a sense of empowerment For this to be possible, the computations must

  1. relate to well-known problems
  2. be shown to be a powerful tool
  3. be understable. Students should write their own code

An introduction to computing trajectories

The goal is to find the velocity and position of an object as functions of time: \(\vec{v}(t)\), \(\vec{r}(t)\)

Equations of motion:

\[\frac{d \vec{r}}{d t} = \vec{v}(t), \qquad \frac{d \vec{v}}{d t} = \vec{a}(t)\]

where \(r\) is position, \(v\) is velocity, \(a\) is acceleration and \(t\) is time.


Newton’s second law of motion:

\[\vec{F} = m\vec{a}\]


\[\vec{F}(r,v,t) = m\vec{a}(r,v,t).\]
Our algorithm is now as follows:
  1. Find the physical forces of the system
  2. Use Newton’s second law to find the acceleration
  3. Calculate ~v(t) and ~r(t) by solving the equations of motion

In this workshop, we will solve step number 3 numerically, using the Euler method.

The Euler Method

A method for solving ordinary differential equations (ODEs)

We can solve the equations of motion numerically using the Euler method.

From the definition of the derivative

\[\frac{d v}{d t} = \lim_{\Delta t \to 0} \frac{v(t+\Delta t) - v(t)}{\Delta t} = a(t)\]

we now remove the limit, making \(\Delta t\) a very small constant.

\[\frac{v(t+\Delta t) - v(t)}{\Delta t} \approx a(t)\]

Solving for \(v(t+\Delta t)\) gives:

\[v(t+\Delta t) \approx v(t) + a(t)\cdot \Delta t\]

We can solve the equations of motion by stepping forward in time

\[v(t+\Delta t) = v(t) + a(t)\cdot \Delta t\]

If \(a(t)\) and \(v(t)\) are known, we can calculate \(v(t+\Delta t)\) plot of eulers method

Our functions are no longer continuous, they have become discretized

We only focus on multiples of our time-step:

\[\begin{split}t \in \{ 0,\ \Delta t,\ 2\Delta t, \ 3\Delta t, \ldots \} \\ t_i \equiv i\cdot\Delta t\end{split}\]

Introduce the shorthand:

\[\begin{split}v(t_i) \equiv v_i \\ r(t_i) \equiv r_i \\\end{split}\]
time discretization

We solve the equations of motion iteratively

\[v_{i+1} = v_i + a_i\cdot\Delta t\]
\[r_{i+1} = r_i + v_i\cdot \Delta t\]

For each time step, we must calculate the acceleration

\[a_i = a(r_i, v_i, t_i).\]

We repeat these steps, starting at our initial conditions \(v_0\) and \(r_0\), until we have reached our end-time \(t_N\)

\[i = 0,1,2,3,\ldots, N.\]

Algorithm for the Euler method

for \(i=0,1,2,3,\ldots, N-1\):
  1. Use the previous results \(x_i\) and \(v_i\) to compute the acceleration: \(a_i =F(x_i, v_i, t_i)/m\)
  2. Compute the new velocity: \(v_{i+1} = v_i +a_i\Delta t\)
  3. Compute the new position: \(r_{i+1} = r_i + v_i\Delta t\)


Moving from physics and math to actual computer code

for \(i=0,1,2,3,\ldots, N-1\):
  1. Use the previous results \(x_i\) and \(v_i\) to compute the acceleration: \(a_i = F(x_i, v_i, t_i)/m\)
  2. Compute the new velocity: \(v_{i+1} = v_i + a_i\Delta t\)
  3. Compute the new position: \(r_{i+1} = r_i + v_i\Delta t\)
for i in range(N):
    a[i] = F(r[i], v[i], t[i])/m
    v[i+1] = v[i] + a[i]*dt
    r[i+1] = r[i] + v[i]*dt

We want the code to look as much as possible like the physics and math we write on paper

\[t_i \Rightarrow \texttt{t[i]} \qquad v_i \Rightarrow \texttt{v[i]} \qquad r_i \Rightarrow \texttt{r[i]}\]

We also need various bookeeping code

Here we define the arrays we will be using

# Import various functions meant for numerical science
import numpy as np
from math import cos,sin,pi

t_0 = 0 # Start time, s
t_end = 10 # End time, s
N = 1000 # Number of time steps

# Create a uniformly spaced time-array
t = np.linspace(t_0, t_end, N+1)

# Calculate the size of a time step
dt = t[1] - t[0]

# Create empty acceleration, velocity and position arrays
a = np.zeros((N+1,2))
v = np.zeros((N+1,2))
r = np.zeros((N+1,2))

# Set initial conditions
v[0] = (100*cos(pi/6), 100*sin(pi/6)) # inital velocity, m/s
r[0] = (0,1)  # initial position, m

Here we define physical constants for our system and define the function that describes the forces

m = 5.5 # mass, kg
g = 9.81 # acceleration of gravity, m/s^2
rho = 1.3 # air density, kg/m^3
C_D = 0.45 # drag coefficient
d = 0.11 # diameter of cannonball, m
A = pi*d**2 # cross-sectional area, m^2

def F(r, v, t):
    return (0, -m*g) - 0.5*rho*C_D*A*abs(v)*v

This example show the forces acting on the cannonball as it sails through the air

\[F(x,v,t) = F_g + F_d(\vec{v}) = -mg\vec{k} - \frac{1}{2}\rho C_D A |\vec{v}|\vec{v}\]

As soon as we have solved the equations of motion, we can plot the result

# Solving equations of motion iteratively
for i in range(N):
    a[i] = F(r[i], v[i], t[i])/m
    v[i+1] = v[i] + a[i]*dt
    r[i+1] = r[i] + v[i]*dt

# Extract x and y coordinates
x = r[:,0]
y = r[:,1]

# Import functionality for plotting
import matplotlib.pyplot as plt

# Plot figure

# Prettify the plot
plt.xlabel('Horizontal distance, [m]')
plt.ylabel('Vertical distance, [m]')
plt.title('Trajectory of a fired cannonball')
plt.axis([0, 900, 0, 250])

# Makes the plot appear on the screen

Numerical Experimentation

Altering parameters let’s us immediately see the consequences

plot_cannonball2 plot_cannonball3

Students can use numerical experimentation to build intuition and knowledge

  • Numerical results can be compared to known analytical solutions. Are numerical results trustworthy?
  • Can study how results are directly changed by parameter choice. Are the parameters chosen reasonable?
  • Can look at systems with and without certain contributions, such as air drag. What is important, and what can be ignored?

Examples of possible projects

Catapults and cannons and sports such as baseball

  • Easy to compare with experimental data, either before or after simulation.
  • Can look into studies of air drag, Reynolds number etc. cannonball

Skydiving and bungeejumping

  • Great study on free fall and terminal velocity
  • Can study how parameters such as cross-sectional area and drag coefficient change as the parachute is opened
  • Can plot the g-forces affecting the jumper. Which sport is more “extreme”?

Pendulum and angular motion


  • Can solve pendulum also for large angles!
  • Energy can beplotted as functions of time
  • Can also simulate double pendulum and chaotic systems

Modelling the solar system

  • Students can gather real data of planetary orbits from NASA webpages
  • Can combine numerical simulation with better graphics