interpolazione polinomiale n.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
Interpolazione polinomiale PowerPoint Presentation
Download Presentation
Interpolazione polinomiale

Loading in 2 Seconds...

play fullscreen
1 / 54

Interpolazione polinomiale - PowerPoint PPT Presentation


  • 177 Views
  • Uploaded on

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

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Interpolazione polinomiale' - amil


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
interpolazione polinomiale1
Interpolazione polinomiale
  • Matrice di Vandermonde
  • Costruzione del polinomio di interpolazione
  • Studio dell’errore
  • Fenomeno di Runge
  • Condizionamento
matrice di vandermonde
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

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

slide6

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

slide7

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

slide8

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

slide9

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

slide12

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

slide14

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

slide17

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.

slide18

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

esempio1
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

slide21

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

esercizio
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
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
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)

slide26

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

esempio2
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
slide28

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

slide29

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

slide31

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

slide32

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

slide33

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

Conclusioni

L’errore dipende:

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

studio della 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.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)

fenomeno di runge
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.

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

%

slide44

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

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

slide47

N = 25, delta = 1e-3

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

rand('state',25)

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

slide50

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.

slide51

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

slide52

I dati hanno questo aspetto:

Interpolando con un polinomio di grado 3 ottengo

slide53

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)

slide54

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: