Numerical Integration – Quadrature, Error, and Applications

Introduction

Numerical integration, or quadrature, is the process of approximating the definite integral
$$I={\int }_a^{b}f(x)dx$$
when an antiderivative of $f$ is unavailable or inconvenient. The idea is to replace the integral by a weighted sum of function values at selected points.

Fundamental Concepts

Definition

A numerical integration rule has the form
$$Q_n(f)=\sum _{i=0}^n{w}_{i}f(x_i)$$
where $x_i$ are nodes in $[a,b]$ and $w_i$ are weights. The rule is exact for all polynomials of degree $\le p$ if the weights and nodes satisfy the corresponding moment equations.

Interpolatory Quadrature

If $f$ is approximated by an interpolating polynomial $P_n(x)$ that matches $f$ at the nodes $x_i$, then
$${\int }_a^{b}f(x)dx\approx {\int }_a^b{P}_n(x)dx$$
The integral of $P_n$ can be computed exactly, yielding the weights $w_i$.

Error Term

For a rule exact for polynomials up to degree $p$, the error for a sufficiently smooth $f$ is
$$E=\frac{f^{(p+1)}(\xi )}{(p+1)!}{\int }_a^b\prod _{i=0}^n(x-x_i)dx,\xi \in (a,b)$$
The magnitude of the error depends on the step size $h$ and the $(p+1{)}^{\text{th}}$ derivative of $f$.

Newton‑Cotes Formulas

Figure 1 – Graph of $f(x)=e^x$ on $[-1,1]$ with nodes for the composite trapezoidal rule (n=4).

Derivation of Simpson’s 3/8 Rule

Using a cubic interpolant through four equally spaced points $x_0,x_1,x_2,x_3$ with spacing $h$, the integral of the interpolant over $[x_0,x_3]$ yields the weights shown above. The derivation follows the standard Lagrange interpolation procedure and integration of each basis polynomial.

Composite Rules

Composite rules apply a basic Newton‑Cotes rule over many subintervals.

Composite Trapezoidal Rule

For $n$ subintervals of width $h=(b-a)/n$:
$$T_n(f)=\frac{h}{2}[f(x_0)+2\sum _{i=1}^{n-1}f(x_i)+f(x_n)]$$

Composite Simpson’s 1/3 Rule

Requires an even number of subintervals $n=2m$:
$$S_n(f)=\frac{h}{3}[f(x_0)+4\sum _{i=1}^{m}f(x_{2i-1})+2\sum _{i=1}^{m-1}f(x_{2i})+f(x_n)]$$

Composite Simpson’s 3/8 Rule

Requires $n$ divisible by 3:
$$S_{3m}(f)=\frac{3h}{8}[f(x_0)+3\sum _{i=1}^{n-1}f(x_i)+2\sum _{i=1}^{n-1}f(x_{3i})+f(x_n)]$$

Figure 2 – Composite trapezoidal rule over $[0,6]$ with $n=6$.

Error Analysis of Composite Rules

For the composite trapezoidal rule:
$$E_T=-\frac{(b-a)h^2}{12}{f}^{\prime \prime }(\xi )$$
For composite Simpson’s 1/3:
$$E_S=-\frac{(b-a)h^4}{180}{f}^{(4)}(\xi )$$
Thus, for smooth functions, Simpson’s rule converges faster (order $h^4$) than the trapezoidal rule (order $h^2$).

Worked Example – Composite Trapezoidal Rule

Compute
$$I={\int }_{-1}^1{e}^{x}dx$$
using the composite trapezoidal rule with $n=4$.

1. Step size $h=\frac{1-(-1)}{4}=0.5$.
2. Nodes: $-1,-0.5,0,0.5,1$.
3. Function values:
   $$f(-1)=e^{-1}=0.3679,f(-0.5)=e^{-0.5}=0.6065,f(0)=1,f(0.5)=1.6487,f(1)=2.7183$$
4. Apply the formula:
   $$T_4=\frac{0.5}{2}[0.3679+2(0.6065+1+1.6487)+2.7183]$$
   $$=0.25[0.3679+5.4114+2.7183]=0.25\times 8.4976=2.1244$$
5. Exact value: $I=e^1-e^{-1}=2.7183-0.3679=2.3504$.
6. Absolute error: $|2.3504-2.1244|=0.2260$.

Figure 3 – Error plot for composite trapezoidal rule on $[-1,1]$ with $n=4$.

Double Integral Example – Simpson’s 1/3 Rule

