Scientific Computing · Differential Equations

Local Truncation Error in the Explicit Euler Method

Where the O(h^2) per-step error in forward Euler comes from, why the global error drops to O(h), and how roundoff limits the step size from below.

Take the simplest initial value problem: $y’ = f(t, y)$ with $y(t_0) = y_0$, on an interval where $f$ is smooth and the solution exists and is unique. Forward Euler evaluates $f$ at the current point and steps along the tangent line:

$$ y_{n+1} = y_n + h,f(t_n, y_n), \qquad t_{n+1} = t_n + h $$

Almost every exposition states that the method is first order, error $O(h)$, and moves on. The distinction that matters for anyone choosing a step size is between the error committed in one step and the error that accumulates across the whole march, because the two differ by a factor of $1/h$ and have different dependencies on the problem.

One step against the exact solution

Define the local truncation error as the residual left when the exact solution is substituted into the recurrence. Let $y(t)$ solve the ODE exactly and expand around $t_n$:

$$ y(t_n + h) = y(t_n) + h,y’(t_n) + \frac{h^2}{2}y”(t_n) + \frac{h^3}{6}y'''(t_n) + \cdots $$

Since $y’ = f$, the first two terms of this expansion are exactly what Euler computes when the starting value is exact. The remainder of the first order is

$$ \tau_{n+1} = \frac{h^2}{2}y”(t_n) + O(h^3) $$

so the local error is $O(h^2)$ with a coefficient $y”/2$. Differentiating $y’ = f(t,y)$ gives $y” = f_t + f_y f$, so the leading coefficient depends on the partial derivatives of the right-hand side along the true trajectory. A problem where $f_y$ is large relative to typical scales of the solution has a large local error coefficient even if $f$ itself is small.

The word “truncation” refers to the fact that Euler keeps the Taylor series up to the linear term and truncates the rest. The truncation error is a property of the formula, not of the floating-point arithmetic that evaluates it, which is a separate contribution.

Why the accumulated error is one power worse

The bound that matters is on the global error $e_n = y_n - y(t_n)$. Each step introduces a local error of size $C h^2$, but there are $N = (t_{\text{end}} - t_0)/h$ steps, so a naive sum suggests the global error is $N \cdot C h^2 = C h (t_{\text{end}} - t_0)$, that is, $O(h)$. The naive sum is the right answer, but it is not a proof: it ignores that earlier errors are themselves propagated forward by the dynamics, and that the propagation can amplify them.

The standard argument uses a Lipschitz condition. Suppose $f$ satisfies $|f(t,u) - f(t,v)| \le L|u - v|$ for all relevant $u, v$. Subtracting the exact relation $y(t_{n+1}) = y(t_n) + h f(t_n, y(t_n)) + \tau_{n+1}$ from the Euler recurrence gives

$$ |e_{n+1}| \le |e_n| + h L |e_n| + |\tau_{n+1}| \le (1 + hL)|e_n| + C h^2 $$

Applying the recurrence repeatedly and using $1 + hL \le e^{hL}$ yields, for $t_n - t_0 = nh$,

$$ |e_n| \le \frac{C h}{L}\left(e^{L(t_n - t_0)} - 1\right) $$

after dropping a decaying term at the origin. Two facts follow. The global error is $O(h)$, one order below the local error, which is why a second-order local error produces a first-order method. And the constant is exponential in $L$ times the length of the interval: a stiff problem with large $L$ can have a modest local error and a useless global one.

The factor $e^{L T}$ is not an artefact of the proof. It is the honest statement that differential equations amplify perturbations, and it is the reason step-size selection is driven by the solution’s sensitivity rather than by the interval length.

A worked step: linear test problem

For $y’ = \lambda y$, Euler gives $y_{n+1} = (1 + h\lambda)y_n$, so $y_n = (1 + h\lambda)^n y_0$. The exact value is $y_0 e^{\lambda t_n}$. Compare the two at fixed $t_n$ as $h \to 0$:

$$ (1 + h\lambda)^{t_n/h} \to e^{\lambda t_n} $$

The convergence is first order. Comparing the logarithms, $\log(1 + x) = x - x^2/2 + O(x^3)$ with $x = h\lambda$, so each step carries a relative error of roughly $-(h\lambda)^2/2$, and $t_n/h$ steps accumulate

$$ \text{relative error} \approx -\frac{\lambda^2 h, t_n}{2} $$

for small $h\lambda$. The relative error is proportional to $h$ and to the curvature of the exponential, which is the general pattern: the error constant encodes how strongly the exact solution bends within a single step.

Implementing the recurrence is a handful of lines, and the structure of the loop determines how the error behaves:

/* Forward Euler on y' = f(t, y), uniform step h. */
for (size_t n = 0; n + 1 < steps; n++) {
    double t = t0 + (double)n * h;
    double k1 = f(t, y[n]);
    y[n+1] = y[n] + h * k1;
}

Two details in the loop are worth naming. Accumulating t by multiplication keeps the time grid free of drift. Repeated addition of h can differ from the nominal grid after many steps. And storing every sample costs $8$ bytes per component per step, which for a long run dwarfs the arithmetic cost. That memory pressure favors methods which carry only a few past states.

Euler also has no built-in error estimate, so any step-size control has to be built from something else: either a second integration at half the step, or a richer method whose stages supply two combinations. The second route leads to the embedded pairs of adaptive step-size control, and the algebraic conditions those methods satisfy are the subject of Butcher tableaux and order conditions. Both depend on the same Taylor expansion used above, extended term by term.

Where the step size stops helping

The truncation analysis says the error goes to zero as $h \to 0$. Arithmetic says otherwise. Each evaluation of the update introduces a rounding error on the order of the machine epsilon times the magnitude of $y$, call it $\varepsilon |y|$ with $\varepsilon \approx 2^{-52}$ for IEEE 754 double precision. Roundoff contributions also accumulate, and they too are amplified by the dynamics, giving a total of roughly

$$ \frac{\varepsilon |y| e^{LT}}{h} $$

since a shorter step means more steps and therefore more rounding decisions. The total error is approximately the sum of two terms moving in opposite directions in $h$. The minimum sits where they cross, at $h \approx \sqrt{\varepsilon}$ scaled by the problem constants, which for double precision is around $10^{-8}$ of the solution’s natural time scale.

This floor is why adaptive controllers stop respecting their own tolerance requests: below a certain step the error estimate stops falling and starts rising. It also explains a common symptom in the tolerance logic of adaptive solvers, where a tolerances below roughly the square root of machine epsilon are unreachable.

Failure modes specific to forward Euler

Euler is not merely low order. It has a structural weakness. On a decaying problem, $y’ = -a y$ with $a > 0$, the iteration is $y_{n+1} = (1 - ah) y_n$. The method is stable only if $|1 - ah| \le 1$, which requires $h \le 2/a$. A step above that threshold does not merely produce a large error. It produces a solution that alternates in sign and grows without bound, diverging from a decaying exact solution. The quadratic approximation of the exponential has infinite negative growth below $-2$, and the iteration maps it onto a growing sequence.

That bound is the first appearance of the stability constraint treated quantitatively in A-stability and the Dahlquist Barrier. For a system of equations the relevant $a$ is the largest-magnitude eigenvalue of the Jacobian, so a system mixing a slow component and a fast decaying component forces a step size set by the fast, short-lived one.

Two further degeneracies occur in practice. A right-hand side that is not Lipschitz continuous, or that is continuous only piecewise, defeats the error bound entirely, and the error can be $O(\sqrt{h})$ or worse near the discontinuity. And allocating storage per component while the integration only needs the previous value wastes memory for long runs. The loop above stores history only to keep the example inspectable, which is not what a real integrator does.

Euler remains the right method in exactly two situations: when the right-hand side is cheap enough that first-order accuracy at a fine step is affordable, and when the method is embedded as the first-order member of a step-doubling or Richardson extrapolation pair. For everything else, the cost of a better method is a small multiple of the f evaluations per step, and the error reduction is a power of $h$.

The stability constraint deserves one further comment. It is stated above for a real, negative eigenvalue, but a problem whose Jacobian eigenvalues lie off the real axis changes the geometry of the region rather than removing it: the admissible step is set by the largest eigenvalue magnitude, and the region’s shape decides how much of the complex plane is covered. That geometry is the subject of A-stability and the Dahlquist barrier, and it is why a scheme with the same order can be usable on one problem and hopeless on another.

A second constraint of the same kind appears when the problem changes character along the solution, which is where stiffness becomes the deciding property and the order of a method ceases to be the leading consideration.

A second structural limit sits alongside the stability constraint: order. The Taylor expansion that gives the local error also gives the conditions a higher-order method must satisfy, and those conditions are what make the per-step cost grow with order. There is no way to buy accuracy without either more stages per step or more steps, and the two currencies are not freely exchangeable when the right-hand side is expensive.

References

  1. Numerical Analysis, 10th edition; Chapter 5, Initial-Value Problems for Ordinary Differential Equations — Burden & Faires; Cengage
  2. Numerical Methods for Ordinary Differential Equations, 3rd edition — J. C. Butcher; Wiley
  3. Solving Ordinary Differential Equations I: Nonstiff Problems — Hairer, Norsett & Wanner; Springer