Ammesso e non concesso che l'implementazione sia corretta, l'ordine dell'errore è legato alle derivate possedute dall'integranda $f=\sqrt(x)$: in particolare il metodo di Simpson dipende da $f^{(4)}(\xi)$, con $\xi$ che sta nell'intervallo $(0,1)$. Probabilmente questo comporta una certa perdita di ordine.
Ho provato in velocità a scrivere qualche riga e confermo quanto dici: l'ordine decade esattamente come $n^{-1.5}$, come mostro nel seguente codice, dove
Cavalieri è una function che prende in ingresso $f$ definite come function handle, gli estremi $a,b$ sopra definiti e il numero di nodi. Per comodità ti metto la function Cavalieri in spoiler, ma l'ho scritta in 5 minuti, senza pensarci neanche ed è la versione naive.
Testo nascosto, fai click qui per vederlo
function I=Cavalieri(f,a,b,m)
if (mod(m,2)~=0)
disp('il numero di intervalli deve essere pari')
else
h=(b-a)/m;
s1=0;
for i=2:2:m
s1=s1+f(a+(i-1)*h);
end
s2=0;
for i=3:2:m-1
s2=s2+f(a+(i-1)*h);
end
I=h/3*(f(a) + 4*s1 + 2*s2 +f(b));
end
- Codice:
clear all
close all
nrange=[100:10:2000];
a=0;b=1;
f=@(x) sqrt(x);
I_ex=quad(f,a,b); %calcolo il risultato di riferimento
err=NaN;
count=0;
for n=nrange
count++;
I=Cavalieri(f,a,b,n);
err(count)=abs(I-I_ex);
end
loglog(nrange,err,'*',nrange,err(end)*(nrange/nrange(end)).^(-4),'r',...
nrange,err(end)*(nrange/nrange(end)).^(-1.5),'k')
legend('err','n^-4','n^-1.5')
xlabel('n')
ylabel('err')
A conferma di quanto detto:
EditHo visto ora la risposta di Raptorista