Interpolazione polinomiale

Gabriella Puppo

Interpolazione polinomiale

• Matrice di Vandermonde

• Costruzione del polinomio di interpolazione

• Studio dell’errore

• Fenomeno di Runge

• Condizionamento

Matrice di Vandermonde

La matrice di Vandermonde è la matrice dei coefficienti di un problema di interpolazione, e dipende solamente dalla disposizione dei nodi di griglia, X0, X1, …, Xm.

E’ la matrice (m+1) per (N+1) definita da:

A(i,j) = xi j

La matrice di Vandermonde quadrata può essere costruita con la function vander(v), dove v è un vettore che contiene i nodi di griglia. Il risultato è una matrice che ha le colonne disposte in modo opposto rispetto alle convenzioni solite

Esempio

Costruisce la matrice 3 X 3, relativa ai nodi:X_0 = 2, X_1=3, X_2=4:

>> vander([2,3,4])ans = 4 2 1 9 3 1 16 4 1

Griglie di interpolazione

1) Griglia equispaziata, sull’intervallo [a,b] con m+1 nodi:Esempio con a=1, b=5, m=10 (numero intervalli):

>> a=1; b=5; m=10;>> h=(b-a)/m;>> x=a:h:bx = Columns 1 through 7 1.0000 1.4000 1.8000 2.2000 2.6000 3.0000 3.4000 Columns 8 through 11 3.8000 4.2000 4.6000 5.0000

Oppure, sempre per creare una griglia equispaziata, posso usare il comando linspace:

LINSPACE Linearly spaced vector. LINSPACE(x1, x2) generates a row vector of 100 linearly equally spaced points between x1 and x2. LINSPACE(x1, x2, N) generates N points between x1 and x2.

Esempio, con a=1, b=5, m=11 (numero nodi):

>> a=1; b=5; m=11;>> x=linspace(a,b,m)x = Columns 1 through 8 1.0000 1.4000 1.8000 2.2000 2.6000 3.0000 3.4000 3.8000 Columns 9 through 11 4.2000 4.6000 5.0000

function x=x_cheb(n,ab)% X_CHEB(N,[A,B]): Calcola la griglia di Chebichev con N+1% punti, sull'intervallo [A,B]% X_CHEB(N) Calcola la griglia di Chebichev su [-1,1]if nargin==1 a=-1; b=1;else a=ab(1); b=ab(2);endfor i=0:n x(i+1)= (a+b)/2 -(b-a)/2*cos(pi*(2*i+1)/(2*(n+1)));end

function x=x_gauss(n,ab)% X_GAUSS(N,[A,B]): Calcola la griglia di Gauss Lobatto con N+1% punti, sull'intervallo [A,B]% X_GAUSS(N) Calcola la griglia di Gauss Lobatto su [-1,1]if nargin==1 a=-1; b=1;else a=ab(1); b=ab(2);endfor j=0:n x(j+1)= (a+b)/2 -(b-a)/2*cos(pi*j/n);end

Esempio: Costruiamo la griglia di Gauss-Lobatto sull’intervallo[1,5] con 11 punti:

>> x_gauss(10,[1,5])ans = Columns 1 through 7 1.0000 1.0979 1.3820 1.8244 2.3820 3.0000 3.6180 Columns 8 through 11 4.1756 4.6180 4.9021 5.0000

Studiamo ora la disposizione dei nodi per la griglia equispaziatae per la griglia di Chebishev, per m=10 e per m=20:

Studio del condizionamento delle matrici di Vandermonde

Costruisco le matrici di Vandermonde, corrispondenti ad una griglia equispaziata contenente N+1 punti, per N che va da 1 a 100 e ne stimo il numero di condizionamento con la function cond per ogni valore di N.Poi stampo un grafico con i risultati

% Studio della stima numerica del condizionamento della % matrice di Vandermonde costruita su un vettore V di nodi% uniformemente distribuiti su [-1,1]%c(1)=1;% La piu' piccola matrice che costruisce e' 2X2for n=1:100 h=2/n; v=-1:h:1; a=vander(v); c(n+1)=cond(a);endplot(log10(c))

Listato dello script con_vander.m

Come si vede, le matrici di Vandermonde sono mal condizionate

Esercizio: Stimare numericamente il rango delle matrici di Vandermonde, su una griglia equispaziata, per N che va da 1 a 100.

Costruzione del polinomio di interpolazione

• Chiamare la function polyfit, per calcolare i coefficienti del polinomio di interpolazione

• Chiamare la function polyval per valutare il polinomio di interpolazione in corrispondenza dei punti richiesti

Per costruire il polinomio di interpolazione utilizzando le functions di Matlab, devo:

Function polyfit.mLa funzione polyfit calcola i coefficienti del polinomio di interpolazione. La chiamata è:

>> p=polyfit(x,fx,n)

dove:x è il vettore che contiene i nodi di griglia;fx è il vettore che contiene i valori da interpolare sui nodi di griglia;n è il grado del polinomio di interpolazione.

N.B. Se n = length(x)-1, allora il polinomio di interpolazione soddisfa esattamente le condizioni di interpolazione.Altrimenti, se n < length(x)-1, polyfit interpola nel senso dei minimi quadrati.

Attenzione! polyfit fornisce i coefficienti in un ordine diverso da quello usuale:Supponiamo che la chiamata di polyfit sia stata

>> a=polyfit(x,fx,n)

Allora il polinomio p(x) è dato da:p(x) = a(1)* x n + a(2)* x n-1 + … + a(n) * x + a(n+1),cioè i coefficienti vengono forniti in ordine inverso rispetto al solito e scalati di 1, in modo da evitare indici uguali a zero,infatti in Matlab tutti i vettori devono avere indici interi e positivi.

Esempio.Considero la funzione f(x) = x^3Costruisco una griglia equispaziata con 4 punti, ed interpolo con il polinomio di grado 3 p(x). In questo caso, l’interpolazione è esatta, cioè f(x) = p(x).I coefficienti del polinomio, usando la notazione solita sarebbero:p = [0 0 0 1].Poiché Matlab li ordina in maniera opposta, troveremo:

>> f=inline('x.^3');>> x=-1:2/3:1;>> fx=f(x);>> p=polyfit(x,fx,3)p = 1.0000 0.0000 0.0000 0

Function polyval.mLa function polyval valuta il polinomio definito dal vettore di coefficienti a sull’insieme dei valori definito dal vettore x:

>> pxx=polyval(a,xx);

Quindi per calcolare il valore che il polinomio:

P(x) = 2*x^2 -2*x -1

assume nei punti x=0 e x=1, posso dare il comando:

>>y=polyval([2 -2 -1],[0 1]) y = -1 -1

Esempio

Calcolare il polinomio di interpolazione di grado N=5 dellafunzione:f(x) = exp(x) * sin(2x),usando una griglia equispaziata sull’intervallo [0,2].

Visualizzare i risultati, riportando sullo stesso grafico:- la funzione f(x)- il polinomio di interpolazione p(x)- i punti di interpolazione

Posso ottenere i risultati richiesti usando i comandi:

>> f=inline('exp(x).*sin(2*x)');>> n=5;>> h=2/n;>> x=0:h:2; % Costruisce la griglia con 6 nodi>> fx=f(x); % Calcola f sui nodi di griglia>> p=polyfit(x,fx,n); % Costruisce il polinomio di grado 5>> xx=0:0.01:2; % Costruisce una griglia grafica>> pxx=polyval(p,xx); % Valuta il polinomio sulla griglia xx>> plot(xx,f(xx))>> hold;>> plot(xx,pxx,'g') % Disegna il grafico del polinomio>> plot(x,fx,'r*') % Disegna i punti di interpolazione

Ottengo:

Esercizio

Interpolare la funzione f(x) = x sin(x) sull’intervallo [0,4] con un polinomio di grado 8 e stampare un grafico che contenga:- Il grafico di f;- Il grafico del polinomio;- I punti di interpolazione.Ottenere i risultati richiesti usando:- La griglia equispaziata;- La griglia di Chebishev

Studio dell’errorePer stimare l’errore di interpolazione fra una funzione f(x) ed il suo polinomio di interpolazione P(x), di grado N, costruito sulla griglia di interpolazione X, devo:

- Calcolare il polinomio di interpolazione P(x)

- Preparare una griglia XX più fitta di quella di interpolazione

- Calcolare sia f(x) che P(x) sulla griglia XX

- Calcolare la norma infinito del vettore f(XX) - P(XX).

Function errore_pol.m

La function errore_pol ha due modalità di funzionamento: - se voglio solo il calcolo dell’errore, chiamo: errore_pol(f,x), dove f deve essere una stringa e x contiene le ascisse dei punti da interpolare

