Numerical Integration - Theory & Concepts

Numerical Integration

Learning Objectives

Describe numerical integration (quadrature) and its geometric meaning.
Apply Newton-Cotes formulas, including trapezoidal and Simpson's rules.
Differentiate between open and closed forms of integration methods.
Compare computational methods for improved accuracy such as composite rules, Romberg integration, and adaptive quadrature.
Utilize Gauss quadrature to integrate functions more efficiently.
Evaluate integrals using Monte Carlo integration for high-dimensional problems.

Numerical integration, also known as quadrature, involves approximating the definite integral of a mathematical function or discrete data. Geometrically, it represents finding the area under a continuous curve.

Numerical Integration (Quadrature)

A broad family of algorithms for calculating the numerical value of a definite integral.

Newton-Cotes Formulas

Newton-Cotes Formulas

A group of numerical integration formulas based on evaluating the integrand at equally spaced points.

The Newton-Cotes formulas are the most common numerical integration schemes. They are based on the strategy of replacing a complicated function or tabulated data with an approximating interpolating polynomial that is easy to integrate analytically.

Open and Closed Forms

Newton-Cotes formulas have two forms: closed forms, where data points at the beginning and end of the limits of integration ( $a$ and $b$ ) are known and utilized, and open forms, where integration limits extend beyond the range of the data points. Open forms are highly useful for improper integrals or for solving initial-value ordinary differential equations.

Trapezoidal Rule (Closed Form)

The simplest closed Newton-Cotes formula. It approximates the area under the curve using a straight line (a first-order polynomial) connecting the two points at the ends of the interval.

Trapezoidal Rule

Approximates the integral of a function using a single trapezoid between two points.

I \approx (b - a)\frac{f(a) + f(b)}{2}

Variables

Symbol	Description	Unit
$I$	Approximate integral value	-
$a$	Lower limit of integration	-
$b$	Upper limit of integration	-
$f(a)$	Function value at lower limit	-
$f(b)$	Function value at upper limit	-

The true error for a single application is $E_t = -\frac{1}{12} f''(\xi)(b-a)^3$ , indicating the rule is exact for linear functions (where $f''=0$ ). For improved accuracy, the multiple-application (or composite) trapezoidal rule divides the interval into $n$ multiple segments of width $h = \frac{b-a}{n}$ and applies the rule to each.

Composite Trapezoidal Rule

Approximates the integral by dividing the domain into multiple segments and summing their individual trapezoidal areas.

I \approx \frac{h}{2} \left[ f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n) \right]

Variables

Symbol	Description	Unit
$I$	Total approximate integral	-
$h$	Segment width ( $\frac{b-a}{n}$ )	-
$x_0$	Lower limit of integration ( $a$ )	-
$x_n$	Upper limit of integration ( $b$ )	-
$x_i$	Intermediate points	-
$n$	Number of segments	-

Global Error Bound for Composite Trapezoidal Rule

Determines the upper bound of the global truncation error.

|E_a| \le \frac{(b-a)^3}{12n^2} \max_{x \in [a,b]} |f''(x)| = \frac{(b-a)h^2}{12} \max_{x \in [a,b]} |f''(x)|

Variables

Symbol	Description	Unit
$E_a$	Global truncation error	-
$a$	Lower limit of integration	-
$b$	Upper limit of integration	-
$n$	Number of segments	-
$h$	Segment width	-
$f''(x)$	Second derivative of the function	-

This demonstrates that the composite rule has an overall truncation error of $O(h^2)$ .

Interactive Simulation

Use the simulation below to explore how the trapezoidal rule approximates areas using linear interpolation.

Composite Trapezoidal Rule

Approximating the integral of $f(x)$ from $a=0$ to $b=0.8$ using multiple segments.

Number of Segments (n): 2

1 Segment20 Segments

Integration Results

Step Size (h)

0.4000

Exact Value

1.640533

Approximate Integral (Area)

1.068800

Formula

I \approx \frac{h}{2} \left[ f(x_0) + 2 \sum_{i=1}^{n-1} f(x_i) + f(x_n) \right]

True Percent Relative Error:

34.8504%

As n increases, the error approaches zero.

Midpoint Rule (Open Form)

The simplest open Newton-Cotes formula is the midpoint rule, which approximates the area using a rectangle based on the function value at the midpoint $x_m = (a+b)/2$ . Since it doesn't use the endpoints $a$ or $b$ , it's highly useful when the integrand has an integrable singularity at a boundary limit.

Midpoint Rule

Approximates the integral using a single rectangle based on the midpoint.

I \approx (b - a) f(x_m)

Variables

Symbol	Description	Unit
$I$	Approximate integral value	-
$a$	Lower limit of integration	-
$b$	Upper limit of integration	-
$x_m$	Midpoint of the interval ( $\frac{a+b}{2}$ )	-
$f(x_m)$	Function value at the midpoint	-

Simpson's Rules

Simpson's rules use higher-order polynomials to connect points, yielding better accuracy for smooth functions than the trapezoidal rule.

