1D Free Fall Simulation with Euler's Method

The Anatomy of a Physics Simulation

At its core, a numerical simulation is a program that steps through time, repeatedly applying the laws of physics to update the state of a system. For our purposes, this process can be broken down into four main parts:

1. Setup: Define the physical constants (like gravity, g) and simulation parameters (like the time step, dt, and total duration).
2. Initialization: Create data structures (like NumPy arrays) to hold the simulation results and set the initial state of the object (its starting position and velocity).
3. The Simulation Loop: This is the heart of the simulation. It's a loop that iterates through time, calculating the new velocity and position at each step using the Euler method update rules.
4. Analysis & Visualization: Once the loop is complete, you process or plot the results to understand what happened.

We will now build a script that implements these four parts.

Implementing the Simulation in Python

Let's translate this structure into Python code. We'll simulate dropping a ball from a height of 100 meters and track its motion for 5 seconds.

1. Setup and Initialization

First, we import the necessary libraries—numpy for numerical operations and matplotlib.pyplot for plotting, both of which you've used in Module 1. Then, we define our constants and parameters.

import numpy as np
import matplotlib.pyplot as plt




# --- 1. Setup ---
g = 9.8      # Acceleration due to gravity (m/s^2)
y0 = 100     # Initial height (m)
v0 = 0       # Initial velocity (m/s)
t_max = 5    # Total simulation time (s)
dt = 0.1     # Time step (s)

Next, we initialize the arrays that will store our results. We need one array for time, one for position (y), and one for velocity (v). We'll use np.linspace for the time array and np.zeros to create empty arrays for position and velocity, which we'll fill in during the simulation.

# --- 2. Initialization ---
# Calculate the number of steps required
n_steps = int(t_max / dt)




# Create arrays. We need n_steps + 1 points to store the initial state plus all the steps.
t = np.linspace(0, t_max, n_steps + 1)
y = np.zeros(n_steps + 1)
v = np.zeros(n_steps + 1)




# Set the initial conditions in the first element of our arrays
y[0] = y0
v[0] = v0

2. The Simulation Loop

Now for the core of the program. We will use a for loop that runs for n_steps. Inside the loop, for each time step, we'll apply the semi-implicit Euler method update rules we discussed:

1. Calculate the new velocity: $v_{i+1}=v_i+a\cdot \Delta t$
2. Calculate the new position using the new velocity: $y_{i+1}=y_i+v_{i+1}\cdot \Delta t$

In our case, the acceleration $a$ is simply the constant $-g$.

The following resource provides excellent, clear Python code for this exact task. Pay close attention to how it sets up the arrays and uses a for loop to iterate and apply the Euler method.

Notice that the resource calculates the new position using the old velocity (x[i] = x[i-1] + dt * v[i-1]). This is known as the explicit Euler method. As we established previously, we will use the semi-implicit version, which updates the position using the newly calculated velocity. This small change often improves the stability and energy conservation of the simulation.

Here is the Python implementation of our simulation loop using the semi-implicit method:

# --- 3. Simulation Loop ---
for i in range(n_steps):



    # The only force is gravity, so acceleration is constant
    a = -g
    



    # Update velocity (v_new = v_old + a * dt)
    v[i+1] = v[i] + a * dt
    



    # Update position using the NEW velocity (y_new = y_old + v_new * dt)
    y[i+1] = y[i] + v[i+1] * dt

This loop takes the values at step i, calculates the new state, and stores it in the i+1 position of our arrays, ready for the next iteration.

3. Visualizing the Results

A simulation isn't complete until you can see the results. We can use matplotlib to plot the height and velocity of the ball as a function of time. This will give us an intuitive understanding of the object's motion.

# --- 4. Visualization ---
plt.figure(figsize=(12, 6))




# Subplot 1: Position vs. Time
plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title("Position vs. Time")
plt.xlabel("Time (s)")
plt.ylabel("Height (m)")
plt.grid(True)




