Metodo dei trapezi per punti non equispaziati

[feddy]

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:

Codice:

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?

[Raptorista]

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?

[feddy]

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?

[Raptorista]

Anche secondo me il problema sta nel fatto che polyfit fa "qualcosa", ma non sai bene che cosa.

[feddy]

Credo di aver trovato il problema. Ho notato che polyfit è particolarmente instabile per punti non equispaziati.

Il comando polyfit¹ 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.

---

Note

1. 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 ↑