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:b
x=
Columns 1 through 7
1.0000 1.4000 1.8000
Columns 8 through 11
3.8000 4.2000 4.6000
2.2000
5.0000
2.6000
3.0000
3.4000
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);
end
for 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);
end
for 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 equispaziata
e 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
Listato dello script con_vander.m
% 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' 2X2
for n=1:100
h=2/n;
v=-1:h:1;
a=vander(v);
c(n+1)=cond(a);
end
plot(log10(c))
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
Per costruire il polinomio di interpolazione utilizzando
le functions di Matlab, devo:
• 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
Function polyfit.m
La 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^3
Costruisco 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.m
La 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 della
funzione:
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’errore
Per 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 P
n=length(x)-1;
fx=feval(f,x);
p=polyfit(x,fx,n);
% Costruisce la griglia su cui calcolare l'errore
a=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
Se risolvo lo stesso problema usando la griglia di Gauss-Lobatto,
ottengo questi risultati:
5 1.22841e+01
10 5.86816e+00
20 1.82158e-03
40 1.13485e-11
Per la funzione f(x) = abs(x-1), con la griglia equispaziata:
5 4.85716e-01
10 1.31240e+00
20 7.88075e+00
40 3.75116e+06
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
L’errore dipende:
• dalla regolarità della funzione
• dalla disposizione dei nodi della griglia
La disposizione dei nodi della griglia influenza l’errore attraverso
la funzione nodale
Studio della funzione nodale
La funzione nodale ha un ruolo importante nell’andamento
dell’errore. Per studiare la funzione nodale, dobbiamo scrivere una
function che:
• 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.
Function nodale.m
function nodale(x)
% NODALE(X) Calcola il grafico della funzione nodale
%
definita su una griglia X
n=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));
end
end
plot(xg,om)
Funzione nodale su x=[-1,1]
Funzione nodale su [-2,2]
Fenomeno di Runge
Interpolo 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
Per studiare il condizionamento di un problema di interpolazione:
• 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
Function condiz_pol.m
function 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 polinomi
a=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
Esempio
Studio 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 equispaziata
e interpolando nel senso dei minimi quadrati con un polinomio di
grado 5, ottengo: