Simpson's Rule, answer inaccurate ?

I was trying to create a Matlab code using the Simpson's rule for a particular function.

%
function [y]= epow(x)
y=exp(-x.^2);
end
%
function [area]=simp(a,b,N)
format long

h=(b-a)/N;
area=0;
x=a:h:b;

for n=2:N,
term1=epow(x(n+1));
term2=epow(x(n));
term3=epow((x(n)+x(n+1))/2);
subarea=(h/6)*(term2+(4*(term3))+term1);

area=subarea+area;
end
end

With boundaries a=0, b=1 and the amount of subintervals N=128. Now when I run this function, I get simp(0,1,128)=0.739011791757018. However,the exact should be 0.842700792949715. So my answer isn't accurate and I can't find what I am doing wrong.

Thanks in advance.

Answer 1

[dpb]

0.842 is solution to the error function, erf(1). BUT, the error function is normalized such that INT(0,inf)=1 so is defined as 2/sqrt(pi)*INT(exp(-x^2) over [0,x].

Hence the "right" answer to the question raised is

>> erf(1)*sqrt(pi)/2
ans =
  0.7468
>>

So, your solution is off only a little. Don't have time to fully debug at the moment where you got off other than the normalization but one thing I note is that using

h=(b-a)/N;
x=a:h:b;

is going to accumulate error in x over the interval. Use instead

x=linspace(a,b,N);

or compute n*h and add over the range to keep the rounding to a minimum. Not sure it that'll quite fixup the result, but it'll probably help noticeably.

>> quad(@(x) exp(-x.^2),0,1)
ans =
  0.7468
>>

for comparison, too...

[Nick de Vries]

Thanks for the dpb. When I increase the number of subintervals the answer is almost the same as the exact (maybe this is the way to increase the accuracy?).

When I replace

a:h:b;

by

linspace(a,b,N);

I get an error: Attempted to access x(129); index out of bounds because numel(x)=128.

Futhermore I still don't quite understand why I should add 2/sqrt(pi) to the error function in case of Simpson's rule. I also used the Trapezoidal rule and there I could directly compare both answers without adding the 2/sqrt(pi). I there a simple way to get erf(1)=0.7468 instead of erf(1)=0.842 in Matlab? Sorry if this sounds obvious, I am new to Matlab.

[dpb]

A) Because your loop runs over 2:N but the addressing covers N+1 at the last element. 0:h:1 generates 129 total points whereas linspace actually creates the number of points requested.

In the particular case of 1/128 the difference isn't material as the delta is representable; I had intended a revision to remove the comment but it didn't seem to get saved...it's a general rule, however, that is true so not a bad thing to keep in mind. Use length(x)-1 as the upper limit in the loop instead of N to solve the indexing issue.

B) erf() is defined such that I([0:inf])==1. The 2/sqrt(pi) is the required normalization constant such as to make that be so.