numerical integration with recursive trapezoid rule

I am pretty new to Matlab and have to use the recursive trapezoid rule in a function to integrate f = (sin(2*pi*x))^2 from 0 to 1. The true result is 0.5 but I with this I get nothing close to it (approx. 3*10^(-32)). I can't figure out where the problem is. Any help is greatly appreciated.

function I = integrate(func,a,b)
         % f: function to integrate
         % a: lower bound
         % b: upper bound
         % I: integral of f from a to b computed by recursive trapezoidal
         %    rule
         
         f = str2func(func);
         n = 1;
         S = 0;
         h = b - a;
         newT = h / 2 * (f(a) + f(b)); 
         oldT = 0; 
         
         while abs(newT - oldT) > 0.00001
                oldT = newT; 
                h = (b - a) / 2^n;
                
                for i = 1 : 2^n
                    S = S + f(a + (2 * i - 1) * h);
                end
                
                newT = oldT / 2 + h * S; 
                n = n + 1;
         end   
         
         I = newT;
end

Answer 1

[James Tursa]

Some issues are immediately apparent.

First, you don't reset S=0 inside the while loop. Isn't S supposed to contain only the accumulated function values at the new points?

Second, your indexing doesn't appear to be correct in the for-loop. When you plug in the largest index you get this:

f(a + (2 * i - 1) * h) = f(a + (2 * 2^n - 1) * (b-a)/2^n) = f(a + 2*(b-a) - (b-a)/2^n) = f(2*b - a - (b-a)/2^n)

As n gets large the argument doesn't approach b, it approaches 2*b-a. It appears you need the top index of the for-loop to be 2^(n-1) for this indexing to work properly.

Are you coding this from an algorithm you were given? If so, can you post the algorithm (even if it is an image)?