Numerical Treatment of ODE

The population models built earlier each came with a closed-form solution: Malthusian growth gave a clean exponential, the logistic model an explicit S-curve, the saturation model an exponential approach to a limit. Being able to write the solution down as a formula is the comfortable case, and also the rare one. The overwhelming majority of ordinary differential equations have no closed-form solution at all — there is simply no finite combination of elementary functions that satisfies them. The moment a model grows past a textbook example, this is what happens, and the only way forward is to compute an approximate solution numerically, step by step. This page is about how that is done, and about the ways it can go wrong.

Initial value problems and boundary value problems

A differential equation on its own does not pin down a single solution — it describes a whole family of curves, one for each way of fixing the constants of integration. To single out the one solution we actually want, the equation has to be paired with side conditions. The population-dynamics models are exactly this: one or several ODEs together with enough extra conditions to make the solution unique. There are two standard ways to supply those conditions, and they split ODE problems into two families.

The first fixes the solution at the start of the time interval — we are told the state at the beginning and asked how the system develops from there.

This is the population-dynamics setting from before. The initial value $y_a$ plays exactly the role the starting population $p_0$ did — the headcount we measure at the outset and then watch evolve. So $y(a) = y_a$ is not just a piece of notation; it is the concrete starting datum, the analog of $p_0$, that turns the whole family of solution curves into the single one we care about.

The second family fixes the solution at both ends of the interval — we are told where the trajectory starts and where it must end, and asked to find the path between.

The optimal trajectory of a space shuttle is the standard example: the craft leaves a known launch point and must arrive at a prescribed target, and the problem is to find the flight path joining the two. Most of what follows concerns initial value problems, which are the natural shape for population dynamics; boundary value problems are taken up later, once the machinery for IVPs is in place.

Reading the prototype

A few features of the prototype $\dot y(t) = f(t, y(t))$ are worth drawing out, because the notation hides some choices.

The right-hand side $f$ is not the $f$ of the population models. There, $f(p, q)$ named the per-capita growth rate — the factor multiplying the population in $\dot p = f(p, q)\, p$. Here $f$ is the entire right-hand side of the ODE: the whole expression for $\dot y$, the rate of change itself. Same letter, a different job.

The unknown $y$ may be a single scalar function of time or a vector of several functions stacked together. The system form above is the vector case written out component by component: $n$ populations evolving at once, where — as is typical in these models — each species’ rate of change depends on all the others. The coupled predator–prey and competition systems are of exactly this shape, $n$ first-order equations sharing one set of unknowns.

And every equation here is first-order: only the first derivative $\dot y$ appears, never $\ddot y$ or anything higher. That is a genuine simplification, and it is not as restrictive as it looks. A higher-order equation — the damped oscillator’s $\ddot p$, for instance — can be folded back into this first-order form by introducing the lower derivatives as extra unknowns, a reduction taken up later in this chapter. So the first-order IVP is not a special case; it is the standard form everything else is converted into.

One last case is worth flagging. If $f$ happens to depend only on $t$ and not on $y$ — that is, $\dot y = f(t)$ — then the ODE is no longer really a differential equation to be solved but a plain integration, $y(t) = y_a + \int_a^t f(s)\,ds$, with nothing to do beyond evaluating an integral. The interesting case, and the one that forces the numerical machinery of this page, is the usual one: $f$ depends on the unknown $y$ itself, tying the rate of change at each instant to the current state, so the solution cannot be had by integrating a known function of $t$ alone.

A short excursus: discretization

Before any method can run, there is a gap to bridge. A model lives in the continuum — real-valued quantities, smooth functions, an unbroken interval of time — but a computer is a finite machine that can hold only finitely many numbers and perform only finitely many operations. Every numerical method therefore begins with the same move.

This happens at several levels, and it is worth seeing each one, because the approximations they introduce are the ultimate source of the errors the rest of this page has to control.

Representing numbers

The first thing to discretize is the real number line itself. A machine cannot store an arbitrary real number — between any two reals lie uncountably many others, and only finitely many fit in memory. There are three standard schemes, differing in how they trade off range (how large and small the representable numbers go) against resolution (how finely spaced they are):

• Integer arithmetic keeps only whole numbers, so the spacing is always exactly one. Nothing between two integers can be stored: a value like $2.5$ has no slot, and a division such as $5 / 2$ silently drops its remainder. For continuous quantities that is hopelessly coarse.
• Fixed-point arithmetic narrows the spacing to let fractions in. It pins the radix point at a fixed position and keeps a fixed number of digits after it, so the machine really stores an integer such as $123$ but reads it as $1.23$ (the point is just understood to sit two places in). Now fractions fit — but the spacing it buys is uniform across the whole range: choose a step of $0.01$ and every representable number sits $0.01$ from its neighbor, near zero and near a million alike. That is the wrong place to spend precision. Out among large values the fixed step is far finer than anyone needs, while near zero that same step is too coarse to resolve anything smaller than itself. What you actually want is fine steps for small numbers and coarse steps for large ones, and a single absolute spacing can never be both.
• Floating-point arithmetic fixes that imbalance by letting the radix point “float” to wherever the digits are needed, so both the range and the resolution vary with the number’s size. The spacing now scales with magnitude — tight near zero, wide out among large values — but the relative precision stays the same everywhere: the same handful of significant digits whether the quantity is tiny or huge. This is what makes it possible to hold both very large and very small numbers sensibly, and it is what essentially all scientific computation uses.

