Numerical integration

Updated 03 Nov 2026

• •

  Quadrature and sampled data

  • •

    Quadrature -> approximate $I=\int_a^b f(x)\,dx$ using weighted function/sample values.

  • •

    Function available -> evaluate it at nodes chosen by the rule.

  • •

    Samples only -> respect available sample locations; new evaluations require a justified interpolant/model.

  • •

    integration accumulates quantities: power -> energy, velocity -> displacement, depth over area -> volume.

  • •

    check units, sample order, smoothness, and known discontinuities before choosing a rule.

• •

  Composite Newton-Cotes rules

  • •

    let $n$ be the number of equal-width subintervals, $h=(b-a)/n$, and $f_i=f(a+ih)$.

    • •

      there are $n+1$ sample points, not $n$.

  • •

    Trapezoidal rule

    • •

      $T_n=h[(f_0+f_n)/2+\sum_{i=1}^{n-1}f_i]$.

    • •

      uses a straight line on each subinterval; global error is $O(h^2)$ for sufficiently smooth functions.

    • •

      Nonuniform data -> $T=\sum_{i=0}^{n-1}(x_{i+1}-x_i)(y_i+y_{i+1})/2$.

  • •

    Simpson 1/3

    • •

      requires even $n$ and equal spacing for this formula.

    • •

      $S_{1/3}=\frac h3[f_0+f_n+4\sum_{i=1,3,\ldots,n-1}f_i+2\sum_{i=2,4,\ldots,n-2}f_i]$.

    • •

      integrates local quadratic interpolants; exact for cubics in exact arithmetic.

  • •

    Simpson 3/8

    • •

      requires $n$ divisible by 3 and equal spacing.

    • •

      $S_{3/8}=\frac{3h}{8}[f_0+f_n+3\sum_{1\le i<n,\,3\nmid i}f_i+2\sum_{1\le i<n,\,3\mid i}f_i]$.

    • •

      integrates local cubic interpolants; exact for cubics in exact arithmetic.

    • •

      for odd $n\ge3$, use 3/8 on three subintervals and 1/3 on the remaining even number.

  • •

    Error bounds

    • •

      with $M_r=\max_{[a,b]}|f^{(r)}|$, sufficient smoothness gives:

      • •

        trapezoidal: $|E_T|\le(b-a)h^2M_2/12$.

      • •

        simpson 1/3: $|E_{1/3}|\le(b-a)h^4M_4/180$.

      • •

        simpson 3/8: $|E_{3/8}|\le(b-a)h^4M_4/80$.

    • •

      error is exact integral minus approximation; the bounds use its magnitude.

    • •

      the derivative bound must be justified; a difference between two estimates is not automatically a rigorous bound.

  • •

    Comparison example

    • •

      $\int_0^1 4/(1+x^2)\,dx=\pi$.

    • •

      with $n=4$: trapezoidal $3.1311764706$, Simpson 1/3 $3.1415686275$.

    • •

      for $\int_1^2\ln x\,dx$, $M_2=1$ gives $|E_T|\le1/(12n^2)$; $n=29$ suffices for error below $10^{-4}$.

• •

  Romberg integration

  • •

    apply Richardson extrapolation to composite trapezoidal estimates on successively halved grids

  • •

    First column -> $R_{k,1}=T_{2^{k-1}}$ for $k=1,2,\ldots$.

  • •

    Recurrence -> $R_{k,j}=[4^{j-1}R_{k,j-1}-R_{k-1,j-1}]/(4^{j-1}-1)$.

    • •

      cancels successive even-power error terms when the smooth error expansion applies

    • •

      the second column equals composite Simpson estimates on the corresponding refined grids.

    • •

      this is quadrature extrapolation, not a recursive least-squares method.

  • •

    Stopping

    • •

      compare successive diagonal estimates using absolute-plus-relative tolerance and a refinement limit.

    • •

      their difference is a convergence indicator, not a certified true-error bound.

    • •

      discontinuities, nonsmoothness, noise, and roundoff can defeat the expected improvement.

  • •

    Checked example: $\int_0^2e^{-x}\,dx=1-e^{-2}\approx0.8646647168$.

    • •

      first column: $1.1353352832$, $0.9355470828$, $0.8826039513$.

    • •

      second column: $0.8689510150$, $0.8649562408$.

    • •

      third diagonal: $0.8646899225$, correcting the PDF's printed value.

• •

  Gauss-Legendre quadrature

  • •

    choose both nodes and weights, rather than requiring equally spaced nodes

  • •

    Map to $[-1,1]$ -> $x(\xi)=(a+b)/2+(b-a)\xi/2$.

  • •

    Rule -> $\int_a^b f(x)\,dx\approx\frac{b-a}{2}\sum_{i=1}^m w_i f(x(\xi_i))$.

  • •

    Two points -> $\xi=\pm1/\sqrt3$, both weights 1.

  • •

    Three points -> $\xi=0,\pm\sqrt{3/5}$; weights $8/9$ at zero and $5/9$ at each outer node.

  • •

    an $m$-point rule is exact for polynomials of degree at most $2m-1$ in exact arithmetic.

    • •

      example: two points integrate $\int_0^2(x^3-2x+1)\,dx=2$ exactly mathematically.

    • •

      accuracy for other functions depends on the integrand; no fixed node count guarantees a universal tolerance.

  • •

    split at known discontinuities or use suitable composite/adaptive quadrature for localized structure.

