6. Numerical Differentiation and Integration

Introduction and Context

Differentiation and integration are central operations in engineering analysis. Differentiation measures local change—such as velocity from displacement, heat flux from temperature gradient, or stress from displacement—whereas integration accumulates a quantity over an interval or region, such as work, mass, energy, volume, or probability. In analytical mathematics these operations are often performed symbolically. Engineering data, however, are commonly available only as discrete measurements, tabulated values, simulation outputs, or experimentally sampled signals. Numerical differentiation and integration provide methods for extracting useful engineering information from such data.

The central difficulty is that numerical methods introduce approximation. Differentiation is particularly sensitive to measurement noise because it amplifies high-frequency errors. Integration is generally more stable, but its accuracy depends on the spacing and quality of data, the smoothness of the function, and the quadrature rule selected. A competent computational engineer must therefore select a method not merely because it is available, but because its accuracy, stability, computational cost, and data requirements suit the physical problem.

This chapter develops finite-difference differentiation, error estimation and Richardson extrapolation, Newton–Cotes quadrature, Gaussian quadrature, multidimensional integration, and adaptive integration. The discussion emphasizes both mathematical formulation and engineering judgement.

Learning Objectives

After studying this chapter, the reader should be able to:

1. Derive and apply forward, backward, and central finite-difference formulas.

2. Distinguish truncation error, round-off error, and data or measurement error.

3. Improve derivative and integral estimates using Richardson extrapolation.

4. Apply composite trapezoidal and Simpson rules to tabulated data.

5. Explain the construction and advantages of Gaussian quadrature.

6. Transform quadrature formulas to arbitrary intervals and multidimensional domains.

7. Design or select adaptive integration procedures.

8. Assess numerical results using error estimates, convergence checks, and physical reasoning.

Discrete Derivatives: From Difference Quotients to Finite Differences

For a differentiable function f(x), the derivative at x is defined by

f'(x) = lim(h→0) [f(x+h) − f(x)]/h.

Replacing the limiting process by a finite spacing h gives the forward-difference approximation

f'(x) ≈ [f(x+h) − f(x)]/h.

Taylor expansion about x gives

f(x+h) = f(x) + hf'(x) + h²f''(x)/2 + h³f'''(x)/6 + ⋯.

Therefore, the forward formula has truncation error O(h), meaning that reducing h by approximately one-half reduces the leading discretization error by one-half, provided round-off and data errors do not dominate.

Using values on the left gives the backward formula,

f'(x) ≈ [f(x) − f(x−h)]/h,

also of order O(h). The central difference uses information on both sides:

f'(x) ≈ [f(x+h) − f(x−h)]/(2h).

Taylor expansions show cancellation of the even-order terms, producing an error O(h²). This usually makes the central difference more accurate than either one-sided formula for the same h. A five-point central formula can achieve fourth-order accuracy:

f'(x) ≈ [f(x−2h) − 8f(x−h) + 8f(x+h) − f(x+2h)]/(12h).

Higher-order formulas require more data points and may be less desirable near boundaries or when the data contain noise.

The second derivative is approximated centrally by

f''(x) ≈ [f(x+h) − 2f(x) + f(x−h)]/h²,

with error O(h²). Although useful in mechanics, heat transfer, and fluid mechanics, this formula is more sensitive to noise than a first-derivative formula because the denominator contains h².

At a boundary, central differences may be unavailable. One-sided formulas are then required. A second-order forward approximation is

f'(x) ≈ [−3f(x) + 4f(x+h) − f(x+2h)]/(2h).

The selection of a formula is therefore influenced by location: central differences are preferred in the interior, while one-sided approximations are necessary at boundaries.

Accuracy, Stability, and the Choice of Step Size

Finite-difference error has several sources. Truncation error results from neglecting higher-order Taylor terms. Round-off error results from finite-precision arithmetic, especially when two nearly equal function values are subtracted. Experimental or simulation uncertainty may be amplified by the derivative operation. For a first derivative, a simplified error model is

E(h) ≈ Ch^p + Dε/h,

where p is the order of the finite-difference formula, ε represents data or machine error, and C and D depend on the function and method. The first term decreases as h decreases, whereas the second increases. Consequently, making h extremely small is not always beneficial.

For noisy measurements, numerical differentiation should generally be preceded by appropriate data processing. Possible approaches include fitting a physically justified polynomial locally, applying a smoothing spline, or using a regression-based derivative estimate. Smoothing must be used cautiously: excessive smoothing can remove real physical features such as shock fronts, resonances, or sharp thermal gradients.