The floating-point idea has a precise definition.

Resolution — accuracy is relative

The single most important property of a floating-point system is how finely it resolves numbers, and the key fact is that it does so relatively, not absolutely.

Here is where it comes from. Take two neighboring machine numbers with the same exponent $E$. Their mantissas differ by one, so the numbers themselves are $M \cdot B^E$ and $(M+1) \cdot B^E$, and subtracting leaves an absolute gap of just $B^E$ between them. To turn that into a relative gap, divide it by the value itself, $M \cdot B^E$. This fraction is largest where the value is smallest, and the smallest a normalized mantissa is allowed to be is $|M| = B^{t-1}$ (the lower bound built into the definition). Putting that worst case in, the $B^E$ cancels top and bottom:

So the relative gap never exceeds $\rho$. The absolute spacing $B^E$ itself grows as the numbers grow, so neighbors far out on the line are spaced far apart while neighbors near zero are packed tightly together, but the relative spacing stays capped at $\rho$ everywhere. This is exactly the behavior one wants: precision is proportional to magnitude, the same handful of significant digits whether the quantity is tiny or huge.

Beyond numbers

Discretization does not stop at the numbers. The same move — replace an infinite or continuous object with a finite one — is applied to functions and to the operations of calculus:

• An infinite series is cut off after finitely many terms, turning it into a polynomial. The sine, for example, is computed in practice not from its defining infinite series but from a truncated polynomial that agrees with it closely on the range of interest.
• A continuous interval of time is replaced by a finite grid of discrete points $t_0, t_1, t_2, \dots$, and the solution is sought only at those points rather than at every instant.
• A derivative — defined as a limit, an inherently continuous notion — is replaced by a difference quotient taken over a small but finite step.

The derivative is the limit of this ratio as $h \to 0$; keeping $h$ small but positive is what makes it computable. This last substitution is the seed of the whole numerical treatment of ODEs — it turns $\dot y = f(t, y)$ from a statement about instantaneous rates into a rule that steps the solution forward one finite $h$ at a time. But every one of these discretizations — rounded numbers, truncated series, finite steps in place of limits — buys that computability at the price of some error. Keeping those errors identified and under control is what the next two sections are about, and what the rest of the page ultimately depends on.

Rounding and round-off error

A numerical algorithm picks up error from several sources. The most basic one is built into floating-point arithmetic itself: a real number that is not a machine number cannot be stored exactly and must be replaced by one that can.

Every real $x$ that is not itself a machine number falls strictly between two that are — its nearest representable neighbors below and above,

If $x$ happens to be a machine number, both neighbors coincide with $x$ itself; otherwise $x$ sits between $f_l(x)$ and $f_r(x)$, and storing it means choosing one of the two. That choice is made by a rounding map.

Which of the two neighbors $\mathrm{rd}(x)$ actually returns is fixed by the rounding mode:

• rounding up always takes the upper neighbor, $\mathrm{rd}(x) = f_r(x)$, i.e. it rounds toward $+\infty$;
• rounding down always takes the lower neighbor, $\mathrm{rd}(x) = f_l(x)$, i.e. it rounds toward $-\infty$;
• correct rounding (round to nearest) takes whichever neighbor is closer, with a fixed tie-breaking rule for the exact midpoint between them;
• truncating (round toward zero) simply cuts the number off, keeping the leading digits and throwing away the rest. Since dropping digits only ever shrinks the magnitude, this lands on the neighbor nearer to zero. That matches rounding down for positive values, but for negative ones it rounds the opposite way: $-2.7$ truncates to $-2$, whereas rounding down would give $-3$.

Whatever the mode, the price of storing $x$ is the gap between it and the machine number chosen for it.

Arithmetic on machine numbers

Rounding is not a one-time cost paid only on input. Even when both operands are machine numbers, the result of an arithmetic operation usually is not — multiply two three-digit numbers and the exact product can need six digits — so it has to be rounded back into $\mathbb{F}$ as well. The idealized model of this is clean:

This is the behavior the IEEE floating-point standard prescribes — though, as a technological matter, not every computer implements it faithfully. For error analysis it is handier to weaken the exact equality into a bound: every machine operation returns the true result perturbed by a small relative amount,

The perturbation $\varepsilon(a, b)$ depends on the operands, but its size is always controlled. The $O(\rho)$ here is big-O notation, read “on the order of $\rho$”: it means the bound $\tilde{\varepsilon}$ is at most some fixed constant times the resolution $\rho$, so it stays proportional to $\rho$ and shrinks right along with it as the precision improves. A single rounding, then, is utterly harmless. The real question is what happens when millions of them compound over a long computation: do the individual errors stay independent and roughly cancel, or do they reinforce and grow? Estimating that accumulated influence is the job of round-off error analysis, and the standard worked example is Horner’s method for evaluating a polynomial — the classic case for tracing how round-off builds up across a long chain of multiply–add steps.

Further sources of errors

Round-off is the most basic source of error, but it is not the only one, and on modern hardware it is rarely the worst. Double-precision arithmetic carries so many significant digits that the resolution $\rho$ is minuscule, so the errors that actually dominate a real computation usually come not from the number format but from the deliberate approximations a method makes. Three of these are worth naming.

This is the direct cost of the discretization that every method on this page relies on, and it is the one these methods spend most of their effort controlling. Once a concrete method is in hand, the discretization error splits into a local part (made in a single step) and a global part (accumulated across all steps), and the central questions become how fast each shrinks as the step size goes to zero.

