Solving ODEs Numerically: Step Size, Accuracy, and Stability

A smooth trajectory on screen may come from a sequence of finite numerical steps. Step size affects the error at each step, but also whether decay turns into growth or an oscillator gains energy. The double pendulum and Rössler attractor notes describe their models and RK4 implementations; the orbital playground uses velocity Verlet. A decaying equation and a harmonic oscillator, both with exact solutions, give a useful starting point for judging these integrators.

What a step approximates

An initial value problem specifies a rate of change and a starting state:

$$y'(t)=f(t,y(t)),\qquad y(t_0)=y_0.$$

The state $y$ may be a scalar or a vector containing quantities such as position and velocity. With step size $h>0$, compute approximations $y_n\approx y(t_n)$ at $t_n=t_0+nh$. An exact step satisfies

$$y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(t,y(t))\,dt.$$

The difficulty is that $y(t)$ inside the integral is still unknown. A numerical method approximates the increment using known states and several slope evaluations. Interpolating the result for display can make a curve smoother without improving the computed states. The error here mainly comes from discretizing time; rounding is a separate issue, covered in Floating-Point Numbers and Stable Computation.

Explicit Euler: local and global error

Section 6.2 of Driscoll and Braun’s Fundamentals of Numerical Computation gives explicit Euler: follow the slope at the start throughout the step.

$$y_{n+1}=y_n+h f(t_n,y_n).$$

Consider an equation whose answer can be checked directly:

$$y'=-2y,\qquad y(0)=1,\qquad y(t)=e^{-2t}.$$

For $h=0.1$, the recurrence is $y_{n+1}=0.8y_n$, so $y_n=0.8^n$.

Step $n$	Time $t_n$	Euler approximation $y_n$	Exact value $e^{-2t_n}$
0	0.0	1.000000000	1.000000000
1	0.1	0.800000000	0.818730753
2	0.2	0.640000000	0.670320046
10	1.0	0.107374182	0.135335283

Local one-step error measures the error in a step starting from the exact state. For a smooth solution, Taylor expansion gives

$$y(t_n+h)-\bigl[y(t_n)+hf(t_n,y(t_n))\bigr] =\frac{h^2}{2}y''(t_n)+O(h^3).$$

Euler’s one-step error is therefore $O(h^2)$. At the initial state, $y''(0)=4$, giving a leading term of $2h^2$. The actual absolute errors are $0.018730753$ for $h=0.1$ and $0.004837418$ for $h=0.05$, a reduction approaching a factor of four. The textbook defines local truncation error by dividing this error by $h$, so its quantity is $O(h)$. Check that convention before comparing error orders.

Global error compares the iterated value $y_n$ with $y(t_n)$, including propagation of earlier errors. The textbook’s convergence theorem assumes uniform bounds on the local error per unit step and on the state derivative of the step’s increment function. On a fixed finite interval, with a sufficiently smooth right-hand side and solution and these bounds satisfied, Euler has global error $O(h)$. Reaching the same endpoint takes about $1/h$ steps, and earlier errors affect later slopes; the last step’s local error cannot represent the whole computation.

At $T=1$, the global absolute errors for $h=0.1$ and $0.05$ are approximately $0.027961101$ and $0.013758629$, a ratio of $2.032$. Halving the step usually approaches a factor-of-two reduction for a first-order method.

Classical RK4: more slopes for greater accuracy

The classical fourth-order Runge–Kutta method in §6.4 evaluates four slopes per step. Here $k_i$ denotes a slope, without a factor of $h$:

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

On a fixed finite interval, with a sufficiently smooth right-hand side and solution and the preceding error-propagation conditions, the one-step error is $O(h^5)$ and the global error is $O(h^4)$. In the asymptotic convergence range, before rounding dominates, halving the step should reduce global error to about $1/16$ of its previous value.

For $y'=-2y$, substituting the four slopes gives the step multiplier

$$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},\qquad z=-2h.$$

This matches the first five terms of $e^z$. The approximation at $T=1$ can thus be checked independently using $R(-2h)^{1/h}$. The program below gives endpoint absolute errors of $4.265\times10^{-6}$, $2.452\times10^{-7}$ and $1.470\times10^{-8}$ for $h=0.1,0.05,0.025$. Successive ratios, $17.396$ and $16.682$, approach 16. A higher-order method still has to satisfy its stability conditions.

Absolute stability: when decay becomes growth

Section 11.3 defines absolute stability using $y'=\lambda y$: with $h$ fixed, does the numerical solution remain bounded as the number of steps tends to infinity? Setting $z=h\lambda$, explicit Euler gives

