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
回答（1 个）
James Tursa
2022-5-16
编辑：James Tursa
2022-5-16
0 个投票
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)?
类别
在 帮助中心 和 File Exchange 中查找有关 Construct and Work with Object Arrays 的更多信息
2022-5-16
2022-5-16
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