This takes a few common shapes. A series computed term by term has to be cut off after $N$ terms, leaving out the infinite tail — the same series-to-polynomial truncation seen earlier, now viewed as a source of error. Root-finding by Newton’s method likewise produces only a sequence of ever-better approximations to a true root; the iteration has to be stopped somewhere, either after a fixed number $N$ of steps or once the change from one step to the next becomes insignificant, and whatever distance remains to the exact root is truncation error.

Data error is different in kind from the other two — it has nothing to do with the algorithm or the machine, but is baked into the problem before any computation starts. It cannot be reduced by a better method or finer precision, only carried along, and how much it grows or shrinks as it passes through a computation turns out to be a property of the problem itself, taken up later on this page.

Round-off, discretization, truncation, data — each is a separate way the computed answer can drift from the true one, and a sound numerical method has to keep an eye on all of them at once.

One-step methods

With the machinery in place, we can build the first actual solvers. They all share one shape: start from the known initial value and march forward, producing the solution at one grid point after another. The simplest family computes each new point using only the single point just before it.

The whole design problem is bottled up in $\Phi$: every one-step method on this page — Euler, Heun, Runge–Kutta — keeps the skeleton new value = old value + step size $\times$ estimated growth and differs only in how it estimates that growth. Keep that template in mind; it is the single thread running through everything below. (Methods that do reuse older points $y_{k-1}, y_{k-2}, \dots$ — the multistep methods — come later in this chapter.)

Two symbols that look alike: $y(t)$ versus $y_k$

One notational point underlies everything that follows. The function $y(t)$, with its time argument in parentheses, is always the exact solution — the true trajectory the ODE and initial value pin down uniquely. We never actually have it, except at the starting point $y(a) = y_a$; in the interesting cases there is no formula to evaluate it with. The subscripted $y_k$, with an index and no parentheses, is always the computed approximation at the grid point $t_k$, and $y_k \approx y(t_k)$ — close, but in general not equal.

Each approximation is built on the previous one. Only the first value is exact, $y_0 = y_a = y(a)$; every step after that feeds the method its own previous output, so $y_{k+1}$ is built from the approximation $y_k$, not from the true $y(t_k)$. That is how error compounds. The exact $y(t)$ never enters the computation — it appears only in the analysis, as the yardstick the computed $y_k$ is judged against.

The Euler method

The starting point is the difference quotient from the discretization excursus — replacing a derivative by a slope over a small finite step is the finite-differences approximation. It says the derivative over that step is approximately the slope of the chord,

The prototype IVP $\dot y(t) = f(t, y(t))$ says that this same derivative equals $f(t, y(t))$. Setting the two equal and solving for the value one step ahead turns a statement about instantaneous rates into a recipe for moving forward:

Read it in the running template: the population a short time $\delta t$ from now equals the population now, $y(a)$, plus the time elapsed times the growth rate $f(a, y(a))$ estimated at the current instant. Here $a$ is the start of the time interval and $y(a)$ the known initial value — the headcount we measure at the outset. We do not know the true rate over the whole step, so we use the only rate we have: the one at the point we are standing on.

Nothing forces us to stop after one step. Having landed on an estimate at $a + \delta t$, we treat it as the new “current” point and repeat, and repeat again. That iteration is the Euler method.

This is the move from continuum to discrete made concrete: one continuous equation, valid at every instant, becomes a sequence of discrete equations, one per grid point. The grid covers a finite range $t \in [a, b]$, the spacing between neighboring points is the step size $\delta t$, and the method hands back the solution only at $t_0, t_1, t_2, \dots$ — never in between.

That the Euler method runs — that the recurrence mechanically produces a value at every grid point — is clear. Whether what it produces is any good is not, and three questions decide it. They organize the rest of this section.

• Does it converge? As the step size shrinks, $\delta t \to 0$, do the computed points close in on the true solution of the ODE? A method that refuses to improve no matter how fine the grid is worthless; this is the bare minimum we demand, and the next section makes it precise as consistency and convergence.
• How fast does it converge? Suppose it does converge. Halving the step size doubles the work — does it buy enough accuracy to be worth it? If halving $\delta t$ only halves the error, progress is painfully slow; if it cuts the error to a quarter or a sixteenth, refinement pays off fast. This is the method’s order, and it is the single number that separates Euler from Heun from Runge–Kutta.
• Is there a largest safe step? Going finer is accurate but expensive; going coarser is cheap but risky. Is there a sharp threshold — a maximal $\delta t$ below which the method behaves and above which it falls apart? The tension is real: forecasting two seconds ahead in one step is fine, leaping two months ahead in one step is meaningless, and somewhere between lies a limit. Pinning that limit down is the subject of stability and stiffness, later on this page.

The first two questions are answered next. The third has to wait — it turns out to depend on the method and the equation together, in ways that take real work to untangle.

Discretization error: local, global, and what “convergence” means

There is a second route to the Euler formula that comes not from the difference quotient but from calculus, and it has the advantage of exposing the error directly. Take the exact solution $y(t)$ and expand it one step ahead with a Taylor series:

The remainder $R$ collects every higher-order term — the ones carrying $\ddot y$ and beyond. Over a short step those terms are small, so we drop them; and we replace $\dot y(t_k)$ with $f(t_k, y_k)$, since the ODE says the rate equals $f$ and $y_k$ is the value we actually have at that step. What is left is exactly the Euler step. The same method, reached two ways: the difference quotient is the numerical route, the truncated Taylor series the analytic one. The Taylor route makes plain what we threw away — the remainder $R$ — and that discarded remainder is the error. Measuring it is the whole game.

That error comes in two flavors, and telling them apart is the key conceptual step.

Reading the formula: the fraction is the scheme’s difference quotient — the slope Euler would use if it sat exactly on the true curve at time $t$ — and $f(t, y(t))$ is the rate the ODE actually prescribes there. Their difference is how far the finite-step slope deviates from the true derivative, sampled on the exact solution. The maximum over $[a, b]$ then takes the largest such deviation anywhere on the interval, so $l(\delta t)$ is a single worst-case figure for the whole run.

The distinction is the difference between the error this one step introduces (local, measured as if the past were perfect) and the error you are actually stuck with (global, the real distance from truth at $t_k$). One would expect the global error to be vastly larger than the local one — after all, it is a running total of many local mistakes, so naively it should pile up step after step without bound. It does not, and the reason is the whole reason numerical methods are usable at all: the per-step errors do not march in lockstep. Some push the estimate above the true curve, some below; they partly cancel rather than purely accumulate. If they really did add up relentlessly, the global error would explode and no amount of computing would give a trustworthy answer. (Exactly when this benign cancellation fails — when errors reinforce instead — is the stability question flagged above.)

Both definitions need the exact solution $y(t)$, which in general we do not have. So to actually measure the error, we supply a stand-in for $y(t)$: either rerun the same problem on a much finer grid and treat that high-resolution result as the “truth,” or test the method on one of the rare problems whose exact solution is known in closed form (the population models from earlier are exactly such cases). Comparing the computed $y_k$ against that reference is how $l(\delta t)$ and $e(\delta t)$ get evaluated in practice.

Now the two minimum demands on a method can be stated cleanly.

Consistency is the weak condition — it only controls one step at a time, and says nothing about what happens once many steps are chained together. It guarantees that each individual step, taken in isolation, can be made as accurate as we like by shrinking $\delta t$. But many small errors can still combine into a large one, so “every step is nearly right” does not by itself force “the final answer is nearly right.” For that we need the accumulated error to vanish — the stronger property, convergence.

Both properties say the error eventually disappears, but neither says how fast. That rate is the practically decisive quantity, and it has a name.

The $O$ here is ordinary big-O notation — the same one met earlier for round-off, read “on the order of.” Writing $l(\delta t) = O(\delta t^k)$ just says the error stays below some fixed constant times $\delta t^k$ once $\delta t$ is small. The two senses of order then line up exactly: the order of a method is the exponent $k$ sitting inside that big-O bound, so it is the same notion as big-O’s “order of magnitude,” not a clashing second use of the word.

Comparing the two errors, step by step

A picture makes the local-versus-global distinction concrete. Plot time along the horizontal axis at the grid points $t_0 = a, t_1, t_2, \dots$ and the solution value $y(t)$ along the vertical. The red curve is the true solution of the original IVP, the one passing through the initial value $y_0 = y_a$. The Euler points $y_1, y_2, y_3, y_4$ are obtained by the calculation rule, each one a tangent-line step: from the current point, follow the local slope straight for one step of width $\delta t$.

The other colored curves are the subtle part, and the source of the “where did those come from?” confusion. Each one is the exact solution of the same ODE, but launched from a computed point rather than from $y_a$. The green curve is the true trajectory that would unfold if the system were restarted exactly at $(t_1, y_1)$; the blue one is the true trajectory restarted at $(t_2, y_2)$; the black one restarted at $(t_3, y_3)$. The label "$y(t),\ y(t_1) = y_1$" means precisely “the solution curve $y(t)$ satisfying $y(t_1) = y_1$.” They are not extra numerical outputs — they are the exact-solution family, one member through each point the method actually produced.

With that, the two errors are read straight off the figure:

• Local discretization error (the short bars): at each step, the vertical gap between the Euler point $y_{k+1}$ and the exact curve that passes through the previous point $y_k$. It is “local” because it measures one Euler step against the true trajectory launched from where that step began — exactly the “assume the past was perfect” definition above.
• Global discretization error (the tall bars): the vertical gap between the Euler point $y_k$ and the red true solution. It is the error you are really carrying — the full distance from the original trajectory, every inherited mistake included.

A last caution about the figure. Drawn this way, the tall global-error bars seem to grow without limit, step after step. That impression is an artifact of one particular example. In general the global error can grow, shrink, or oscillate as you march along — the partial cancellation of local errors discussed above means it need not balloon. If it genuinely grew without bound, numerical integration would be pointless; that it usually does not is what makes the whole enterprise work.

Order, and getting more of it: the Heun method

For the Euler method, as long as the second derivative $\ddot y$ and the partial derivative $f_y = \partial f / \partial y$ (how the rate $f$ responds to a change in the state $y$) stay bounded — never running off to infinity anywhere on $[a, b]$ — the analysis above pins down the order exactly: it is consistent of first order, $l(\delta t) = O(\delta t)$, and convergent of first order, $e(\delta t) = O(\delta t)$. In the language just built, $O(\delta t)$ means halving the step size only halves the error — honest, but slow. And the warning attached to the convergence definition bites here too: there exist methods that are perfectly consistent yet do not converge, so first-order consistency is reassuring but not the end of the story.

Slow convergence is the practical problem. To reach a target accuracy with Euler you may need a punishingly small step, and small steps mean many of them — high computational cost. The cure is a higher order: a method with $e(\delta t) = O(\delta t^2)$ reaches the same accuracy with far larger steps, because halving now quarters the error. So we want order, and the question is how to manufacture it.

One route is staring at us. The Taylor derivation of Euler dropped the remainder $R$ after the first-derivative term; keep the next term instead and you get a higher-order method directly. The trouble is that the next term carries $\ddot y$, which by the chain rule unpacks into derivatives of $f$ itself — and the term after that into still more. The formulas grow ugly fast, and they demand that you supply derivatives of $f$ by hand. For a real model that is often impractical.

The alternative sidesteps derivatives entirely: instead of evaluating higher derivatives of $f$ at one point, evaluate $f$ itself at several points within the step and combine the results. Methods built this way are the Runge–Kutta methods, and the simplest one above Euler is the method of Heun.

The idea follows the running template — new value = old value + step size $\times$ estimated growth — and improves only the growth estimate. Euler uses a single slope, $f(t_k, y_k)$, sampled at the start of the step. But the true rate changes across the step, and trusting the starting slope alone systematically lags. The natural fix is to average the slope at the start with the slope at the end of the step — an arithmetic mean of the two derivatives, which tracks the curve far better than either endpoint alone.

There is one obstruction: the slope at the end is $f(t_{k+1}, y_{k+1})$, and $y_{k+1}$ is the very thing we are trying to compute — we do not know it yet. Heun’s resolution is to predict the endpoint with a cheap Euler step, then use that prediction to evaluate the ending slope:

The payoff is concrete: where Euler’s error halves when you halve the step, Heun’s quarters. The cost is two evaluations of $f$ per step instead of one — a trade almost always worth making.

Runge–Kutta methods

Heun is the first rung of a ladder. Push the same idea further — sample $f$ at more points inside the step and take a cleverly weighted average — and you climb to higher and higher order. The whole family is the Runge–Kutta methods, and its most famous member, often meant when someone says simply “the Runge–Kutta method,” uses four slope samples.

The four samples read as a sequence of ever-better guesses at the slope across the step: $T_1$ is the slope at the start (plain Euler); $T_2$ is the slope at the midpoint, reached by stepping there with $T_1$; $T_3$ is the midpoint slope again, but reached using the improved $T_2$; and $T_4$ is the slope at the endpoint, reached using $T_3$. The weighting $1 : 2 : 2 : 1$ leans on the two midpoint estimates, which carry the most information about the step’s average behavior. Fourth order is the headline: RK4 reaches accuracies with modest steps that Euler could only dream of — which is why it is the workhorse default for non-stiff problems (what makes a problem stiff is taken up later on this page, alongside stability). The price is four evaluations of $f$ per step, and the same template still holds: new value = old value + step size $\times$ (weighted) estimated growth.

Multistep methods

The Runge–Kutta ladder buys order by sampling $f$ at more points inside each step — RK4 calls $f$ four times per step. That accounting hides where the real cost lives. The arithmetic that combines the samples $T_1, \dots, T_4$ is trivial; the expense is the evaluations of $f$ themselves. When $f$ is a tidy formula — a polynomial, a sine — an evaluation is nearly free and four of them cost nothing to worry about. But $f$ is the entire right-hand side of the model, and for a serious model it can be enormous: if a single evaluation of $f$ is itself a multi-hour supercomputer run (a weather forecast, say), then RK4’s four-evaluations-per-step is not the difference between 10 and 40 milliseconds but between six hours and a full day. At that scale the natural question is whether order can be bought without paying for extra evaluations of $f$.

It can — by looking the other way in time. Runge–Kutta manufactures its extra slope samples by probing forward into the step, and every probe is a fresh evaluation. But the method has already evaluated $f$ at the previous grid points $t_{k-1}, t_{k-2}, \dots$ on earlier steps, and those numbers are just sitting there. Reusing that history instead of computing new samples buys higher order essentially for free.

Before writing one down, a shorthand the rest of this page leans on: write

for the slope the method computes at grid point $t_k$. Once $y_k$ is known, $f_k$ is one evaluation of $f$, and a method that reuses history is one that keeps the old $f_k$ values around rather than recomputing them.

The prominent representative is the second-order Adams–Bashforth method, which reaches one step into the past:

It still fits the running template — new value = old value + step size $\times$ estimated growth — and differs only in the growth estimate: it extrapolates from the last two slopes $f_k$ and $f_{k-1}$ instead of probing new ones. The only fresh slope a step needs is $f_k$ — the other one, $f_{k-1}$, is reused from the previous step — so a whole step costs one evaluation of $f$, the same as Euler — yet the method is second order, the same as Heun, which costs two.

Where does the $3 f_k - f_{k-1}$ come from? The integral form from the remark above is the key. Producing $y_{k+1}$ means approximating $\int_{t_k}^{t_{k+1}} f(t, y(t))\, dt$, and the multistep idea is to replace the integrand $f$ by a polynomial that is cheap to integrate exactly.

The polynomial is fixed by data we already have, and integrating a polynomial is exact and instant, so the only cost is the interpolation itself — no new evaluation of $f$. The number $p$ of points sets the order:

