Simulating Projectile motion with Python — I

The most common physical phenomenon in our day to day life is the motion of an object in air. It may be observed in any game or similar other activities. We may observe an object moving up in air and coming down due to gravity or it may be thrown in air making an angle with horizontal in which case it travels in air for some time and finally, falls on the ground. We may simulate the motion in computer with the knowledge of basic numerical methods and coding. This helps in understanding and viewing the effect of various physical parameters like angle of projection or air resistance etc. on the prjectile motion.

In this section, I shall demonstrate the simulation of a simple vertical motion of a projectile without considering air resistance. Consider throwing a projectile vertically in the upward direction and expect it to come back and hit the ground after ten second. The motion of the projectile can be described by \begin{align} \frac{d² y}{dt²} = -g, \hspace{0.5cm}y(0)=0 \end{align} Here, $g$ is the acceleration due to gravity. This is a second order differential equation. To solve this type equation, we need two initial conditions — initial position and initial velocity. If we know the initial position and if we do not have prior knowledge of initial velocity which together progress the simulation, then we can apply a numerical technique, called “shooting method” to obtain correct value of velocity. So, for simulating this projectile motion, we have to do basically two operations.

1. Solving the differential equation.
2. Estimating correct initial velocity of projection by “shooting method”.

The differential equation can be solved by any standard method like Euler method or Runge Kutta method. But, we can make use of SciPy package to solve this initial value problem to save time and space.

First of all, we shall import necessary library packages like numpy for creating arrays, SciPy package for solving initial value problem and matplotlib for graphical demonstration of the simulation.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

To solve the second order differential equation numerically, we have to split up into two first order differential equations. \begin{align} &\frac{dy}{dt} = v \ & \frac{dv}{dt} = — g \end{align}

def model(t, z):
  y, v = z
  dy_dt = v
  dv_dt = - g
  dz_dt = [dy_dt, dv_dt]
  return dz_dt

def yFunc(v):
  v0 = v
  t = np.linspace(0, 10, 100)
  t_span = (t[0], t[-1])
  y0 = 0
  y0 = [y0, v0]
  sol = solve_ivp(model, t_span, y0, t_eval=t)
  return sol

g  = 9.8
v_guess = [35, 55, 75]
color = ['red', 'blue', 'green']
for i in range(len(v_guess)):
  sol = yFunc(v_guess[i])
  t = sol.t
  y = sol.y[0, :]
  plt.plot(t, y, ls=':', color=color[i], label='v = %0.2f'%v_guess[i])
plt.xlabel('$t$', fontsize=14)
plt.ylabel('$y(t)$', fontsize=14)
plt.axhline(ls='--', color='black')
plt.legend(loc='best')
plt.grid()
plt.show()

Since, we did not have have any prior knowledge of initial velocity, let us start with a guess value of it. It is seen that the solver function has simulated the projectile motion for three different guess values of initial velocities. But, it failed to arrive at the initial point after ten second. The value $v=55 m/s$ is little close to the level of expectation. Other two values are far away from the actual one. This means our guess on initial velocity is not correct. So, we try the simulation with a range of initial velocities.

@np.vectorize
def y(v0):
  t = np.linspace(0, 10, 100)
  t_span = (t[0], t[-1])
  y0 = 0
  y0 = [y0, v0]
  sol = solve_ivp(model, t_span, y0, t_eval=t)
  t = sol.t
  y, v = sol.y
  return y[-1]

v0 = np.linspace(0, 100, 100)

plt.plot(v0, y(v0), 'b--')
plt.axhline(c='r')
plt.xlabel('$v_{guess}$', fontsize=14)
plt.ylabel('$y$', fontsize=14)
plt.grid()
plt.show()

So, we see that for the velocity around 50, the projectile exactly returns to the initial point. For finding the exact value, we can apply bisection method to loacte the root. We can use SciPy package to apply bisection method.

from scipy.optimize import bisect

exact_velocity = bisect(y, 40, 60)
print("Exact velocity = %0.2f"%exact_velocity)
Exact velocity = 49.00

So, if we throw the projectile with an initial velocity 49 m/s, it will come back exactly after 10 sec. to the initial point. The way to arrive at correct value of initial velocty is called “shooting method”.

In the next post, I shall come up with the case of inclined projectile, i.e. when the projectile is fired making some angle with the horizontal and also, consider the effect of air resistance on the projectile motion.