Approssimazione Di Funzioni

Calcolo numerico

1.

Seconda tesina – Approssimazione di funzioni

Problema

Si consideri la funzione 

f (x ) =

3x 2 + 2 1+ x 2 ⋅ (1+ 2/ 3x 2 )

 

x  X [a,b] = [-2,1]. Si vuole:  a.

Approssimare  f  in [a,b] con i polinomi interpolatori  p k ,  k  = 1,2,3,4, costruiti su  3,4,6,9 nodi di Chebychev(implementando la formula di Lagrange. 

b.

Approssimare  f   in  [a,b]  con  una  spline  cubica  interpolante  f   negli  stessi  nodi  di  p k  (utilizzando i moduli di matlab) 

Confrontare  i  grafici  dei  polinomi  considerati,  discutere  l’errore  totale  che  si  commette  nei vari casi. 

2.

Introduzione

Per  svolgere  le  operazioni  richieste  sono  stati  implementati  vari  programmi  in  Matlab, elencati di seguito, in appendice è presente il codice completo:  Funzione

Descrizione

fcalc(x)

Calcola la funzione data nei punti del vettore x 

f1calc(x)

Implementa la derivata della f(x) che calcola nei punti x 

chebspace(a,b,n) Restituisce un vettore di n nodi distribuiti in [a,b] secondo Chebyschev  lpoly(x)

Restituisce i polinomi di base di Lagrange, dati i nodi x 

Lag=lfitn(L,y)

Restituisce  la  matrice  Lag  dei  coefficienti  del  polinomio  interpolatore  di Lagrange, a partire dai polinomi di base di Lagrange L ed i valori y  assunti dalla funzione nei nodi. 

V=lebesgue(L) 

Restituisce il valore della funzione di Lebesgue 

Il  computer  utilizzato  ha  una  precisione  di  17  cifre,  pertanto  ciascun  valore  sarà  soggetto in primo luogo ad un errore di troncamento. In oltre la distanza minima tra due  -16 numeri dell’ordine dell’unità ottenuta con il comando eps, è pari a 2.2204·10 , per cui  l’errore  per  la  rappresentazione  di  qualsiasi  valore  dell’ordine  dell’unità  è  di  circa     -16 1,1·10 .  Le  ordinate  dei  nodi  saranno  perciò  affette  almeno  da  tale  errore,  se  non  sono  subentrati errore di ordine maggiore nel calcolo della funzione.  L’utilizzo di nodi approssimati comporta in primo luogo la presenza di un errore di propagazione, la cui entità è proporzionale alla funzione di Lebesgue Λ.  Poiché  n

Λ = ∑ L k ( x ) ,  k =0

Roberto Patrizi, 09114657

pag 1 di 8 

Calcolo numerico

Seconda tesina – Approssimazione di funzioni

è funzione esclusivamente dei nodi, è possibile minimizzare tale errore operando con nodi  distribuiti  secondo  Chebyshev.La  funzione  chebspace(a,b,n)  implementa  proprio  la  formula per il calcolo delle ascisse dei nodi, le rispettive ordinate si ottengono si ottengono con la funzione fcalc(x). 

2.1. Polinomi interpolatori di Lagrange In definitiva la sequenza di comandi:  %% Nodi di Chebichev nnodi=2; X=chebspace(-2,1,nnodi); Y=fcalc(X); %% Polinomi interpolatori di Lagrange L=lpoly(X); Lag=lfit(L,Y); %% Grafico curva interpolata yy=polyval(Lag,xx); xx=[-2:0.05:1]; yf=fcalc(xx); plot(xx,yf,'lineWidth',2,'color',[0,0,0]) plot(xx,yy,'k');

Fornisce i seguenti grafici al variare di nnodi: 

  Le interpolazioni con 3 e 4 punti si discostano fortemente dalla funzione, com’era  ovvio prevedere. Con 3 punti in particolare ho un polinomio interpolatore di grado 2, che  non può seguire la funzione nei suoi 3 cambi di direzione, presenti solamente per polinomi di grado almeno pari a 4. Infatti anche con 4 nodi la curva interpolata, di grado 3,  cresce  per  x>1,  mentre  la  funzione  decresce  per  x>1.  Con  6  nodi  ho  un  polinomio  di  Roberto Patrizi, 09114657

pag 2 di 8 

Calcolo numerico

Seconda tesina – Approssimazione di funzioni

grado  5  che  riesce  a  seguire  la  funzione,  mentre  con  9  nodi  si  comincia  ad  avere  una  buona riproduzione della funzione.  L’espressione analitica del polinomio interpolatore si può ottenere nei vari casi dai  due blocchi di istruzioni riportati nella pagina precedente, identificati con %%nodi di Chebichev e  %%Polinomi interpolatori di Lagrange. Si ottiene una variabile identificata con Lag ,che rappresenta la matrice dei coefficienti del polinomio interpolatore, da intendersi come: 

p( x ) = Lag (1)x n + Lag ( 2)x n −1 + ... + Lag ( n − 1)x + Lag ( n )   2.2. Interpolazione tramite Spline Matlab permette di costruire una spline e di specificarne i valori delle derivate agli  estremi. La derivata della funzione data è: 

f '(x ) = −

3x ( 6x 4 + 3x 2 − 4 ) ( x 2 + 1)3 ⋅ ( 2x 2 + 3)

 

Questa espressione è stata implementata nel file f1calc.m, che calcola il valora assunto dall’espressione data di f’(x).  La sequenza di istruzioni per ottenere le spline è la seguente  %% Nodi di Chebichev nnodi=2; X=chebspace(-2,1,nnodi); Y=fcalc(X); f1a=f1calc(X(nnodi+1)); %I nodi sono in ordine inverso f1b=f1calc(X(1)); %% Calcolo spline e funzione xx=[-2:0.05:1]; yf=fcalc(xx); ps=spline(X,[f1a,Y,f1b]); ys=ppval(ps,xx); %% Grafico funzione plot(xx,yf); plot(X,Y,'o'); plot(xx,yy);

Si ottengono i seguenti grafici 

 

Roberto Patrizi, 09114657

pag 3 di 8 

Calcolo numerico

Seconda tesina – Approssimazione di funzioni

  Si osserva una migliore accuratezza della riproduzione già con 4 nodi, grazie alla  presenza delle condizioni sulle derivate, che tuttavia non riescono comunque a riprodurre  la funzione con soli 3 nodi.  Anche la Spline con 3 nodi infatti non è altro che una parabola se non vengono  specificate le condizioni sulle derivate. I tre punti sono sufficienti a definire univocamente  i coefficienti dell’equazione di secondo grado che è identica sia per la spline che per i polinomi interpolatori di Lagrange. 

3.

Analisi dell’errore di interpolazione

L’errore di propagazione è trascurabile per il problema posto. Si ha infatti un nu-16 mero esiguo  di  nodi,  la  cui precisione  è  dell’ordine  di  10 .  L’errore  di  propagazione,  si  può stimare come prodotto tra la precisione dei nodi e la funzione di Lebesgue,   Si  osserva  che  il  Valore  della  funzione  di  Lebesgue  nell’intervallo  di  interesse,  non  supera  mai  le  due  unità,  per  cui  l’errore  di  propagazione  è  trascurabile,  come  era  ovvio  del  resto:  si  dispone  infatti  di  pochi  dati  ed  accurati.  L’errore  maggiore  è  dovuto  alla  differenza  tra  la  funzione  data  e  la  sua  approssimazione  polinomiale,  tanto  più  grande  quanto  minore  è  il  numero  di  nodi.  L’errore di troncamento è ottenibile, come  da  definizione,  dalla  differenza  tra  la  funzione e la sua approssimazione polinomiale:  e(x ) = f (x ) - p n (x )  Il calcolo dell’errore di troncamento in Matlab si effettua, dopo aver calcolato spline e polinomi interpolatori, con le seguenti righe di codice:  N=zeros(1,length(X)); xx=[-2:0.05:1]; yf=fcalc(xx); yl=polyval(Lag,xx); a=f1calc(X(nnodi+1)); b=f1calc(X(1)); pp=spline(X,[a,Y,b]); ysd=ppval(pp,xx); ys=spline(X,Y,xx);

Roberto Patrizi, 09114657

pag 4 di 8 

Calcolo numerico

Seconda tesina – Approssimazione di funzioni

%% Errore es=yf-ys; esd=yf-ysd; el=yf-yl;

Si ottengono i seguenti grafici: 

  Per  la  spline  è  stato  graficato  (la  linea  sottile  continua  indicata  con  Spline  Der)  l’errore  che si commette anche  nel  caso  in  cui  si specificano  le  condizioni  sulle  derivate  agli  estremi  dell’intervallo,  che,  come  si  osserva,  si  mantiene  generalmente  inferiore  rispetto all’errore sia per l’approssimazione tramite polinomi di Lagrange che per le spline  senza condizioni sulle derivate.  Per 3 nodi, solo la condizione sulla derivata garantisce una migliore approssimazione della funzione con la spline rispetto ai polinomi con base di Lagrange. Per 4 nodi il  miglioramento è notevole, infine aumentando il numero di nodi, la presenza delle condizioni  sulle  derivate  agli  estremi  dell’intervallo  provoca  modifiche  soltanto  agli  estremi  dell’intervallo,  mentre  all’interno  le  due  spline  sono  identiche.  Bastano  9  punti  -3 nell’intervallo considerato affinché l’errore si mantenga al di sotto dei 6 · 10 . 

Roberto Patrizi, 09114657

pag 5 di 8 

Calcolo numerico

4.

Seconda tesina – Approssimazione di funzioni

Funzioni Matlab

function Y=fcalc(X) % % function Y=fcalc(X) % f(x) = (3x^2 + 2) / ( sqrt(1+x^2) * (1+2/3 x^2) ) Y = (3*X.^2 + 2) ./ ( sqrt(1+X.^2) .* (1+2/3*X.^2) );

function Y=f1calc(X) % % Y=f1calc(x) % %Tale funzione è la derivata della f implementata in fcalc Y=-(3*X.*(6*X.^4 + 3*X.^2 - 4)) ./ ( ((X.^2+1).^(3/2)) .* ((2*X.^2+3).^2) );

function X=chebspace(a,b,n) % % X=chebspace(a,b,n) % %Restituisce un vettore di punti distribuiti secondo Chebyshev. %a,b sono gli estremi dell'intervallo, n è il numero di nodi -1. X=0:n; %è un vettore temporaneo che porta gli indici X= ((b-a)/2) * cos( (2*X+1) * pi ./ (2*(n+1)) ) + (b+a)/2; %NB: i nodi sono ordinati da b verso a!

Roberto Patrizi, 09114657

pag 6 di 8 

Calcolo numerico

Seconda tesina – Approssimazione di funzioni

function L=lpoly(X) % % L=lpoly(X) % %Funzione per il calcolo numerico dei polinomi interpolatori di Lagrange. In ingresso ho i nodi X. %In uscita il vettore L contiene i coefficienti dei polinomi di Lagrange % L(k)= [x-X(j)] / [X(k)-X(j)] per j=1..n+1, j diverso da k % n=length(X); L=[]; for k=1:n %Ciclo del calcolo dei vari L (L1,L2,L3...) temp=1; for j=1:n %Ciclo per il calcolo della base k-esima Lk if j~=k temp=conv(temp,[1,-X(j)]/(X(k)-X(j))); end end %temp contiene la base k-esima che voglio assegnare a L(k). Il polinomio è al %massimo di grado n, =>temp contiene al massimo n+1elementi, L deve avere n righe %ed n+1 colonne con i vettori temp %allineati a destra. L=accoda(L,temp); end function Lag=lfit(L,Y) % % Lag=lfit(L,Y) % %Restituisce il vettore Lag dei coefficienti del polinomio di Lagrange %interpolante i punti X,Y %In ingresso ho i polinomi di base di Lagrange L calcolati sui nodi X. n=length(Y); Lag=0; for k=1:n Lag=Lag+Y(k)*L(k,:); end

function V=lebesgue(L,xx) % % V=lebesgue(L,xx) % %Calcola la funzione di Lebesgue a partire dal vettore dei coefficienti dei %polinomi di base di Lagrange L, (ottenibili con la funzione L=lpolyn(nod)), %e ne calcola il valore nei punti indicati dal vettore xx. Il risultato è %un vettore V della stessa lunghezza di xx che contiene il valore assunto %dalla funzione di lebesgue in ciascun punto di xx. % n=length(xx); for i=1:n V(i)=0; for j=1:length(L) V(i)=V(i)+abs(polyval(L(j,:),xx(i))); end end

Roberto Patrizi, 09114657

pag 7 di 8 

Calcolo numerico

Seconda tesina – Approssimazione di funzioni

function C=accoda(A,B) % % C=accoda(A,B) % %A e B sono matrici non necessariamente con lo stesso numero di colonne. %La matrice C risultante è data dalla concatenazione di A,B, in modo tale %che se le matrici siano allineate a destra, ed eventualmente gli spazi %vuoti sulla sinistra sono colmati da zeri % [n,m]=size(A); j=length(B); if m==0

C=B; elseif m==j C=vertcat(A,B); elseif m>j B(:,j+1:m)=0; B=circshift(B,[0,m-j]); C=vertcat(A,B); elseif m<j A(:,m+1:j)=0; A=circshift(A,[0,j-m]); C=vertcat(A,B); else error('impossibile concatenare le matrici' ); end

Roberto Patrizi, 09114657

pag 8 di 8