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.
Example 17.5.1. Adaptive Gauss integration.
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{.}\)
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]\)
Question 17.5.2. Gauss-Kronrod error estimate.
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?