Numerical Calculus

Real scientific problems rarely have closed-form integrals or analytic derivatives. The temperature profile in a reactor, the drag on an airfoil, the cumulative dose of a drug — all are defined by data or by integrals nobody has solved on paper. Numerical calculus turns those problems into arithmetic.

This chapter covers the two pillars: numerical differentiation and numerical integration (quadrature), together with how to think about their accuracy.

The big picture

Differentiation in finite precision is ill-conditioned — you are subtracting nearby numbers. Integration is well-conditioned — errors average out. This asymmetry shapes every algorithm below.

Finite differences

The definition $f'(x) = \lim_{h\to 0} (f(x+h)-f(x))/h$ suggests we just pick a small $h$. In floating point, two things fight:

• Truncation error decreases as $h$ shrinks (the formula gets more accurate)
• Round-off error increases as $h$ shrinks (subtraction of nearby numbers loses digits)

The optimal $h$ is a balance. For the forward difference $f'(x)\approx (f(x+h)-f(x))/h$, the sweet spot is around $h \approx \sqrt{\epsilon_\text{mach}} \approx 10^{-8}$.

Two lessons:

• The centered difference is one order more accurate ($O(h^2)$ vs $O(h)$) for the same cost.
• The optimal $h$ is not "as small as possible". Going below $10^{-8}$ makes things worse.

For most production work, use scipy.differentiate.derivative (modern SciPy) or numpy.gradient for sampled data.

np.gradient uses second-order centered differences in the interior and second-order one-sided at the boundaries. This is the right tool when all you have is a table of values.

Quadrature: integration as a weighted sum

Every integration rule has the form $\int_a^b f(x)\,dx \approx \sum_i w_i\, f(x_i)$ — pick nodes $x_i$ and weights $w_i$, multiply, add. The art is in choosing them.

Notice quad returns both the value and an error estimate. This is the right tool for general one-dimensional integrals: it adaptively refines around peaks and singularities until the estimated error is below the requested tolerance.

Handling singular and oscillatory integrands

quad is robust enough to handle integrable singularities and moderately oscillatory integrands. For truly oscillatory cases, SciPy offers specialized routines like quad with the weight argument set to "sin" or "cos".

Multidimensional integration

For 2-D and 3-D, deterministic rules still work:

But cost explodes: a rule that uses $n$ nodes per axis costs $n^d$ total. Beyond 4-5 dimensions, switch to Monte Carlo, whose error scales as $1/\sqrt{N}$ regardless of dimension (we devote a whole chapter to it later).

A multi-file calculus toolkit

Adaptive Simpson is a beautiful 20-line algorithm: it estimates its own error by comparing one Simpson's rule against the sum of two half-interval Simpson rules, and recurses only where the disagreement exceeds the tolerance. SciPy's quad is a much more sophisticated relative, but the idea is the same.

Challenge: estimate $\pi$ by integration

Check your understanding

QuestionSelect one

You compute $f'(x_0)$ by forward differences and try ever-smaller step sizes $h$. The error first decreases but eventually starts increasing as $h \to 0$. Why?

The CPU clock skew accumulates

Subtracting $f(x_0+h)$ and $f(x_0)$ loses significant digits when the values are nearly equal, so round-off error dominates for tiny $h$

The compiler reorders operations

Python's float type uses arbitrary precision so the rounding only happens at the end

QuestionSelect one

Why does scipy.integrate.quad return both a value and an estimated error bound?

It needs the error to log it

It uses an adaptive algorithm that already computes two different approximations on every sub-interval and uses their difference as a built-in error estimator

Returning the error is required by the IEEE standard

The error is the variance of the integrand