Skip to main content

Section 17.3 Recursive subdivision algorithm

The previous section showed how the accuracy of integration can be automatically estimated. Next question is, what can be done when it is insufficient? Using more evaluation points is logical. But if the algorithm keeps the interval of integration \([a, b]\) the same and just increases the number of points on that interval (e.g., in Simpson's rule), this is not the best way to adapt. Typically, the accuracy is poor because of insufficient smoothness or large derivative at some parts of the interval, while other parts are fine.

A better way to adapt is recursive subdivision.

  1. Given an interval \([a, b]\text{,}\) compute the integral over it and estimate the error.
  2. If the error is small enough, return the result of computation.
  3. If the error is too large, divide the interval into two halves \([a, c]\) and \([c, b]\) where \(c = (a+b)/2\text{,}\) and restart the process from step 1, for each half separately.

This process involves recursion, because at step 3 the function evaluating the integral has to call itself. Schematically it looks like this.

function I = adaptive(f, a, b)
    I1 = % first computation
    I2 = % second computation, more accurate
    if abs(I1 - I2) % "small enough"
        I = I2;
    else 
        c = (a+b)/2; 
        I = adaptive(f, a, c) + adaptive(f, c, b);
    end
end

One has to carefully consider when the difference abs(I1 - I2) is “small enough”. One approach is to ask for small relative error, like abs(I1 - I2) < 1e-6 * abs(I2). However this risks infinite recursion for functions like \(f(x) = \sqrt{x}\) near \(0\text{,}\) which is both small and poorly behaved. No matter how small an interval \([0, \epsilon]\) we use, dividing it further will result in a change that is not very small in relative terms. The process terminates with Matlab reporting: “Out of memory. The likely cause is an infinite recursion within the program.”

If one uses an absolute error bound like abs(I1 - I2) < 1e-6, there is another issue. We do not know how many subintervals will be used in the process, so even if each of them contributes an error of less than \(10^{-6}\text{,}\) the total may be significantly higher.

A third approach is to enforce enforce error bound relative to the size of interval, for example abs(I1 - I2) < 1e-6 * (b-a). This avoids the issue with functions that are problematic and small, such as \(f(x) = \sqrt{x}\) near \(0\text{.}\) But it runs into the same issue with functions that are problematic and large, such as \(f(x) = 1/\sqrt{x}\) near \(0\text{.}\)

A combination approach, such as abs(I1 - I2) < 1e-6 * (b-a) + 1e-9, allows the process to stop if the difference became small relative to interval size or just very small in absolute terms. This allows the subdivision to end even for functions like \(f(x) = 1/\sqrt{x}\text{.}\)