Metodo di Simpson per radice quadrata
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
Grazie in anticipo
Risposte
Vediamo come hai calcolato quel risultato.
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.
A conferma di quanto detto:
Edit
Ho visto ora la risposta di Raptorista
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.
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
"Raptorista":
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.
E come si spiega questo fatto?
Ho scritto la mia ipotesi nelle prime righe. Per sicurezza aspetterei la risposta di Raptorista.
[ot]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}$[/ot]
[ot]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}$[/ot]
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
EDIT ho visto pra la tua risposta scusa
Non preoccuparti, piuttosto evita di citare tutto il mio messaggio quando rispondi 
Succede la stessa cosa se fai l'integrale in un intervallo non problematico, per esempio \([0.5,1]\)?
Ciao Raptorista,
in questo caso (col mio codice, che suppongo essere lo stesso dell'OP) trovo l'ordine di convergenza aspettato.
in questo caso (col mio codice, che suppongo essere lo stesso dell'OP) trovo l'ordine di convergenza aspettato.
La mia domanda era per magicfillo, ma grazie per la conferma. Sicuramente è colpa della derivata che esplode.
Sì immaginavo, ma essendo una cosa che mi interessa non sono riuscito ad aspettare 
Appena arrivo a casa provo! Quindi come posso spiegare la situazione? È per il fatto che la derivata non è definita in zero?
Ho provato, ed effettivamente torna il normale esponente di -4 per Simpson.
Se la derivata non è limitata allora la stima dell'errore non è utile e non c'è motivo per cui l'errore debba decadere con una potenza piuttosto che un'altra. Nella pratica, come calcoli la derivata in zero?
"Raptorista":
Se la derivata non è limitata allora la stima dell'errore non è utile e non c'è motivo per cui l'errore debba decadere con una potenza piuttosto che un'altra. Nella pratica, come calcoli la derivata in zero?
Non ho ben capito la domanda, io non calcolo la derivata in zero, la consegna è di trovare la legge di potenza e spiegare come mai differisce dalla teoria, non saprei come dirlo bene.
Ah già, colpa mia che stavo pensando a un'altra cosa. Non calcoli la derivata, correttamente. La spiegazione è quella che ho dato sulla stima dell'errore.
"Raptorista":
Ah già, colpa mia che stavo pensando a un'altra cosa. Non calcoli la derivata, correttamente. La spiegazione è quella che ho dato sulla stima dell'errore.
Cioè quindi che la derivatà non è limitata, corretto?
Sì. Adesso non ricordo se ci sono ipotesi sulla derivata della funzione nel teorema della stima dell'errore, ma se non ci sono limitazioni allora la stima è comunque valida, ma la costante potrebbe essere non limitata o cose del genere.
Si può quindi dire che l'errore non segue la legge di potenza teorica in quanto la derivata della funzione non è definita in zero?
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo