Metodo dei trapezi per punti non equispaziati
Ciao a tutti,
come da titolo, volevo implementare in MatLab il metodo dei trapezi per un insieme di punti non equispaziati.
Per la regola dei trapezi, $ int_(a)^(b) f(x) dx = sum_(i = 1) ^(m) int_(x_i)^(x_(i+1)) f(x) dx ~~ sum_(i = 1) ^(m) ((x_(i+1)-x_i)/2)*(f(x_i)+f(x_i+1)) $
L'ho implementato con la seguente function:
Ora, per testarlo, ho provato tale routine sull'insieme di dati
$x=[0,2,5,6,7,8,10]$ e $y=[0,4,5,8,8.5,9,12]$.
Il risultato che trovo è che l'integrale vale $61.5.$
Provando a interpolare questo insieme di punti col comando polyfit, e trovando il polinomio interpolante, mi risulta che il valore esatto calcolato da MatLab è $70.6$.
Cambiando i valori trovo risultati che si differiscono di molto.
Eppure, testando su un insieme di punti non equispaziati ma tali che il polinomio interpolante sia una retta (esempio:$x=[0,2,5,6,7,8,10]$ e $y=[0,2,5,6,7,8,10]$) trovo il risultato esatto.
La domanda che mi sorge è: il metodo dei trapezi per punti non equispaziati è corretto (dal punto di vista del coedice) ma alquanto impreciso proprio per come è fatto, oppure ho sbagliato nell'implementarlo?
come da titolo, volevo implementare in MatLab il metodo dei trapezi per un insieme di punti non equispaziati.
Per la regola dei trapezi, $ int_(a)^(b) f(x) dx = sum_(i = 1) ^(m) int_(x_i)^(x_(i+1)) f(x) dx ~~ sum_(i = 1) ^(m) ((x_(i+1)-x_i)/2)*(f(x_i)+f(x_i+1)) $
L'ho implementato con la seguente function:
function I=trapeziNonEq(x,y)
%INPUT
%x->vettore contenente i punti xi
%y->vettore delle f(xi)
%OUTPUT
%valore dell'integrale tra a=x(1) e b=x(end)
m=length(x);
I=0;
for i=1:1:m-1
I=I+((x(i+1)-x(i))/2)*(y(i)+y(i+1))
end
Ora, per testarlo, ho provato tale routine sull'insieme di dati
$x=[0,2,5,6,7,8,10]$ e $y=[0,4,5,8,8.5,9,12]$.
Il risultato che trovo è che l'integrale vale $61.5.$
Provando a interpolare questo insieme di punti col comando polyfit, e trovando il polinomio interpolante, mi risulta che il valore esatto calcolato da MatLab è $70.6$.
Cambiando i valori trovo risultati che si differiscono di molto.
Eppure, testando su un insieme di punti non equispaziati ma tali che il polinomio interpolante sia una retta (esempio:$x=[0,2,5,6,7,8,10]$ e $y=[0,2,5,6,7,8,10]$) trovo il risultato esatto.
La domanda che mi sorge è: il metodo dei trapezi per punti non equispaziati è corretto (dal punto di vista del coedice) ma alquanto impreciso proprio per come è fatto, oppure ho sbagliato nell'implementarlo?
Risposte
Prova a confrontare con un risultato calcolato a mano: per esempio per i punti \([(0,1),(1,2),(1.5,2)]\) il risultato dovrebbe essere, se non ho contato male i quadretti, \(I = 2.5\). Ti viene?
Sì, viene corretto.
Penso che il mio risultato non combacia perché il comando polyfit con cui vado a costruire il polinomio interpolante può generare un polinomio interpolante che non rispecchia l'esatto andamento della funzione. Dovrebbe essere quello il motivo. Anche perché la formula da implementare è molto semplice.
Che dici, puo' essere?
Penso che il mio risultato non combacia perché il comando polyfit con cui vado a costruire il polinomio interpolante può generare un polinomio interpolante che non rispecchia l'esatto andamento della funzione. Dovrebbe essere quello il motivo. Anche perché la formula da implementare è molto semplice.
Che dici, puo' essere?
Anche secondo me il problema sta nel fatto che [inline]polyfit[/inline] fa "qualcosa", ma non sai bene che cosa.
Credo di aver trovato il problema. Ho notato che polyfit è particolarmente instabile per punti non equispaziati.
Il comando polyfit[nota]https://it.mathworks.com/help/matlab/ref/polyfit.html. Documentazione del comando. Ho notato che negli esempi vengono infatti sempre presi punti "equally spaced". Non c'è scritto però se questa scelta sia la migliore o meno[/nota] restiuisce i coefficienti del polinomio interpolante un insieme di punti.
Ho testato su più funzioni il seguente schema:
-Definisco una funzione (per esempio $sin(x)*cos(x)$).
-Su un insieme di punti (volontariamente non equispaziati), trovo le immagini:
e.g $xi=[0,1,3,7,10], yi=[f(0),f(1),f(3),f(7),f(10)]$.
-tramite polyfit trovo i coefficienti e costruisco con un ciclo opportuno il polinomio interpolante $q(x)$
-plotto i seguenti due grafici.
1.Quello ottenuto dalla funzione che ho definito in partenza sull'intervallo $[0,10]$ col comando plot(t,f(t),'b')
2.Quello ottenuto dalla valutazione del polinomio interpolante sullo stesso intervallo: plot(t,q(t),'r')
Noto che il polinomio interpolante passa sì per i nodi, ma evidentemente in modo differente da quello corretto.
E' questo il motivo dell'errore nel calcolo dell'area sottesa.
Il comando polyfit[nota]https://it.mathworks.com/help/matlab/ref/polyfit.html. Documentazione del comando. Ho notato che negli esempi vengono infatti sempre presi punti "equally spaced". Non c'è scritto però se questa scelta sia la migliore o meno[/nota] restiuisce i coefficienti del polinomio interpolante un insieme di punti.
Ho testato su più funzioni il seguente schema:
-Definisco una funzione (per esempio $sin(x)*cos(x)$).
-Su un insieme di punti (volontariamente non equispaziati), trovo le immagini:
e.g $xi=[0,1,3,7,10], yi=[f(0),f(1),f(3),f(7),f(10)]$.
-tramite polyfit trovo i coefficienti e costruisco con un ciclo opportuno il polinomio interpolante $q(x)$
-plotto i seguenti due grafici.
1.Quello ottenuto dalla funzione che ho definito in partenza sull'intervallo $[0,10]$ col comando plot(t,f(t),'b')
2.Quello ottenuto dalla valutazione del polinomio interpolante sullo stesso intervallo: plot(t,q(t),'r')
Noto che il polinomio interpolante passa sì per i nodi, ma evidentemente in modo differente da quello corretto.
E' questo il motivo dell'errore nel calcolo dell'area sottesa.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo