Numerical Differentiation & Integration, Revisited

Numerical Methods & Floating-Point Computation

Chapter 9 · Numerical Differentiation & Integration, Revisited

Calculus & Optimization Chapters 4 and 9 introduced numerical differentiation and integration as basic techniques and this course's own Chapter 4 explained why the simplest version breaks down. This chapter goes further: three genuinely better techniques for differentiation, and one genuinely better strategy for integration — each directly informed by what Chapters 1-8 already established about cancellation, stability, and conditioning.

Central Differences: A Free Upgrade Over Forward Differences

The forward-difference formula, (f(x+h)-f(x))/h, uses one point ahead of x. The central-difference formula uses one point on each side: (f(x+h) − f(x−h)) / (2h). Both approximate the same derivative, but their error behaves very differently as h shrinks — forward error shrinks proportionally to h itself, while central error shrinks proportionally to h².

Verified directly — for f(x)=sin(x) at x=1 (true derivative cos(1)≈0.5403023058681398)

At h=10⁻⁴: forward error ≈4.21×10⁻⁵, central error ≈9.00×10⁻¹⁰ — central is already almost five orders of magnitude more accurate for the same step size and the same number of extra function evaluations. Both formulas eventually hit the same cancellation wall this course's Chapter 4 already explained: at h=10⁻¹⁶, forward error balloons to 0.54 (completely wrong) while central error reaches 0.0148 — degraded, but noticeably more resistant to the collapse than forward differences, since central differencing's symmetric structure cancels out more of the leading error terms before cancellation error takes over.

Complex-Step Differentiation: Sidestepping Cancellation Entirely

Both formulas above eventually fail because they subtract two nearly-equal real numbers. Complex-step differentiation is a genuinely different trick: evaluate the function at a small imaginary step, f(x + ih), and take the imaginary part: f'(x) ≈ Im(f(x+ih)) / h. For a function that's analytic (most ordinary functions — polynomials, sin, exp, etc. all qualify), this formula has essentially no cancellation error, because it never subtracts two real, nearly-equal quantities at all — the real and imaginary parts of a complex number are stored and manipulated independently.

import cmath
derivative = cmath.sin(x + 1j*h).imag / h
    

Verified directly — accurate to the exact same double-precision value all the way down to h=10⁻¹⁰⁰

Computing f'(1) for f(x)=sin(x) via complex-step differentiation at h=10⁻⁸, 10⁻¹⁶, 10⁻²⁰, 10⁻³⁰, and even 10⁻¹⁰⁰: every single one returned exactly 0.5403023058681398 — the true value to full double precision, with a measured error of 0.000 at every tested h. There is no "too small" step size for this method, unlike forward or central differences, both of which collapsed catastrophically once h got small enough.

Why this isn't a free lunch in general

Complex-step differentiation requires the function to support complex-number arithmetic internally (which ordinary sin, cos, exp, and polynomial code usually does, but a function containing abs(), comparisons, or other non-analytic operations generally doesn't handle correctly), and it only computes first derivatives cleanly — it's a specialized tool for a specific situation, not a universal replacement for finite differences.

Richardson Extrapolation: Cancel the Leading Error Term

Central differences have error that behaves predictably: D(h) = f'(x) + C·h² + O(h⁴) for some constant C. Richardson extrapolation exploits this directly — compute the same central-difference estimate at two step sizes, h and h/2, and combine them to cancel out the h² term algebraically: R = (4·D(h/2) − D(h)) / 3. What's left over is O(h⁴) — a dramatically better estimate, built entirely from two ordinary central-difference calculations already in hand.

Verified directly — four orders of magnitude of accuracy, for free

At h=0.1: D(h) has error 9.00×10⁻⁴, D(h/2) has error 2.25×10⁻⁴ — but the Richardson combination has error just 1.13×10⁻⁷, nearly four orders of magnitude better than either individual estimate. At h=0.001, the Richardson result's error is 3.5×10⁻¹⁴ — right at the edge of what double precision can represent at all.

Adaptive Quadrature: Spending Effort Where It's Needed

Calculus & Optimization Chapter 9 covered fixed-grid numerical integration — evaluating the function at evenly-spaced points across the whole interval, regardless of how the function actually behaves. Adaptive quadrature instead estimates the local error on each subinterval and only subdivides further where that error is too large — concentrating function evaluations where the function is actually changing quickly, and spending almost none where it's nearly flat.

Verified directly — half the function evaluations, for comparable accuracy

Integrating f(x) = 1/(1 + 10000(x−0.5)²) over [0,1] — a function that's nearly zero almost everywhere except for a sharp, narrow spike right at x=0.5 — against the exact value 0.031015979856434922: a fixed-grid Simpson's rule needs n=500 (501 function evaluations) to reach an error of ≈3.16×10⁻⁹. Adaptive Simpson's rule reaches a comparable error, ≈3.53×10⁻⁹, using only 265 function evaluations — roughly half as many, because most of the fixed grid's points were wasted evaluating the function where it's already almost exactly zero.

Why this matters more as functions get more expensive

The saving here (roughly 2×) is on a function that costs almost nothing to evaluate. For a real-world integrand that might take milliseconds or seconds per evaluation — a physics simulation step, a Monte Carlo sample, a call to an external service — the difference between 265 and 501 evaluations is the difference between a computation finishing twice as fast, for the same accuracy, with no extra implementation cost once the adaptive logic is written once.

Where This Connects

This chapter's finding	What it resolves or sets up
Complex-step differentiation avoids cancellation entirely	A genuinely different resolution to Chapter 4's cancellation problem than the "restructure the algebra" fix used for the quadratic formula
Richardson extrapolation cancels the leading error term algebraically	The same pattern underlies higher-order Runge-Kutta methods and other advanced numerical techniques beyond this course's scope
Adaptive methods concentrate effort where local error is largest	Directly informs how Chapter 10's capstone audit should treat any fixed-step-size code it finds — a specific, checkable red flag

Hands-On Exercises

Exercise 1

Using this chapter's own verified numbers, explain why central differences beat forward differences by nearly five orders of magnitude at h=10⁻⁴, but both still eventually collapse at very small h — what does central differencing fix, and what does it not fix?

📄 View solution

Exercise 2

Using this chapter's own verified complex-step results, explain specifically why the formula Im(f(x+ih))/h has no cancellation error, even though it still involves a division by a small h — your answer should identify what operation is genuinely missing compared to the real-valued finite-difference formulas.

📄 View solution

Exercise 3

A colleague argues that adaptive quadrature is strictly better than fixed-grid integration and should always be used. Using this chapter's own verified example and its own tip box about evaluation cost, describe a realistic situation where the extra implementation complexity of adaptive quadrature might not be worth it.

📄 View solution

Chapter 9 Quick Reference

Central difference (f(x+h)-f(x-h))/(2h): O(h²) error, verified nearly 5 orders of magnitude more accurate than forward differences at h=10⁻⁴ — but still eventually collapses from cancellation at extreme small h
Complex-step differentiation Im(f(x+ih))/h: verified exact to full double precision at every tested h down to 10⁻¹⁰⁰ — completely sidesteps cancellation by never subtracting two real numbers, at the cost of requiring an analytic, complex-compatible function
Richardson extrapolation (4D(h/2)-D(h))/3: verified nearly 4 orders of magnitude accuracy improvement over either individual central-difference estimate, using calculations already computed
Adaptive quadrature: verified roughly half the function evaluations of a fixed grid for comparable accuracy, by concentrating evaluations where a function's local behavior actually demands them
Next chapter: Capstone — auditing a real, ordinary-looking codebase for exactly these nine chapters' worth of numerical bugs