erf(3) evaluation with recursive trapezoid rule

6 vues (au cours des 30 derniers jours)
redundant
redundant le 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
  1 commentaire
Sulaymon Eshkabilov
Sulaymon Eshkabilov le 27 Nov 2020
by the way, don't forget to include:
f = @(t)(exp(-t^2));

Connectez-vous pour commenter.

Réponses (1)

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

Catégories

En savoir plus sur Discrete Math dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by