$$y_{n+1}=(1+z)y_n,\qquad y_n=(1+z)^n y_0.$$

Its stability region is $|1+z|\le1$, the closed disk of radius 1 centered at $-1$ in the complex plane. Actual decay to zero requires the strict inequality $|1+z|<1$. For real $\lambda<0$, this becomes

$$0<h<\frac{2}{|\lambda|}.$$

For $\lambda=-2$, the boundary is $h=1$. There the multiplier is $-1$: a nonzero solution alternates sign without decaying. Beyond the boundary its amplitude grows. With $h=1.1$ and $y_0=1$, the values are

$$1,\;-1.2,\;1.44,\;-1.728,\;2.0736.$$

The exact solution stays positive and decays. Even inside the stability region, Euler alternates sign when $0.5<h<1$, so stability alone ensures neither quantitative accuracy nor positivity. For complex $\lambda$, check $|1+h\lambda|$ rather than substituting $|\lambda|$ into the real-valued restriction.

RK4’s region on this test equation is $|R(z)|\le1$ and also has a boundary. Stability and global error order concern different limits: convergence at a fixed endpoint as $h\to0$ does not allow an indefinitely long simulation at a fixed $h$.

Stiffness and implicit Euler

Section 11.4 explains stiffness through multiple time scales: fast decaying modes restrict an explicit method’s step while the slow evolution of interest requires a long observation. For example,

$$x'=-x,\qquad w'=-1000w,\qquad x(0)=w(0)=1.$$

The decay time scales are 1 and $0.001$. The slow component continues evolving after the fast component has become tiny, yet explicit Euler still requires $h<0.002$ to damp the fast mode. Rounding and discretization errors can excite that mode again. Stiffness depends on the equation, method and required accuracy; the apparent steepness of a curve is insufficient to identify it.

Implicit Euler evaluates the slope at the unknown next state:

$$y_{n+1}=y_n+h f(t_{n+1},y_{n+1}).$$

Rearranging for the test equation gives

$$y_{n+1}=\frac{y_n}{1-h\lambda}.$$

The multiplier yields the stability region $|1-z|\ge1$, containing the entire left half-plane, so the method is A-stable. For real $\lambda<0$, any $h>0$ damps the mode. The fast component with $h=0.1$, for instance, has multiplier $1/101$. This permits stepping over a fast transient; resolving that transient still requires smaller steps. Implicit Euler retains first-order global error. Stability does not increase its order of accuracy.

The textbook’s implementation of implicit methods in §6.7 explains why the unknown next state generally requires rootfinding. For a nonlinear equation, each implicit Euler step requires solving $G(u)=0$, where

$$G(u)=u-y_n-hf(t_{n+1},u).$$

One Newton correction solves

$$\bigl[I-hJ_f(t_{n+1},u^{(k)})\bigr]\Delta u=-G(u^{(k)}), \qquad u^{(k+1)}=u^{(k)}+\Delta u.$$

Here $J_f$ is the state Jacobian of $f$ and $I$ is the identity matrix. Evaluating or approximating the Jacobian, solving linear systems and repeating iterations all cost work. Check convergence and keep the solve error below the time-discretization error budget; a failed implicit solve is not an acceptable step.

Long-term oscillator energy: Euler and symplectic methods

Take a harmonic oscillator with unit mass and unit angular frequency:

$$q'=p,\qquad p'=-q,\qquad (q(0),p(0))=(1,0).$$

Its exact solution is $q(t)=\cos t$, $p(t)=-\sin t$, with constant energy $E=(q^2+p^2)/2=1/2$. Explicit Euler updates both variables from their old values:

$$q_{n+1}=q_n+hp_n,\qquad p_{n+1}=p_n-hq_n.$$

Adding the two squares cancels the cross terms, giving

$$E_{n+1}=(1+h^2)E_n.$$

Any fixed positive step increases energy repeatedly. With $h=0.1$ and $T=100$, there are 1000 steps and $E=\tfrac12(1.01)^{1000}\approx10479.577819$. Implicit Euler does the opposite: it divides energy by $1+h^2$ each step, leaving approximately $0.000024$ at the same endpoint. Its bounded trajectory has acquired artificial damping.

Lecture 2 of Hairer’s Geometric Numerical Integration course gives symplectic Euler. For this separable Hamiltonian system, update momentum first, then position using the new momentum:

$$p_{n+1}=p_n-hq_n,\qquad q_{n+1}=q_n+hp_{n+1}.$$

