Numerical Differentiation & Integration

Compute derivatives and integrals from data or formulas without antiderivatives. Learn forward, backward and central differences and their error orders, the trapezoidal rule, Simpson's 1/3 and 3/8 rules, error behaviour, Richardson extrapolation and Gaussian quadrature — with a Python comparison.

Velocity from GPS positions, heat flux from thermocouple readings, energy from sampled power, and areas of irregular sections all need numerical calculus. FEM stiffness matrices are built with Gaussian quadrature.

An analogy

Numerical integration is like estimating the area of a lake by laying down trapezoids or curved strips and adding their areas.

Finite differences

Derived from Taylor series. Central differences are more accurate for the same h. Differentiating noisy data amplifies noise — smooth or fit first.

Newton–Cotes rules

Trapezoidal (straight lines), Simpson 1/3 (parabolas, even panels), Simpson 3/8 (cubics, panels in multiples of 3). Composite rules split the interval into many panels.

Error control

Halving h reduces trapezoidal error 4× and Simpson 16×. Richardson extrapolation combines two step sizes to cancel the leading error (Romberg integration). Adaptive quadrature refines where needed.

Gaussian quadrature

Chooses optimal points and weights: n points integrate polynomials of degree 2n − 1 exactly. Standard in FEM element integration.

Formulas

Worked example: Central vs forward difference

Given: f = sin x at x = 1, h = 0.1

Find: Errors.

Answer: Central ≈ 48× more accurate

Worked example: Gauss 2-point

Given: ∫_{−1}^{1} x² dx (exact 2/3)

Find: Gauss estimate.

Answer: 2/3 exactly

Trapezoidal vs Simpson for ∫₀¹ e^(−x²) dx

import math

f = lambda x: math.exp(-x * x)          # no elementary antiderivative
a, b, n = 0.0, 1.0, 10
h = (b - a) / n
x = [a + i * h for i in range(n + 1)]
y = [f(xi) for xi in x]

trap = h * (y[0] / 2 + sum(y[1:-1]) + y[-1] / 2)
simp = h / 3 * (y[0] + 4 * sum(y[1:-1:2]) + 2 * sum(y[2:-1:2]) + y[-1])
exact = math.sqrt(math.pi) / 2 * math.erf(1)

print(f"trapezoidal: {trap:.10f}  error {abs(trap - exact):.2e}")
print(f"simpson 1/3: {simp:.10f}  error {abs(simp - exact):.2e}")
print(f"exact:       {exact:.10f}")
# trapezoidal: 0.7462107961  error 6.13e-04
# simpson 1/3: 0.7468249483  error 8.15e-07
# exact:       0.7468241328

With the same 11 function values, Simpson's rule is ~750× more accurate because it fits parabolas instead of straight lines (O(h⁴) vs O(h²)).

Common mistakes

Applications

Key takeaways

Going further