• $p = 1$: $P$ is a constant through the single point $(t_k, f_k)$, and integrating it recovers $y_{k+1} = y_k + \delta t\, f_k$ — the Euler method.
• $p = 2$: $P$ is the line through $(t_{k-1}, f_{k-1})$ and $(t_k, f_k)$, and integrating it over the next step gives exactly the $3 f_k - f_{k-1}$ rule above.
• general $p$: order $p$, at the fixed price of one evaluation of $f$ per step.

That is the payoff Runge–Kutta could not offer: more order without more evaluations.

Starting up: the first few steps

The cheapness comes with a debt at the very beginning. A $p$-step method needs $p$ past points before it can take a step, but at the start of the integration only the initial value $y_0$ is on hand — there is no history yet to reuse.

The fix is to bootstrap. Run a one-step method (or a smaller multistep method) for the first few steps, just long enough to generate the missing $y_1, \dots, y_{p-1}$ and their slopes, and then switch to the full multistep method once enough history has accumulated. With a sufficiently accurate starter this preserves the method’s order: the startup phase is a handful of steps and does not spoil the asymptotics.

There is also a disadvantage shared by every method so far, explicit one-step and multistep alike: the step size $\delta t$ often has to be kept very small — not for accuracy, but to stop the computation from going unstable. Push $\delta t$ too large and the computed solution sprouts spurious oscillations that bear no resemblance to the true trajectory. Why this happens, and why it is a property of the method rather than the problem, is the subject of the stability discussion later on this page; for now the consequence is what matters. Small steps mean many steps, and many steps mean high computational cost. The remedy — methods that stay well-behaved at large step sizes — is implicit methods.

Implicit methods

Every method up to here shares a structural feature that is easy to overlook because it is so convenient: each computes $y_{k+1}$ from a formula whose right-hand side contains only things already known. Plug in $y_k$ (and any history) and read off $y_{k+1}$ directly.

The alternative is to let $y_{k+1}$ appear on the right-hand side as well — to use the not-yet-known new slope $f_{k+1} = f(t_{k+1}, y_{k+1})$ in the very formula meant to produce $y_{k+1}$.

Applying this to the interpolation idea — interpolating through the previous points and the new grid point $t_{k+1}$ — gives the implicit multistep family.

The first-order member is important enough to have its own name.

The second-order Adams–Moulton member, averaging the start and end slopes, is the implicit trapezoidal method — the same trapezoidal average the Heun method used, but with the exact end slope $f_{k+1}$ rather than a predicted one.

That “solving an equation” is the catch. Since $f$ is generally nonlinear, $y_{k+1}$ is the unknown of a nonlinear equation, and there are two ways to extract it:

• Brute force: solve the nonlinear equation directly each step, e.g. by a Newton iteration for the root $y_{k+1}$. Accurate, but a full root-find per step is costly.
• Predictor–corrector: the easier and far more common route, taken up next.

The trade-off is the whole reason to go implicit. A single implicit step costs much more than an explicit one — you are solving an equation, not just “stepping further.” But, as the stability and stiffness sections will show, implicit methods stay well-behaved at far larger step sizes (sometimes a hundred times larger) with no loss of accuracy, so they need far fewer steps. Whether that is a net win depends on the problem: many cheap steps versus a few expensive ones. There is no universal rule, but for the hard cases coming up — stiff problems — implicit methods turn from merely attractive into necessary.

Predictor–corrector

The cheap way around the implicit equation is almost obvious once stated. If the obstruction is that we do not yet know $y_{k+1}$ to put into the right-hand side, then estimate it first with an explicit method and use that estimate wherever the true $y_{k+1}$ was needed. An implicit step becomes two explicit ones — approximating the approximation.

This is not a new trick — it is exactly what the Heun method already did: its predictor is a plain Euler step and its corrector is the implicit trapezoidal rule evaluated at that prediction,

Seen this way, Heun is the predictor–corrector pairing of explicit Euler with the second-order Adams–Moulton rule.

The accounting comes out well. A predictor–corrector step costs about twice an explicit step — one evaluation for the predictor, one for the corrector — instead of the open-ended cost of a full Newton solve, while keeping the large-step-size advantage of the implicit method it imitates. So the assessment of implicit methods, fairly stated: a single step is more expensive (solving an equation, or running two explicit stages, rather than just “going further”), but the number of steps is generally much smaller because larger step sizes are admissible, and the total cost can come out lower.

Higher-order equations and systems

Two loose ends remain before leaving the methods behind, and both reduce to one move: everything so far was written for a scalar first-order equation, and both the vector case and the higher-order case fold back into it.

Systems first, because they need almost nothing new. The population models with several interacting species are systems of first-order ODEs, and the entire machinery above — Euler, Heun, RK4, multistep, implicit — carries over verbatim. The recipe is identical; only the objects change. Wherever a scalar $y_k$ appeared, read a vector $\mathbf{y}_k$; wherever the scalar slope appeared, read the vector-valued right-hand side. The step formulas stay the same symbol for symbol, and everything provable for one ODE holds for a system by reading the scalars as vectors.

Higher-order equations take one extra step. The population models only ever involved the first derivative, but mechanics routinely produces the second — acceleration, hence force — as in the damped oscillator $\ddot p + \mu \dot p + \omega^2 (p - p_\infty) = 0$ derived among the population models. The prototype on this page is first-order by design, and the first-order paragraph earlier promised that higher-order equations fold back into it. Here is how.

The reduction trades order for size: an $n$-th order scalar equation becomes a first-order system of $n$ equations. Nothing is lost and no new solver is needed — the order is gone, absorbed into a slightly larger system, which the methods above already handle. This is why the first-order IVP is not a narrow special case but the standard form every ODE is converted into.

