Skip to main content

Section 17.5 Examples and questions

These are additional examples for reviewing the topic we have covered. When reading each example, try to find your own solution before clicking “Answer”. There are also questions for reviewing the concepts of this section.

Implement adaptive integration as in Section 17.3 using the midpoint rule for \(I_1\) and Gauss rule with \(n=2\) for \(I_2\text{.}\) Test it on the integral \(\int_0^9 x^{-1/2}\,dx = 6\text{.}\)

Answer

When implementing Gauss rules on a general interval \([a, b]\text{,}\) it helps to introduce the midpoint \(c = (a+b)/2\) and half-length \(k = (b-a)/2\text{.}\) This is because the rule was originally designed for \([-1, 1]\) and the linear function from \([-1, 1]\) to \([a, b]\) is \(kx+c\text{.}\) The coefficient \(k\) will appear in the formulas due to the change of variable formula.

f = @ (x) x.^(-1/2);
fprintf('The integral is %.6f\n', adaptive(f, 0, 9));

function I = adaptive(f, a, b)
    c = (a+b)/2;
    k = (b-a)/2;
    I1 = k*2*f(c);    % midpoint, same as Gauss n=1
    I2 = k*(f(c - k/sqrt(3)) + f(c + k/sqrt(3)));  % Gauss n=2
    if abs(I1 - I2) < 1e-6 * (b-a) + 1e-9 
        I = I2;
    else
        I = adaptive(f, a, c) + adaptive(f, c, b);
    end
end

Note that the function is not defined at 0 which is an endpoints of the interval of integration. This would be a problem if we used, for example, trapezoidal or Simpson's rules. Not a problem for midpoint or Gauss rules. By the way, the midpoint rule can be viewed as the Gauss rule with 1 point, because the degree 1 Legendre polynomial is \(P_1(x)=x\text{,}\) with the zero at the middle of interval \([-1,1]\)

The current version of Wikipedia article Gauss–Kronrod quadrature formula says: “The integral is then estimated by the Kronrod rule K15 and the error can be estimated as |G7-K15|”. This estimate, while reasonable, is quite conservative: the actual error will be much less. Can you explain why?