Problem with matlab, simpsonrule

Well I have a script and function that is supposed to approximate the integral of cos between 0 to 1 with the simpsonrule. the problem I have is that it works great for 10,100, and 1000 itterations but when I go to 10000 it kinda fails. I have a log plot that shows the change in the error linear. The error to the interval is all nice and linear until I reach 10000. at 10000 it kinda "breaks" and the change in error is less than it "should" be. I have no idea why it fails at 10000. I have checked everything but I fail to see whats wrong, especially since it works for the other amount of itterations.

This is the simpson

summa = f(x(1))+2*sum(f(x(3:2:end-2)))+4*sum(f(x(2:2:end)))+f(x(end));

I divide it by h/3 later.

I have obviously done something wrong but I can't figure out what. I have gone trhough the code multiple times and even tried a different type of simpson that uses a forloop but gets basically the same result.

help would be appriciated thanks

Here is all of the code I use

SCRIPT

exakt_varde =sin(1); %Exakta värdet av integralen av sin(1)

for i=1:4

    n(i)=10^i

    fel(i)=abs(simpson (@cos,0,1,n(i))-exakt_varde);

end

h=(1-0) ./n; %

[h;fel] %skriv ut felet för olika värden på h

k=(log(fel) ./log(h)) % Sriv ut värdet på k i de olika delintervallen

plot(log(h), log(fel), 'o-')

title('Fel för Simpsonmetoden av cos(x) mellan 0 och 1')

xlabel('ln(h)')

ylabel('ln(fel)')

FUNCTION

function integral=simpson(f,a ,b ,n) %f är funktionen, a och b start-slut

format long % Vi vill ha fler decimaler än 4 i utskriften

h=(b-a)/n; %intervalllängden beräknas

x=linspace(a,b,n+1)%x innehåller de n+1 ändpunkterna för

delintervallen

summa=0; % nollställning av summan

%for i=1:2:n

%summa=summa+(f(x(i)))+(4*(f(x(i+1))))+f(x(i+2));

%end

summa = f(x(1))+2*sum(f(x(3:2:end-2)))+4*sum(f(x(2:2:end)))+f(x(end));

%Simpsonsuträkningen föruotom det sista som vi lagt nedanför för att göra

%föregående formel kortare

integral=h/3*summa; %slutligen multipliceras med h/3 enligt formeln

Answer 1

[Mike Hosea]

Looks to me like you're just running out of floating point precision. Double precision has about 16 decimal digits max, and you'll probably lose at least one or two to round-off error.

[Jaksan]

Thanks for the answer. That sounds plausible as the error with 1000 itterations max out at the 15th decimal and with 10000 itteration it just goes as low as it can with 1*10^-14.

Is there some way to adress this?

[Mike Hosea]

No, not really. You can kick the proverbial can down the road by computing your quadrature rule in a symbolic environment that provides multiple precision. Then you can set your floating point numbers to have as much precision as you want, but how much precision do you really need and why?