- se voglio sia il calcolo dell’errore che il grafico di f e del polinomio (con in più i punti di interpolazione), chiamo: errore_pol(f,x,1)

function err= errore_pol(f,x,dummy)% ERR=ERRORE_POL(F,X) Calcola l'errore nella norma infinito% fra F ed il polinomio P di interpolazione di F,% costruito sulla griglia X% Se ci sono tre argomenti in input, stampa anche% il grafico di F e di Pn=length(x)-1;fx=feval(f,x);p=polyfit(x,fx,n);% Costruisce la griglia su cui calcolare l'errorea=x(1); b=x(n+1); h=(b-a)/200; xx=a:h:b;fxx=feval(f,xx);pxx=polyval(p,xx);err = norm(fxx-pxx,inf);if nargin==3 % Stampa il grafico di F(XX), P(XX) e dei punti di % interpolazione plot(xx,fxx), hold on; plot(xx,pxx,'g') plot(x,fx,'r*')end

Esempio

• Calcolare l’errore di interpolazione per la funzione f(x)=exp(x)*cos(4x), su [-3,3], usando una griglia equispaziata con N=5, 10, 20 e 40 nodi.

• Stessa cosa per la funzione f(x) = abs(x-1)

• Rifare i conti per la griglia di Gauss-Lobatto

Per calcolare l’errore di interpolazione per la funzione f(x)=exp(x)*cos(4x), su [-3,3], usando una griglia equispaziata con N=5, 10, 20 e 40 nodi, e visualizzare i risultati posso usare questo script:

% Stampa i grafici di confronto fra funzione e% polinomio di interpolazione su [a,b]a=-3; b=3;f=inline('exp(x).*cos(4*x)');k=0;for n=[5,10,20,40] x=linspace(a,b,n); k=k+1; subplot(2,2,k) deg=sprintf('Polinomio di grado %d',n-1) err(k)=errore_pol(f,x,1) title(deg)end

…continua...

k=0;for n=[5,10,20,40] k=k+1; fprintf('%2.0f %10.5e \n',n,err(k))end

Ottengo questi risultati:

5 2.13962e+01 10 3.13726e+01 20 2.78234e-01 40 1.83374e-06

Quindi in questo caso l’errore di interpolazione diminuisce, all’aumentare di n, anche se ricevo dei warnings, perché la matrice di Vandermonde è mal condizionata

Griglia equispaziata

5 1.22841e+01 10 5.86816e+00 20 1.82158e-03 40 1.13485e-11

Se risolvo lo stesso problema usando la griglia di Gauss-Lobatto, ottengo questi risultati:

5 4.85716e-01 10 1.31240e+00 20 7.88075e+00 40 3.75116e+06

Per la funzione f(x) = abs(x-1), con la griglia equispaziata:

Se uso la griglia di Gauss-Lobatto, i risultati migliorano, ma la velocità di convergenza è molto più bassa di prima

5 6.81282e-01 10 3.03245e-01 20 1.37322e-01 40 4.78853e-02

Conclusioni

• dalla regolarità della funzione

• dalla disposizione dei nodi della griglia

L’errore dipende:

La disposizione dei nodi della griglia influenza l’errore attraverso la funzione nodale

Studio della funzione nodale

• Accetti in input il vettore X contenente le ascisse della griglia.

• Imposti una griglia grafica, sulla quale disegnare il grafico della funzione nodale.

• Calcoli il valore della funzione nodale su ogni punto della griglia grafica.

• Disegni il grafico.

La funzione nodale ha un ruolo importante nell’andamento dell’errore. Per studiare la funzione nodale, dobbiamo scrivere una function che:

Function nodale.m

function nodale(x)% NODALE(X) Calcola il grafico della funzione nodale% definita su una griglia Xn=length(x);a=x(1); b=x(n); h=(b-a)/200% Costruisce la griglia grafica XG xg=a:h:b; np=length(xg);for jp=1:np om(jp)=1; for i=1:n om(jp)=om(jp)*(xg(jp)-x(i)); endendplot(xg,om)

Funzione nodale su x=[-1,1]

Funzione nodale su [-2,2]

Fenomeno di RungeInterpolo la funzione f(x) = 1 / (1 +x^2)su intervalli del tipo [-a,a], con a=1, a=3, a=5,usando polinomi di grado N = 10 e N = 20.

Considero:1) griglie equispaziate;2) griglie di Gauss-Lobatto.

Griglia equispaziata

