Lecture 6: January 29

Quadrature and Integration

The extended trapezoidal approximation and Romberg integration are very useful to compute most one-dimensional integrals with well-behaved integrands. These algorithms are discussed in more detail in Chapter 4 of "Numerical Recipes":

• §4.2 Elementary Algorithms discusses the extended trapezoid and Simpson rules, and

• §4.3 Romberg Integration discusses how to improve the trapezoid rule using Neville extrapolation to zero step size.

Programming the Rectangular Approximation

// Integ.java

import comphys.*;

public class Integ {

    // compute integral of f(x) form x = a to x = b
    // using rectangular approximation

    double f (double x) {
        return Math.cos(x);
    }

    double a;                   // lower limit
    double b;                   // upper limit
    int n;                      // number of intervals
    double dx;
    double delta;
    double sum;

    void initial () {
        a = Console.readDouble("lower limit a = ");
        b = Console.readDouble("upper limit b = ");
        n = 2;                  // initial number of intervals
        dx = (b - a) / n;
        delta = dx;             // delta = dx only for n = 2
        sum = f(a);             // -> (f(a) + f(b))/2 for trapezoid rule
        System.out.println();
        System.out.println("  n\trectangular approximation");
        System.out.println();
    }

    void sumf () {
        double x = a + dx;
        while (x < b) {
            sum += f(x);
            x += delta;
        }
    }

    void doubleIntervals () {   // can't use double which is a keyword!
        // double the number of intervals
        delta = dx;
        n *= 2;
        dx /= 2;
    }

    void output () {
        System.out.println(n + "\t" + Easy.format(sum * dx, 8));
    }
    
    public static void main (String[] args) {
        int doublings = 4;
        if (args.length > 0)
            doublings = Integer.parseInt(args[0]);
        Integ integ = new Integ();
        integ.initial();
        do {
            for (int i = 0; i < doublings; i++) {
                integ.sumf();
                integ.output();
                integ.doubleIntervals();
            }
        } while (Console.getYes("continue " + doublings + " more doublings"));
    }
}

If you are worried that Java might be too slow for your numerical analysis applications, try a Fortran version of PROGRAM integ translated as Integ.f on your computer system!

Multidimensional Quadrature

An algorithm, such as the trapezoid rule, which uses a regular lattice of points cannot easily be extended to multidimensional integrals for two reasons:

• First, if the total number of points in the lattice is n, the number of points available in each of d dimensions is of order n^1/d. Thus the trapezoid rule in one dimension has a global error of order 1/n². In d dimensions the global error will be of order 1/n^2/d.

• Second, when the dimensionality of the integral is larger than 5 - 10, using a regular lattice becomes impractical because the number of lattice points increases exponentially with dimension. For example, to evaluate the classical equilibrium partition function for a gas of 100 molecules, one would have to perform a 600-dimensional phase space integral over coordinates and momenta. If one applies the elementary trapezoid rule (only 2 points!!) in each dimension, the lattice has 2⁶⁰⁰ points, and the partition function is a sum of 2⁶⁰⁰ = 4 x 10¹⁸⁰ terms!!