A convergence study is more reliable than selecting h by intuition alone. Compute the derivative using h, h/2, and perhaps h/4. If the results approach a stable value at the expected rate, the calculation is likely in the truncation-error-dominated region. If the estimates fluctuate, round-off or measurement noise may dominate.

Richardson Extrapolation

Richardson extrapolation improves an approximation when its leading error has a known power of h. Suppose an approximation A(h) satisfies

A(h) = A + Ch^p + higher-order terms,

where A is the exact value. Calculating the approximation again with h/2 gives

A(h/2) = A + C(h/2)^p + ⋯.

Eliminating the leading error produces

A_R = [2^p A(h/2) − A(h)]/(2^p − 1).

For the central difference, p = 2. Thus,

A_R = [4A(h/2) − A(h)]/3.

The extrapolated value often has error O(h^(p+1)) or, more precisely, the next noncancelled even power in a symmetric scheme. Richardson extrapolation can also be applied to numerical integration, where it forms the foundation of Romberg integration.

Extrapolation is not a cure for noisy data. If the values A(h) and A(h/2) contain independent measurement errors, the subtraction involved in extrapolation may magnify uncertainty. It is most effective when the function is smooth and discretization error dominates.

Numerical Integration and Quadrature Concepts

For a definite integral,

I = ∫_a^b f(x) dx,

quadrature replaces the continuous function by a weighted combination of function values. A quadrature rule has the form

I ≈ Σ w_i f(x_i),

where x_i are nodes and w_i are weights. The quality of a rule depends on its degree of exactness, convergence order, node placement, smoothness assumptions, and computational cost.

The simplest rule is the trapezoidal rule,

I ≈ (b−a)[f(a)+f(b)]/2.

Its error is proportional to (b−a)^3 f''(ξ), for some ξ in the interval. The composite trapezoidal rule divides [a,b] into n equal panels of width h=(b−a)/n:

I ≈ h[f_0/2 + f_1 + f_2 + ⋯ + f_{n−1} + f_n/2].

It is second-order accurate for smooth functions. It is especially useful for tabulated data because it requires equally spaced values and is simple to implement.

The midpoint rule evaluates the function at the center of each panel and is also second-order accurate. Simpson’s one-third rule fits a quadratic polynomial through three equally spaced points:

∫_{x_0}^{x_2} f(x)dx ≈ (h/3)[f_0 + 4f_1 + f_2].

The composite Simpson rule is

I ≈ (h/3)[f_0 + f_n + 4(f_1+f_3+⋯+f_{n−1}) + 2(f_2+f_4+⋯+f_{n−2})],

where n must be even. For smooth functions its error is O(h^4), making it substantially more accurate than the composite trapezoidal rule at comparable spacing. Simpson’s rule is exact for polynomials of degree three or less, although its practical performance depends on smoothness and data quality.

Newton–Cotes rules are obtained by interpolating f(x) at equally spaced nodes and integrating the interpolating polynomial. Closed rules include both endpoints; open rules omit them. Increasing the degree does not guarantee improved numerical behaviour. High-degree Newton–Cotes formulas can suffer from oscillatory interpolation and large alternating weights. Composite low-order rules are usually safer than a single high-order rule over a wide interval.

Gaussian Quadrature

Gaussian quadrature differs fundamentally from Newton–Cotes integration because both nodes and weights are selected to maximize accuracy. An n-point Gauss–Legendre rule on [−1,1] integrates every polynomial of degree up to 2n−1 exactly:

∫_{−1}^{1} f(ξ)dξ ≈ Σ_{i=1}^n w_i f(ξ_i).

The nodes ξ_i are the roots of the Legendre polynomial P_n(ξ). For two-point Gauss–Legendre quadrature,

ξ_1 = −1/√3, ξ_2 = 1/√3,

and both weights equal 1. Hence,

∫_{−1}^{1} f(ξ)dξ ≈ f(−1/√3)+f(1/√3).

To transform the rule to [a,b], use

x = (a+b)/2 + (b−a)ξ/2,

dx = (b−a)dξ/2.

Therefore,

∫_a^b f(x)dx ≈ (b−a)/2 Σ w_i f((a+b)/2 +(b−a)ξ_i/2).

Gaussian quadrature is highly efficient when function evaluations are expensive, as in finite-element material models or computational fluid dynamics. It does not require function values at the interval endpoints, which can be advantageous when the integrand is singular or difficult to evaluate there. However, a standard Gauss rule is not automatically appropriate for singular, discontinuous, or sharply peaked integrands; specialized weighted rules or interval subdivision may be required.

Multidimensional Integration