Griglia di Gauss-Lobatto

Condizionamento

• Definisco una funzione f(x) e una griglia X

• Costruisco il polinomio P(x) che interpola i dati (X, f(X))

• Disegno il grafico di f(x) e di P(x)

• Perturbo i dati f(X), per esempio F(X)=f(X) +10K * r K=-3, -4, ecc. (dove r è un vettore di numeri casuali)

• Calcolo il polinomio di interpolazione dei dati perturbati e ne disegno il grafico

Per studiare il condizionamento di un problema di interpolazione:

Function condiz_pol.mfunction deltap = condiz_pol(f,x,delta,dummy)% CONDIZ_POL(F,X,DELTA,DUMMY) Studia la variazione del polinomio% di interpolazione rispetto a variazioni nei dati in ingresso% F e' la funzione esatta (stringa di caratteri);% X e' la griglia di interpolazione;% DELTA e' l'ampiezza della perturbazione (che viene costruita% con deltaF = DELTA*(RAND(1,N+1)-0.5)% DUMMY e' un segnaposto: se c'e' vengono prodotti i grafici% di F, del polinomio con i dati esatti e del polinomio% con i dati perturbati% DELTAP (in output) e' la norma della differenza relativa fra il% polinomio esatto e quello perturbato%

n=length(x)-1;fx=feval(f,x);p=polyfit(x,fx,n);% Costruisce la griglia su cui calcola la differenza fra i% due polinomia=x(1); b=x(n+1); h=(b-a)/200;xx=a:h:b;pxx=polyval(p,xx);% perturba i dati fx e costruisce il polinomio perturbato:fxp=fx+delta*rand(1,n+1);p=polyfit(x,fxp,n);ppxx=polyval(p,xx);deltap = norm(pxx-ppxx,inf)/norm(pxx,inf);if nargin == 4 fxx=feval(f,xx); plot(xx,fxx) hold on plot(xx,pxx,'g') plot(x,fx,'r*') plot(xx,ppxx,'m')end

EsempioStudio la propagazione di errori nei dati per il polinomio che interpola la funzione

f(x) = exp(x) * cos( 4*x) su [0,3]

L’ampiezza della perturbazione nei dati è delta. Il grado del polinomio di interpolazione è N.Nei grafici che seguono, la funzione è disegnata in blu, il polinomio che interpola i dati esatti è in verde, mentre il polinomio che interpola i dati perturbati è in magenta.

N = 20; delta = 1e-3

N = 25, delta = 1e-3

Qui e nel grafico seguente, la routine rand è stata inizializzata con

rand('state',25)

Con la griglia di Gauss-Lobatto: N = 25, delta = 1e-3

Minimi quadrati

Per interpolare dei dati con un polinomio nel senso dei minimi quadrati, uso un numero M+1 di dati, con M > N, dove N è il grado del polinomio di interpolazione.Ottengo un sistema sovradeterminato, che si risolve con il metodo QR.Notare che ora le condizioni di interpolazione non sono soddisfatte esattamente.Per studiare l’interpolazione nel senso dei minimi quadrati, ci appoggeremo alle functions di Matlab.

Il polinomio che si ottiene con i minimi quadrati puo’ essere parecchio lontano dalla funzione che si vuole approssimare.

Qui la funzione f(x) = exp(x)*cos(4x) è stata interpolata con un polinomio di grado 5, utilizzando 11 nodi equispaziati.

L’ interpolazione nel senso dei minimi quadrati è utile quando si ha un gran numero di dati, o quando la funzione da interpolare non è regolare.

Supponiamo di avere assegnato un fenomeno che è dato da una legge incognita. Abbiamo a disposizione un elevato numero di valori sperimentali, soggetti ad un errore.

Posso simulare una situazione di questo tipo sovrapponendo una perturbazione casuale ad una funzione, per esempio:

>> f=inline('x.^3 -3.*x +2');>> fx=f(x)+0.2*rand(1,100);

I dati hanno questo aspetto:

Interpolando con un polinomio di grado 3 ottengo

Interpolando la funzione f(x) = abs(x-1), con un polinomio di grado 10, che interpola f esattamente, ottengo:

Qui le oscillazioni sono dovute alla non regolarità di f(x)

Il risultato migliora se uso l’interpolazione nel senso dei minimi quadrati.Usando 11 punti di interpolazione su una griglia equispaziatae interpolando nel senso dei minimi quadrati con un polinomio di grado 5, ottengo: