CODEWITHSUNDEEP
matlabnumerical-methods
Djamillah
Djamillah

Reputation: 289

What is wrong with my Simpson algorithm?

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.

Upvotes: 0

Views: 112

Answers (2)

David
David

Reputation: 8459

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.

Upvotes: 1

Geoff
Geoff

Reputation: 1603

Djamillah - your code looks fine though the initialization of h is probably valid only for the case where a==0, so you may want to change this line of code to

h = (b-a)/(2^n+1);

I wonder if x = a:h:b; will always be valid - sometimes b may be included in the list, and sometimes it might not be, depending upon h. You may want to reconsider and use linspace instead

x = linspace(a,b,2^n+1);

which will guarantee that x has 2^n+1 points distributed evenly in the interval [a,b]. h could then be initialized as

h = x(2)-x(1);

Also, when determining the even and odd indices, we need to ignore the last element of x for both even and odd. So instead of

xodd  = x(2:2:end-1);
xeven = x(3:2:end);

do

xodd  = x(2:2:end-1);
xeven = x(3:2:end-1);

Finally, rather than using a vector y (how is this set?) I might just use the function handle to the function that I'm integrating instead and replace the calculation above as

Isimp = delta/3*(func(x(1)) + 4*sum(func(xodd)) + 2*sum(func(xeven)) + ...
        func(x(end)));

Other than these tiny things (which are probably negligible), there is nothing in your algorithm to indicate a problem. It produced similar results to a version that I have.

As for the order of convergence, should it be O(n^4) or O(h^4)?

Upvotes: 3

Related Questions