Chapter 7 develops the theory and practice of approximating definite integrals using finite computational methods. Students who have studied calculus know that the definite integral of a function represents accumulated change, area under a curve, total work, total probability, or any quantity built up continuously over an interval. They also know, from the Fundamental Theorem of Calculus, that definite integrals can sometimes be evaluated exactly by finding antiderivatives. But many integrals that arise in science, engineering, finance, and data analysis cannot be evaluated symbolically — either because no closed-form antiderivative exists, because the function is known only from measured data, or because exact evaluation would be impractical. Numerical integration, sometimes called numerical quadrature, replaces the continuous accumulation process with a structured finite sum that approximates the integral to within a controllable error.

This chapter introduces the major families of numerical integration rules: the rectangle rules, the midpoint rule, the trapezoidal rule, and Simpson's rule. It develops composite versions of these rules for greater accuracy, derives error bounds so students can judge reliability, introduces adaptive integration as an idea, and applies numerical integration to problems from physics, finance, engineering, and data science. Throughout, the chapter connects numerical integration to the Riemann sum interpretation of the definite integral from calculus and to the error analysis framework from Chapter 2.

The ability to evaluate integrals numerically is one of the most broadly useful tools in applied mathematics. Probability distributions that cannot be integrated symbolically — including the normal distribution — are evaluated numerically. Engineering simulations compute forces, energies, and fluxes by numerical quadrature. Financial models compute option prices, expected values, and risk measures using integrals over probability distributions. Signal processing computes transforms and convolutions numerically. Scientific computing packages from astrophysics to climate modeling rely on efficient, reliable numerical integration at every step.

Numerical integration is also foundational to the numerical solution of differential equations. Euler's method and Runge-Kutta methods, which students will study in Chapter 12, are built from the same ideas as numerical integration rules. Understanding how and why numerical integration works prepares students for the much more complex algorithms that appear in scientific computing.

The Fundamental Theorem of Calculus gives an elegant formula for the definite integral:

where \( F \) is any antiderivative of \( f \). This formula is powerful, but it requires finding \( F \) in closed form. Not every function has a closed-form antiderivative. Some functions are defined only by data. Others have antiderivatives that exist theoretically but cannot be expressed in terms of elementary functions.

Consider the following examples.

In all of these cases, numerical integration provides a way to approximate the integral to any desired accuracy, given sufficient function evaluations.

Before developing the approximation methods, it is worth reviewing the geometric and physical meaning of the definite integral, because these meanings guide the intuition behind the approximation rules.

The definite integral \( \int_a^b f(x)\, dx \) can be interpreted in at least three ways that support numerical thinking.

When \( f(x) \geq 0 \), the definite integral equals the area of the region bounded by the curve \( y = f(x) \), the \( x \)-axis, and the vertical lines \( x = a \) and \( x = b \). When \( f \) changes sign, the integral gives the signed area, counting regions below the axis as negative.

The integral accumulates the continuous quantity \( f(x)\,dx \) over the interval \( [a, b] \). If \( f(x) \) represents a rate of change — velocity, power, density, population growth rate — then \( \int_a^b f(x)\, dx \) gives the total accumulated quantity: displacement, energy, mass, total population change.

The definite integral is defined as the limit of Riemann sums. Partition \( [a, b] \) into \( n \) subintervals of equal width \( h = (b-a)/n \). Choose a sample point \( x_i^* \) in each subinterval \( [x_{i-1}, x_i] \). Then

Numerical integration methods are finite Riemann sums. They stop before taking the limit, using a finite number of subintervals. The key insight is that for smooth functions, a well-chosen finite sum gives an excellent approximation.

The simplest numerical integration rules approximate the integrand by a constant on each subinterval. This is the same as using a single-term Riemann sum on each piece.

Partition \( [a, b] \) into \( n \) subintervals of equal width:

The partition points are \( x_0 = a \), \( x_1 = a + h \), \( x_2 = a + 2h \), \( \ldots \), \( x_n = b \).

