Numerical Differentiation and Integration
Differentiation and integration are fundamental operations in applied mathematics. When analytical forms are unavailable, numerical approximations based on interpolating polynomials and finite differences provide practical alternatives.
Finite Difference Formulas
Forward difference f'(x) ≈ (f(x + h) − f(x))/h has O(h) accuracy. Backward difference f'(x) ≈ (f(x) − f(x − h))/h is similar. Central difference f'(x) ≈ (f(x + h) − f(x − h))/(2h) is O(h2) and usually preferred. Second derivatives use f''(x) ≈ (f(x + h) − 2f(x) + f(x − h))/h2.
Choice of Step Size
Numerical differentiation suffers from a fundamental tension. Smaller h reduces truncation error but increases round-off error from cancellation. An optimal h balances the two; for most central differences, h proportional to the square root of machine epsilon minimizes total error.
Newton–Cotes Integration Rules
The Newton–Cotes rules integrate a polynomial interpolant through equally spaced nodes. Closed rules include the trapezoidal rule, Simpson's 1/3 rule, and Simpson's 3/8 rule; open rules (like the midpoint rule) use interior nodes only.
Trapezoidal Rule
The trapezoidal rule approximates ∫ab f(x) dx ≈ (b − a)/2 × (f(a) + f(b)). Composite trapezoidal rule subdivides [a, b] into n intervals and sums trapezoids; error is O(h2). Simple and robust.
Simpson's Rules
Simpson's 1/3 rule fits a parabola through three points: ∫ab f(x) dx ≈ h/3 × (f(x0) + 4f(x1) + f(x2)). Composite Simpson's has O(h4) accuracy for smooth integrands. Simpson's 3/8 rule uses four points and is useful when the number of intervals is not even.
Romberg Integration
Romberg integration applies Richardson extrapolation to trapezoidal estimates at successively halved step sizes, systematically eliminating leading error terms. It produces very accurate integrals for smooth functions with modest work.
Gaussian Quadrature
Gaussian quadrature chooses both the nodes and the weights to maximize accuracy. An n-point Gauss rule integrates polynomials of degree up to 2n − 1 exactly. Gauss–Legendre quadrature is standard for finite intervals; Gauss–Chebyshev, Gauss–Hermite, and Gauss–Laguerre address weighted or infinite-interval integrals.
Adaptive Quadrature
Adaptive quadrature automatically refines the subdivision of [a, b] where the integrand varies rapidly. The algorithm compares estimates from two rules and subdivides until a tolerance is met. This is the standard method in scientific software.
Multiple Integrals
Multi-dimensional integrals can be computed by tensor products of 1D rules for low dimension, but the cost grows exponentially (curse of dimensionality). In high dimensions, Monte Carlo and quasi-Monte Carlo methods become competitive with costs independent of dimension.
Summary
Numerical differentiation is delicate but useful; numerical integration is robust and widely used. Trapezoidal, Simpson's, Romberg, and Gaussian rules, combined with adaptive strategies, solve nearly all practical one-dimensional integration problems, while Monte Carlo methods handle high-dimensional cases.