“Symplectic” refers to preserving the geometric structure of Hamiltonian flow. In this two-dimensional example, the update matrix has determinant 1 and preserves phase-space area. The method is still first-order and does not preserve the original energy at every step. Substitution into the update equations shows that, for this oscillator, it exactly preserves a modified energy:

$$H_h=\frac12(q^2+p^2-hqp).$$

For fixed $0<h<2$, this quadratic form is positive definite. Using $|qp|\le E$ and initial $H_h=1/2$ gives

$$\frac{1/2}{1+h/2}\le E_n\le\frac{1/2}{1-h/2}.$$

Thus for $h=0.1$, energy oscillates between bounds of approximately $0.476190$ and $0.526316$. These bounds are specific to this linear oscillator and cannot simply be transferred to every Hamiltonian system. The invariant also requires constant $h$; varying the step invalidates this argument. Bounded energy can coexist with accumulated phase error.

Theorem 2 in the same lecture gives second-order velocity Verlet: a half momentum step, a full position step and another half momentum step. For $q''=a(q)$,

$$\begin{aligned} p_{n+1/2}&=p_n+\frac h2a(q_n),\\ q_{n+1}&=q_n+hp_{n+1/2},\\ p_{n+1}&=p_{n+1/2}+\frac h2a(q_{n+1}). \end{aligned}$$

It is also symplectic for separable Hamiltonian systems; the orbital playground uses it for central gravity. The example below sets $a(q)=-q$ to compare long-term energy behavior and endpoint state error directly.

A runnable comparison

Save the following as ode_compare.py and run python3 ode_compare.py. It uses only the Python standard library. The Euler and RK4 functions accept state sequences of any length. f(t, y) must return the same number of derivative components, or the step raises ValueError. The two implicit functions are closed-form updates specifically for these linear examples; the symplectic Euler and Verlet functions are specific to this oscillator.

from math import cos, exp, hypot, sin


def slope(f, t, y):
    values = tuple(f(t, y))
    if len(values) != len(y):
        raise ValueError("f(t, y) must have the same length as y")
    return values


def euler(f, t, y, h):
    return tuple(a + h * b for a, b in zip(y, slope(f, t, y)))


def rk4(f, t, y, h):
    def shift(k, scale):
        return tuple(a + scale * b for a, b in zip(y, k))
    k1 = slope(f, t, y)
    k2 = slope(f, t + h / 2, shift(k1, h / 2))
    k3 = slope(f, t + h / 2, shift(k2, h / 2))
    k4 = slope(f, t + h, shift(k3, h))
    return tuple(a + h * (b + 2 * c + 2 * d + e) / 6
                 for a, b, c, d, e in zip(y, k1, k2, k3, k4))


def decay(t, y):
    return (-2 * y[0],)


def decay_implicit(f, t, y, h):
    return (y[0] / (1 + 2 * h),)


def oscillator(t, y):
    q, p = y
    return (p, -q)


def oscillator_implicit(f, t, y, h):
    q, p = y
    return ((q + h * p) / (1 + h * h),
            (p - h * q) / (1 + h * h))


def symplectic_euler(f, t, y, h):
    q, p = y
    p = p - h * q
    return (q + h * p, p)


def verlet(f, t, y, h):
    q, p = y
    half_p = p - h * q / 2
    q = q + h * half_p
    return (q, half_p - h * q / 2)


print("local Euler: h, one-step error")
for h in (0.1, 0.05):
    print(f"{h:.3f} {abs(1 - 2 * h - exp(-2 * h)):.9f}")

print("decay at T=1: method, h, y, absolute error")
for name, step in (("Euler", euler), ("implicit", decay_implicit), ("RK4", rk4)):
    errors = []
    for n in (10, 20, 40):
        h, y = 1 / n, (1.0,)
        for j in range(n):
            y = step(decay, j * h, y, h)
        error = abs(y[0] - exp(-2))
        errors.append(error)
        print(f"{name:8s} {h:.3f} {y[0]:.9f} {error:.3e}")
    print("error ratios: " + " ".join(f"{a / b:.3f}"
          for a, b in zip(errors, errors[1:])))

print("Euler with h=1.1: y0 through y4")
y, values = (1.0,), [1.0]
for j in range(4):
    y = euler(decay, j * 1.1, y, 1.1)
    values.append(y[0])
print(" ".join(f"{value:.4f}" for value in values))

