Interpolazione polinomiale

10 04 2011

Cominciamo a spiegare cosa è l’interpolazione polinomiale. Immaginate di avere una curva analitica da studiare y=”qualcosa di complicato”, oppure una serie di dati in coppie (x1,y1), (x2,y2)…. con tutti gli xi diversi.
Come potete ben immaginare non sarà sempre un oggetto semplice da studiare.

Un tipo di curve particolarmente semplici da studiare sono i polinomi, e quindi ci piacerebbe trovare dei polinomi che assomiglino alle funzioni che vogliamo studiare, in modo da avere un lavoro più semplice da svolgere.
Per chi conosce relativamente bene la matematica avrà magari pensato ai polinomi di Taylor o al caso particolare McLaurin, ma noi vogliamo prendere polinomi che siano anche molto semplici da costruire. Fare le derivate non è sempre una operazione semplice e poi questi due polinomi sono adatti per fare una analisi locale in un punto, non sempre approssimano una funzione su un intervallo.

Il metodo è utilizziamo è il seguente (supponendo di dover sudiare una funzione in una variabile o qualcosa di analogo)
Scegliamo sull’intervallo che vogliamo studiare n+1 punti. Fissati questi punti sull’ascissa delle x li valutiamo e otteniamo la coppie (x1,y1),(x2,y2),…

A questo punto si può dimostrare che esiste un solo polinomio di grado n (ad esempio per due punti distinti passa una sola retta, per tre punti distinti una sola parabola…). Il problema (ma neanche tanto complicato) è come costruirlo.

Il primo metodo è quello di Vandermonde. Si considera un generico polinomio di grado n e si impongono tutte le condizioni. Dovremo quindi poi risolvere un sistema lineare, che coincide con una matrica, guarda a caso, chiamata matrice di Vandermonde…

Il secondo metodo, forse più “diretto” e che uso più spesso (ma siete liberi di usare il metodo che volete) è quello di Lagrange.
Si costruiscono n polinomi L_i(x) di grado n tali che $latexL_i(x_j)=0$ per ogni i diverso da j, mentre L_i(x_i)=1. Il polinomio interpolante sarà quindi semplicemente \sum_i(L_i(x)*y_i). Provare per credere se inserite nel polinomio il valore x_i otterrete il valore y_i.

Facciamo un esempio, dati i nostri punti da interpolare:
(1,2), (2,3), (3,4)
Costruiamo
L_1(x)=\frac{(x-2)(x-3)}{(1-2)(1-3)}\ \ \ L_2(x)=\frac{(x-1)(x-3)}{(2-1)(2-3)}\ \ \ \L_3(x)=\frac{(x-1)(x-2)}{(3-1)(3-2)}
E quindi abbiamo P_3(x)=2*L_1(x)+3*L_2(x)+4*L_3(x) il nostro polinomio interpolante!

Il terzo metodo, concettualmente più complicato ma da un punto di vista operativo molto più conveniente.
Infatti se adesso volessimo il punto (4,5) ai nostri punti precedenti, dovremmo rifare tutti i conti per ottenere il polinomio interpolante. Con il metodo delle differenze divise i conti si riducono drasticamente se vogliamo aggiungere dei punti.

Sono un poco più complicate da spiegare e mi porterebbero via tanto spazio, rischiando di rendere poco leggibile il tutto, pertanto vi rimando a Wikipedia, purtroppo l’articolo in italiano non esiste, e purtroppo la spiegazione in inglese non è (a mio parere) molto buona, trovo quella in tedesco più comprensibile, qui l’algoritmo viene anche chiamato di Neville-Aitken ma temo che non tutti conosciate abbastanza questa lingua.

In ogni caso, il metodo di Lagrange sarà quello che useremo per interpolare un polinomio, vediamo un esempio!