What makes an ODE hard to solve

Convergence and order describe a method working as intended: refine the grid and the answer improves, at a predictable rate. But several things can still derail a real computation, and they are worth separating because they have different owners. Some trouble belongs to the problem — baked in before any method is chosen — and some belongs to the method. The rest of this page works through three such phenomena: ill-conditioning, instability, and stiffness. Keeping straight which is a property of the problem and which of the method is the whole point, and a summary table at the end lays them side by side.

Ill-conditioning

The first phenomenon has nothing to do with the algorithm at all. It is a property of the problem.

Ill-conditioning is the close quantitative cousin of an ill-posed problem: an ill-posed problem fails outright to depend continuously on its data, while an ill-conditioned one depends on them continuously but so steeply that, in finite precision, it behaves almost as badly. Such problems are very hard to treat numerically, because exact input is almost never available — and an example makes the danger concrete.

Now blur the input by the tiniest amount. Suppose the first initial condition is off by an infinitesimal $\varepsilon$ — say $\varepsilon = 10^{-100}$, a perturbation so small the problem is, for every practical purpose, unchanged:

Solving again gives

At first glance this looks reassuringly close to the original. The first term is $(1 + \text{tiny})\, e^{-t} \approx e^{-t}$, and the second carries the microscopic coefficient $\varepsilon/(N + 2)$. Surely negligible.

It is not. The second term hides a positive exponential, $e^{(N + 1)t}$, and in the long run that changes everything. No matter how minuscule its coefficient, an exponentially growing factor eventually overwhelms an exponentially decaying one: there is always a time beyond which $\frac{\varepsilon}{N + 2} e^{(N + 1)t}$ dwarfs the $e^{-t}$ term, however small $\varepsilon$ is (and the larger $N$ is, the faster that runaway term grows). So the two solutions do not merely differ slightly — they part ways completely:

The true solution decays to zero; the infinitesimally perturbed one blows up to infinity.

This is ill-conditioning in its purest form, and it is exactly why numerical methods have to worry. The exact input is rarely available — initial data come from measurements (the data error), and intermediate results pick up round-off along the way — so some perturbation is unavoidable. For a well-conditioned problem that is harmless; for an ill-conditioned one it is fatal, and no choice of method or precision can rescue it, because the trouble lives in the problem, not the algorithm. It is also why careful people working from slightly different data can reach completely opposite conclusions — and why ill-conditioning is such a convenient excuse: when nothing can be done, nothing is anyone’s fault.

Stability

Ill-conditioning lets the problem off the hook for the method. The next phenomenon is the reverse: a well-conditioned problem, no excuses available, that a method still gets badly wrong. This is where the gap left open earlier finally gets filled.

Recall the warning attached to convergence: a method can be consistent — each isolated step accurate in the limit — without being convergent, because consistency says nothing about whether inherited errors stay under control as steps are chained. The property that controls them, and so bridges consistency and convergence, is stability. The sharpest way to see why it is needed is a method that is consistent and yet not convergent.

To pin the blame on the method, take a problem that is impeccably well-conditioned:

whose exact solution

starts at $1$ and decays smoothly to the limit $\tfrac12$ as $t \to \infty$, the $e^{-2t}$ term dying away. The problem is well-conditioned: perturbing the start to $y_\varepsilon(0) = 1 + \varepsilon$ shifts the solution by only $|y_\varepsilon(t) - y(t)| = \varepsilon\, e^{-2t}$, an error that itself shrinks with time. Small input error, small and shrinking output error — the good case.

Now solve it with the midpoint rule.

Substituting $f(t_k, y_k) = -2 y_k + 1$ turns it into the explicit recurrence

started from $y_0$ and the exact second value $y(\delta t)$. Mechanically it runs without complaint. But watch what it produces deep into the integration:

The true solution sits placidly at $\tfrac12$. The computed one swings between large negative and large positive values, growing in amplitude — and, crucially, refining the step size does not cure it. Shrinking $\delta t$ from $1.0$ to $0.01$ only postpones the blow-up; the oscillation still takes over and the values still explode. At every step size, eventually, it goes wrong.

The picture makes the failure unmistakable:

The exact solution is the flat curve resting at $\tfrac12$. The midpoint approximations start out tracking it, then each peels away into growing oscillations — the larger the step size, the sooner the departure. Even the finest step ($\delta t = 0.01$) only holds the line longer before a swelling oscillation envelope takes over. A method for which this kind of blow-up simply cannot happen is called stable; one for which it can is unstable. And the lesson of the figure is that instability is not a flaw in the model or the data — the problem here is flawless — but purely a matter of having chosen the wrong algorithm.

The midpoint rule is the cautionary case laid bare. It is consistent — its symmetric difference quotient genuinely approximates the first derivative as $\delta t \to 0$ — yet it is not convergent, as the blow-up showed, and on an interval unbounded to the right the failure persists no matter how small $\delta t$ is. Consistency without stability is not enough.

What is enough is the two together. This is the relation that closes the gap consistency left open:

Stability usually comes at a price — frequently a condition that the step size be small enough, sometimes very small (this is the source of the small-step requirement flagged back at the multistep methods). Which methods are stable is partly settled wholesale:

• all explicit one-step methods (Euler, Heun, RK4) are stable;
• the midpoint rule is not;
• the Adams–Bashforth and Adams–Moulton multistep families are stable as $s$-step methods for $s > 1$.

