I was trying to write an algorithm to approximate integrals with Simpson's method. When I try to plot it in a loglog plot, however, I don't get the correct order of accuracy which is O(h^4) (I get O(n)). I can't find any errors though. This is my code:
%Reference solution with Simpson's method (the reference solution works well)
yk = 0;
yj = 0;
href = 0.0001;
mref = (b-a)/href;
for k=2:2:mref-1
yk=yk+y(k*href+a);
end
for j=1:2:mref
yj=yj+y(href*j+a);
end
Iref=(href/3)*(y(a)+y(b)+2*yk+4*yj);
%Simpson's method
iter = 1;
Ehmatrix = 0;
for n = 0:nmax
h = b/(2^n+1);
x = a:h:b;
xodd = x(2:2:end-1);
xeven = x(3:2:end);
yodd = y(xodd);
yeven = y(xeven);
Isimp = (h/3)*(y(x(1))+4*sum(yodd)+2*sum(yeven)+y(b));
E = abs(Isimp-Iref);
Ehmatrix([iter],1) = [E];
Ehmatrix([iter],2) = [h];
iter = iter + 1;
end
figure
loglog(Ehmatrix(:,2),Ehmatrix(:,1))
a and b are the integration limits and y is the integrand that we want to approximate.
Taking into account Geoff's suggestions, and making a few other changes, it all works as expected.
%Reference solution with Simpson's method (the reference solution works well)
a=0;
b=1;
y=@(x) cos(x);
nmax=10;
%Simpson's method
Ehmatrix = [];
for n = 0:nmax
x = linspace(a,b,2^n+1);
h = x(2)-x(1);
xodd = x(2:2:end-1);
xeven = x(3:2:end-1);
yodd = y(xodd);
yeven = y(xeven);
Isimp = (h/3)*(y(x(1))+4*sum(yodd)+2*sum(yeven)+y(b));
E = abs(Isimp-integral(y,0,1));
Ehmatrix(n+1,:) = [E h];
end
loglog(Ehmatrix(:,2),Ehmatrix(:,1))
P=polyfit(log(Ehmatrix(:,2)),log(Ehmatrix(:,1)),1);
OrderofAccuracy=P(1)
You were getting O(h)
accuracy because xeven=x(3:2:end)
was wrong. Replacing it by xeven=x(3:e:end-1)
fixes the code, and thus the accuracy.