Interpolazione polinomiale a tratti

1. Interpolazione polinomiale a tratti Gabriella Puppo

2. Interpolazione polinomiale a tratti • Interpolazione polinomiale a tratti funzione interp1 • Calcolo dell’errore • Interpolazione di funzioni non regolari • Esercizi

3. 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

4. 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

5. 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*')

6. Polinomio costante a tratti, N=10.

7. Polinomio lineare a tratti, N = 10;

8. Interpolazione spline, N = 10.

9. 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.

10. 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

11. 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

12. 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:

13. 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);

14. Funzione f(x) = exp(x)*cos(4x): >> err=andamento(f,[0,3],'linear'); >> [p,res]=min_quad(err) p = 1.9424 res = 0.2801 Per l’interpolante lineare a tratti, trovo: >> err=andamento(f,[0,3],'spline'); >> [p,res]=min_quad(err) p = 4.5473 res = 1.2200 mentre, per l’interpolante spline trovo:

15. 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)?

16. Interpolazione spline di una funzione a gradino Come posso migliorare l’interpolazione?

17. 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).

18. 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 valoridiversi

19. 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)