Can You Prove Sensitivity to Initial Conditions Numerically?

A numerical experiment can demonstrate how nearby solutions of an ordinary differential equation separate. It cannot, by itself, prove sensitivity to initial conditions for the exact system: computations have finite precision, and a plot only covers a finite time interval.

The useful goal is therefore to make the evidence reproducible and check that the observed separation is not a numerical artifact. The Lorenz system below provides a concrete example and a way to estimate the growth rate.

What sensitivity means

For a flow, write $\phi_t(x)$ for the state at time $t$ that starts at $x$. For a flow on a metric space $X$, one common definition of sensitivity is that there is a fixed distance $\delta>0$ such that

$$\forall x\in X,\ \forall r>0,\ \exists y\in X,\ \exists t\ge0:\quad 0<d(x,y)<r,\quad d(\phi_t(x),\phi_t(y))>\delta.$$

The exact quantifiers and setting depend on the definition being used. The key idea is that arbitrarily small changes in initial state can eventually produce a difference greater than a fixed distance. A single pair of trajectories separating does not establish this full property.

In practice, researchers often look for an interval where the distance between nearby trajectories grows approximately exponentially. If $d(t)$ is their separation, then

$$d(t)\approx d(0)e^{\lambda t}, \qquad \log d(t)\approx \log d(0)+\lambda t.$$

A positive slope $\lambda$ over an appropriate interval is evidence of local exponential separation. It is a finite-time estimate, not automatically a proof that the system has a positive asymptotic Lyapunov exponent.

A numerical test with the Lorenz system

Consider the classical Lorenz equations:

$$\begin{aligned} \dot{x}&=10(y-x),\\ \dot{y}&=x(28-z)-y,\\ \dot{z}&=xy-\frac{8}{3}z. \end{aligned}$$

Integrate once from $(1,1,1)$ and again from a nearby state $(1+\varepsilon,1,1)$. Keep all parameters the same. Then calculate the Euclidean distance between the two computed states at each time.

This Python example uses SciPy’s adaptive solve_ivp solver. It fits a line to $\log d(t)$ while the distance is between $10^{-6}$ and $10^{-2}$; that interval is meant to capture early growth before the trajectories become widely separated.

import numpy as np
from scipy.integrate import solve_ivp

sigma, rho, beta = 10.0, 28.0, 8.0 / 3.0

def lorenz(t, state):
    x, y, z = state
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

times = np.linspace(0.0, 40.0, 4001)
epsilon = 1e-8
initial_a = [1.0, 1.0, 1.0]
initial_b = [1.0 + epsilon, 1.0, 1.0]

options = dict(method="DOP853", rtol=1e-12, atol=1e-14, t_eval=times)
sol_a = solve_ivp(lorenz, (times[0], times[-1]), initial_a, **options)
sol_b = solve_ivp(lorenz, (times[0], times[-1]), initial_b, **options)

if not (sol_a.success and sol_b.success):
    raise RuntimeError("An integration failed")

distance = np.linalg.norm(sol_a.y - sol_b.y, axis=0)
mask = (distance > 1e-6) & (distance < 1e-2)

if np.count_nonzero(mask) >= 2:
    slope, intercept = np.polyfit(times[mask], np.log(distance[mask]), 1)
    print(f"Fitted log-distance growth rate: {slope:.3f} per time unit")
else:
    print("Too few points in the selected growth interval; adjust the run or interval.")

The fitted slope is an estimate from that run; it is not a universal constant guaranteed by this particular computation. Plotting times against distance and against np.log(distance) can help you check whether the selected portion is actually close to a straight line. Do not fit through the whole run: once the trajectories are far apart, their distance typically stops growing exponentially and fluctuates as they move around the attractor.

How to judge the result

Repeat the experiment with several values of $\varepsilon$, for example $10^{-6}$, $10^{-8}$, and $10^{-10}$. Also tighten the solver tolerances and, if possible, compare with another suitable solver. A credible early growth rate should not depend dramatically on one arbitrary choice of perturbation or numerical settings.

The perturbation must be larger than the effective numerical error. If it is too small, rounding and integration error can dominate the difference you intend to measure. If it is too large, you may skip the close-trajectory regime altogether. Tight tolerances help, but they do not remove all numerical error or certify the exact solution.

What numerical evidence can and cannot establish

A graph showing two trajectories eventually far apart establishes only that those two computed trajectories separated during the simulated time. It does not establish the quantifiers in the formal definition, nor does it prove that the same behavior continues indefinitely. For a rigorous result, use an analytical argument or a validated numerical method, such as interval arithmetic with error bounds, appropriate to the system and claim.

Numerics are still valuable: they can reveal likely sensitivity, estimate a finite-time growth rate, and help identify the behavior a proof should explain. State the equations, initial states, perturbation size, solver, tolerances, time interval, and how you calculated the rate so another reader can assess the evidence.

Common mistakes

• Calling divergence a proof of chaos. Sensitivity is one feature associated with chaos, but a plot alone does not establish a full definition of chaos.
• Interpreting every separation as exponential growth. Look for an interval where $\log d(t)$ is approximately linear, not merely a late-time gap.
• Ignoring numerical error. A computed difference may reflect rounding or solver error, especially when the initial perturbation is extremely small.
• Extrapolating beyond the run. A finite simulation cannot establish indefinite behavior or prove a long-term property for every nearby initial condition.