Implementing Constant Drag in 2D Projectiles

Welcome to the next stage of our simulation journey. In our last lesson, you built an impressive "dashboard" to visualize the flight of a projectile in a perfect, frictionless world. You saw how the plots of position, velocity, and acceleration confirmed the fundamental laws of physics for ideal motion.

Today, we'll take a significant step toward reality by adding a new force to our simulation: air resistance, also known as drag. This force is what prevents a feather from falling as fast as a stone and what makes a real-world baseball trajectory differ from the perfect parabolas of introductory physics. Our goal is to incorporate a simple model of this drag force into our 2D simulation and observe how it transforms the projectile's flight.

What is Drag?

Drag is a resistive force exerted by a fluid (like air) on an object moving through it. Unlike gravity, which always pulls straight down, the drag force always acts in the direction opposite to the object's velocity.

This image compares a projectile's path in a vacuum (dashed line, our previous simulation) with its path when subject to air resistance (solid line). The drag force, $F_{d r a g}$ , always opposes the velocity vector, $v_{0}$ , resulting in a shorter range and a lower maximum height.

As you can see from the diagram, because the velocity vector changes direction throughout the flight, the drag force vector must also change direction. This makes the problem much more complex than the constant force of gravity.

To begin, we will use a simplified model that captures the essential nature of drag. We'll assume the drag force has a constant magnitude, but its direction will always oppose the velocity. This is a powerful first step that introduces the core concepts needed for more advanced aerodynamic models later on.

The Vector Nature of Drag

To correctly model a force that opposes velocity, we need to think in terms of vectors. The "direction of velocity" is represented by a unit vector, often written as $v^$ . A unit vector has a magnitude of 1 and simply points in a certain direction.

This video from the "Let's Code Physics" channel provides an excellent explanation of why drag opposes velocity and introduces the concept of using a unit vector to define its direction.

Projectile Motion 11 - Projectiles with Drag

This video explains the physics of drag and how to represent it as a vector for a simulation. While it eventually builds up to a more complex drag formula, we'll focus on the initial conceptual parts.

Please watch these two segments: The drag force: This part explains that drag always points in the opposite direction of the velocity vector and how that direction changes during the flight. Unit vectors: This section explains what a unit vector is and how to calculate it by dividing a vector by its own magnitude. This is the key mathematical tool we need.

From the video, we can establish our simple drag model. We'll represent the drag force vector, $F_{d r a g}$ , as:

$F_{d r a g} = - k \cdot v^$

Where:

$k$ is a constant that represents the magnitude of our drag force (in Newtons). We'll define this as a new parameter in our simulation.
$v^$ is the unit vector in the direction of the projectile's velocity.

To implement this, we first need to find the components of the unit vector. If the velocity vector is $v = (v_{x}, v_{y})$ , its magnitude (the speed) is $v = v_{x}^{2} + v_{y}^{2}$ . The unit vector is then:

$v^= \frac{v}{v} = (\frac{v _{x}}{v}, \frac{v _{y}}{v})$

Now we can write the components of the drag force:

$F_{d r a g, x} = - k \cdot \frac{v _{x}}{v}$
$F_{d r a g, y} = - k \cdot \frac{v _{y}}{v}$

Finally, we update Newton's Second Law ( $Σ F = m a$ ) to include both gravity and drag.

x-direction: $Σ F_{x} = F_{d r a g, x} ⟹ m a_{x} = - k \frac{v _{x}}{v} ⟹ a_{x} = - \frac{k}{m} \frac{v _{x}}{v}$
y-direction: $Σ F_{y} = F_{g, y} + F_{d r a g, y} ⟹ m a_{y} = - m g - k \frac{v _{y}}{v} ⟹ a_{y} = - g - \frac{k}{m} \frac{v _{y}}{v}$

Notice the crucial difference: the acceleration components $a_{x}$ and $a_{y}$ now depend on the velocity components $v_{x}$ and $v_{y}$ . This creates a "coupled" system where you cannot solve for the motion in one dimension without knowing about the other. This is precisely the kind of problem that numerical methods like the Euler algorithm are designed to solve!

Updating the Python Simulation

Let's modify our script from the previous lesson to include this new drag force.

Add New Parameters: We need to define the projectile's mass (m) and our new constant drag magnitude (k).
Modify the Simulation Loop: The calculation for ax and ay is where the new physics comes in. Inside the loop, for each time step, we must first calculate the current speed v, then use it to find the drag force components and the resulting accelerations.

A small but important detail: what happens if the speed v is zero? We would get a division-by-zero error. We must add a check to ensure that if the speed is very close to zero, we simply set the drag force to zero.

Here is the complete, updated script incorporating the simple drag model. Replace your previous code with this version.

import numpy as np
import matplotlib.pyplot as plt

# --- 1. Setup ---
g = 9.81      # Acceleration due to gravity (m/s^2)
m = 1.0       # Mass of the projectile (kg)
k = 0.5       # Constant drag force magnitude (N)

v0 = 50       # Initial velocity (m/s)
theta_deg = 60  # Launch angle (degrees)
theta_rad = np.radians(theta_deg)

dt = 0.01     # Time step (s)
t_max = 12    # Max simulation time (s)

# --- 2. Initialization ---
n_steps = int(t_max / dt)
t = np.zeros(n_steps + 1)
x = np.zeros(n_steps + 1)
y = np.zeros(n_steps + 1)
vx = np.zeros(n_steps + 1)
vy = np.zeros(n_steps + 1)
ax_arr = np.zeros(n_steps + 1)
ay_arr = np.zeros(n_steps + 1)