Evaluate
$${\iint }_R\frac{dxdy}{xy},R=[2,2.6]\times [4,4.4]$$
using Simpson’s 1/3 rule in both directions.

1. Partition $x$ into two subintervals: $x_0=2,x_1=2.3,x_2=2.6$.
2. Partition $y$ into two subintervals: $y_0=4,y_1=4.2,y_2=4.4$.
3. Compute function values at the 9 grid points:
   $$f(x_i,y_j)=\frac{1}{x_i{y}_j}$$
   (values omitted for brevity).
4. Apply the 2‑D Simpson’s 1/3 rule:
   $$I\approx \frac{h_x{h}_y}{9}[f_{00}+f_{02}+f_{20}+f_{22}+4(f_{01}+f_{10}+f_{12}+f_{21})+16f_{11}]$$
5. Numerical evaluation gives $I\approx 0.1045$.

Figure 4 – Grid and weights for 2‑D Simpson’s 1/3 rule over $R$.

Error Sources in Numerical Integration

Figure 5 – Bar chart of error contributions for composite Simpson’s rule on a smooth function.

Advantages and Disadvantages

Applications

1. Finance – Computing total interest on a loan with a time‑varying interest rate $r(t)$:
   $$\text{Interest}={\int }_0^{T}Pr(t)dt$$
   Numerical integration approximates the area under the rate curve.

2. Traffic Engineering – Estimating average travel time on a route where speed $v(s)$ varies with distance $s$:
   $$\text{Time}={\int }_0^L\frac{ds}{v(s)}$$

3. Data Analytics – Calculating the cumulative distribution function (CDF) of a continuous variable from its probability density function (PDF).

In the real world

• eSewa – The mobile payment app calculates the total transaction fee for a batch of payments. The fee is a function of transaction amount $f(x)=0.02x+5$. eSewa uses the composite trapezoidal rule to integrate this fee function over the range of amounts processed in a day, yielding the total fee revenue.

• Google Maps – When a user requests a route that follows a curved road, Google Maps approximates the route length by integrating the speed profile along the path. The integral
  $${\int }_{\text{route}}ds$$
  is evaluated numerically using Simpson’s 1/3 rule on discretized GPS points.

• NEPSE (National Stock Exchange of Nepal) – The exchange calculates the daily trading volume index by integrating the trading volume density over the trading hours. A composite trapezoidal rule with hourly data points provides a quick estimate of the total volume.

Google Maps route visualization (Image: OathOn, CC BY-SA 4.0, via Wikimedia Commons)

Algorithm – Composite Trapezoidal Rule (Mermaid)

flowchart TD
  "Start"["Start"] --> "Input: f, a, b, n"
  "Compute h"["h = (b-a)/n"] --> "Initialize sum = f(a)+f(b)"
  "Loop i=1 to n-1"["For i=1 to n-1"] --> "x_i = a + i*h"
  "x_i" --> "sum = sum + 2*f(x_i)"
  "sum" --> "End Loop"
  "Result"["I = (h/2)*sum"] --> "Output I"
  "Output I" --> "End"

C Program – Composite Trapezoidal Rule

#include <stdio.h>
#include <math.h>

/* Function to integrate */
double f(double x) {
    return exp(x);   /* Example: e^x */
}

/* Composite trapezoidal rule */
double trapezoidal(double (*func)(double), double a, double b, int n) {
    double h = (b - a) / n;
    double sum = func(a) + func(b);
    for (int i = 1; i < n; ++i) {
        double x = a + i * h;
        sum += 2.0 * func(x);
    }
    return (h / 2.0) * sum;
}

int main(void) {
    double a = -1.0, b = 1.0;
    int n = 4;
    double I = trapezoidal(f, a, b, n);
    printf("Approximate integral = %.6f\n", I);
    return 0;
}

Exam tip

• Derivation questions: Show the Lagrange interpolant and integrate term‑by‑term to obtain weights.
• Error questions: Write the error formula for the requested rule and identify the derivative that appears.
• Worked example: Follow the step‑by‑step procedure: compute $h$, list nodes, evaluate $f$, apply the formula, and compute the error if the exact integral is known.
• Algorithm questions: Sketch a flowchart or write pseudocode; emphasize initialization, loop over subintervals, and final multiplication by $h/2$ or $h/3$.
• Multiple‑choice: Pay attention to subtle differences between Simpson’s 1/3 and 3/8 rules (order, step size divisibility).