Simpson's rule - Ciallo Wiki

1 Introduction

1.1 Problem of numerical integration

Numerical integration aims to approximate a definite integral when the integrand is difficult to integrate analytically or when only function values at discrete points are available. The central task is to replace an exact area computation with a computable sum based on sampled values.

1.2 Where Simpson’s rule fits among integration methods

Simpson’s rule belongs to the family of Newton–Cotes quadrature methods, which approximate an integral by fitting low-degree polynomials through data points. In practice, it is often chosen when the function is sufficiently smooth and an accuracy level better than the trapezoidal rule is desired without moving to more complex schemes.

1.3 Key idea: quadratic approximation on subintervals

The rule approximates the integrand over small sections using parabolic arcs. Instead of using straight-line segments between points, it fits a quadratic polynomial on every pair of adjacent subintervals, then computes the area under that quadratic. These local areas are summed to approximate the overall integral.

2 Mathematical Formulation

2.1 Setup of the interval and partition

Let \(f\) be the function to integrate on \([a,b]\). Choose a positive integer \(n\) and partition the interval into \(n\) equal subintervals of width \[ h=\frac{b-a}{n}. \] For the composite Simpson’s rule, \(n\) must be even. Define nodes \[ x_i=a+ih,\quad i=0,1,\dots,n. \]

2.2 The composite Simpson’s rule formula

With \(n\) even, Simpson’s rule approximates \[ \int_a^b f(x)\,dx \] by \[ \frac{h}{3}\Big[f(x_0)+f(x_n)+4\sum_{i=1,\, i\ \text{odd}}^{n-1} f(x_i)+2\sum_{i=2,\, i\ \text{even}}^{n-2} f(x_i)\Big]. \]

2.2.1 Weighting pattern (1–4–2–4–…–2–4–1)

The coefficients applied to function values follow the alternating weight sequence \[ 1,\ 4,\ 2,\ 4,\ 2,\ \ldots,\ 4,\ 1, \] reflecting the fact that quadratic fits are applied over overlapping pairs of subintervals. Endpoints receive weight \(1\), odd-indexed interior points receive weight \(4\), and even-indexed interior points (excluding endpoints) receive weight \(2\).

2.3 Single-interval Simpson’s rule

For a single pair of subintervals, take \([a,b]\) with midpoint \(m=(a+b)/2\) and step \(h=(b-a)/2\). The Simpson approximation over \([a,b]\) is \[ \int_a^b f(x)\,dx \approx \frac{b-a}{6}\Big[f(a)+4f(m)+f(b)\Big]. \] This is the fundamental block from which the composite rule is built by applying it repeatedly across the partition.

2.4 Relationship to polynomial interpolation

The composite Simpson’s rule can be viewed as interpolatory: on each pair of subintervals, a quadratic polynomial is constructed that matches \(f\) at three nodes. The integral of that polynomial over the local interval becomes the local approximation. Summation across the whole partition yields the composite formula.

3 Derivation and Intuition

3.1 Deriving from Lagrange interpolation

On an interval \([x_{2k},x_{2k+2}]\) with midpoint \(x_{2k+1}\), construct the Lagrange interpolating polynomial \(p_2(x)\) of degree at most two satisfying \[ p_2(x_{2k})=f(x_{2k}),\quad p_2(x_{2k+1})=f(x_{2k+1}),\quad p_2(x_{2k+2})=f(x_{2k+2}). \] Then approximate the integral over that interval by \(\int_{x_{2k}}^{x_{2k+2}} p_2(x)\,dx\).

3.2 Area under interpolating quadratic

Because the interpolant is quadratic and matches the function at the three chosen points, the integral of the interpolant can be written explicitly in terms of \(f(x_{2k})\), \(f(x_{2k+1})\), and \(f(x_{2k+2})\). Carrying out that polynomial integration produces the local form \[ \int_{x_{2k}}^{x_{2k+2}} f(x)\,dx \approx \frac{h}{3}\Big[f(x_{2k})+4f(x_{2k+1})+f(x_{2k+2})\Big], \] with \(h\) the subinterval width of the full partition.

3.3 Combining local approximations over the partition

Apply the local approximation on each pair \((x_{2k},x_{2k+2})\) for \(k=0,1,\dots,\frac{n}{2}-1\), and sum the results. Interior points appear in multiple local contributions, which is why their overall coefficients become \(4\) for odd indices and \(2\) for even indices, while the endpoints appear only once and therefore carry weight \(1\). This yields the composite Simpson’s rule.

4 Accuracy and Error Analysis

4.1 Error term for Simpson’s rule

For sufficiently smooth functions, Simpson’s rule has an error that depends on higher derivatives of \(f\). In the composite setting on \([a,b]\) with \(n\) even, the truncation error scales like a constant multiple of \(h^4\) times the maximum (or an appropriate bound) of the fourth derivative over the interval. This \(h^4\) behavior is a key reason it often performs well on smooth problems.

4.2 Conditions under which accuracy is improved

Accuracy typically improves when:

the integrand is smooth (at least four derivatives exist and are not too large),
the step size \(h\) is small,
the function does not oscillate rapidly within subintervals.

In such cases, the local quadratic approximation better tracks the function’s curvature, leading to a smaller cumulative error.

4.3 Comparison with the trapezoidal rule

The trapezoidal rule can be regarded as using linear interpolation on each subinterval. For smooth functions, its error decreases like \(h^2\), which is generally slower than Simpson’s \(h^4\). As a result, Simpson’s rule often achieves the same accuracy with fewer function evaluations, especially for well-behaved integrands.

4.4 Comparison with higher-order alternatives