• •

  Multiple integrals and engineering use

  • •

    Rectangle

    • •

      map each coordinate to $[-1,1]$ separately, including the scale factors

    • •

      tensor rule: $I\approx\frac{(b-a)(d-c)}4\sum_i\sum_jw_iw_jf(x(\xi_i),y(\eta_j))$.

    • •

      for $\int_{[0,1]^2}(x^2+y^2)\,dA$, one-cell tensor trapezoidal gives 1; $2\times2$ Gauss gives the exact $2/3$.

    • •

      tensor Gauss is exact up to degree $2m-1$ in each coordinate, not just in total degree.

  • •

    Unit right triangle

    • •

      map the unit square: $x=u$, $y=(1-u)v$ for $0\le u,v\le1$

    • •

      include the Jacobian $1-u$: integrate $f(u,(1-u)v)(1-u)$ over that square.

  • •

    Gridded volumes -> repeated 1-D integration can integrate depth over a rectangular grid.

  • •

    Higher dimensions

    • •

      tensor rules need $m^D$ evaluations for $m$ nodes in each of $D$ dimensions.

    • •

      ordinary Monte Carlo has sampling-error scale $O(N^{-1/2})$ under finite-variance assumptions.

    • •

      the variance constant still depends on the problem; it is not always better above a particular dimension.

  • •

    Unit-sensitive examples

    • •

      $E=\int V(t)I(t)\,dt$: volts times amperes times seconds gives joules.

    • •

      integrating kW over seconds gives kJ; divide by 3600 to obtain kWh.

    • •

      power samples $(t,P)=(0,0),(10,40),(30,50),(60,10)$ in seconds/kW give trapezoidal energy $2000$ kJ, or $0.555556$ kWh.

    • •

      heat flux $q''=-k\,dT/dx$ requires spatial distances in units compatible with conductivity.

    • •

      integrating acceleration requires an initial velocity; integrating velocity requires an initial position.

    • •

      integration of sensor bias can drift; smoothing does not by itself remove bias or supply initial conditions.

  • •

    

• •

  Python: quadrature checks

  • •

    standalone NumPy/SciPy example; trapezoid handles nonuniform samples with supplied coordinates.

  • •

    the composite Simpson implementation below deliberately requires equal spacing and even $n$.

  • •

    SciPy's sampled-data simpson(y, x=x) supports other cases, but its usual uniform-grid accuracy claims need qualification.

  • •

    import math
    import numpy as np
    from scipy.integrate import trapezoid
    from numpy.polynomial.legendre import leggauss
    
    def simpson_uniform(f, a, b, n):
          if not isinstance(n, int) or n < 2 or n % 2:
              raise ValueError("even number of subintervals required")
          if not np.isfinite([a, b]).all() or a >= b:
              raise ValueError("finite increasing bounds required")
          x = np.linspace(a, b, n+1)
          y = np.asarray(f(x), dtype=float)
          if y.shape != x.shape or not np.isfinite(y).all():
              raise ValueError("one finite value per node required")
          return (b-a)/(3*n)*(y[0]+y[-1]+4*y[1:-1:2].sum()+2*y[2:-1:2].sum())
    
    def romberg_table(f, a, b, levels=3):
          if not isinstance(levels, int) or not 1 <= levels <= 15:
              raise ValueError("levels must be an integer in 1..15")
          if not math.isfinite(a) or not math.isfinite(b) or a >= b:
              raise ValueError("finite increasing bounds required")
          def evaluate(x):
              y = f(x)
              if not math.isfinite(y):
                  raise ValueError("nonfinite integrand")
              return y
          R = np.full((levels, levels), np.nan)
          R[0, 0] = (b-a)*(evaluate(a)+evaluate(b))/2
          for k in range(1, levels):
              h = (b-a)/(2**k)
              new_sum = math.fsum(evaluate(a+(2*i-1)*h) for i in range(1, 2**(k-1)+1))
              R[k, 0] = R[k-1, 0]/2 + h*new_sum
              for j in range(1, k+1):
                  R[k, j] = R[k, j-1]+(R[k, j-1]-R[k-1, j-1])/(4**j-1)
          return R
    
    x = np.linspace(0., 1., 5)
    assert math.isclose(trapezoid(4/(1+x*x), x=x), 3.131176470588235)
    assert math.isclose(simpson_uniform(lambda x: 4/(1+x*x), 0., 1., 4), 3.1415686274509804)
    R = romberg_table(lambda x: math.exp(-x), 0., 2.)
    assert math.isclose(R[2, 2], 0.8646899224948368)
    nodes, weights = leggauss(2)
    mapped = 1 + nodes
    assert math.isclose(weights @ (mapped**3 - 2*mapped + 1), 2.)
    t = np.array([0., 10., 30., 60.])
    power_kw = np.array([0., 40., 50., 10.])
    energy_kj = trapezoid(power_kw, x=t)
    assert math.isclose(energy_kj, 2000.)
    print("Romberg table:", R, "energy (kWh):", energy_kj/3600)
    

◌ Explore connections in Graph view