Numerical simulation of projectile motion with the explicit Euler method

In this post I build a numerical simulation of projectile motion using the explicit Euler method.
The purpose is not only to draw a curve, but to show the complete numerical workflow:
mathematical model, time discretization, implementation, comparison with the analytical solution,
error analysis and convergence.

1. Physical model

I consider the simplest possible projectile model. The motion takes place in the vertical plane,
the projectile is treated as a point mass, air resistance is neglected, gravity is constant,
and the ground is located at y = 0.

Under these assumptions the only force is gravity. Therefore there is no horizontal acceleration,
while the vertical acceleration is constant and directed downward:

g \simeq 9.81\,\mathrm{m/s^2}

2. Mathematical model

The position is described by $x(t)$ and $y(t)$ , while the velocity components are
$v_x(t)$ and $v_y(t)$ . The equations of motion are:

\begin{cases}\dfrac{dx}{dt}=v_x\\[4pt]\dfrac{dy}{dt}=v_y\\[4pt]\dfrac{dv_x}{dt}=0\\[4pt]\dfrac{dv_y}{dt}=-g\end{cases}

The state vector is:

\mathbf{u}(t)=\begin{pmatrix}x(t)\\y(t)\\v_x(t)\\v_y(t)\end{pmatrix}

So the system can be written compactly as:

\dfrac{d\mathbf{u}}{dt}=f(\mathbf{u},t)=\begin{pmatrix}v_x\\v_y\\0\\-g\end{pmatrix}

3. Time discretization

This is an ordinary differential equation problem, not a partial differential equation problem.
Therefore I discretize time, not space. The coordinates x and y
are state variables updated at each time step; they are not points of a fixed spatial grid.

t_{n+1}=t_n+\Delta t,\qquad \mathbf{u}_n\approx \mathbf{u}(t_n)

4. Python model function

The following function implements the right-hand side of the ODE system.

import numpy as np

def f(u, t):
    x, y, vx, vy = u

    dxdt = vx
    dydt = vy
    dvxdt = 0
    dvydt = -9.81

    return [dxdt, dydt, dvxdt, dvydt]

5. Explicit Euler method

The explicit Euler method advances the state using the derivative evaluated at the current step:

\mathbf{u}_{n+1}=\mathbf{u}_n+\Delta t\,f(\mathbf{u}_n,t_n)

This method is simple and transparent. It is useful for learning because every update is visible,
but it is only first-order accurate.

def euler_step(u, rhs, t, dt):
    """Single explicit Euler step."""
    return np.asarray(u, dtype=float) + dt * np.asarray(rhs(u, t), dtype=float)


def integrate_euler(u0, rhs, t0, t_end, dt):
    """Time integration loop using the explicit Euler method."""
    n = int(np.ceil((t_end - t0) / dt))

    t = t0 + dt * np.arange(n + 1)
    u = np.zeros((n + 1,) + np.shape(u0), dtype=float)

    u[0] = u0

    for k in range(n):
        u[k + 1] = euler_step(u[k], rhs, t[k], dt)

    return t, u

6. Initial conditions

The projectile starts from the origin with initial speed $v_0=30\,\mathrm{m/s}$ and launch angle
$\theta=45^\circ$ . The initial velocity components are:

v_{x0}=v_0\cos(\theta)

v_{y0}=v_0\sin(\theta)

x0 = 0.0
y0 = 0.0

v0 = 30.0
theta_deg = 45.0
theta = np.deg2rad(theta_deg)

vx0 = v0 * np.cos(theta)
vy0 = v0 * np.sin(theta)

u0 = np.array([x0, y0, vx0, vy0], dtype=float)

t0 = 0.0
t_end = 5.0
dt = 0.01

t, u = integrate_euler(u0, f, t0, t_end, dt)

x = u[:, 0]
y = u[:, 1]
vx = u[:, 2]
vy = u[:, 3]

7. Numerical trajectory

The following code plots the part of the trajectory above the ground.

import matplotlib.pyplot as plt

mask = y >= 0

plt.figure(figsize=(8, 5))
plt.plot(x[mask], y[mask], label="Explicit Euler")
plt.axhline(0, linestyle="--", label="Ground")

plt.xlabel("Horizontal position x [m]")
plt.ylabel("Vertical position y [m]")
plt.title("Projectile trajectory - Explicit Euler method")
plt.grid(True)
plt.legend()
plt.show()

Projectile trajectory computed with the explicit Euler method.

8. Analytical solution

For the ideal model without air resistance, the analytical solution is known:

x(t)=x_0+v_{x0}t

y(t)=y_0+v_{y0}t-\dfrac{1}{2}gt^2

This makes the example especially useful, because the numerical result can be compared with an exact reference.

def analytical_solution(t, x0, y0, vx0, vy0, g=9.81):
    x_exact = x0 + vx0 * t
    y_exact = y0 + vy0 * t - 0.5 * g * t**2
    vx_exact = vx0 * np.ones_like(t)
    vy_exact = vy0 - g * t

    return x_exact, y_exact, vx_exact, vy_exact


x_exact, y_exact, vx_exact, vy_exact = analytical_solution(
    t, x0, y0, vx0, vy0
)

error_y = np.abs(y - y_exact)
max_error_y = np.max(error_y)

Comparison between explicit Euler and the analytical solution.

9. Error analysis

The vertical absolute error is:

E_y(t_n)=\left|y_{\mathrm{num}}(t_n)-y_{\mathrm{exact}}(t_n)\right|

For this simulation, using Δt = 0.01 s, the maximum vertical error is approximately:

E_{y,\max}\simeq 0.2453\,\mathrm{m}

Absolute error on the vertical coordinate y.

10. Convergence test

To test convergence, I repeat the simulation using smaller and smaller time steps:

\Delta t=0.2,\;0.1,\;0.05,\;0.025,\;0.0125

The estimated convergence order is computed as:

p=\dfrac{\log(E_{\Delta t}/E_{\Delta t/2})}{\log(2)}

For the explicit Euler method, the expected global order is:

p\simeq 1

dt_values = [0.2, 0.1, 0.05, 0.025, 0.0125]
errors = []

for dt_test in dt_values:
    t_test, u_test = integrate_euler(u0, f, t0, t_end, dt_test)

    y_num = u_test[:, 1]
    _, y_ex, _, _ = analytical_solution(t_test, x0, y0, vx0, vy0)

    error = np.max(np.abs(y_num - y_ex))
    errors.append(error)

for i in range(len(errors) - 1):
    p = np.log(errors[i] / errors[i + 1]) / np.log(2)
    print(f"Estimated order between {dt_values[i]} and {dt_values[i+1]}: {p}")

Convergence test for the explicit Euler method.

11. Stability and accuracy

For this particular model, stability is not the main difficulty. The acceleration is constant and the system is not stiff.
A large time step does not usually make the solution explode, but it makes the trajectory visibly less accurate.

Effect of the time step on the projectile trajectory.

12. Numerical results

Theoretical flight time: 4.325 s
Theoretical range: 91.743 m
Theoretical maximum height: 22.936 m
Maximum vertical error with Δt = 0.01 s: 0.2453 m
Estimated convergence order: approximately 1

13. Conclusion

This example shows the basic structure of a numerical simulation.

First, the physical problem is translated into a system of ordinary differential equations.
Then time is discretized, and the explicit Euler method is used to advance the state step by step.
Finally, the numerical solution is checked through comparison with the analytical solution.

The explicit Euler method is not the most accurate method, but it is excellent for learning because it makes the role of the time step,
the numerical error and the convergence behaviour very clear.