Higher-order Newton–Cotes formulas and other quadrature schemes can sometimes surpass Simpson’s rule, particularly for extremely smooth functions or when specialized nodes are used. However, higher-order methods may introduce additional oscillatory behavior or greater sensitivity to endpoint effects. Gaussian quadrature, for example, can achieve high precision with fewer evaluations by choosing optimal nodes, though it may require more sophisticated setup.

5 Practical Computation

5.1 Choosing the number of subintervals (even requirement)

Because the composite rule relies on pairs of subintervals, \(n\) must be even. If an application suggests an initially chosen step count, it should be adjusted upward to the nearest even integer to use Simpson’s rule directly.

5.2 Handling function evaluations efficiently

The method evaluates \(f\) at equally spaced nodes \(x_i\). When refining the partition, previously computed values at shared nodes can often be reused. This is especially beneficial when each function evaluation is expensive (e.g., involves simulation or complex models).

5.3 Implementation tips and common pitfalls

Common issues include:

forgetting the even requirement for \(n\),
misapplying weights (especially mixing up odd and even indices),
off-by-one errors in loops that assemble the sums,
inconsistent handling of endpoints when computing the composite total.

A robust implementation explicitly constructs the node list or uses carefully indexed loops that match the weighting pattern.

A practical way to gauge accuracy is to compute Simpson’s approximation with step size \(h\) (using \(n\) even) and then recompute with \(h/2\) (using \(2n\), still even). Because the leading error term scales with \(h^4\), the difference between the two approximations can provide an empirical error estimate and a sanity check that refinement is improving results.

6 Examples and Applications

6.1 Integrating a smooth polynomial

If \(f(x)\) is a polynomial of degree at most three, Simpson’s rule yields the exact integral over \([a,b]\) (for the composite version, when the partition aligns with the method’s assumptions). For higher-degree polynomials, the method remains often accurate, with error governed by the next higher derivatives.

6.2 Integrating a non-polynomial function (benchmark example)

For functions such as \(\sin(x)\), \(\exp(x)\), or \(\ln(1+x)\) on a moderate interval, Simpson’s rule typically delivers good results because these functions are smooth and well approximated by quadratics over small spans. In benchmarking, error often drops rapidly when the number of subintervals increases, reflecting the \(h^4\) trend.

6.3 Approximation of definite integrals in applied problems

In applied contexts—such as computing expected values in probability (as integrals), cumulative quantities in physics (e.g., work as an integral), or areas under measured curves—Simpson’s rule can serve as a reliable general-purpose tool. It is particularly useful when the integrand can be evaluated at many points but symbolic integration is unavailable.

6.4 Visualizing the parabolic arcs and resulting area

A helpful visualization is to draw the integrand curve and then superimpose parabolic arcs across each pair of subintervals. The approximation takes the area under these arcs rather than the exact curve area. Seeing how the arcs follow curvature near the nodes can clarify why the technique performs strongly for smooth functions.

7 Variants of Simpson’s Rule

7.1 Composite Simpson’s rule

The composite version applies the single-interval rule repeatedly across a uniform partition with an even number of subintervals. It increases stability for long intervals and allows systematic refinement by reducing \(h\).

7.2 Simpson’s 3/8 rule (overview and when used)

Simpson’s 3/8 rule is another Newton–Cotes method based on a cubic interpolant over three subintervals. It has a different weighting scheme and can be advantageous for particular smoothness patterns or when the number of subintervals is naturally a multiple of three. Like Simpson’s rule, it approximates the integral using the area under a polynomial interpolant.

7.3 Adaptive Simpson’s rule concept

Adaptive Simpson’s rule aims to concentrate computational effort where the integrand is harder to approximate. Instead of using one global partition, it recursively subdivides subintervals and compares results at different refinement levels to decide where further splitting is needed.

7.3.1 Recursion based on local error estimation

A common strategy computes Simpson’s approximation on an interval and then recomputes it after splitting the interval into two halves. The difference between these estimates provides an indicator of local error. If the error proxy exceeds a tolerance, the method recurses on the subintervals; otherwise it accepts the refined approximation and aggregates the results.

8 Connections and Extensions

8.1 Numerical methods perspective (Newton–Cotes connection)

From a broader viewpoint, Simpson’s rule is a special case of Newton–Cotes quadrature. The underlying idea—integrate a fitted polynomial over the same interval—unifies many quadrature rules. The choice of polynomial degree and sampling pattern determines the resulting weights and accuracy properties.

8.2 Relation to Gaussian quadrature (high-level comparison)

Gaussian quadrature chooses nodes and weights to optimize accuracy for polynomials of a given degree, often requiring fewer evaluations than uniform-node methods. Simpson’s rule instead uses equally spaced nodes and fixed weights derived from quadratic interpolation. Both are interpolatory quadrature ideas, but Gaussian quadrature is typically more efficient per evaluation for smooth problems.

8.3 Extension to integrating sampled data

In many real applications, the integrand is not available as a closed-form expression but as discrete samples. Simpson’s rule can still be applied when samples correspond to equally spaced nodes. The method then acts directly on the provided values, producing an approximate integral consistent with a piecewise-quadratic fit to the data.

9 Summary

9.1 When to use Simpson’s rule

Simpson’s rule is well suited for integrating smooth functions over an interval when:

function values can be evaluated at evenly spaced points,
an even number of subintervals is acceptable,
improved accuracy over the trapezoidal rule is desired without large computational overhead.

9.2 Strengths, limitations, and best practices

Strengths include strong accuracy for sufficiently smooth integrands and straightforward implementation using a clear weighting pattern. Limitations arise when the function is not smooth or has sharp features within the interval, where a fixed partition may not capture behavior effectively. Best practices include using step refinement or adaptive strategies, verifying the even subinterval condition, and checking sensitivity by increasing the number of subintervals to confirm convergence.