Interpolazione polinomiale a tratti

Interpolazione polinomiale
a tratti
Gabriella Puppo
Interpolazione polinomiale
a tratti
• Interpolazione polinomiale a tratti
funzione interp1
• Calcolo dell’errore
• Interpolazione di funzioni non regolari
• Esercizi
Interpolazione polinomiale a tratti
Per costruire l’interpolante polinomiale a tratti di una funzione f(x)
su un intervallo [a,b], devo fare i seguenti passi:
• Costruire una griglia X su [a,b]
• Valutare f sulla griglia X
• Calcolare i punti XX sui quali si vuole conoscere il
valore dell’interpolante
• Chiamare la function di interpolazione
polinomiale a tratti
Function interp1.m
La routine di interpolazione polinomiale a tratti di Matlab si
chiama interp1. La sua sintassi è la seguente:
yy = interp1( X, FX, XX, TIPO), dove
X è la griglia di interpolazione;
FX sono i valori da interpolare;
XX sono i punti sui quali si vuole calcolare l’interpolante
TIPO è il tipo di interpolazione richiesta:
- ‘nearest’ : interpolazione con polinomio costante a tratti
- ‘linear’ : (default) polinomio lineare a tratti
- ‘cubic’ : polinomio cubico a tratti (con derivate continue)
- ‘spline’ : interpolazione con spline cubica
Esempio
Consideriamo la funzione f(x) = exp(x) * cos(4x) su [0,3].
Calcoliamo l’interpolante polinomiale a tratti I(x) su una griglia
uniforme di 11 punti.
Disegniamo un grafico di f(x) e di I(x).
Per esempio, per il polinomio costante a tratti:
>> f=inline('exp(x).*cos(4*x)');
>> x=0:3/10:3;
>> fx=f(x);
>> xx=0:0.01:3;
>> fxx=f(xx);
>> i0=interp1(x,fx,xx,'nearest');
>> plot(xx,fxx); hold on, plot(xx,i0,'g')
>> plot(x,fx,'r*')
Polinomio costante a tratti, N=10.
Polinomio lineare a tratti, N = 10;
Interpolazione spline, N = 10.
Errore
Per studiare l’andamento dell’errore per l’interpolazione
polinomiale a tratti di una funzione f , devo:
• Costruire una griglia di interpolazione uniforme X, con
intervalli di ampiezza h;
• Calcolare l’interpolante Ih polinomiale a tratti basato
sulla griglia X;
• Valutare l’interpolante Ih(x) appena trovato su una griglia
XX, molto più fitta di X;
• Calcolare l’errore: ||f(XX) - Ih (XX)|| ;
• Ripetere la procedura dimezzando ogni volta h e stampare
i risultati
• Cercare di ricostruire l’andamento dell’errore in funzione
di h.
function err=andamento(f,ab,s)
% ERR=ANDAMENTO(F,AB,S) Calcola l'errore dovuto all'interpolazione
%
polinomiale a tratti per la funzione F sull'intervallo
%
[A,B]=AB, per l'interpolazione di tipo S
% S è una stringa e può assumere i valori: 'nearest', 'linear'
%
'spline' o 'cubic'
%
La function parte da una griglia con 5 intervalli,
%
e calcola l'errore in norma inf, dimezzando l'ampiezza
%
degli intervalli
n=5;
a=ab(1); b=ab(2); dx=(b-a)/200;
xx=a:dx:b; fxx=feval(f,xx);
h=(b-a)/n;
for k=1:7
x=a:h:b; fx=feval(f,x);
pi=interp1(x,fx,xx,s);
err(k) = norm(pi-fxx,inf);
h=h/2;
end
Funzione f(x) = exp(x)*cos( 4x), su [0,3]:
>> err=andamento(f,[0,3],'nearest')
err =
19.8408 10.1467 4.8818 2.1173 1.0578 0.5291 0.2649
>> err=andamento(f,[0,3],'linear')
err =
6.1800 2.2747 0.5641 0.1401 0.0353 0.0088 0.0022
>> format short e
>> err=andamento(f,[0,3],'spline')
err =
Columns 1 through 5
5.4764e+000 8.7790e-001 3.9363e-002 1.3420e-003 4.0357e-005
Columns 6 through 7
1.1925e-006 7.4154e-008
Calcolo l’andamento dell’errore
Posso verificare che l’andamento dell’errore è del tipo:
e(h) = C h 
Per far questo, passo ai logaritmi:
ln(e) = c +  ln(h)
Ottengo un problema ai minimi quadrati, nel quale  è il
coefficiente della retta incognita:
[, c] = polyfit( ln(h), ln(e), 1)
Adattando questa idea all’output della function andamento.m,
ottengo il programma seguente:
function min_quad
function [p,res]=min_quad(err)
% [P,RES]=MIN_QUAD(ERR) calcola la pendenza (P) e il residuo
% (RES) nella retta di regressione lineare per il problema
% log(ERR) = C + P*log(H)
% Parte con un valore iniziale H=1 (il valore iniziale non e' importante):
h(1)=1; k=length(err);
for i = 2:k
h(i) = h(i-1)/2;
end
x=log(h); fx=log(err);
a=polyfit(x,fx,1);
% Interpola con un polinomio di grado 1
p=a(1);
retta=polyval(a,x);
% Calcola i valori della retta di regressione in x
res=norm(fx-retta);
Funzione f(x) = exp(x)*cos(4x):
Per l’interpolante
lineare a tratti, trovo:
mentre, per
l’interpolante spline
trovo:
>> err=andamento(f,[0,3],'linear');
>> [p,res]=min_quad(err)
p=
1.9424
res =
0.2801
>> err=andamento(f,[0,3],'spline');
>> [p,res]=min_quad(err)
p=
4.5473
res =
1.2200
Esercizio
Valutare l’andamento dell’errore rispetto ad H per
l’interpolazione
polinomiale a tratti sull’intervallo [0,3][, nel caso delle funzioni:
1) f(x) = abs(x - 1.1)
2) f(x) = 0, se x >= sqrt(2); f(x) = 1, se x < sqrt(2) (Questa
function è disponibile nell’archivio come function step(x))
Che risultati si osservano rispetto alla funzione
f(x)=exp(x)*cos(4x)?
Interpolazione spline di una funzione a gradino
Come posso migliorare l’interpolazione?
Interpolazione nel senso dei minimi quadrati con un polinomio di
grado 10:
Devo usare un interpolante diverso su ogni intervallo di regolarità
di f(x).
Esercizio
Considerare la funzione
f(x) = x^2 * sin(3x) su [0,3]
Definire su [0,3] una griglia equispaziata X di 21 punti.
Interpolare la funzione data nel senso dei minimi quadrati sulla
griglia assegnata con polinomi di grado 5, 10, 15, 20.
Calcolare ogni volta la norma 2 del residuo (f(X) - P(X)).
Discutere i risultati ottenuti
Suggerimento: La ||f(X)-P(X)|| = 0 solo nel caso di interpolazione
polinomiale esatta: negli altri casi sui punti della griglia X di
interpolazione, il polinomio e la funzione assumono valori diversi
Esercizio
Considerare la funzione
f(x) = x^2 * sin(3x) su [0,3]
Definire su [0,3] una griglia equispaziata X di 11 punti.
Interpolare la funzione data nel senso dei minimi quadrati sulla
griglia assegnata con un polinomio di grado 4.
Verificare che si trova lo stesso polinomio impostando il sistema
lineare con il metodo delle equazioni normali
Suggerimento: Per il metodo delle equazioni normali, calcolo
esplicitamente la matrice V di interpolazione 11 per 5 data da V(i,j)
= x(i)^(5-j), con i = 1…11, j=1,…5, dove ho usato lo stesso
ordinamento di Matlab per le colonne di V, poi risolvo
V’Vp=V’f(X)