Interpolazione polinomiale

1. Interpolazione polinomiale Gabriella Puppo

2. Interpolazione polinomiale • Matrice di Vandermonde • Costruzione del polinomio di interpolazione • Studio dell’errore • Fenomeno di Runge • Condizionamento

3. 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) = xij 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

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

5. 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 2.2000 2.6000 3.0000 3.4000 Columns 8 through 11 3.8000 4.2000 4.6000 5.0000

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

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

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

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

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

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

13. Come si vede, le matrici di Vandermonde sono mal condizionate

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

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

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

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

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

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

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

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

22. Ottengo:

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

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

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

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

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

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

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

30. Griglia equispaziata

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

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

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

34. dalla regolarità della funzione dalla disposizione dei nodi della griglia Conclusioni L’errore dipende: La disposizione dei nodi della griglia influenza l’errore attraverso la funzione nodale

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

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

37. Funzione nodale su x=[-1,1]

38. Funzione nodale su [-2,2]

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

40. Griglia equispaziata

41. Griglia di Gauss-Lobatto

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

43. 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 %

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

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

46. N = 20; delta = 1e-3

47. N = 25, delta = 1e-3 Qui e nel grafico seguente, la routine rand è stata inizializzata con rand('state',25)

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

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

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