Simpson's 1/3 Rule

Uses a parabola (second-order polynomial) connecting three points, exact for polynomials up to the third degree.

I \approx \frac{h}{3} [f(x_0) + 4f(x_1) + f(x_2)]

Variables

Symbol	Description	Unit
$I$	Approximate integral value	-
$h$	Step size ( $\frac{b-a}{2}$ )	-
$x_0$	First point ( $a$ )	-
$x_1$	Middle point	-
$x_2$	End point ( $b$ )	-

Composite Simpson's 1/3 Rule

Applies Simpson's 1/3 Rule over multiple segments. Requires an even number of segments.

I \approx \frac{h}{3} \left[ f(x_0) + 4\sum_{i=1,3,5}^{n-1} f(x_i) + 2\sum_{j=2,4,6}^{n-2} f(x_j) + f(x_n) \right]

Variables

Symbol	Description	Unit
$I$	Total approximate integral	-
$h$	Segment width	-
$n$	Even number of segments	-
$x_i$	Points with odd indices	-
$x_j$	Points with even indices	-

Global Error Bound for Composite Simpson's 1/3 Rule

Error is strictly bounded by the fourth derivative, showing an order of $O(h^4)$ .

|E_a| \le \frac{(b-a)^5}{180n^4} \max_{x \in [a,b]} |f^{(4)}(x)| = \frac{(b-a)h^4}{180} \max_{x \in [a,b]} |f^{(4)}(x)|

Variables

Symbol	Description	Unit
$E_a$	Global truncation error	-
$a$	Lower limit of integration	-
$b$	Upper limit of integration	-
$n$	Number of segments	-
$h$	Segment width	-
$f^{(4)}(x)$	Fourth derivative of the function	-

Thus, Simpson's 1/3 rule has a high accuracy of order $O(h^4)$ .

Simpson's 3/8 Rule

Uses a cubic curve (third-order polynomial) connecting four points.

I \approx \frac{3h}{8} [f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)]

Variables

Symbol	Description	Unit
$I$	Approximate integral value	-
$h$	Step size ( $\frac{b-a}{3}$ )	-
$x_i$	Evaluation points	-

The 3/8 rule also has a truncation error of $O(h^4)$ globally, and is primarily used when there is an odd number of segments, often combined with the 1/3 rule for the remaining even segments.

Integration with Unequal Segments

Experimental data often involves unequal spatial or temporal intervals. In this case, standard composite rules fail, and the basic single-segment trapezoidal rule is typically applied individually to each sub-interval and summed.

Romberg Integration

Romberg Integration

A numerical integration method that generates a triangular array of integral estimates by repeatedly applying Richardson extrapolation to results from the composite trapezoidal rule.

Romberg integration applies Richardson extrapolation recursively to the composite trapezoidal rule to systematically generate highly accurate integral estimates. By halving the step size and combining the results, it eliminates consecutive leading error terms ( $O(h^2)$ , $O(h^4)$ , $O(h^6)$ ). It requires the function to be evaluated at evenly spaced points, doubling the points at each refinement level.

Adaptive Quadrature

Standard Newton-Cotes formulas use a constant step size $h$ across the entire domain. However, if a function has regions of sharp variation (e.g., a steep spike) alongside regions of slow variation, a constant $h$ must be inefficiently small everywhere simply to capture the spike accurately.

Adaptive Quadrature

Numerical integration methods that dynamically alter the step size depending on the function's local behavior, optimizing accuracy and efficiency.

How Adaptive Quadrature Works

Adaptive methods dynamically adjust the step size $h$ based on an estimated local error. The algorithm computes the integral over a specific interval using two methods (e.g., a basic Simpson's rule and a refined Simpson's rule with twice the segments). If the absolute difference exceeds a specified tolerance, the interval is recursively halved in the problematic region. This concentrates computational effort only where the integrand is "difficult," drastically improving overall efficiency.

Gauss Quadrature

Gauss Quadrature

An integration rule that computes the integral as a weighted sum of function values at specific points within the domain of integration, treating both the evaluation points and their weights as variables to maximize accuracy.

While Newton-Cotes formulas use fixed, equally spaced evaluation points, Gauss quadrature treats both the evaluation points (roots or nodes) and their corresponding weights as unknowns, selecting them optimally to maximize the polynomial degree that can be integrated exactly.

Gauss-Legendre Form & Explicit Values

The most common form uses the roots of Legendre polynomials. An $n$ -point formula perfectly integrates polynomials up to degree $2n-1$ , making it exceptionally efficient for smooth functions.

Explicit values for low-order formulas:

Two-point formula ( $n=2$ ): Exact for 3rd degree polynomials.
- Roots: $t_1 = -1/\sqrt{3} \approx -0.57735$ , $t_2 = 1/\sqrt{3} \approx 0.57735$
- Weights: $c_1 = 1$ , $c_2 = 1$
Three-point formula ( $n=3$ ): Exact for 5th degree polynomials.
- Roots: $t_1 = -\sqrt{3/5} \approx -0.77459$ , $t_2 = 0$ , $t_3 = \sqrt{3/5} \approx 0.77459$
- Weights: $c_1 = 5/9$ , $c_2 = 8/9$ , $c_3 = 5/9$