On each subinterval \( [x_{i-1}, x_i] \), approximate the integrand by the constant \( f(x_{i-1}) \), the function value at the left endpoint. This gives the rectangle with base \( h \) and height \( f(x_{i-1}) \). Summing over all \( n \) subintervals:

On each subinterval \( [x_{i-1}, x_i] \), approximate the integrand by the constant \( f(x_i) \), the value at the right endpoint:

For an increasing function, the left rule underestimates because each rectangle lies below the curve. The right rule overestimates for the same reason. For a decreasing function, the roles reverse. For functions that neither consistently increase nor decrease, the errors partially cancel.

The midpoint rule improves on the rectangle rules by evaluating the integrand at the midpoint of each subinterval rather than at the endpoints. The midpoint value better represents the average of the integrand over a smooth subinterval.

Let \( m_i = (x_{i-1}+x_i)/2 = x_{i-1} + h/2 \) denote the midpoint of the \( i \)-th subinterval. The midpoint rule is:

The midpoint rule is a second-order method: its error decreases proportionally to \( h^2 \). This is a substantial improvement over the \( O(h) \) convergence of the left and right rules. The reason is geometric: the midpoint approximation happens to match a tangent-line approximation at the midpoint, and by a symmetry argument, the overestimate on one side of the midpoint cancels the underestimate on the other side more effectively than either endpoint.

This can be understood from Taylor series. If \( f \) is smooth over \( [x_{i-1}, x_i] \), then

The error on each subinterval is proportional to \( h^3 \). Summing over \( n = (b-a)/h \) subintervals gives a global error proportional to \( h^2 \).

The trapezoidal rule replaces the constant approximation of the rectangle rules with a linear approximation on each subinterval. On each subinterval \( [x_{i-1}, x_i] \), it draws the trapezoid whose two parallel sides have heights \( f(x_{i-1}) \) and \( f(x_i) \).

On a single interval \( [a, b] \) with \( n = 1 \):

This is the area of the trapezoid with base \( (b-a) \) and two heights \( f(a) \) and \( f(b) \).

For \( n \) subintervals of width \( h = (b-a)/n \):

Notice the pattern: the endpoints \( x_0 \) and \( x_n \) each receive a weight of 1, while all interior points receive a weight of 2. This formula arises naturally by applying the single-interval trapezoidal rule to each subinterval and summing:

Adjacent subintervals share an interior point, so interior function values appear in two trapezoids and are counted twice.

The trapezoidal rule is also a second-order method: its error decreases proportionally to \( h^2 \), the same rate as the midpoint rule. For a fixed \( n \), however, the midpoint rule is typically about twice as accurate as the trapezoidal rule for smooth functions (see Section 7.8 for error bounds). The trapezoidal rule has the practical advantage that it requires only the partition points, making it ideal when data are available only at equally spaced points.

Simpson's rule achieves higher accuracy than the trapezoidal rule by replacing the linear approximation on each pair of subintervals with a quadratic approximation. Instead of connecting two consecutive points by a line, it fits a parabola through three consecutive points.

On the interval \( [x_0, x_2] \) with midpoint \( x_1 = x_0 + h \), fit a parabola \( p(x) \) through the three points \( (x_0, f(x_0)) \), \( (x_1, f(x_1)) \), \( (x_2, f(x_2)) \). The integral of this parabola over \( [x_0, x_2] \) equals:

This formula is Simpson's rule applied to a single pair of subintervals. The weight pattern is \( 1, 4, 1 \), which can be remembered as the weights assigned to the left endpoint, midpoint, and right endpoint of each double-subinterval.

For \( n \) subintervals (where \( n \) must be even) of width \( h = (b-a)/n \):

The weight pattern across the full interval is: \( 1, 4, 2, 4, 2, 4, \ldots, 4, 2, 4, 1 \). The endpoints receive weight 1. Interior points at odd indices receive weight 4. Interior points at even indices receive weight 2.

Simpson's rule is a fourth-order method: its error decreases proportionally to \( h^4 \). Halving the step size reduces the error by a factor of approximately \( 2^4 = 16 \). This is a dramatic improvement over the second-order trapezoidal and midpoint rules.

