numerical integration with recursive trapezoid rule

32 views (last 30 days)
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

Answers (1)

James Tursa
James Tursa on 16 May 2022
Edited: James Tursa on 16 May 2022
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)?

Categories

Find more on Simulink in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!