Ordinary Differential Equations

Most laws of nature are written as differential equations. A pendulum, a chemical reaction, a satellite orbit, a virus spreading through a population — all of them say "the rate of change depends on the current state". A solver turns those local laws into a trajectory you can plot, analyze, and reason about.

This chapter takes you from the simplest ODE solver you can write by hand to the production-grade adaptive solvers in SciPy, and shows you how to recognize when a problem is stiff.

The initial-value problem

$$\frac{dy}{dt} = f(t, y), \qquad y(t_0) = y_0$$

Given the law $f$ and a starting state $y_0$, find $y(t)$ for $t > t_0$. Every solver does the same dance: it takes a step from $t_n$ to $t_{n+1} = t_n + h$, predicting $y_{n+1}$ from $y_n$. The differences are how the prediction is made and how the step size $h$ is chosen.

Euler's method — the simplest possible solver

Pretend the slope is constant over the step: $y_{n+1} = y_n + h \, f(t_n, y_n)$. One function evaluation per step, easy to code, and almost always wrong.

With $h = 0.5$ Euler is wildly off — it can even oscillate. Halve $h$ and the error halves (first-order accuracy). To get one extra correct digit, you need ten times as many steps. Unacceptable for serious work.

Runge–Kutta 4 (RK4)

Probe the slope at four cleverly-chosen points within the step and combine them with the right weights:

$$\begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h/2,\, y_n + h k_1/2) \\ k_3 &= f(t_n + h/2,\, y_n + h k_2/2) \\ k_4 &= f(t_n + h,\, y_n + h k_3) \\ y_{n+1} &= y_n + \tfrac{h}{6}(k_1 + 2 k_2 + 2 k_3 + k_4) \end{aligned}$$

Cost: four evaluations per step. Reward: fourth-order accuracy — halving $h$ shrinks the error by sixteen.

With ten steps of size 0.3, RK4 nails the answer to seven decimals. Euler with the same step size would have produced garbage.

scipy.integrate.solve_ivp — what you should actually use

In production, never hand-roll RK4. SciPy gives you adaptive-step solvers that monitor their own error and adjust $h$ on the fly. The default is RK45 (Dormand–Prince), which is RK4 with a built-in error estimator.

Key arguments:

• rtol, atol — relative and absolute error tolerances per step
• method — RK45 (default), RK23 (cheaper), DOP853 (very high accuracy), Radau/BDF/LSODA (stiff)
• dense_output=True — returns a callable sol.sol(t) you can query at any time, not just at the solver's internal steps
• events — detect zero-crossings of arbitrary functions

Event detection

Often you do not want the trajectory; you want to know when something happens. events lets you catch those moments exactly.

solve_ivp solves the equation $g(t, y) = 0$ to root-finding precision, not just to the nearest time step. The reported event time is essentially exact.

Systems of equations

ODE systems are written as a vector field. The solver doesn't care how big the state is.

You should see the classic oscillatory boom-bust cycle: prey grow, predators boom on the abundant food, prey crash, predators starve, prey recover, and the cycle repeats.

Stiff equations — when explicit methods fail

A stiff ODE has time scales that differ by many orders of magnitude. Explicit methods like RK4 are forced to take absurdly tiny steps for stability even though the solution is changing slowly. The classic example is the chemical kinetics problem Robertson's system.

RK45 either fails or takes thousands of evaluations. BDF and LSODA — implicit methods designed for stiff problems — handle it in tens of steps.

Rule of thumb: if RK45 gets unbearably slow or refuses to converge, switch to BDF or LSODA.

A multi-file simulation: the Lorenz attractor

That butterfly-shaped curve is the visual signature of deterministic chaos. Two trajectories starting $10^{-9}$ apart diverge to completely different parts of the attractor within a few time units — yet the system is fully deterministic.

Challenge: detect the apex of a projectile

Check your understanding

QuestionSelect one

You halve the step size $h$ in your explicit RK4 integrator. By roughly what factor does the global error shrink?

2

4

8

16

QuestionSelect one

Your reaction-network ODE has fast equilibria on timescales of $10^{-6}$ seconds and slow dynamics on timescales of $10^{3}$ seconds. solve_ivp with the default RK45 method either crawls or returns "step size too small". What's the right next step?

Lower the absolute tolerance

Switch to an implicit stiff solver such as BDF or LSODA

Switch to DOP853 for higher order

Add a perturbation to break the symmetry