The high accuracy of Simpson's rule comes from the fact that a parabola exactly integrates any cubic polynomial (not just any quadratic). This is because the \( 1, 4, 1 \) weighting pattern happens to cancel the third-order error term as well as the second-order one.

The methods in the previous sections are most powerful when applied as composite rules: the same basic rule is applied repeatedly over \( n \) small subintervals. This section makes the structure of composite rules explicit and discusses how to choose \( n \).

The key observation is:

If you have a good rule for approximating the integral over a single small interval, applying it to every subinterval and summing the results gives a good approximation for the whole integral. The more subintervals, the smaller the error — but also the more function evaluations required.

For \( n \) subintervals of width \( h = (b-a)/n \):

Choosing \( n \) requires balancing accuracy against computational cost. More subintervals give more accuracy but require more function evaluations. For hand computation, \( n \) between 4 and 20 is typical. For computer computation, \( n \) may be hundreds or thousands.

Every numerical integration method introduces an error called the truncation error: the error that would remain even if arithmetic were performed exactly. Understanding and bounding this error is essential for judging the reliability of an approximation.

There are two sources of error in numerical integration. Truncation error arises from replacing the true integrand with an approximating function (a constant, line, or parabola). This error depends on how smooth the integrand is and how small the subintervals are. Roundoff error arises from the finite precision of arithmetic. In practice, roundoff error is usually small relative to truncation error unless \( n \) is extremely large. For most applications, truncation error dominates, and the analysis below focuses on it.

For a function \( f \) with bounded derivatives on \( [a, b] \):

Error bounds tell you the worst-case error. In practice, the actual error is usually much smaller. The bounds are most useful for: (1) estimating \( n \) needed to achieve a given accuracy before computing; (2) certifying that the computed approximation is within a required tolerance; (3) comparing methods to understand which is more efficient.

This section consolidates the error bound formulas and discusses their practical interpretation.

All of the methods developed so far use equally spaced subintervals. But many integrands behave very differently in different parts of the interval: nearly constant in some regions, rapidly changing in others. Using equally spaced subintervals wastes computational effort where the integrand is smooth and risks inaccuracy where it changes rapidly.

Adaptive integration allocates subintervals based on the behavior of the integrand. Where the integrand is smooth and nearly linear or quadratic, fewer subintervals suffice. Where the integrand is rapidly changing or nearly singular, more subintervals are needed.

Most adaptive integration algorithms work as follows:

