Numerical Solutions of ODEs
Almost every law of nature is a statement about rates. A population grows at a rate
proportional to its size; a hot coffee cools at a rate proportional to how much hotter it is than the
room; a pendulum's angle accelerates in proportion to its own sine. Each of these is an
differential
equation — a rule for the slope y' of an unknown function in
terms of where you are. The initial value problem packages this up as
y'(t) = f(t, y), \qquad y(t_0) = y_0.
You are told the slope f(t,y) at every point of the plane, and one
starting point (t_0, y_0). From there the solution is a single curve threading
through that field of slopes. The trouble is that for all but a lucky handful of
f, there is no closed-form formula for that curve — no
combination of \exp, \sin, powers and logs that
satisfies the equation. The three-body problem, the weather, the flow over a wing: none has a tidy
answer you can write down.
So we do the only thing we can. We stand at the start, read off the slope, and march the
solution forward in tiny steps of size h, trading the unattainable
exact curve for a chain of short straight guesses. This page is about how to take those steps well: the
humble Euler method, why halving the step only halves its error, and the mighty
Runge–Kutta 4, which turns that same halving into a sixteen-fold gain.
Follow the tangent: the forward Euler method
Here is the whole idea in one sentence: if you know where you are and you know your slope,
you know roughly where you will be a moment later. Suppose you sit at
(t_n, y_n). The differential equation hands you the slope there,
f(t_n, y_n). Pretend the solution is a straight line with that slope — its
tangent — and walk along it for one step of width h. You land at
-
The update.
t_{n+1} = t_n + h, \qquad y_{n+1} = y_n + h\,f(t_n, y_n).
-
Local error is O(h^2) — the mistake made on a
single step, from dropping every term past the linear one in the Taylor expansion of
y.
-
Global error is O(h) — the error accumulated over the
whole interval. Euler is a first-order method: halve h
and you halve the final error.
Why does the global order drop from the local O(h^2) to
O(h)? Because to cross a fixed interval you need
N = (t_{\text{end}} - t_0)/h steps — and as you shrink
h you take proportionally more of them. Each step contributes an
error of size O(h^2), and there are O(1/h) of
them, so the errors pile up to O(1/h)\cdot O(h^2) = O(h). One power of
h is eaten by the growing number of steps. This bookkeeping — local order
minus one equals global order — is the master key to every method on this page.
Let's march. The cleanest test problem in all of numerical analysis is
y' = y with y(0)=1, whose exact solution is
y(t) = e^{t} — the function that is its own slope. Take
h = 0.1 and step from t=0 to
t=1, where the truth is e \approx 2.71828:
// Forward Euler on y' = f(t,y) = y, y(0) = 1. Exact answer: y(t) = e^t.
function f(t: number, y: number): number {
return y; // the slope equals the height
}
const h: number = 0.1; // step size
let t: number = 0;
let y: number = 1; // initial value y(0) = 1
while (t < 1 - 1e-9) {
y = y + h * f(t, y); // walk one step along the tangent
t = t + h;
console.log(`t = ${t.toFixed(1)} y_euler = ${y.toFixed(6)} e^t = ${Math.exp(t).toFixed(6)}`);
}
console.log(`\nEuler y(1) = ${y.toFixed(6)}`);
console.log(`true e = ${Math.E.toFixed(6)}`);
console.log(`error = ${Math.abs(y - Math.E).toExponential(3)}`);
Euler limps in at about 2.5937 against the true
2.71828 — an error of roughly 0.125, a full
4.6\% low. Look at why: because
y' = y is curving upward, every tangent line sits below the true
curve, so each step lands a little short and the shortfalls accumulate. Euler always undershoots a
convex solution. To do better we could take a smaller h — but there is a
far cleverer route.
The chart: tangent steps chasing a curving truth
Below is exactly what Euler does to y'=y from
y(0)=1. The smooth curve is the true solution
e^{t}; the second curve is the Euler polyline with
h=0.25 — a chain of straight tangent segments. Because the truth curves
upward and each segment flies off straight, the polyline slips a little further below the true curve at
every step. That growing gap is the global error, and it is what shrinking
h (or raising the order) closes.
Sample the slope more cleverly: midpoint and Heun (RK2)
Euler's flaw is that it trusts a single slope — the one at the start of the step — for the
entire crossing. But the slope changes as you go. The fix is to sample the slope more than
once and take a smarter average, which is the whole family idea behind
Runge–Kutta methods.
The midpoint method does this in two taps. First take a tentative half-step
with Euler to peek at the slope in the middle of the interval; then throw that trial away and take the
real full step using that better, middle slope:
k_1 = f(t_n, y_n), \qquad k_2 = f\!\left(t_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2}k_1\right), \qquad y_{n+1} = y_n + h\,k_2.
Heun's method is the sibling that averages the slope at the two ends instead
(a predictor–corrector): predict with Euler, evaluate the slope there, and average it with the starting
slope. Both are second-order: local error O(h^3), global
error O(h^2). That extra power is worth a lot — now halving
h quarters the error instead of merely halving it — and it costs just
two slope evaluations per step instead of one. There is a pattern forming: pay for more slope
samples, buy more accuracy. Runge–Kutta 4 pushes it to its famous sweet spot.
The workhorse: Runge–Kutta 4
RK4 is the method you reach for by default. It samples the slope four
times across each step — once at the start, twice at the middle, once at a trial end — and blends them
with the weights 1:2:2:1, leaning hardest on the two midpoint samples:
-
The four slopes.
\begin{aligned}
k_1 &= f(t_n,\, y_n), \\
k_2 &= f\!\left(t_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2}k_1\right), \\
k_3 &= f\!\left(t_n + \tfrac{h}{2},\; y_n + \tfrac{h}{2}k_2\right), \\
k_4 &= f\!\left(t_n + h,\; y_n + h\,k_3\right).
\end{aligned}
-
The weighted update.
y_{n+1} = y_n + \frac{h}{6}\bigl(k_1 + 2k_2 + 2k_3 + k_4\bigr).
-
Local error O(h^5), global error
O(h^4) — a fourth-order method. Four slope evaluations
per step, roughly four times the work of Euler, for a wildly disproportionate accuracy gain.
Run RK4 on the very same problem y'=y, with the very same step
h=0.1, and compare to Euler's 4.6\% miss:
// Runge–Kutta 4 on y' = y, y(0) = 1, with the SAME h = 0.1 as Euler.
function f(t: number, y: number): number {
return y;
}
const h: number = 0.1;
let t: number = 0;
let y: number = 1;
while (t < 1 - 1e-9) {
const k1: number = f(t, y);
const k2: number = f(t + h / 2, y + (h / 2) * k1);
const k3: number = f(t + h / 2, y + (h / 2) * k2);
const k4: number = f(t + h, y + h * k3);
y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4); // weighted blend of the four slopes
t = t + h;
}
console.log(`RK4 y(1) = ${y.toFixed(12)}`);
console.log(`true e = ${Math.E.toFixed(12)}`);
console.log(`error = ${Math.abs(y - Math.E).toExponential(3)}`);
Euler was off by about 1.3\times10^{-1}; RK4, using the identical step size,
is off by around 2\times10^{-6} — five orders of magnitude
better for only four times the arithmetic. That astonishing bargain is exactly why RK4 has been the
default general-purpose ODE integrator since Runge and Kutta devised it around 1900, and why it still
powers everything from game physics to spacecraft trajectories.
The order made visible: halve the step, watch the error fall
The order of a method is a promise about what happens when you halve
h. A first-order method (Euler) should see its global error fall by a
factor of 2; a fourth-order method (RK4) by a factor of
2^4 = 16. Let's not take that on faith — let's measure it. This table solves
y'=y to t=1 at a sequence of ever-smaller steps
and prints the error and the ratio to the previous row:
// Global error at t = 1 for y' = y, halving h each row.
// Watch the "ratio" column: ~2 for Euler (order 1), ~16 for RK4 (order 4).
function f(t: number, y: number): number { return y; }
function euler(h: number): number {
let t = 0, y = 1;
while (t < 1 - 1e-9) { y = y + h * f(t, y); t += h; }
return Math.abs(y - Math.E);
}
function rk4(h: number): number {
let t = 0, y = 1;
while (t < 1 - 1e-9) {
const k1 = f(t, y);
const k2 = f(t + h / 2, y + (h / 2) * k1);
const k3 = f(t + h / 2, y + (h / 2) * k2);
const k4 = f(t + h, y + h * k3);
y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
t += h;
}
return Math.abs(y - Math.E);
}
const hs: number[] = [0.2, 0.1, 0.05, 0.025, 0.0125];
let prevE: number | null = null;
let prevR: number | null = null;
console.log(" h Euler err ratio RK4 err ratio");
for (const h of hs) {
const eE = euler(h);
const eR = rk4(h);
const rE = prevE === null ? NaN : prevE / eE;
const rR = prevR === null ? NaN : prevR / eR;
const f2 = (x: number) => (Number.isNaN(x) ? " - " : x.toFixed(1).padStart(5));
console.log(
`${h.toFixed(4)} ${eE.toExponential(2)} ${f2(rE)} ${eR.toExponential(2)} ${f2(rR)}`,
);
prevE = eE;
prevR = eR;
}
The ratio columns tell the whole story. Euler's error ratio hovers around
2 — every halving of h buys you one factor of
two, one bit. RK4's ratio hovers around 16 — every halving buys you a factor
of sixteen, better than a full decimal digit per halving. This is what "order" means, made
concrete: it is the exponent on the payoff you get for spending on a smaller step.
Higher order is not free, and the returns diminish sharply. Each extra order demands more slope
evaluations per step (RK4 needs 4; classical high-order schemes need many more), and a Taylor
expansion only helps if the solution is that smooth — a function with a kink or a sharp
feature gains nothing from extra derivatives it doesn't possess. In practice the smart move is not a
higher fixed order but adaptive step size: run two methods of nearby order (the
famous Runge–Kutta–Fehlberg and Dormand–Prince pairs do this), use
their difference to estimate the local error, and automatically shrink
h where the solution is wild and grow it where the solution is calm. That
is what ode45 in MATLAB and solve_ivp in SciPy actually run under the hood
— RK4-ish accuracy, with the step size steering itself.
The catch: accuracy is not stability
So far, smaller h has meant less error, full stop. But there is a second,
sneakier failure mode that has nothing to do with accuracy and everything to do with the step size
being too big for the problem. It is called stiffness, and it is where
explicit methods like Euler and RK4 can catastrophically explode.
Consider the decay y' = -\lambda y with a large rate
\lambda > 0. The true solution
y = y_0 e^{-\lambda t} races quietly to zero. But one Euler step multiplies
the value by (1 - \lambda h). If \lambda h > 2,
that factor has magnitude greater than 1, so instead of decaying the numbers
grow — and flip sign — at every step, oscillating out to infinity while the true
answer sits placidly near zero. The method is unstable.
A higher-order method is not an unconditionally stable method. Take the stiff decay
y' = -15y from y(0)=1 — the true solution
plummets harmlessly to zero. With h = 0.2 the Euler amplification factor is
1 - 15\cdot 0.2 = -2, so each step doubles the magnitude and flips
the sign. Run it and watch the numbers detonate:
// Forward Euler on the STIFF decay y' = -15 y, y(0) = 1.
// True solution y = e^(-15 t) → 0. But h is too big for stability.
function f(t: number, y: number): number {
return -15 * y;
}
const h: number = 0.2; // stability needs h < 2/15 ≈ 0.133 — this is too big!
let t: number = 0;
let y: number = 1;
for (let n = 1; n <= 8; n++) {
y = y + h * f(t, y); // multiplies y by (1 - 15h) = -2 each step
t = t + h;
console.log(`step ${n}: t = ${t.toFixed(1)} y_euler = ${y.toFixed(2)} true = ${Math.exp(-15 * t).toExponential(2)}`);
}
console.log(`\nEuler has exploded to ${y.toFixed(1)} while the truth is essentially 0.`);
The true answer never exceeds 1; Euler screams past
\pm 100 within a few steps. The cure is not a fancier explicit
method — RK4 has its own (larger, but still finite) stability limit and will blow up too if pushed.
The cure is an implicit method, next.
Backward (implicit) Euler evaluates the slope at the arrival point instead of
the departure point:
y_{n+1} = y_n + h\,f(t_{n+1},\, y_{n+1}).
Notice y_{n+1} appears on both sides — you must solve an equation
(often a nonlinear one, via Newton's
method) at every step, which is why implicit methods cost more per step. The payoff is
enormous stability. For y'=-\lambda y the backward-Euler factor is
1/(1+\lambda h), whose magnitude is below 1 for
every positive h: it is unconditionally stable
and simply cannot blow up, no matter how stiff the problem or how large the step. Stiff systems
— chemical kinetics with fast and slow reactions, electrical circuits with wildly separated time
constants — are exactly where you pay for an implicit solver and are glad you did.