# Subplot 2: Velocity vs. Time
plt.subplot(1, 2, 2)
plt.plot(t, v, 'r-') # 'r-' for a red line
plt.title("Velocity vs. Time")
plt.xlabel("Time (s)")
plt.ylabel("Velocity (m/s)")
plt.grid(True)

plt.tight_layout()
plt.show()

The Complete Script and Your First Simulation

Here is the complete script. I encourage you to copy this into your Python environment, run it, and see the output for yourself.

import numpy as np
import matplotlib.pyplot as plt




# --- 1. Setup ---
g = 9.8      # Acceleration due to gravity (m/s^2)
y0 = 100     # Initial height (m)
v0 = 0       # Initial velocity (m/s)
t_max = 5    # Total simulation time (s)
dt = 0.1     # Time step (s)




# --- 2. Initialization ---
n_steps = int(t_max / dt)
t = np.linspace(0, t_max, n_steps + 1)
y = np.zeros(n_steps + 1)
v = np.zeros(n_steps + 1)

y[0] = y0
v[0] = v0




# --- 3. Simulation Loop ---
for i in range(n_steps):
    a = -g
    v[i+1] = v[i] + a * dt
    y[i+1] = y[i] + v[i+1] * dt




    # Optional: Stop the simulation if the ball hits the ground
    if y[i+1] < 0:
        y[i+1] = 0 # Set final height to 0 to prevent it from going "underground"



        # Truncate arrays to the point of impact for cleaner plots
        t = t[:i+2]
        y = y[:i+2]
        v = v[:i+2]
        break




# --- 4. Visualization ---
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(t, y)
plt.title("Position vs. Time")
plt.xlabel("Time (s)")
plt.ylabel("Height (m)")
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(t, v, 'r-')
plt.title("Velocity vs. Time")
plt.xlabel("Time (s)")
plt.ylabel("Velocity (m/s)")
plt.grid(True)

plt.tight_layout()
plt.show()

When you run the code, observe the plots. The "Position vs. Time" graph should show a downward-curving parabola, which is the classic trajectory of an object under constant acceleration. The "Velocity vs. Time" graph should show a straight line with a constant negative slope. The slope of this line is, of course, the acceleration due to gravity, $-g$.

This video provides a great visual and conceptual walkthrough of what our code is doing—approximating the motion by taking small, linear steps.

An Alternative: The while Loop

Our simulation used a for loop, which is perfect when you want to run the simulation for a fixed amount of time. An alternative approach, especially useful when your simulation should end based on a condition (like a rocket running out of fuel or an object hitting the ground), is to use a while loop.

For example, we could have structured our loop like this:

#
# ... (Setup and initialization as before) ...
#




# Using a while loop
i = 0
while y[i] > 0 and i < n_steps:
    a = -g
    v[i+1] = v[i] + a * dt
    y[i+1] = y[i] + v[i+1] * dt
    i += 1




# Truncate arrays to the actual number of steps taken
t = t[:i+1]
y = y[:i+1]
v = v[:i+1]

This version stops automatically once the height y becomes zero or less. This is a common and powerful pattern in simulation programming that you will find useful later on. The document below uses this while loop approach.

Conclusion

In this lesson, you have successfully built your first physics simulation. You translated the physical model of 1D free fall into a complete, working Python program using Euler's method. This process is the template we will expand upon throughout the rest of the course.

Key takeaways from this lesson are:

• A simulation consists of setup, initialization, a main loop, and visualization.
• The simulation loop iteratively applies physics rules (the Euler update equations) over small time steps (dt) to predict the future state of a system.
• Both for loops (for fixed time) and while loops (for state-based conditions) are effective ways to structure a simulation.
• Visualizing the results with matplotlib is essential for interpreting the simulation's output.

You now have a solid foundation in 1D simulation. In the next lesson, we will expand this code to handle two dimensions, allowing us to simulate the full arc of a projectile.