homework and exercises - Add air resistance to projectile motion

I am given an initial x and y position and initial velocity and I was asked to graph the trajectory in 1 second intervals.

This is what I have so far:

If $x_0 = 1, v_{0x} = 70, y_0 = 0, v_{0y} = 80, a_x = 0, a_y = -9.8$, and time will be 0,1,2,3... and so on.

Using these equations on every second you can find the plot will be a bell shaped with the highest point being ~ 325 m at about 600 seconds:

$$ x = x0 + (v_{0x})t + 1/2((a_x)t^2) $$ $$ y = y0 + (v_{0y})t + 1/2((a_y)t^2) $$

Usually in physics, we are taught in perfect condition with no air resistance. But what if there was air resistance? How would it affect this problem? How can i add it to my calculations and see what the difference is?

Answer

With drag force $- \alpha \left|{\dot{\bf r}}\right| {\dot{\bf r}}$ and gravitational force $-mg {\bf {\hat{y}}}$, the equations of motion are (see my answer to this question) $$ \begin{align} {\ddot{x}} &= - \beta {\dot{x}} \sqrt{{\dot x}^2+{\dot y}^2} \\ {\ddot{y}} &= - g - \beta {\dot{y}} \sqrt{{\dot x}^2+{\dot y}^2} \end{align} $$ where $\beta = \alpha / m$, $\alpha = \rho C_d A / 2$, $\rho$ is the air density, $C_d$ is the drag coefficient, and $A$ is the cross-sectional area of the projectile.

I've never seen the above equations solved analytically. Here's some Python code that plots the solution with and without air resistance. To run, copy into a text file myFile.py and run

python myFile.py

from the command line.

from pylab import *
import math

# Physical constants
g = 9.8
m = 1.0
rho = 1.0

Cd = 1.0
A = math.pi * pow(0.01,2.0)
alpha = rho * Cd * A / 2.0
beta = alpha / m

# Initial conditions
X0 = 1.0
Y0 = 0.0
Vx0 = 70.0
Vy0 = 80.0


# Time steps
steps = 1000
t_HIT = 2.0*Vy0/g
dt = t_HIT / steps

# No drag
X_ND = list()
Y_ND = list()


for i in xrange(steps+1):
  X_ND.append(X0 + Vx0 * dt * i)
  Y_ND.append(Y0 + Vy0 * dt * i - 0.5 * g * pow(dt * i,2.0))

# With drag
X_WD = list()
Y_WD = list()
Vx_WD = list()
Vy_WD = list()


for i in xrange(steps+1):
  X_ND.append(X0 + Vx0 * dt * i)
  Y_ND.append(Y0 + Vy0 * dt * i - 0.5 * g * pow(dt * i,2.0))

# With drag
X_WD = list()
Y_WD = list()
Vx_WD = list()
Vy_WD = list()


X_WD.append(X0)
Y_WD.append(Y0)
Vx_WD.append(Vx0)
Vy_WD.append(Vy0)

stop = 0
for i in range(1,steps+1):
  if stop != 1:
    speed = pow(pow(Vx_WD[i-1],2.0)+pow(Vy_WD[i-1],2.0),0.5)


    # First calculate velocity
    Vx_WD.append(Vx_WD[i-1] * (1.0 - beta * speed * dt))
    Vy_WD.append(Vy_WD[i-1] + ( - g - beta * Vy_WD[i-1] * speed) * dt)

    # Now calculate position
    X_WD.append(X_WD[i-1] + Vx_WD[i-1] * dt)
    Y_WD.append(Y_WD[i-1] + Vy_WD[i-1] * dt)

    # Stop if hits ground
    if Y_WD[i] <= 0.0:

      stop = 1

# Plot results
plot(X_ND, Y_ND)
plot(X_WD, Y_WD)
show()