Many engineering quantities are defined over areas, volumes, or parameter spaces. A rectangular-domain integral can be approximated by applying one-dimensional quadrature successively:

I = ∫_a^b ∫_c^d f(x,y)dy dx.

For a tensor-product Gaussian rule,

I ≈ Σ_i Σ_j w_i w_j f(x_i,y_j),

after mapping each coordinate to its physical interval. In two dimensions, an n-point rule in each direction requires n² function evaluations; in three dimensions it requires n³. This rapid growth is called the curse of dimensionality.

For nonrectangular domains, one may divide the region into subregions, transform the region to a reference element, or use triangular and tetrahedral quadrature rules. In finite-element analysis, numerical integration on reference elements is fundamental for constructing stiffness matrices, mass matrices, and internal force vectors. The Jacobian of the coordinate transformation must be included:

∫_Ω f(x,y)dΩ = ∫_{\hat{Ω}} f(x(ξ,η),y(ξ,η)) |J| dξdη.

A negative or nearly zero Jacobian signals an invalid or severely distorted element and should not be treated as a mere quadrature issue.

Adaptive Integration

Uniform spacing is inefficient when the integrand is smooth in one region and rapidly varying in another. Adaptive integration addresses this by estimating the error locally and subdividing only those intervals requiring greater resolution.

A common procedure is adaptive Simpson integration. First compute Simpson’s estimate S(a,b) over the full interval. Then compute S(a,m) and S(m,b), where m=(a+b)/2. The difference

|S(a,m)+S(m,b)−S(a,b)|

provides an error indicator. If the indicator is below a prescribed tolerance, accept the refined estimate; otherwise, recursively subdivide one or both subintervals. A corrected estimate often uses

S_adaptive = S(a,m)+S(m,b) + [S(a,m)+S(m,b)−S(a,b)]/15.

A robust implementation also imposes a maximum recursion depth, avoids repeated function evaluations by reusing previously computed values, and distinguishes absolute from relative tolerance. Absolute tolerance is important when the integral is near zero; relative tolerance is useful when its magnitude is large.

Adaptive methods are valuable for boundary layers, localized peaks, rapidly changing material laws, and integrands with piecewise smooth behaviour. They can perform poorly when discontinuities are not detected, when the function is highly oscillatory, or when the error estimator is misleading. Known discontinuities should normally be supplied as subdivision points rather than left for the algorithm to discover.

Worked Engineering Examples

Example 1: Estimating a velocity from displacement data

Suppose displacement is represented by s(t)=t³−2t²+3t, and the velocity is required at t=2. The exact velocity is s'(t)=3t²−4t+3, so v(2)=7. With h=0.1, the central difference gives

[s(2.1)−s(1.9)]/0.2,

which is close to 7 and has O(h²) error. A forward difference uses [s(2.1)−s(2)]/0.1 and has O(h) error. Repeating the central calculation with h=0.05 and applying Richardson extrapolation gives a still more accurate result, provided the displacement values are exact or sufficiently precise.

Example 2: Heat absorbed during a transient process

If a measured heat-transfer rate q(t) is available at equally spaced times, the accumulated heat is Q=∫q(t)dt. The composite trapezoidal rule is appropriate when the data are tabulated and may not conform exactly to a simple polynomial. If q(t) varies smoothly and an even number of panels is available, composite Simpson integration can provide higher accuracy. A comparison of both estimates, together with halving the time step where possible, provides practical evidence of discretization error.

Example 3: Work over a displacement path

For a force F(x), the work W=∫_{x_0}^{x_1}F(x)dx. If F is evaluated through an expensive nonlinear constitutive model, a two- or four-point Gaussian rule may achieve the required accuracy with fewer evaluations than a uniformly spaced rule. If the force changes abruptly because of contact or frictional transitions, the path should be divided at those transitions and quadrature applied separately.

Illustrative Case Study: Adaptive Integration in a Thermal Component

Background

Consider an illustrative electronic cooling component whose heat-generation rate along a one-dimensional path is modeled by q(x). Most of the component produces heat gradually, but a narrow region near a resistance hotspot produces a sharp peak. The total heat generation is required for a thermal balance calculation.

Problem

A coarse uniform trapezoidal rule underestimates the peak because no sampling point lies close to its maximum. Increasing the number of points over the whole component would improve the estimate but would also increase computational cost.

Analysis and Approach

The engineer first evaluates q(x) at the endpoints and midpoint and forms a Simpson estimate. The interval is then bisected, and the two subinterval estimates are compared with the original estimate. The region containing the hotspot produces a large discrepancy and is recursively subdivided. Smooth regions satisfy the tolerance after only one or two subdivisions. Previously calculated values at shared interval boundaries are reused.

The calculation is repeated with a tighter tolerance. If the total heat changes negligibly, the result is considered converged. The engineer also compares the result against a high-resolution reference calculation and checks that the integrated heat is nonnegative and consistent with independent power estimates.

Outcome

Adaptive integration concentrates function evaluations around the hotspot and uses relatively few evaluations in smooth regions. It provides a more reliable result than a coarse uniform rule and a lower computational cost than applying the finest spacing everywhere.

Learning Points

The case illustrates that quadrature selection is a modelling decision, not merely an algebraic exercise. Local error behaviour matters; known physical transitions should guide interval subdivision; and convergence and physical checks are essential before an integral is used in a design decision.

Limitations, Risks, and Common Mistakes

A frequent mistake is to assume that a smaller step always improves a derivative. Subtractive cancellation and measurement noise may make the result worse. Another is to apply a central formula at a boundary without recognizing that the required data do not exist. For integration, common errors include using Simpson’s rule with an odd number of panels, omitting the Jacobian in a coordinate transformation, and confusing the weights or nodes of a Gaussian rule after mapping intervals.

Numerical differentiation of noisy experimental data should not be reported with unjustified precision. Integration across a discontinuity should be split at the discontinuity. Oscillatory functions require a spacing that resolves the oscillation; otherwise, aliasing can produce a plausible but incorrect result. In multidimensional problems, tensor-product rules may become computationally prohibitive, motivating sparse grids, quasi-Monte Carlo methods, or problem-specific dimensional reduction.

Professional and Computational Practice

A reliable implementation should document the formula, step size, node locations, tolerance, stopping criterion, and estimated error. It should test the method on functions with known derivatives or integrals and include limiting cases. Software libraries are useful, but their output should be interpreted rather than accepted blindly. In safety-critical engineering, independent calculations, dimensional checks, bounds, and reproducibility of numerical settings are important forms of computational governance.

Automatic differentiation can compute derivatives of differentiable computer programs without finite-difference step-size selection. It is valuable in optimization and machine learning, but it does not replace numerical quadrature and may not handle discontinuous logic or externally measured data without special treatment. Modern engineering workflows frequently combine automatic differentiation, adaptive quadrature, surrogate models, and uncertainty quantification. The fundamental requirements remain unchanged: understand the mathematical assumptions, monitor error, and relate numerical output to physical behaviour.

Best Practices

Use central differences in smooth interior regions and one-sided formulas at boundaries. Select the step size through convergence studies rather than arbitrary rules. Use Richardson extrapolation only when the expected error structure is valid. Prefer composite low-order Newton–Cotes rules for general tabulated data and Gaussian quadrature for smooth, expensive-to-evaluate functions. Subdivide at discontinuities and use adaptive methods for localized variation. In multidimensional integration, include the Jacobian and inspect element quality. Finally, report estimates with precision justified by numerical and physical uncertainty.

Chapter Summary

Numerical differentiation estimates local rates of change from discrete values, with forward and backward formulas of first order and central formulas commonly of second order. Higher-order formulas and Richardson extrapolation can improve accuracy, but noise and round-off impose practical limits. Numerical integration replaces an integral by weighted function evaluations. Composite trapezoidal and Simpson rules are convenient for tabulated data, while Gaussian quadrature achieves high polynomial exactness with optimally selected nodes and weights. Multidimensional integration requires tensor products or reference-element transformations and must include the Jacobian. Adaptive integration allocates computational effort according to local error indicators. Across all methods, convergence studies, error estimates, physical checks, and awareness of data quality are essential for credible engineering computation.

Key Terms

Finite difference: An approximation to a derivative formed from function values at discrete points.

Forward difference: A one-sided derivative approximation using the current and a succeeding point.

Backward difference: A one-sided derivative approximation using the current and a preceding point.

Central difference: A derivative approximation using points symmetrically located around the evaluation point.

Truncation error: Error caused by neglecting higher-order terms in a numerical approximation.

Round-off error: Error arising from finite-precision arithmetic.

Richardson extrapolation: An error-cancellation technique using approximations at two step sizes.

Quadrature: Numerical approximation of a definite integral.

Newton–Cotes rule: A quadrature rule based on interpolation at equally spaced points.

Trapezoidal rule: A quadrature method that approximates the integrand by straight-line segments.

Simpson’s rule: A quadrature method based on quadratic interpolation over pairs of panels.

Gaussian quadrature: A rule with specially selected nodes and weights that maximizes polynomial exactness.

