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]\)