print("oscillator at T=100, h=0.1: method, E_min, E_max, E_final, state_error")
for name, step in (("Euler", euler), ("implicit", oscillator_implicit),
                   ("RK4", rk4), ("symplectic", symplectic_euler), ("Verlet", verlet)):
    y, energies = (1.0, 0.0), [0.5]
    for j in range(1000):
        y = step(oscillator, j * 0.1, y, 0.1)
        energies.append((y[0] ** 2 + y[1] ** 2) / 2)
    error = hypot(y[0] - cos(100), y[1] + sin(100))
    print(f"{name:10s} {min(energies):.6f} {max(energies):.6f} "
          f"{energies[-1]:.6f} {error:.3e}")

Output:

local Euler: h, one-step error
0.100 0.018730753
0.050 0.004837418
decay at T=1: method, h, y, absolute error
Euler    0.100 0.107374182 2.796e-02
Euler    0.050 0.121576655 1.376e-02
Euler    0.025 0.128512157 6.823e-03
error ratios: 2.032 2.016
implicit 0.100 0.161505583 2.617e-02
implicit 0.050 0.148643628 1.331e-02
implicit 0.025 0.142045682 6.710e-03
error ratios: 1.966 1.983
RK4      0.100 0.135339548 4.265e-06
RK4      0.050 0.135335528 2.452e-07
RK4      0.025 0.135335298 1.470e-08
error ratios: 17.396 16.682
Euler with h=1.1: y0 through y4
1.0000 -1.2000 1.4400 -1.7280 2.0736
oscillator at T=100, h=0.1: method, E_min, E_max, E_final, state_error
Euler      0.500000 10479.577819 10479.577819 1.438e+02
implicit   0.000024 0.500000 0.000024 9.935e-01
RK4        0.499993 0.500000 0.499993 8.332e-05
symplectic 0.476190 0.526316 0.521321 5.665e-02
Verlet     0.498750 0.500000 0.499724 4.222e-02

Energy minima and maxima include the initial state and every grid point. state_error is the Euclidean distance from the endpoint state to $(\cos100,-\sin100)$. RK4 has a small error on this finite interval; the symplectic methods have bounded energy variation but still differ from the exact state. Compare endpoint accuracy, phase and long-term invariants separately when choosing a method, rather than ranking methods solely by order.

Adaptive step control and checking a solution

Section 6.5 introduces embedded error estimation: compute two approximations of different orders in the same step and use their difference to estimate the lower-order formula’s local error. Reject a step exceeding the target and retry with smaller $h$; accept it only when the estimate is small enough. If an estimate $\eta$ scales approximately as $h^{r+1}$, a possible adjustment is

$$h_{\mathrm{new}}=s h\left(\frac{\mathrm{tol}}{\eta}\right)^{1/(r+1)},\qquad 0<s<1,$$

with limits on step growth and shrinkage. Here $r$ is the order of the formula whose error is estimated, and $s$ is a safety factor. The exponent must match the error estimator’s order; this is a local-control model, not a universal update rule for every solver.

The SciPy solve_ivp reference describes RK45 as an embedded 5(4) method: estimate error with the fourth-order formula and advance with the fifth-order one. It differs from classical RK4 above. The local error scale is atol + rtol * abs(y); choose an appropriate atol near zero, with separate values for components of different scales. Use RK45 for nonstiff problems and Radau or BDF for stiff ones. These local-control settings do not guarantee the endpoint error.

Keep the comparison target fixed when checking a solution:

1. Halve the step at the same endpoint. For an order-$r$ method, compute $y_h,y_{h/2},y_{h/4}$ in the asymptotic range. The difference ratio $\|y_h-y_{h/2}\|/\|y_{h/2}-y_{h/4}\|$ should approach $2^r$; estimate the error in $y_{h/2}$ by $\|y_{h/2}-y_h\|/(2^r-1)$. For the finest grid, estimate the error in $y_{h/4}$ by $\|y_{h/4}-y_{h/2}\|/(2^r-1)$. This requires stability, sufficient smoothness, a nonzero leading error term and rounding that does not yet dominate. For whole trajectories, compare common times.
2. Check the model’s invariants. For this oscillator, check energy alongside phase or state error; for central gravity, check energy and angular momentum. Close invariants do not establish an accurate trajectory, and a dissipative model should not be forced to conserve energy.
3. Use an independent reference. Prefer an exact solution when available. Otherwise choose a reference solver suited to the equation, tighten its tolerances and tighten them again to check the reference is stable. When comparing solve_ivp with a fixed-step program, t_eval selects stored output times, not internal step sizes. Use identical equations, parameters, initial conditions and comparison times.

Chaos amplifies small state differences. For the double pendulum and Rössler attractor, first check step convergence on a finite short interval, then compare relevant invariants or statistical features over longer times. A diverging long trajectory alone cannot distinguish sensitivity to initial conditions from integration error.

My reading · Opened notes are remembered in this browser.