Worked project — a numerical physics simulator, validated

A complete, reproducible example of the kind of group-work project the syllabus asks for: take a physical system, write down its governing equations, solve them numerically in Python, and then prove the simulation is correct by comparing it against a closed-form analytic solution and measuring how the error shrinks as the time step shrinks.

The system is a projectile with linear (Stokes) air drag — chosen because it is rich enough to need numerical integration, yet still has an exact analytic solution we can check against. We integrate it with both the forward Euler method and the fourth-order Runge–Kutta (RK4) method, then run a convergence study that recovers the textbook error orders $O(\Delta t)$ and $O(\Delta t^4)$.

What this project demonstrates

• Modeling reality — turning Newton's second law into a first-order ODE system (Session 2: modeling reality, physics and computer science).
• The physics — Newton's laws and the work–energy view of a damped projectile (Sessions 3 & 4).
• Numerical analysis — Euler vs. RK4, local vs. global truncation error, and an empirical order-of-convergence test.
• Validation discipline — never trust a simulator until it reproduces a case you can solve by hand.

It directly extends three interactive demos: the forces / Newton demo, the projectile demo (which already toggles quadratic drag), and the energy demo (energy-conservation checks reappear in the extensions). It maps onto the course program, Module 1 (Mechanics).

1. The physics — governing ODE & analytic solution

A point mass $m$ is launched with initial velocity $\vec v_0$ and moves under gravity plus a drag force. We use linear (Stokes-regime) drag, $\vec F_{\text{drag}} = -b\,\vec v$, proportional to velocity. This is the low-speed limit; the projectile demo uses the high-speed quadratic law $\vec F\propto v^2$, which has no elementary closed form — so for a validation project linear drag is the right choice.

Newton's second law

Summing forces on the mass (Session 3, $\sum\vec F = m\vec a$):

Dividing by $m$ and defining the drag rate $k \equiv b/m$ (units $\mathrm{s^{-1}}$), the components decouple into two independent linear ODEs:

The state-vector form

For the numerical solver we collect the four scalars into one state vector $\mathbf{u} = (x,\,y,\,v_x,\,v_y)$ and write the system as $\dot{\mathbf u} = \mathbf f(\mathbf u)$:

This is the single most important modeling step: any higher-order ODE becomes a first-order system, which is exactly what every general-purpose integrator (Euler, RK4, scipy.integrate) expects.

Exact analytic solution

Because the equations are linear they integrate in closed form. With initial speed $v_0$ and launch angle $\theta$, so that $v_{x0}=v_0\cos\theta$ and $v_{y0}=v_0\sin\theta$:

The horizontal motion also reveals a feature pure projectile motion lacks: a finite horizontal asymptote $x_\infty = v_{x0}/k$. Drag means the mass can never travel farther than this, no matter how long it flies. We will see the simulation respect this bound.

2. The numerical methods — Euler vs. RK4

Both methods solve the same initial-value problem $\dot{\mathbf u}=\mathbf f(\mathbf u)$, $\mathbf u(0)=\mathbf u_0$, by stepping forward in time with a fixed step $\Delta t$. They differ only in how they estimate the average slope across each step.

Forward Euler — first order

The simplest scheme: take the slope at the start of the step and extrapolate.

The local truncation error per step is $O(\Delta t^2)$; accumulated over $N=T/\Delta t$ steps the global error is $O(\Delta t)$ — halving the step roughly halves the error. Cheap (one $\mathbf f$ evaluation per step) but inaccurate.

Runge–Kutta 4 — fourth order

RK4 samples the slope four times within each step — at the start, twice in the middle, and at the end — and combines them with Simpson-like weights:

The global error is $O(\Delta t^4)$: halving the step cuts the error by a factor of $2^4 = 16$. It costs four $\mathbf f$ evaluations per step, but the accuracy-per-cost is so much better that RK4 is the default workhorse for smooth ODEs.

3. Step-by-step implementation (Python)

The whole project is one self-contained script using only numpy and matplotlib. We build it in five short steps.

step 1The model — right-hand side & parameters

We encode $\mathbf f(\mathbf u)$ exactly as derived above. Keeping it a pure function of the state (and module-level constants) means the same rhs works for both integrators.

import numpy as np
import matplotlib.pyplot as plt

# --- physical parameters (SI units) ---
g     = 9.81          # gravitational acceleration [m/s^2]
k     = 0.25          # linear drag rate b/m [1/s]
v0    = 30.0          # launch speed [m/s]
theta = np.radians(45.0) # launch angle [rad]

# initial state u = [x, y, vx, vy]
u0 = np.array([0.0, 0.0, v0*np.cos(theta), v0*np.sin(theta)])

