Using Romberg integration with improper integral

928 Views Asked by At

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.

1

There are 1 best solutions below

1
On BEST ANSWER

@ThomasAndrews lent the solution. I changed my function implementation so that it's manually defined at the undefined point. A good work around!

function y = f2(t)
if(t == 0)
    y = 1;
    return;
end
y = sin(t) ./ t;
end