Case Study 1: The Derivative That Got Worse

DataField.Dev

Case Study 1: The Derivative That Got Worse

"The first principle is that you must not fool yourself — and you are the easiest person to fool." — Richard Feynman

Executive Summary

You have inherited a small Fortran module that computes forces in a molecular-dynamics-style code by numerically differentiating a potential-energy function, $F = -\,dU/dx$. It works — until a colleague "improves its accuracy" by shrinking the finite-difference step, after which the forces come back subtly noisier and the simulation's energy stops being conserved. Your job is to read the routine, reproduce the regression, diagnose it, and fix it without changing what it computes. The culprit is the round-off floor of §22.4: past an optimal step $h^{*}\!\sim\!\sqrt\varepsilon$, catastrophic cancellation (Chapter 20) makes a forward difference less accurate, not more. The fix is a better formula (central difference), a principled step size, and — where possible — an analytic derivative.

Skills applied - Reading a one-sided finite difference and identifying its order (§22.1). - Recognizing the truncation-vs-round-off trade-off and the optimal step (§22.4). - Connecting a "smaller $h$ made it worse" bug to catastrophic cancellation (Chapter 20). - Replacing a first-order formula with a second-order one and measuring the improvement (§22.1, §22.4).

Background

The code integrates particles under a potential $U(x)$; the force on a particle is $F = -\,dU/dx$. For most potentials the derivative is known analytically, but this module was written to accept an arbitrary user-supplied $U$ as a procedure, so it differentiates numerically. That is a legitimate design choice — but it makes the module a live specimen of everything §22.4 warned about, and the "accuracy improvement" that broke it is one of the most common mistakes in scientific programming.

For a reproducible audit we use a test potential with a known answer, $U(x) = x^2$, whose exact derivative is $U'(x) = 2x$, so at $x = 2$ the true force magnitude is $|U'(2)| = 4$. Any deviation from $4$ is pure numerical error, and we can see it to the digit.

Phase 1 — Read the Inherited Routine

Here is the differentiator, lightly cleaned up. Read it before reading on.

module forces
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  private
  public :: dp, dudx_forward

  abstract interface
    pure function potential(x) result(u)
      import :: dp
      real(dp), intent(in) :: x
      real(dp)             :: u
    end function potential
  end interface

  real(dp), parameter :: h_step = 1.0e-13_dp     ! "smaller = more accurate", allegedly

contains

  pure function dudx_forward(u, x) result(dudx)
    procedure(potential)  :: u
    real(dp), intent(in)  :: x
    real(dp)              :: dudx
    dudx = (u(x + h_step) - u(x)) / h_step        ! forward difference, O(h)
  end function dudx_forward

end module forces

Two facts jump out once you have read §22.1 and §22.4. First, it is a forward difference — the least accurate of the three, only $O(h)$. Second, and fatally, the step is hard-coded to $h = 10^{-13}$, on the folk theory that a smaller step always means a better derivative. It is exactly the theory the chapter demolished.

Phase 2 — Reproduce the Regression

The module's git history shows the step was once $10^{-2}$ and someone changed it to $10^{-13}$ "for accuracy." Let us see what each choice actually delivers for $U'(2) = 4$. With the original sane step, the forward difference is honestly first-order and fully hand-checkable:

! forward difference of x^2 at x=2, h=1e-2:
! (U(2.01) - U(2)) / 0.01 = (4.0401 - 4) / 0.01 = 0.0401/0.01 = 4.01

h = 1.0e-2 : dU/dx = 4.0100000000 , error = 1.00e-02   (= h; a clean O(h) result)

Error $0.01 = h$, precisely the forward difference's truncation error for $f'' = 2$ (leading term $\tfrac{h}{2}f'' = h$). Now the "improved" step. The following table is illustrative — the exact digits in the round-off-dominated regime depend on the machine, so we do not certify them — but the shape is real and reproducible, and it is the whole story:

  forward dU/dx for U(x)=x^2 at x=2  (true value 4), double precision (illustrative)
    h          dU/dx           |error|      what dominates
  1e-2         4.0100          1e-2         truncation  (error ~ h)
  1e-4         4.000100        1e-4         truncation
  1e-8         4.00000002      ~2e-8        *** near the optimal h* ***
  1e-11        4.0000027       ~3e-6        round-off creeping in
  1e-13        4.00036         ~4e-4        round-off dominates (cancellation)
  1e-15        3.6 or 4.4      ~4e-1        round-off catastrophe

The "improvement" to $10^{-13}$ moved the routine from the truncation-dominated left side of the valley, past the optimal step near $10^{-8}$, and onto the round-off-dominated right side — trading a controllable $O(h)$ error for an uncontrollable cancellation error thousands of times larger. The forces are now noisy in their fourth digit, and that noise, injected every timestep, is what breaks energy conservation.

Phase 3 — Diagnose: The Round-Off Floor

The mechanism is precisely §22.4's pitfall. The numerator $u(x+h) - u(x)$ subtracts two values that, as $h \to 0$, agree in ever more leading digits; the subtraction is catastrophic cancellation (Chapter 20), leaving mostly round-off, which the division by the tiny $h$ then amplifies. Truncation error falls like $\tfrac{h}{2}|U''|$; round-off error grows like $\varepsilon|U|/h$. Their sum is minimized at

$$ h^{*} \approx 2\sqrt{\frac{\varepsilon\,|U|}{|U''|}} \sim \sqrt{\varepsilon} \approx 10^{-8}, $$

with a best-possible error near $10^{-8}$ — half the machine's digits, gone, no matter how carefully the code is written, simply because a forward difference subtracts near-equal numbers. The colleague chose $h$ five orders of magnitude below the optimum.

The diagnosis in one sentence: the routine is not buggy in the sense of a typo — it is numerically mis-designed, using the least accurate formula at a step size on the wrong side of the round-off floor.

Phase 4 — Fix It

Three changes, in order of impact.

(1) Use a central difference. It is $O(h^2)$ instead of $O(h)$, so it reaches a better floor ($\sim\!\varepsilon^{2/3}\approx 10^{-11}$) at a larger, safer step ($h^{}\sim\varepsilon^{1/3}\approx 10^{-5}$). For our test potential it is even better than that: $U = x^2$ has $U''' = 0$, so the central difference is exact* for any step, which makes the fix vivid.

(2) Choose $h$ from $\varepsilon$, not from folklore. Scale the step to the point and to machine epsilon rather than hard-coding it: $h = \varepsilon^{1/3}\max(|x|, 1)$ is the standard rule for a central difference.

(3) Prefer the analytic derivative when you have it. If the user can supply $U'$, take it; numerical differentiation is a fallback, not a default.

Here is the repaired differentiator and its certified output on the test potential:

pure function dudx_central(u, x) result(dudx)
  procedure(potential)  :: u
  real(dp), intent(in)  :: x
  real(dp)              :: dudx, h
  h = epsilon(1.0_dp)**(1.0_dp/3.0_dp) * max(abs(x), 1.0_dp)   ! ~ 6e-6 near x=2
  dudx = (u(x + h) - u(x - h)) / (2.0_dp * h)                  ! central difference, O(h^2)
end function dudx_central

For $U(x) = x^2$ at $x = 2$, the central difference returns the true $2x = 4$ to full precision, because a quadratic has no third derivative to leave a truncation error behind:

central dU/dx at x=2 : 4.0000000000 , error ~ 1e-15   (exact for a quadratic, up to round-off)

To confirm the fix is not an artefact of the too-easy $U = x^2$, apply the same central difference to $U(x) = x^4$ at $x = 1$ (true $U'(1) = 4$), where the central formula has a genuine $O(h^2)$ error. With a sane $h = 0.01$ the estimate is the exactly hand-computable $4 + 4h^2 = 4.0004$, error $4\times10^{-4}$ — and it shrinks as you refine toward $h^{*}$, the opposite of the broken routine.

Phase 5 — Sanity Check and Hand-Off

Before shipping, run the §22.4 convergence check as a regression test: differentiate $U = x^4$ at $x = 1$ on $h = 0.1, 0.05, 0.025$ and confirm the error ratio is $\approx 4$ (second order). If a future edit accidentally reintroduces a one-sided formula, the ratio will fall to $\approx 2$ and the test will catch it. Document the step-size rule in a comment citing this chapter, and note in the module header that the analytic derivative should be supplied when available. The forces are now smooth, energy conservation is restored, and the module carries a test that encodes why.

Discussion Questions

The broken routine used $h = 10^{-13}$. Estimate its error for $U = x^2$ from the round-off model $\varepsilon|U|/h$, and compare with the illustrative table.
Why is $U = x^2$ a dangerously reassuring test case for a central-difference fix? What does adding the $U = x^4$ test buy you?
The fix scales $h$ as $\varepsilon^{1/3}\max(|x|,1)$. Why the $\max(|x|,1)$ factor — what goes wrong near $x = 0$ without it, and what goes wrong for very large $x$?
Energy drift is a systematic symptom of a random-looking force error. Why does noise in the force, rather than a constant bias, break energy conservation over many steps?

Your Turn: Extensions

Option A (analyze). Instrument the original routine to print $u(x+h) - u(x)$ (the raw numerator) for $h$ from $10^{-2}$ to $10^{-15}$ and watch the number of significant digits collapse. At what $h$ does the numerator lose half its digits?
Option B (build). Add a dudx_richardson that Richardson-extrapolates two central differences at $h$ and $h/2$ (order $p = 2$) to get an $O(h^4)$ derivative, and measure the new convergence ratio ($\approx 16$).
Option C (extend). Generalize the module to a gradient of a multivariable $U(\mathbf{x})$ by central-differencing one coordinate at a time, and estimate the cost: how many function evaluations for a gradient in $d$ dimensions, and how does that motivate analytic or automatic differentiation?

Key Takeaways

"Smaller $h$ is more accurate" is false for finite-difference derivatives; there is an optimal step near $\sqrt\varepsilon$ (forward) or $\varepsilon^{1/3}$ (central), and going below it loses accuracy to cancellation.
A "smaller $h$ made it worse" report is a fingerprint of the round-off floor — reach for Chapter 20, not a debugger.
The cheapest robustness upgrade is a central difference at a $\varepsilon$-scaled step; the best is an analytic derivative when one exists.
Encode the fix as a convergence-ratio regression test so the mistake cannot silently return.