1. Apply a basic quadrature rule (often Simpson's rule or a related pair of rules) to a given subinterval.
2. Estimate the error on that subinterval by comparing two approximations of different order or step size.
3. If the error estimate is below the desired tolerance, accept the approximation for that subinterval.
4. If the error estimate exceeds the tolerance, split the subinterval in half and apply the procedure recursively to each half.

This process continues until the estimated error over every subinterval is below the tolerance. The total integral is the sum of the accepted approximations.

Adaptive integration is the strategy used by most practical scientific computing libraries. Functions like quad in Python's SciPy, integral in MATLAB, and analogous functions in R, Julia, and Mathematica all use adaptive strategies. When you call such a function, the algorithm is automatically concentrating computational effort where the integrand is hardest to resolve.

Adaptive integration is particularly valuable for integrands with near-singularities or sharp peaks; integrands that are nearly constant in most of the interval but rapidly changing in a small region; integrals where the desired accuracy is not known in advance; and situations where function evaluations are inexpensive relative to the cost of programming decisions.

In many applications, the integrand is not given by a formula. It is measured at discrete points: experimental data, sensor readings, financial returns, weather observations. Numerical integration of tabular data requires adapting the rules developed above to situations where the function is known only at the given data points.

When data are given at equally spaced points \( x_0, x_1, \ldots, x_n \) with spacing \( h \), the trapezoidal rule and Simpson's rule apply directly. The trapezoidal rule requires function values at all \( n+1 \) partition points. Since the data are already at these points, no additional function evaluations are needed.

Simpson's rule requires an even number of subintervals (an odd number of data points). If the data have an even number of points, apply Simpson's rule to all but the last pair and use the trapezoidal rule on the final interval — or apply a modified endpoint correction.

If the data are not equally spaced, the standard composite formulas do not apply directly. Options include:

1. Apply the trapezoidal rule on each subinterval with its own width: \( T = \sum_{i=1}^n \frac{x_i - x_{i-1}}{2}[f(x_{i-1}) + f(x_i)] \).
2. Fit an interpolating polynomial or spline through the data and integrate the interpolant (see Chapter 4 and Chapter 5).
3. Use a specialized rule for non-uniform data.

Numerical integration appears throughout quantitative science. The following examples illustrate its breadth.

In mechanics, the work done by a force \( F(x) \) acting over a displacement from \( a \) to \( b \) is

If \( F(x) \) is measured experimentally (for example, from a strain gauge), the integral must be approximated numerically. In thermodynamics, the heat absorbed or released by a substance as its temperature changes from \( T_1 \) to \( T_2 \) is

where \( c_p(T) \) is the specific heat capacity, often given by empirical data.

The present value of a continuous income stream \( r(t) \) over a period \( [0, T] \) at discount rate \( \delta \) is

When \( r(t) \) depends on economic conditions and is not given by a simple formula, this integral must be computed numerically.

In the Black-Scholes framework for option pricing, the price of a European call option requires evaluating the cumulative normal distribution:

Since \( e^{-t^2/2} \) has no closed-form antiderivative, this integral is approximated numerically to high precision using composite Simpson's rule or adaptive quadrature.

In signal processing, the energy of a signal \( x(t) \) over a time window \( [0, T] \) is

Measured signals are discrete sequences, so the integral is approximated by the trapezoidal rule or Simpson's rule over the sample points.

In feedback control systems, the integral of the error signal determines the "I" term in a PID controller. Numerical integration of the error at each time step is essential for real-time control.

In probability and statistics, the probability that a random variable \( X \) with density \( f(x) \) falls between \( a \) and \( b \) is

For distributions like the normal, beta, gamma, or Student's \( t \), these integrals must be approximated numerically. Statistical tables, confidence intervals, hypothesis test p-values, and credible intervals in Bayesian analysis all depend on numerical integration of probability densities.

Students learning numerical integration frequently make errors of two kinds: computational errors in applying the formulas, and conceptual errors in interpreting or choosing methods. The following list addresses the most common problems.

Chapter 7 closes the Part III development of numerical calculus. Students have now studied numerical differentiation (Chapter 6) and numerical integration (Chapter 7). These two tools are the computational foundation for the most important application of numerical calculus: solving differential equations.

The initial value problem \( y' = f(t, y) \), \( y(t_0) = y_0 \) can be rewritten in integral form:

To advance from \( t_n \) to \( t_{n+1} = t_n + h \), a numerical ODE solver approximates:

Euler's method uses the rectangle rule to approximate this integral:

The improved Euler method and Runge-Kutta methods use higher-order approximations — analogous to the midpoint rule and Simpson's rule — to achieve better accuracy per step. The entire theory of ODE solver accuracy, step size selection, and stability grows directly from the numerical integration ideas developed in this chapter.

Chapter 8 will develop Taylor polynomials and series as another family of approximation tools. Taylor-based ODE methods explicitly use derivative information to achieve high-order accuracy. These methods connect the local polynomial approximation ideas of Chapters 6–8 to the ODE solver framework of Chapter 12.

The local truncation error and global error of ODE solvers parallel the subinterval error and total error of composite integration rules. The stability concepts developed for ODE solvers are direct extensions of the stability ideas introduced in Chapter 2.

Students who understand how and why composite trapezoidal and Simpson's rules work — and how to bound and control their errors — are well prepared for the much richer theory of ODE solvers.

This chapter introduced numerical integration: the approximation of definite integrals by finite weighted sums of function values.

Numerical integration is essential throughout physics, engineering, finance, statistics, and data science, and it provides the conceptual foundation for the numerical ODE methods studied in Chapter 12.