# Set initial conditions
x[0] = 0
y[0] = 0
vx[0] = v0 * np.cos(theta_rad)
vy[0] = v0 * np.sin(theta_rad)

# --- 3. Simulation Loop ---
i = 0
while y[i] >= 0 and i < n_steps:
    # --- Calculate forces and acceleration ---
    
    # Calculate current speed
    speed = np.sqrt(vx[i]**2 + vy[i]**2)
    
    # Calculate acceleration components including drag
    if speed > 1e-8: # Avoid division by zero if speed is negligible
        drag_force_x = -k * (vx[i] / speed)
        drag_force_y = -k * (vy[i] / speed)
    else:
        drag_force_x = 0
        drag_force_y = 0
        
    ax = drag_force_x / m
    ay = -g + (drag_force_y / m)
    
    # Store acceleration values
    ax_arr[i] = ax
    ay_arr[i] = ay
    
    # --- Update state using Euler's method ---
    
    # Update velocities
    vx[i+1] = vx[i] + ax * dt
    vy[i+1] = vy[i] + ay * dt
    
    # Update positions
    x[i+1] = x[i] + vx[i+1] * dt
    y[i+1] = y[i] + vy[i+1] * dt
    
    # Update time
    t[i+1] = t[i] + dt
    
    # Increment counter
    i += 1

# Store acceleration for the last point
ax_arr[i] = ax
ay_arr[i] = ay
    
# Trim arrays to actual flight time
num_points = i + 1
t = t[:num_points]
x = x[:num_points]
y = y[:num_points]
vx = vx[:num_points]
vy = vy[:num_points]
ax_arr = ax_arr[:num_points]
ay_arr = ay_arr[:num_points]

# --- 4. Visualization ---
# (The visualization code is identical to the previous lesson)
fig, axs = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Projectile Motion with Simple Drag', fontsize=16)

# Plot 1: Trajectory
axs[0, 0].plot(x, y)
axs[0, 0].set_title("Trajectory")
axs[0, 0].set_xlabel("Range (m)")
axs[0, 0].set_ylabel("Height (m)")
axs[0, 0].grid(True)
axs[0, 0].axis('equal')

# Plot 2: Position vs Time
axs[0, 1].plot(t, x, label='x (horizontal)')
axs[0, 1].plot(t, y, label='y (vertical)')
axs[0, 1].set_title("Position vs. Time")
axs[0, 1].set_xlabel("Time (s)")
axs[0, 1].set_ylabel("Position (m)")
axs[0, 1].grid(True)
axs[0, 1].legend()

# Plot 3: Velocity vs Time
axs[1, 0].plot(t, vx, label='vx (horizontal)')
axs[1, 0].plot(t, vy, label='vy (vertical)')
axs[1, 0].set_title("Velocity vs. Time")
axs[1, 0].set_xlabel("Time (s)")
axs[0, 1].set_ylabel("Velocity (m/s)")
axs[1, 0].grid(True)
axs[1, 0].legend()

# Plot 4: Acceleration vs Time
axs[1, 1].plot(t, ax_arr, label='ax (horizontal)')
axs[1, 1].plot(t, ay_arr, label='ay (vertical)')
axs[1, 1].set_title("Acceleration vs. Time")
axs[1, 1].set_xlabel("Time (s)")
axs[1, 1].set_ylabel("Acceleration (m/s^2)")
axs[1, 1].grid(True)
axs[1, 1].legend()
axs[1, 1].set_ylim(-g*1.5, g*0.5) # Adjust y-limits for better view

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

print(f"Simulation finished with drag force magnitude k = {k} N.")
print(f"Range: {x[-1]:.2f} m")
print(f"Max Height: {np.max(y):.2f} m")
print(f"Time of Flight: {t[-1]:.2f} s")

Analysis and Interpretation

Run the new script. First, try running it with k = 0. The output should be identical to your previous simulation. This is an important verification step.

Now, set k to a value like 0.5 or 1.0 and run it again. Examine your dashboard of plots. What has changed?

Trajectory: The path is no longer a perfect, symmetric parabola. The range and maximum height are both reduced. The descent after the peak is noticeably steeper than the ascent.
Position vs. Time: Both x and y curves are altered. The rate of increase of x is no longer constant.
Velocity vs. Time: The vx plot is no longer a flat line; the horizontal velocity now decreases over time. The vy plot is no longer a straight line but a curve, as the vertical drag component changes.
Acceleration vs. Time: This is perhaps the most telling plot. ax and ay are no longer constant! They change at every time step because they depend on the changing velocity.

Experiment by changing the value of k. A larger k represents more drag. You will see the effects become much more pronounced, with the trajectory becoming shorter and more distorted.

Conclusion

You have successfully added a non-trivial new force to your physics engine! This is a major step up in complexity and realism. You've seen firsthand how a velocity-dependent force fundamentally changes the motion of an object.

Key takeaways from this lesson:

Drag is a resistive force that always opposes an object's velocity.
We can model this behavior using a unit vector, $v^$ , to separate the direction of motion from its magnitude (speed).
A simple drag model can be created with a constant force magnitude, $k$ , applied in the direction of $- v^$ .
Adding a velocity-dependent force makes the acceleration components non-constant, demonstrating the power of numerical simulation for problems without simple analytical solutions.

In the next lesson, we will explore another fundamental aspect of simulation: comparing our numerical results to known analytical solutions for the cases where they exist (like our original, no-drag model) to validate the accuracy of our code.

Recommended resources

What is Drag?

The Vector Nature of Drag

Updating the Python Simulation

Analysis and Interpretation

Conclusion