I'm trying to write a MATLAB function that implements Romberg's method of integration.
Problem is that I'm trying to approximate:$$\int_0^1\frac{\sin{x}}{x}dx$$ but the function is not defined at $x = 0$.
This is my MATLAB code:
function r = romberg(f,a,b,n)
h = (b - a) ./ (2.^(0:n-1));
r(1,1) = (b - a) * (feval(f, a) + feval(f, b)) / 2;
for j = 2:n
subtotal = 0;
for i = 1:2^(j-2)
subtotal = subtotal + feval(f, a + (2 * i - 1) * h(j));
end
r(j,1) = r(j-1,1) / 2 + h(j) * subtotal;
for k = 2:j
r(j,k) = (4^(k-1) * r(j,k-1) - r(j-1,k-1)) / (4^(k-1) - 1);
end
end;
It replies with NaN when run with the integrand above. I have a tolerance of 10^-9 set, thinking I might use a+TOL as my lower bound, but that might be wrong.
@ThomasAndrews lent the solution. I changed my function implementation so that it's manually defined at the undefined point. A good work around!