erf(3) evaluation with recursive trapezoid rule

51 views (last 30 days)
redundant on 20 Nov 2013
I get an answer of 0.995231 with my code but I was wondering if I'm doing something wrong since when I do erf(3) the anwser is 1.
% Compute erf(3) by the recursive trapezoid rule
error = 0.0001;
a = 0;
b = 2;
coeff = 2 / sqrt(pi);
count = 0;
%- the two initial values (of integration part only)
n = 1;
h = (b - a) / (2^n);
I1 = 0.5 * h * (f(a) + f(b)) + h * f((a+b)/2);
count = count + 3;
n = n + 1;
h = h / 2;
I2 = 0.5 * I1 + h * (f((a+b)/4) + f((a+b)*3/4));
I3 = 0.5 * I2 + h * (f((a+b)/4) + f(a+b));
count = count + 2;
%- initial diff
diff = abs(I3 - I1);
%- loop
while (diff > error)
n = n + 1;
h = h / 2;
tmp_sum = 0;
for k = 1:(2^(n-1))
tmp_sum = tmp_sum + f(a + (2 * k - 1) * h);
count = count + 1;
end
I1 = I2;
I2 = I3;
I3 = (1/2) * I1 + h * tmp_sum;
diff = abs(I3 - I1);
end
%- Calculate the final result of erf(3)
result = I3 * coeff;
fprintf('Result: %f\n', result);
fprintf('Number of function evaluations: %d\n', count);
end
Sulaymon Eshkabilov on 27 Nov 2020
by the way, don't forget to include:
f = @(t)(exp(-t^2));

Sulaymon Eshkabilov on 27 Nov 2020
You need to set:
b = 3;