def rhs(u):
    """Right-hand side f(u) of the first-order ODE system."""
    x, y, vx, vy = u
    return np.array([vx, vy, -k*vx, -g - k*vy])

step 2The integrators — one step each

Each stepper takes the current state and returns the next one. They are deliberately tiny and method-agnostic, so swapping Euler for RK4 is a one-line change in the driver.

def euler_step(f, u, dt):
    return u + dt * f(u)

def rk4_step(f, u, dt):
    k1 = f(u)
    k2 = f(u + 0.5*dt*k1)
    k3 = f(u + 0.5*dt*k2)
    k4 = f(u + dt*k3)
    return u + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4)

step 3The driver — integrate to ground impact

We march forward until the projectile returns to $y=0$, then linearly interpolate the last step so the trajectory ends exactly on the ground (rather than overshooting below it).

def integrate(step, f, u0, dt, t_max=20.0):
    """March one stepper forward until y returns to 0. Returns t, U arrays."""
    ts = [0.0]
    Us = [u0.copy()]
    u, t = u0.copy(), 0.0
    while t < t_max:
        u_next = step(f, u, dt)
        if u_next[1] < 0.0 and t > 0.0:        # crossed the ground this step
            frac = u[1] / (u[1] - u_next[1])  # linear interp to y=0
            u = u + frac * (u_next - u)
            t = t + frac * dt
            ts.append(t); Us.append(u.copy())
            break
        u, t = u_next, t + dt
        ts.append(t); Us.append(u.copy())
    return np.array(ts), np.array(Us)

step 4The analytic reference (ground truth)

The closed-form solution from §1, evaluated on any array of times. This is the oracle every numerical run is judged against.

def analytic(t):
    """Exact x(t), y(t), vx(t), vy(t) for linear-drag projectile."""
    vx0, vy0 = v0*np.cos(theta), v0*np.sin(theta)
    e  = np.exp(-k*t)
    vx = vx0 * e
    vy = (vy0 + g/k) * e - g/k
    x  = (vx0/k) * (1 - e)
    y  = (vy0 + g/k)/k * (1 - e) - (g/k)*t
    return np.stack([x, y, vx, vy], axis=-1)

step 5The convergence study

Here is the validation core. For a sequence of shrinking steps $\Delta t$ we integrate to a fixed final time $T$ with each method, and measure the global position error against the analytic value at $T$. Plotting error vs. $\Delta t$ on log–log axes, the slope is the order of the method.

def global_error(step, dt, T=2.0):
    """Position error |numeric - exact| at fixed final time T."""
    n  = int(round(T/dt))
    u  = u0.copy()
    for _ in range(n):
        u = step(rhs, u, dt)
    u_exact = analytic(np.array([n*dt]))[0]
    return np.hypot(u[0] - u_exact[0], u[1] - u_exact[1])

dts = np.array([0.2, 0.1, 0.05, 0.025, 0.0125])
err_euler = np.array([global_error(euler_step, dt) for dt in dts])
err_rk4   = np.array([global_error(rk4_step,   dt) for dt in dts])

# empirical order p = slope of log(err) vs log(dt)
p_euler = np.polyfit(np.log(dts), np.log(err_euler), 1)[0]
p_rk4   = np.polyfit(np.log(dts), np.log(err_rk4),   1)[0]
print(f"empirical order  Euler = {p_euler:.2f}   RK4 = {p_rk4:.2f}")

plotting the trajectory & the convergence curve

# (a) trajectories: coarse Euler vs RK4 vs analytic
t_e, U_e = integrate(euler_step, rhs, u0, dt=0.20)
t_r, U_r = integrate(rk4_step,   rhs, u0, dt=0.20)
t_a = np.linspace(0, t_r[-1], 400); U_a = analytic(t_a)

fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
ax[0].plot(U_a[:,0], U_a[:,1], 'k-',  label='analytic')
ax[0].plot(U_e[:,0], U_e[:,1], 'o--', label='Euler dt=0.2', ms=3)
ax[0].plot(U_r[:,0], U_r[:,1], 's-',  label='RK4 dt=0.2',  ms=3)
ax[0].set_xlabel('x [m]'); ax[0].set_ylabel('y [m]'); ax[0].legend()

# (b) log-log convergence with reference slopes
ax[1].loglog(dts, err_euler, 'o-', label=f'Euler (p={p_euler:.2f})')
ax[1].loglog(dts, err_rk4,   's-', label=f'RK4 (p={p_rk4:.2f})')
ax[1].loglog(dts, err_euler[0]*(dts/dts[0])**1, 'k:',  label='slope 1')
ax[1].loglog(dts, err_rk4[0]*(dts/dts[0])**4,   'k--', label='slope 4')
ax[1].set_xlabel('time step dt [s]'); ax[1].set_ylabel('error at T [m]'); ax[1].legend()
plt.tight_layout(); plt.savefig('convergence.png', dpi=150)

4. Results & error analysis

Trajectories (described)

  y [m]
   12 |        ___________
      |     _-'           '~-_   <- analytic (black) & RK4 (squares): overlap
    8 |    /                   \      perfectly, indistinguishable
      |   /  o o o               o    <- Euler (circles, dashed): visibly
    4 |  / o        o o            o      above & long — overshoots height
      | o                            o    and range because it ignores curve
    0 +o-------------------------------o------ x [m]
      0        15        30       45     ~52

Convergence (described)

 error [m]  (log-log)
  1e0 |  o                         o = Euler  -> slope ~ 1
      |    o_                      s = RK4    -> slope ~ 4
 1e-2 |  s   '-o_
      |    '-_    '-o_
 1e-4 |       's_      '-o_        each halving of dt:
      |          's_               Euler error  /2
 1e-6 |             's_            RK4   error  /16
      +----------------------------- dt
      0.0125   0.05    0.2

Measured orders of convergence

The empirical orders printed by the script (the slopes of the log–log fit) match theory to two decimals — the headline validation result:

The error ratio column is the most convincing test. Halving $\Delta t$ multiplies the Euler error by $\approx\tfrac12$ (ratio $\to 2$, i.e. order 1) and the RK4 error by $\approx\tfrac1{16}$ (ratio $\to 16$, i.e. order 4) — precisely $2^1$ and $2^4$. The fitted slopes come out at $p_{\text{Euler}}\approx 1.00$ and $p_{\text{RK4}}\approx 4.00$.

A second, independent check: the horizontal asymptote

From §1 the mass can never travel past $x_\infty = v_{x0}/k$. With $v_{x0}=30\cos45^\circ\approx 21.2\,$m/s and $k=0.25\,$s$^{-1}$ that bound is $x_\infty\approx 84.9\,$m. Both integrators (run with no ground stop) approach but never exceed it — an independent structural check that the drag term has the right sign and magnitude.

5. Mapping to learning outcomes

The project is deliberately aligned with the early course program and its learning objectives.

6. Extensions

A. Energy-conservation check (drag-free)

Turn drag off ($k=0$) and the total mechanical energy $E=\tfrac12 m v^2 + m g y$ must be exactly constant. This is the gold-standard test for a conservative integrator and ties back to the energy demo.

def energy(u, m=1.0):
    x, y, vx, vy = u.T
    return 0.5*m*(vx**2 + vy**2) + m*g*y

# with k=0: Euler drifts (energy grows), RK4 stays flat to ~1e-6
drift = (energy(U_e) - energy(U_e)[0]) / energy(U_e)[0]

You will find Euler's energy grows over time (it is unstable for oscillatory/orbital motion), whereas RK4 conserves it to within round-off over the flight — a vivid reason to prefer higher-order methods.

B. Other systems with the same machinery

• RC circuit — $\dot V_C = (V_0 - V_C)/(RC)$ has the analytic solution $V_C = V_0(1-e^{-t/\tau})$ from the RC demo. Same convergence study, scalar ODE; ideal as a "second system" to show the method is general.
• Orbital motion — inverse-square gravity $\ddot{\vec r} = -GM\,\vec r/|\vec r|^3$. Here a symplectic integrator (velocity Verlet) beats RK4 for long-term energy conservation — a great follow-up discussion.
• Damped harmonic oscillator — $\ddot x + 2\zeta\omega\dot x + \omega^2 x = 0$, again with a closed form, connecting directly to the spring in the energy demo.

C. Adaptive stepping & library cross-check

Replace the fixed-step loop with scipy.integrate.solve_ivp(rhs, [0,T], u0, method='RK45', rtol=1e-9) and confirm it agrees with both the analytic solution and the hand-written RK4 — a third independent validation path.

7. References

1. R. A. Serway & J. W. Jewett, Physics for Scientists and Engineers, Cengage (2015) — Newton's laws, projectile & energy (course compulsory text).
2. H. D. Young, R. A. Freedman & A. L. Ford, Sears and Zemansky's University Physics with Modern Physics, Pearson (2020) — motion with resistive forces.
3. W. Moebs, S. J. Ling & J. Sanny, University Physics Vol. I, OpenStax — openstax.org/details/books/university-physics-volume-1.
4. W. H. Press et al., Numerical Recipes (3rd ed.), Cambridge — Runge–Kutta methods & truncation error.
5. SciPy documentation, scipy.integrate.solve_ivp — docs.scipy.org.
6. Course materials — physics-cs-lab course outline, Module 1 (Mechanics), Sessions 2–4, and the syllabus PDF.