Sensitive dependence, and what a solver can promise

6 September 2026

A differential equation promises less than it appears to. It fixes the future completely — and yet, for some systems, knowing the present to forty decimal places buys you only a few dozen units of time. This post works through where that boundary sits, and it doubles as the page every part of this site is measured against: every figure below is drawn in your browser, and every one of them can be poked at.

What determinism actually gives you

The classical existence theorem costs almost nothing to state, and it is worth being precise about what it does not say.

Uniqueness is a statement about exact real numbers. It says two solutions from the same point agree; it says nothing about two solutions from points a billionth apart. That gap is the whole subject.

Start with a case where the gap stays closed. The logistic equation

has every solution converging to $x=K$, whatever it starts from. Nearby starts stay nearby, forever: the flow is contracting, and a solver’s early rounding error gets forgotten rather than amplified.

The same picture in the phase plane, for a damped oscillator rather than a logistic growth, shows the other classical ending: everything spirals into a single point.

Where it breaks

Lorenz’s system is the standard counterexample, and it is worth writing out because its nonlinearity is so slight — two quadratic terms, nothing more:

with the usual $\sigma =10$, $\rho =28$, $\beta =8/3$.

Sampling one long orbit densely enough turns the trajectory into a picture of the attractor itself — the set the system settles onto, of which any particular trajectory is only one thread.

Compare that with the surface swept out by the logistic solutions of $\text{(1)}$ as the initial condition varies. Every slice at fixed $x_0$ is one solution curve, and the whole sheet is smooth, monotone and dull — which is exactly the point.

The cost of a billionth

Take two Lorenz trajectories whose starting points differ by ${10}^{-9}$ in the $x$ coordinate alone, and watch the distance between them.

The straight section is the content of the figure: on a log axis, exponential growth is a line, and its slope is the leading Lyapunov exponent $\lambda \approx 0.906$. So

The logarithm in $\text{(3)}$ is the bad news. Buying ten more decimal places of initial accuracy buys about $10\log 10/\lambda \approx 25$ more time units of prediction — and then you are back where you started.

Where the exponential comes from

Nothing so far explains why the separation should grow exponentially rather than, say, quadratically. The answer is that a small perturbation obeys its own linear equation — the variational equation — whose coefficient matrix is the Jacobian of the field, evaluated along the trajectory:

Discretise $\text{(4)}$ on the same grid the solver uses, write $M_k$ for the one-step amplification $I+hJ(x_k)+O(h^2)$, and stack the perturbations at every step into one vector. The whole history then satisfies a single linear system:

The blanks are genuine zeros, and there are a great many of them: of the $n^2$ blocks in $\text{(5)}$ only $2n-1$ are non-zero, two per row at worst. A matrix this sparse is never assembled, let alone inverted. It is lower triangular, so forward substitution solves it in one sweep — which is precisely what running the solver forward in time is:

So the amplification of an initial error is a product of matrices, one per step. Products grow or shrink geometrically, and the average logarithmic growth rate of $\Vert M_{n-1}\cdots M_0\Vert$ is the definition of the leading Lyapunov exponent. The straight line in the separation plot is $\text{(6)}$ seen on a log axis, and $\lambda \approx 0.906$ is its slope.

This also settles what a better integrator can and cannot do. A higher-order method computes each $M_k$ more accurately; it does not make the product smaller. The exponent belongs to $J$, and $J$ belongs to the equation.

What the integrator controls

None of this excuses a sloppy solver. Over the interval where prediction is meaningful, the method’s order is exactly what decides how much work a given accuracy costs.

A method of order $p$ has global error $O(h^p)$, so on log–log axes each method should trace a straight line of slope $p$. That is a claim worth checking rather than repeating:

The three orders come out of a Taylor expansion that each method truncates at a different place:

// One RK4 step. The weights are chosen so the error terms through h⁴ cancel.
export function rk4Step(fn, t, state, h) {
  const k1 = fn(t, state);
  const k2 = fn(t + h / 2, add(state, k1, h / 2));
  const k3 = fn(t + h / 2, add(state, k2, h / 2));
  const k4 = fn(t + h, add(state, k3, h));
  return state.map((v, i) => v + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]));
}

What to take away

1. Uniqueness is not predictability. Theorem 1.2 holds for the Lorenz system as much as for the logistic one; it simply does not say what you want.
2. The horizon grows logarithmically. By $\text{(3)}$, precision is the expensive way to buy time, and it is the only way.
3. Order still matters. Within the horizon, the difference between Euler and RK4 is seven orders of magnitude for four times the work.

A theory is the more impressive the greater the simplicity of its premises, the more different kinds of things it relates, and the more extended its area of applicability.

— Einstein, on why two quadratic terms deserve this much attention

Everything above renders from about three kilobytes of figure recipes; the trajectories are integrated in your browser rather than shipped as data. See the spike notes for the measurements behind that choice.