To apply the method, the limits of integration must be algebraically transformed from the original domain $[a, b]$ to the standardized interval $[-1, 1]$ using a linear change of variables.

Gauss-Legendre Change of Limits Formula

Transforms integration limits to the standard interval $[-1, 1]$ to apply Gauss quadrature.

\int_{a}^{b} f(x) dx \approx \frac{b-a}{2} \sum_{i=1}^n c_i f\left( \frac{b+a}{2} + \frac{b-a}{2}t_i \right)

Variables

Symbol	Description	Unit
$a$	Original lower limit	-
$b$	Original upper limit	-
$n$	Number of evaluation points	-
$c_i$	Optimal weight for the $i$ -th node	-
$t_i$	Root of the $n$ -th degree Legendre polynomial	-

Improper Integrals

Improper integrals feature limits of $\pm \infty$ or integrands with singularities at a boundary. They can be evaluated numerically by using open Newton-Cotes formulas (which inherently avoid evaluating the boundary points) or by applying an algebraic change of variable (e.g., $x = 1/t$ , $dx = -1/t^2 dt$ ) to transform an infinite limit into a finite limit (like $t \to 0$ ).

Multiple Integrals

Double or triple integrals (e.g., calculating volumes, moments of inertia, or centroids) can be evaluated by successive, nested applications of 1D numerical integration rules. For instance, a double integral $\iint f(x,y) \,dx \,dy$ is solved numerically by holding $y$ constant, integrating over $x$ at each discrete $y$ -node, and then integrating those resulting 1D areas over the $y$ -axis.

Monte Carlo Integration

For low-dimensional integrals (1D, 2D, 3D), Gauss quadrature or Romberg integration are highly efficient. However, for high-dimensional multiple integrals (e.g., 10D or more), the number of required function evaluations grows exponentially for grid-based methods (the "curse of dimensionality").

Monte Carlo Method Conceptual Background

Monte Carlo integration uses random uniform sampling to approximate the integral. While the convergence rate is slow ( $O(1/\sqrt{N})$ ), it is fundamentally independent of the number of dimensions. Therefore, for very high-dimensional problems (common in statistical mechanics, quantum physics, and quantitative finance), Monte Carlo methods remain the only computationally feasible approach.

Monte Carlo Integration

Approximates an integral by randomly sampling points within the domain volume.

\int_V f(\mathbf{x}) d\mathbf{x} \approx V \cdot \frac{1}{N} \sum_{i=1}^N f(\mathbf{x}_i)

Variables

Symbol	Description	Unit
$V$	Total multidimensional volume of the domain	-
$N$	Number of random samples	-
$\mathbf{x}_i$	Randomly sampled point vector	-
$f(\mathbf{x}_i)$	Function evaluated at the sampled point	-

Civil Engineering Applications

Numerical integration is a foundational tool in civil engineering analysis:

Section Properties: Computing irregular cross-sectional areas, centroids, and moments of inertia by numerically integrating geometry data.
Earthwork Calculations: Using the trapezoidal or Simpson's rule to estimate soil excavation or fill volumes between discrete survey stations along a roadway.
Structural Analysis: Integrating distributed load profiles to determine shear force diagrams, and subsequently integrating shear to determine bending moment diagrams.

Key Takeaways

Newton-Cotes formulas replace the exact integrand with an easily integrable polynomial. They exist in both closed (uses boundaries) and open (avoids boundaries) forms.
The Trapezoidal rule uses a linear approximation with a global composite error strictly proportional to $O(h^2)$ .
Simpson's rules use parabolic (1/3 rule) and cubic (3/8 rule) approximations, offering much higher accuracy with a global composite error proportional to $O(h^4)$ .
Composite rules divide the domain into multiple smaller segments to drastically increase accuracy while maintaining simple low-order polynomial evaluations.
Romberg integration elegantly uses Richardson extrapolation to efficiently improve simple trapezoidal rule estimates to high orders of accuracy.
Adaptive Quadrature dynamically halves the local step size $h$ only in regions where the integrand varies rapidly, optimizing computational effort.
Gauss quadrature treats both nodes and weights as unknowns to maximize polynomial degree accuracy ( $2n-1$ ). It uses the Change of Limits formula to map any interval $[a, b]$ onto the standardized interval $[-1, 1]$ .
Improper integrals can be evaluated using open formulas (like the midpoint rule) or by transforming limits algebraically (e.g., mapping infinity to zero).
Multiple integrals are solved by successive nested applications of single-variable integration rules.
Monte Carlo Integration relies on random sampling and becomes the method of choice for extremely high-dimensional integrals, bypassing the curse of dimensionality.

Previous TopicNumerical Differentiation - Examples & Applications

Quiz Me

Next TopicNumerical Integration - Examples & Applications

Prev Next

Quiz Me