Projectile Motion Simulation

The Story Told by Graphs

Before we write the code, let's establish what we expect to see. In physics, graphs of position, velocity, and acceleration versus time tell a complete story of an object's motion. The relationship between them is fundamental.

• The slope of a position-time graph at any point gives the instantaneous velocity.
• The slope of a velocity-time graph at any point gives the instantaneous acceleration.

This video from the "Professor Dave Explains" channel provides an excellent conceptual overview of this graphical analysis.

Based on this and your physics knowledge, what should our plots for simple projectile motion (no air drag) look like?

1. Acceleration vs. Time: Since gravity is the only force, horizontal acceleration ($a_x$) should be zero, and vertical acceleration ($a_y$) should be a constant $-g$. The plots should be flat horizontal lines.
2. Velocity vs. Time: Since $a_x=0$, the horizontal velocity ($v_x$) should be constant. Since $a_y$ is constant, the vertical velocity ($v_y$) should decrease in a straight line.
3. Position vs. Time: Since $v_x$ is constant, the horizontal position ($x$) should increase in a straight line. Since $v_y$ is changing linearly, the vertical position ($y$) should follow a parabolic curve.

Our goal is to write code that generates plots to verify these expectations.

Modifying the Code to Plot Everything

Let's modify the script from the previous lesson. We need to make two main changes:

1. Store the acceleration values during the simulation so we can plot them later.
2. Use Matplotlib's subplots feature to create a grid of several plots in a single figure.

1. Storing Acceleration

In our previous simulation loop, we calculated ax and ay on the fly but didn't save them. Let's add arrays to store them, just like we did for position and velocity.

First, add ax_arr and ay_arr to your initialization section:

# Initialization
...
vx = np.zeros(n_steps + 1)
vy = np.zeros(n_steps + 1)
ax_arr = np.zeros(n_steps + 1) # Add this
ay_arr = np.zeros(n_steps + 1) # Add this
...

Next, inside the while loop, store the acceleration values for the current step i.

# Simulation Loop
while y[i] >= 0 and i < n_steps:



    # Acceleration components
    ax = 0
    ay = -g
    



    # Store acceleration values
    ax_arr[i] = ax # Add this
    ay_arr[i] = ay # Add this




    # Update velocities
    ...

Finally, we must trim these new arrays at the end, along with the others.

# Trim arrays
...
vx = vx[:num_points]
vy = vy[:num_points]
ax_arr = ax_arr[:num_points] # Add this
ay_arr = ay_arr[:num_points] # Add this

2. Creating a Plotting Dashboard

Instead of one plot, we'll create a 2x2 grid of plots. This gives us a comprehensive "dashboard" view of the flight. The plt.subplots() function is perfect for this. It returns both a figure object (the whole window) and an axes object, which in this case will be a 2D array corresponding to our grid. We can then access each individual plot using indices like axs[0, 0].

This tutorial on plotting with Python provides a great example of creating subplots.

By adapting that example, we can create a much more informative visualization section for our script.

The Complete Script for Visualization

Here is the complete, updated script. It includes the modifications to store acceleration and the new visualization section that generates a 2x2 dashboard of plots. Replace the previous script with this one and run it.

import numpy as np
import matplotlib.pyplot as plt




# --- 1. Setup ---
g = 9.81     # Acceleration due to gravity (m/s^2)
v0 = 50      # Initial velocity (m/s)
theta_deg = 60 # Launch angle (degrees)
theta_rad = np.radians(theta_deg) # Convert to radians

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:



    # Acceleration components
    ax = 0
    ay = -g
    



    # Store acceleration values
    ax_arr[i] = ax
    ay_arr[i] = ay
    



    # 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] = 0
ay_arr[i] = -g
    



# 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 ---
# Create a 2x2 grid of subplots
fig, axs = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Projectile Motion Analysis', fontsize=16)




# Plot 1: Trajectory (y vs x)
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[1, 0].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(-10, 1) # Set y-axis limits for better viewing




# Adjust layout to prevent titles from overlapping
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

print(f"Simulation finished.")
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")

When you run this script, examine the output plots. Do they match the theoretical expectations we discussed? You should see the flat lines for acceleration, the constant vx and linear vy, and the linear x and parabolic y positions. This visual confirmation is powerful—it proves that our simple numerical simulation correctly reproduces the fundamental physics of projectile motion.

Conclusion

In this lesson, you've moved beyond simply running a simulation to analyzing its output in detail. Visualizing the data is as important as generating it, as it allows for verification, debugging, and deeper physical insight.

Here are the key takeaways:

• The state of a moving object can be fully described by plotting its position, velocity, and acceleration against time.
• These plots are physically related: velocity is the slope of position, and acceleration is the slope of velocity.
• Matplotlib's subplots feature is a powerful tool for creating a dashboard of related graphs to analyze simulation data.
• By plotting the kinematic variables, we can visually confirm that our simulation behaves according to the laws of physics.

You now have a robust framework for simulating and analyzing 2D motion. In our next lesson, we will make our simulation more realistic by introducing a new force: air resistance. You will see how adding a drag force, which depends on velocity, changes these "perfect" graphs into something much more complex and interesting.