Metodo di Simpson per radice quadrata

[magicfillo]

Ciao a tutti, in un problema mi è richiesto di calcolare l'integrale della radice di x tra 0 e 1 con il metodo del trapezio e di Simpson. Dalla teoria si sa che l'errore dei due metodi va con $n^-2$ per il trapezio e con $n^-4$ per simpson, ma andando a calcolare la legge di potenza in questo caso viene che l'errore va con $n^-(3/2)$. Come si spiega? Come può essere coerente con la teoria? (n è il numero dei sottointervalli)

Grazie in anticipo

[Raptorista]

Vediamo come hai calcolato quel risultato.

[feddy]

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:



Edit
Ho visto ora la risposta di Raptorista

[magicfillo]

Vediamo come hai calcolato quel risultato.

L'ho ottenuto calcolando l'errore per entrambi i metodi e facendo un grafico log-log tra numero di sottointervalli ed errore. Per entrambi i metodi la retta ha una pendenza di -1.5, che è l''esponente di n. Non mi spiego come questo possa essere coerente con la teoria.

[magicfillo]

E come si spiega questo fatto?

[feddy]

Ho scritto la mia ipotesi nelle prime righe. Per sicurezza aspetterei la risposta di Raptorista.

Testo nascosto, perché contrassegnato dall'autore come fuori tema. Fai click in quest'area per vederlo.

Mi è rcentemente capitata una cosa simile quando stavo risolvendo con differenze finite al secondo ordine un problema differenziale che ho poi notato avere soluzione nemmeno $C^(1)$, e a causa di questo la pendenza era scesa a circa $- \frac{3}{2}$

[magicfillo]

Se andiamo a prendere l'errore $E(n)=kn^(-4)$ con k che contiene la lunghezza dell'intervallo, la derivata quarta eccetera, e ne facciamo il logaritmo decimale: $log(E)=log(k)-4log(n)$, la derivata quarta va a modificare quindi solo l'intercetta della retta ma l'esponente della n dovrebbe essere invariato.

EDIT ho visto pra la tua risposta scusa

[feddy]

Non preoccuparti, piuttosto evita di citare tutto il mio messaggio quando rispondi

[Raptorista]

Succede la stessa cosa se fai l'integrale in un intervallo non problematico, per esempio \([0.5,1]\)?

[feddy]

Ciao Raptorista,

in questo caso (col mio codice, che suppongo essere lo stesso dell'OP) trovo l'ordine di convergenza aspettato.