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?