Solving ODEs

Why this matters

Most physical models are written as rates of change. A capacitor's voltage changes at a rate set by the current; a drug's concentration falls at a rate proportional to how much is left; a mass on a spring accelerates according to the spring force; a reactor's temperature rises with the heat released. These are ordinary differential equations (ODEs). Only the simplest have neat formula solutions. For everything else, including the realistic versions with friction, non-linear springs and coupled reactions, you solve them numerically by stepping forward in time.

Core idea: state, rate, step

Write the system in first-order form:

dy/dt = f(t, y),    y(t0) = y0

y is the state (a number or a NumPy array), and f returns its rate of change. A second-order equation becomes two first-order ones by adding the velocity to the state. For a mass on a spring, m x'' = -k x - c x':

import numpy as np

def spring(t, y, k=4.0, c=0.0, m=1.0):
    x, v = y
    return np.array([v, (-k * x - c * v) / m])     # [dx/dt, dv/dt]

Euler's method

The simplest step: follow the current slope for a time h.

y_{n+1} = y_n + h f(t_n, y_n)

def euler(f, y0, t):
    y = np.zeros((len(t), np.size(y0)))
    y[0] = y0
    for n in range(len(t) - 1):
        h = t[n + 1] - t[n]
        y[n + 1] = y[n] + h * f(t[n], y[n])
    return y

The loop over time steps is unavoidable (each step needs the previous one), but each step is vectorised over the state. Euler is easy to understand and badly inaccurate: its error is proportional to h, so halving the step only halves the error. On the undamped spring it slowly gains energy, and the oscillation grows when it should not.

Fourth-order Runge-Kutta (RK4)

RK4 samples the slope four times per step and averages them:

k1 = f(t,       y)
k2 = f(t + h/2, y + h k1 / 2)
k3 = f(t + h/2, y + h k2 / 2)
k4 = f(t + h,   y + h k3)
y_next = y + h (k1 + 2 k2 + 2 k3 + k4) / 6

Its error scales like h^4: halving the step cuts the error by 16. For the same accuracy it is usually far cheaper than Euler, despite four function calls per step.

Worked example: comparing the two

Exponential decay dy/dt = -0.5 y, y(0) = 1, has the exact solution e^(-0.5 t). Over 10 s with h = 0.1:

method   y(10)        error
exact    0.0067379
Euler    0.0059205    8.2e-4   (12% relative error)
RK4      0.0067379    1.8e-9

RK4 is accurate to about eight significant digits with the same 100 steps.

The real tool: solve_ivp

In practice you use SciPy, which picks step sizes automatically to meet an error tolerance:

from scipy.integrate import solve_ivp

sol = solve_ivp(spring, t_span=(0, 10), y0=[1.0, 0.0],
                method="RK45", rtol=1e-8, atol=1e-10,
                t_eval=np.linspace(0, 10, 501), args=(4.0, 0.3, 1.0))
x = sol.y[0]          # sol.y has shape (n_states, n_times)
print(sol.success, sol.message)

Points to know:

• The function signature is f(t, y): time first. (The older odeint uses f(y, t), a common source of bugs.)
• rtol and atol control accuracy. The defaults (1e-3, 1e-6) are loose; tighten them and check the answer does not change.
• t_eval chooses where you get output; the solver's internal steps are separate.
• events= stops the integration when something happens, such as a projectile hitting the ground.
• Stiff problems (chemical kinetics with very fast and very slow reactions, circuits with tiny time constants) make explicit methods take absurdly small steps. Switch to method="LSODA", "Radau" or "BDF".

Checking a numerical solution

1. Convergence: halve the step (or tighten rtol) and confirm the answer stops changing.
2. Conserved quantities: for an undamped spring, energy 0.5 k x² + 0.5 m v² should stay constant. Drift exposes a poor method.
3. Known limits: set damping to zero, or compare with an analytic case.
4. Units and signs: a decay that grows usually means a sign error in f.

Common mistakes

• Wrong argument order (f(y, t) versus f(t, y)).
• Forgetting to convert a second-order equation to first-order form.
• Trusting Euler with a step that "looks small". Always do a convergence check.
• Returning a list from f and then adding it to arrays: return a NumPy array.
• Loose default tolerances in solve_ivp for long simulations.

Quick recap

• Write dy/dt = f(t, y) with a state vector; second-order becomes first-order by adding velocities.
• Euler: y + h f, error ~ h. RK4: four slopes, error ~ h^4.
• Use solve_ivp(f, (t0, t1), y0, rtol=..., t_eval=...) in real work; use a stiff method when steps collapse.
• Check convergence and conserved quantities.

Try it

In the Do tab you will implement Euler and RK4 steppers in NumPy and test them on exponential decay and a spring.