Stiffness

The third phenomenon is the subtlest, because it strikes a method that has every good property and still fails. Consider another well-conditioned problem:

with exact solution

This solution is almost boring. The $e^{-1000t}$ term collapses to nothing almost instantly — by $t = 0.01$ it is already about $e^{-10} \approx 0.00005$ — so the curve drops from $2$ to $1$ in a razor-thin initial layer and then sits flat at $1$ for the entire rest of the domain. Practically a horizontal line.

Apply the explicit Euler method:

The behavior is governed entirely by the factor $q = 1 - 1000\delta t$ raised to the $(k+1)$-th power. If $|q| < 1$ the term decays and $y_k \to 1$ as it should (at $\delta t = 0.002$ exactly, $q = -1$ and the term oscillates forever without decaying); but $|q| > 1$ — which happens as soon as $\delta t > 0.002$ — makes $q^{k+1}$ grow in magnitude and flip sign each step, so the computed solution oscillates and diverges, instead of resting at the flat line $1$.

Here is the unsettling part. Explicit Euler is consistent; it is stable, like every explicit one-step method; and consistency plus stability give convergence. It has the full set of guarantees — and it diverges anyway.

That is exactly what happens here: the negligible, short-lived term $e^{-1000t}$ — utterly insignificant after the first sliver of the interval — dictates a step size $\delta t < 0.002$ over the whole domain, just to keep the method from exploding. To reproduce what is essentially the constant $1$, the method is forced into thousands of tiny steps. Enormous effort, trivial answer.

The resolution of the paradox is that consistency, stability, and convergence are asymptotic properties: they all carry the silent qualifier “for sufficiently small $\delta t$.” Convergence is a statement in big-O form — like a complexity bound, it describes the limit $\delta t \to 0$ and promises nothing for a finite $\delta t$ above the threshold. None of the guarantees is violated by the blow-up; the computation has simply not yet entered the asymptotic regime where they take effect, and for this problem the entry ticket — $\delta t < 0.002$ — is punishingly small. This is the case for all explicit methods, which makes them unsuitable for stiff ODEs.

The remedy is implicit methods. Apply the implicit Euler method to the same problem, evaluating the slope at the new point:

Now the controlling factor is $\frac{1}{1 + 1000\delta t}$, which lies strictly between $0$ and $1$ for every positive $\delta t$. The term always decays, the solution always settles to $1$ — no oscillations, convergent at any step size. The tiny-step tyranny is gone, and one can take large steps to trace the flat line cheaply.

The reason is structural. An explicit method builds its approximation out of polynomials in $\delta t$ — and a polynomial (other than a constant) always runs off to $\pm\infty$ as its argument grows, so it cannot mimic a bounded, decaying solution once the step is large. An implicit method instead produces rational functions, and a rational function like $\frac{1}{1 + x}$ stays bounded — it tends to $0$ as $x \to \infty$ rather than blowing up — so it can track confined behavior at large step sizes. That is the deeper reason large steps work for implicit methods and fail for explicit ones. Ergo: for stiff ODEs, always use implicit methods.

The three phenomena at a glance

Three separate things, then, can make the numerical solution of an ODE hard, and the practical first question for each is always the same: whose fault is it — the problem’s or the method’s?

• Bad conditioning is a threatening property of the underlying problem, with nothing to do with the method at all. In the extreme it leaves very few numerical options, because no algorithm can undo a sensitivity baked into the problem.
• Instability is a threatening property of the method — it can, for instance, force tiny steps or blow up outright. Here implicit methods are often superior to explicit ones.
• Stiffness is a threatening property of the problem, but a tractable one: implicit methods handle it, and in fact for a stiff problem they are not just preferable but absolutely necessary.

The same split organizes the method properties met along the way. The table below gathers everything, sorted by what it is a property of — because that is what tells you whether a better algorithm can fix it (method properties) or whether you are stuck coping with it (problem properties).

The one relation worth carrying away from all of this binds the method properties together: consistency + stability = convergence. Consistency alone (the midpoint rule had it) is not enough; stability is the second ingredient that turns locally-accurate steps into a globally-correct answer.

Boundary value problems

Everything so far has been about initial value problems — the natural shape for population dynamics, where the state is known at the start and propagated forward. The boundary value problem (BVP) set aside at the very beginning is a different animal, and it needs a different machine.

The defining contrast is informational. An IVP hands you everything at one end and lets you march: from the state at $t_a$ you read off the slope and step forward, point by point. A BVP fixes the solution at both ends and asks for the path between, so there is nothing to march from — you cannot step left to right when the constraint that pins the trajectory down lives at the far end you have not reached yet. The fix is to stop marching and instead discretize the whole interval at once, turning the differential equation into a system of equations solved simultaneously.

A useful rule of thumb sets the stage: an ODE needs as many side conditions as its order to have a unique solution. A BVP supplies two conditions (one at each end), so it generically pairs with a second-order equation. The general second-order BVP is

Rather than the general $f$, take the linear special case — linear because, although it has products like $a(t)\dot y$, the unknown $y$ and its derivatives appear only to the first power and never multiplied together:

The simplest sub-case, $a(t) = 0$ with $b(t) > 0$, is the one to build the method on; under exactly those conditions the BVP is guaranteed a unique solution.

Discretizing the second derivative

Lay a grid over the interval. With $n$ subintervals the step size is $h = \delta t = (t_b - t_a)/n$, and the grid points are