Example Task¶
Projectile Motion with Air Resistance¶
Goal: investigate the motion of real projectiles (e.g. batted baseballs) compared to the ideal case
This notebook is an example of the format that will be used in the real task sessions. It starts with a short background, then come a handful of numbered tasks and a couple of guiding questions for the 2-minute presentation. All notebooks contain some starting code that will help you begin working on the tasks. Following the "Tasks" section below is one possible solution, included here so you can see the level of depth and the kind of write-up expected.
Basic background¶
Frictionless projectile motion has an exact solution. Launched from the origin with speed $V_0$ at angle $\theta$ above the horizontal, ignoring air resistance:
$$x(t) = V_0\cos\theta\, t, \qquad y(t) = V_0\sin\theta\, t - \tfrac12 g t^2,$$
which gives a parabola, a range $R = V_0^2\sin(2\theta)/g$, and a maximum height $H = V_0^2\sin^2\theta/(2g)$. This is the trajectory shown as a dashed curve in most textbook figures. It's perfectly symmetric, as the ball takes as long to climb as it does to fall, and it lands at the same angle $\theta$ it was launched at, just flipped below the horizontal.
However, in many real-life scenarios, air resistance plays a crucial part and cannot be ignored. For example, baseball and golf players sometimes claim that hit balls appear to fall almost straight down at the end of their flight, at least noticeably steeper than the angle they went up at. Is this a real effect of air resistance, or an optical illusion? A common model for the drag force is
$$\vec F^{(f)} = -k\,m\,|\vec v|^{\,n}\,\frac{\vec v}{|\vec v|},$$
i.e. a force opposing the velocity, with magnitude growing as some power $n$ of the speed. Three values of $n$ are commonly used depending on the speed regime: $n=1$ (low speed, viscous drag), $n=3/2$ (an intermediate regime), and $n=2$ (high speed, the most common choice for a thrown or batted ball, and the regime we'll focus on). Newton's second law in components gives
$$\ddot x = -k\,|\vec v|^{\,n-1}\dot x, \qquad \ddot y = -g - k\,|\vec v|^{\,n-1}\dot y, \qquad |\vec v| = \sqrt{\dot x^2+\dot y^2}.$$
Unlike the frictionless case, there is no closed-form solution once drag is included, the trajectory has to be found numerically, by turning this into four first-order ODEs ($x,\dot x,y,\dot y$) and integrating with (for example) RK4.
Tasks¶
Task #1: The code below is an RK4 integrator for frictionless projectile motion, launched at $V_0=35\,\mathrm{m/s}$, $\theta=45°$. Check that the numerically integrated range and maximum height agree with the analytic formulas above.
Task #2: Add air resistance to the equations of motion, using the quadratic drag law ($n=2$).
Task #3: For the same launch speed (35 m/s) and angle (45°), sweep the drag strength $k$ and determine the impact angle (the angle below the horizontal at which the ball lands). A realistic value for a baseball is around $k\approx0.005\,\mathrm{m^{-1}}$ (from $k=\tfrac12 C_d\rho_{\rm air}A/m$ with $C_d\approx0.3$). Include that value in your sweep, but also look at a wider range so you can see the trend. At a realistic drag strength, how much steeper is the impact angle than the launch angle? Is there a drag strength strong enough to make the ball land "almost straight down" (impact angle $\gtrsim80°$)?
Task #4: Repeat Task #3 for $n=1$ and $n=3/2$ as well, and see whether your conclusion about "falling almost straight down" depends on which drag law you use.
Guiding questions for the presentation (2 minutes only!)¶
- What impact angle did you find at a realistic drag strength, and how does it compare to the $45°$ launch angle?
- Does the time spent ascending vs. descending tell the same story as the impact angle?
- What went wrong or surprised you?
- If you had another hour, what would you check next?
import numpy as np
import matplotlib.pyplot as plt
g = 9.8 # m/s^2
def projectile_deriv(state):
"""
state = [x, vx, y, vy]
Frictionless projectile motion: no air resistance.
"""
x, vx, y, vy = state
return np.array([vx, 0.0, vy, -g])
def rk4_step(f, state, dt, *args):
k1 = f(state, *args)
k2 = f(state + 0.5*dt*k1, *args)
k3 = f(state + 0.5*dt*k2, *args)
k4 = f(state + dt*k3, *args)
return state + (dt/6.0)*(k1 + 2*k2 + 2*k3 + k4)
def simulate_until_ground(f, state0, dt, *args, t_max=30.0):
"""Integrate until y drops back below 0, then linearly interpolate the landing point."""
state = state0.copy()
traj = [state.copy()]
t = 0.0
while t < t_max:
state_new = rk4_step(f, state, dt, *args)
t += dt
if state_new[2] < 0:
frac = state[2] / (state[2] - state_new[2])
traj.append(state + frac*(state_new - state))
return np.array(traj)
state = state_new
traj.append(state.copy())
return np.array(traj)
# Example run: V0=35 m/s, theta=45 degrees
V0 = 35.0
theta = np.deg2rad(45)
state0 = np.array([0.0, V0*np.cos(theta), 0.0, V0*np.sin(theta)])
traj = simulate_until_ground(projectile_deriv, state0, 0.01)
# plot the example trajectory
plt.plot(traj[:, 0], traj[:, 2])
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.title("Frictionless trajectory")
plt.hlines(0, 0, traj[-1, 0], color='k', ls=':')
plt.axis('equal')
plt.show()
R_analytic = V0**2 * np.sin(2*theta) / g
H_analytic = V0**2 * np.sin(theta)**2 / (2*g)
R_numeric = traj[-1, 0]
H_numeric = traj[:, 2].max()
print(f"Range: numeric={R_numeric:.3f} m, analytic={R_analytic:.3f} m, "
f"relative difference={abs(R_numeric-R_analytic)/R_analytic:.2%}")
print(f"Height: numeric={H_numeric:.3f} m, analytic={H_analytic:.3f} m, "
f"relative difference={abs(H_numeric-H_analytic)/H_analytic:.2%}")
Range: numeric=125.000 m, analytic=125.000 m, relative difference=0.00% Height: numeric=31.250 m, analytic=31.250 m, relative difference=0.00%
Task #2¶
adding air resistance
We redefine projectile_deriv here to include drag. n=2 by default gives the quadratic
drag law we'll focus on first, but setting it to other values will come in handy for task #4.
def projectile_deriv(state, k=0.0, n=2):
"""
state = [x, vx, y, vy]
Projectile motion with a drag force -k*m*|v|^n * v/|v| (per unit mass, so k already
has the mass divided out). k=0 recovers the frictionless case from Task #1.
"""
x, vx, y, vy = state
speed = np.hypot(vx, vy)
if speed > 0 and k > 0:
drag = k * speed**(n - 1)
ax = -drag * vx
ay = -g - drag * vy
else:
ax, ay = 0.0, -g
return np.array([vx, ax, vy, ay])
# Quick check: k=0 should reproduce the Task #1 trajectory exactly
traj_check = simulate_until_ground(projectile_deriv, state0, 0.01, 0.0, 2)
print("Matches Task #1 (frictionless) run:", np.allclose(traj_check, traj))
# Comparing frictionless vs. a realistic(ish) k
k_baseball = 0.5 * 0.3 * 1.2 * (np.pi * 0.0369**2) / 0.145 # 0.5*Cd*rho_air*A/m
print(f"k for a real baseball: {k_baseball:.4f} 1/m")
traj_drag = simulate_until_ground(projectile_deriv, state0, 0.005, k_baseball, 2)
plt.figure(figsize=(7, 4))
plt.plot(traj[:, 0], traj[:, 2], '--', label="frictionless")
plt.plot(traj_drag[:, 0], traj_drag[:, 2], '-', label=f"with drag, k={k_baseball:.4f}")
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.legend()
plt.axis('equal')
plt.title("Frictionless vs. realistic baseball drag")
plt.show()
Matches Task #1 (frictionless) run: True k for a real baseball: 0.0053 1/m
Task #3¶
sweeping the drag strength and finding the impact angle and the fraction of time spent ascending
def flight_summary(k, n=2, dt=0.005):
traj = simulate_until_ground(projectile_deriv, state0, dt, k, n)
vx_land, vy_land = traj[-1, 1], traj[-1, 3]
impact_angle = np.degrees(np.arctan2(-vy_land, vx_land))
t_flight = (len(traj) - 1) * dt # slightly approximate due to the final interpolated point
t_apex = np.argmax(traj[:, 2]) * dt
return traj, impact_angle, t_apex, t_flight
k_values = np.concatenate([[0.0], np.geomspace(0.001, 0.05, 12)])
results = [flight_summary(k) for k in k_values]
impact_angles = np.array([r[1] for r in results])
frac_ascending = np.array([r[2] / r[3] for r in results])
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(k_values, impact_angles, 'o-')
axes[0].axhline(45, color='k', ls=':', label="launch angle (45 deg)")
axes[0].axvline(k_baseball, color='r', ls='--', label=f"realistic baseball k={k_baseball:.4f}")
axes[0].set_xlabel("k (1/m)")
axes[0].set_ylabel("impact angle (deg)")
axes[0].legend()
axes[0].set_title("Impact angle vs. drag strength")
axes[1].plot(k_values, frac_ascending, 'o-', color='C1')
axes[1].axhline(0.5, color='k', ls=':', label="frictionless (symmetric)")
axes[1].axvline(k_baseball, color='r', ls='--')
axes[1].set_xlabel("k (1/m)")
axes[1].set_ylabel("fraction of flight time ascending")
axes[1].legend()
axes[1].set_title("Time asymmetry vs. drag strength")
plt.tight_layout()
plt.show()
_, angle_baseball, t_apex_b, t_flight_b = flight_summary(k_baseball)
print(f"At the realistic baseball k={k_baseball:.4f}: "
f"impact angle = {angle_baseball:.1f} deg (vs. 45 deg launch), "
f"ascending fraction = {t_apex_b/t_flight_b:.3f}")
# Is there a k that gets us to a near-vertical (>=80 deg) landing?
steep_idx = np.argmax(impact_angles >= 80)
if impact_angles[steep_idx] >= 80:
print(f"Impact angle first reaches >=80 deg at k={k_values[steep_idx]:.4f}")
else:
print("No k in our tested range reaches an 80 deg impact angle.")
At the realistic baseball k=0.0053: impact angle = 54.4 deg (vs. 45 deg launch), ascending fraction = 0.477 No k in our tested range reaches an 80 deg impact angle.
In the tested range, no drag strength caused the ball to fall "almost vertically" to the floor. However, by extending the tested range (to non-realistic values) can result in a 80° impact angle:
k_extreme = np.array([0.05, 0.1, 0.2, 0.5, 1.0, 2.0])
for k in k_extreme:
_, angle, _, _ = flight_summary(k, n=2)
print(f"k={k:5.2f} impact_angle={angle:5.1f} deg "
f"(={k/k_baseball:5.0f}x the realistic baseball k)")
k= 0.05 impact_angle= 71.2 deg (= 9x the realistic baseball k) k= 0.10 impact_angle= 75.4 deg (= 19x the realistic baseball k) k= 0.20 impact_angle= 78.7 deg (= 38x the realistic baseball k) k= 0.50 impact_angle= 82.0 deg (= 94x the realistic baseball k) k= 1.00 impact_angle= 83.8 deg (= 188x the realistic baseball k) k= 2.00 impact_angle= 85.1 deg (= 377x the realistic baseball k)
Task #4¶
changing the drag exponent n
fig, ax = plt.subplots(figsize=(7, 4.5))
for n in [1, 1.5, 2]:
angles_n = [flight_summary(k, n=n)[1] for k in k_values]
ax.plot(k_values, angles_n, 'o-', label=f"n={n}")
ax.axhline(45, color='k', ls=':', label="launch angle")
ax.set_xlabel("k (1/m)")
ax.set_ylabel("impact angle (deg)")
ax.legend()
ax.set_title("Impact angle vs. drag strength, for different drag laws")
plt.show()
Summary of results¶
- What impact angle did you find at a realistic drag strength, and how does it compare to the launch angle? At $k\approx0.0053\,\mathrm{m^{-1}}$ (a real baseball, $n=2$), the impact angle is about $54$-$55°$, which is noticeably steeper than the $45°$ launch angle. However, it's only a roughly 10° steepening, not the near-vertical drop ($\gtrsim80°$) that the "falls straight down" claim would suggest. Reaching 80° requires $k$ roughly 40–90 times larger than a real baseball's drag coefficient (crossing 80° somewhere between $k\approx0.2$ and $k\approx0.5$), which isn't physically realistic for a ball moving through air.
- Does the time asymmetry tell the same story? Yes, at the realistic $k$, only about 47.7% of the flight time is spent ascending (vs. exactly 50% with no drag), so the ball also spends measurably longer falling than it took to rise.
- What went wrong or surprised you? It was surprising to see that the folklore saying the ball seems to be falling straight down is so largely exaggerated. We assumed that professional athletes have very precise estimates of these effects through continued experimenting. Although it is possible that their experience is real, but not only explained by drag.
- With another hour: It would have been interesting to check whether a spinning ball behaves differently than a non-spinning one. We could also have tested a range of initial velocity values to see how this affects the results.