Node: A point at which the integrand is evaluated.

Weight: A coefficient multiplying a function value in a quadrature sum.

Degree of exactness: The highest polynomial degree integrated exactly by a quadrature rule.

Adaptive integration: Integration that refines intervals according to estimated local error.

Tensor-product quadrature: Multidimensional quadrature formed by combining one-dimensional rules.

Jacobian: The determinant accounting for local scaling in a coordinate transformation.

Convergence study: Repeated computation at refined resolution to assess numerical stability.

Curse of dimensionality: Rapid growth in computational cost as the number of dimensions increases.

Review Questions

Short Answer Questions

1. Why is a central-difference formula generally more accurate than a forward-difference formula of the same spacing?

2. Distinguish truncation error from round-off error.

3. State the condition on the number of panels for composite Simpson’s one-third rule.

4. What is meant by the degree of exactness of a quadrature rule?

5. Why must a Jacobian be included in a transformed multidimensional integral?

Long Answer Questions

1. Derive the central-difference approximation for the first derivative and identify its leading error order.

2. Explain Richardson extrapolation and discuss the circumstances in which it may fail to improve a result.

3. Compare composite trapezoidal, composite Simpson, and Gaussian quadrature for engineering applications.

4. Describe adaptive Simpson integration, including its error indicator and stopping strategy.

5. Explain how multidimensional Gaussian quadrature is constructed on a rectangular and a transformed domain.

Analytical / Application Questions

1. A temperature signal is measured at equal time intervals but contains moderate random noise. Which derivative-estimation strategy would you use, and how would you validate its result?

2. A force–displacement curve contains a sharp change caused by contact. Select an integration strategy and explain how the transition should be treated.

3. An engineer reports that reducing the finite-difference step from 10⁻³ to 10⁻⁶ made the derivative less stable. Explain the likely cause and propose a diagnostic procedure.

22. Multiple Choice Questions with Answers

Question 1

Which finite-difference approximation has truncation error of order O(h²) for a smooth function?

A. Forward difference

B. Backward difference

C. Central difference

D. First-order one-sided difference

Correct Answer: C. Central difference

Question 2

For a central-difference approximation with leading error Ch², Richardson extrapolation using h and h/2 is given by:

A. [A(h/2)−A(h)]/2

B. [4A(h/2)−A(h)]/3

C. [A(h)+A(h/2)]/2

D. [2A(h)−A(h/2)]/3

Correct Answer: B. [4A(h/2)−A(h)]/3

Question 3

What is the principal reason that an extremely small finite-difference step can produce an unstable derivative?

A. The function becomes discontinuous

B. Gaussian weights become invalid

C. Subtractive cancellation and round-off error increase

D. The derivative necessarily becomes zero

Correct Answer: C. Subtractive cancellation and round-off error increase

Question 4

Composite Simpson’s one-third rule requires:

A. An odd number of panels

B. An even number of panels

C. Exactly two panels

D. Unequally spaced nodes only

Correct Answer: B. An even number of panels

Question 5

An n-point Gauss–Legendre rule on [−1,1] integrates exactly every polynomial of degree up to:

A. n

B. n−1

C. 2n−1

D. 2n

Correct Answer: C. 2n−1

Question 6

When transforming a Gaussian quadrature rule from [−1,1] to [a,b], the integral must be multiplied by:

A. a+b

B. b−a

C. (b−a)/2

D. 2/(b−a)

Correct Answer: C. (b−a)/2

Question 7

In a two-dimensional tensor-product quadrature rule, using n nodes in each coordinate direction requires how many function evaluations?

A. 2n

B. n+2

C. n²

D. 2ⁿ

Correct Answer: C. n²

Question 8

Adaptive integration primarily improves efficiency by:

A. Using the smallest possible step everywhere

B. Refining intervals where the estimated local error is large

C. Ignoring discontinuities

D. Replacing all function values with endpoint values

Correct Answer: B. Refining intervals where the estimated local error is large

Question 9

If an integral is evaluated after a coordinate transformation, omission of the Jacobian generally causes:

A. A change in the integration domain’s scaling being ignored

B. Automatic error cancellation

C. The integrand to become analytic

D. The quadrature order to double

Correct Answer: A. A change in the integration domain’s scaling being ignored

Question 10

Which situation most strongly favours splitting an integration interval before applying quadrature?

A. A smooth polynomial integrand

B. A known discontinuity or sharp physical transition

C. A constant integrand

D. A function whose value is zero everywhere

Correct Answer: B. A known discontinuity or sharp physical transition