Ho scritto le function lagr.m per il metodo di lagrange, e per chi volesse anche diffdiv_coeff.m, che servono per creare i polinomi interpolanti.
Vogliamo, ad esempio, interpolare con un polinomio di grado 5 la funzione sin(2*pi*x) nell’intervallo [-1,1].
Creiamo quindi uno script-file (considerate che da ora in poi qualsiasi esercizio lo dovete fare su uno script-file, per poterlo correggere, rileggere, salvare…) e cominciamo a scrivere:

f=inline('sin(2*pi*x)');
xnodi=linspace(-1,1,5+1); %ci servono n+1 punti per avere un polinomio di grado n!
ynodi=f(xnodi); %qui stiamo valutando la nostra funzione nei nodi d'interpolazione
x=linspace(-1,1);
p_lag = lagr(xnodi,ynodi,x);
y=f(x);
plot(x,y,x,p_lag)

Dopo aver eseguito lo script-file verrà visualizzato un grafico contenente le due funzioni. È auspicabile pensare che all’aumentare del grado n del polinomio aumenti la somiglianza tra le due funzioni. Provate infatti ad aumentare il numero di nodi a 11.

La funzione lagr.m si occupa di valutare il polinomio interpolanti sull’asse delle x, in sostanza quello che ci restituisce è un vettore con la valutazioni del polinomio interpolante.

Purtroppo però non è vero che prendendo nodi equispaziati (come li abbiamo scelti nell’esempio precedente) la funzione interpolante aumenti il grado di precisione all’aumentare del grado del polinomio stesso. Esiste infatti un famoso controesempio, il controesempio di Runge. Per visualizzarlo potete scaricare il file runge.m e provare ad eseguirlo con Octave.

Possiamo notare che al crescere di n le due funzioni non sono poi tanto simili e si può dimostrare che queste differenze diventano sempre più evidenti.

Risolto quindi il problema di come studiare una funzione generica attraverso un polinomio ci dobbiamo porre il problema di come scegliere i nodi d’interpolazione, in modo che al crescere del loro numero il polinomio interpolante assomigli alla funzione da studiare. Dal controesempio di Runge appare ovvio che scegliendo più nodi ai bordi si costringerebbe il polinomio interpolante a stare più vicino alla funzione da approssimare.

Si può dimostrare, che scegliendo i nodi di Tschebyscheff (valori che sono soluzione di un particolare polinomio), la funzione interpolante tenderà sempre al polinomio da studiare. Per fortuna questi nodi non sono difficili da calcolare, se stiamo lavorando nell’intervallo [-1,1] questi sono tutti nella forma: x_j=\cos((2_j-1)*\frac{\pi}{2(n+1)}), dove n è il grado di interpolazione e j assume i valori da 1 a n.

Se lavoriamo su un generico intervallo [a,b] basterà effettuare la trasposizione t_j=(b-a)*x_j/2+(b+a)/2 e valutare il problema con i nodi t_j.
Proviamo quindi a eseguire runge2.m. Avremo però bisogno di mettere tra le nostre function nodicheb.m che si occupa di creare i giusti nodi.
Possiamo osservare che all’aumentare dei nodi il polinomio interpolante diventa sempre più simile alla nostra funzione.

Ovviamente non è strettamente necessario usare le mie function per interpolare un polinomio, le ho messe a disposizione soltanto per scopo didattico. Octave integra già questa possibilità attraverso il comando polyfit(xnodi, ynodi, n), dove n è il grado del polinomio che vogliamo ottenere e xnodi, ynodi sono due vettori di lunghezza n+1.

A questo punto potete passare agli esercizi. Ovviamente potete interpolare anche altri tipi di funzioni!


Azioni

Informazione

Lascia un Commento

Fill in your details below or click an icon to log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Log Out / Modifica )

Foto Twitter

You are commenting using your Twitter account. Log Out / Modifica )

Foto di Facebook

You are commenting using your Facebook account. Log Out / Modifica )

Connecting to %s