Interpolazione polinomiale

11 dicembre 2007

1 Interpolazione polinomiale

Siano dati 1 punti _{0} _{}a due a due distinti e in ordine crescente (cioé 1), e i valori _{0} _{}ivi assunti da una funzione ( ).

Il problema dell’interpolazione polinomiale (cf. [1, p.131], [8, p.289]) consiste nel calcolare il polinomio _{}di grado tale che

( _{}) _{}0  (1)

Perché questo problema abbia senso, necessita garantire l’esistenza ed unicitá di tale polinomio _{}, questione che puó essere provata in vari modi (ad esempio in termini della non singolaritá della matrice di Vandermonde [18]).

Nel caso 1, il problema diventa relativamente semplice. Non è altro che il calcolo della retta che passa per due punti assegnati _{0} ( _{0} _{0}) e _{1} ( _{1} _{1}) cioė da

\frac{_{0}}{_{1} _{0}} \frac{_{0}}{_{1} _{0}}  (2)

abbiamo facilmente

0 \frac{_{1}}{_{0} _{1}}  _{1} \frac{_{0}}{_{1} _{0}}  (3)

Notiamo subito che i polinomi

0( ) : \frac{_{1}}{_{0} _{1}}  (4)

1( ) : \frac{_{0}}{_{1} _{0}}  (5)

hanno la particolarità che ( _{}) _{}

dove al solito _{}è il delta di Kronecker , cioè

1 _{}0 per

Da (3), (4), (5) si ha ovviamente

( ) _{0} _{0}( ) _{1} _{0}( )  (6)

L’idea dei polinomi di Lagrange è di estendere questa proprietà al caso 1. Lagrange (cf. [12], [5]) osservò che il polinomio di grado

( ) : _{0}^{}\frac{_{}}{_{}_{}}\frac{( _{0}) ( _{1} )( _{1} ) ( _{})}{( _{}_{0}) ( _{}_{1} )( _{}_{1} ) ( _{}_{})}  (7)

è tale che _{}( _{}) _{}. Quindi si vede facilmente che

( ) _{0}^{}_{}_{}( )  (8)

Figura 1: Carl Runge (1856-1927) e Giuseppe Luigi Lagrange (1736-1813).

Esempio Calcoliamo il polinomio di grado 2 che assume nei nodi _{0} 2, _{1} 1, _{2} 3 rispettivamente i valori _{0} 2, _{1} 11, _{2} 17. Come si vede dalla definizione, i polinomi di Lagrange sono

0( )  \frac{( _{1})( _{2})}{( _{0} _{1})( _{0} _{2})} \frac{( 1)( 3)}{( 2 1)( 2 3)} \frac{( 1)( 3)}{15}
1( )  \frac{( _{0})( _{2})}{( _{1} _{0})( _{1} _{2})} \frac{( ( 2))( 3)}{(1 ( 2))(1 3)} \frac{( 2)( 3)}{6}
2( )  \frac{( _{0})( _{1})}{( _{2} _{0})( _{2} _{1})} \frac{( ( 2))( 1)}{(3 ( 2))(3 1)} \frac{( 2)( 1)}{10}

Si vede subito che i polinomi _{0}, _{1}, _{2} sono di secondo grado (e quindi tali le loro combinazioni lineari), che

0( _{0}) _{0}( 2) 1 _{0}( _{1}) _{0}(1) 0 _{0}( _{2}) _{0}(3) 0
1( _{1}) _{1}(1) 1 _{1}( _{0}) _{1}( 2) 0 _{1}( _{2}) _{1}(3) 0
2( _{2}) _{2}(3) 1 _{2}( _{0}) _{2}( 2) 0 _{2}( _{1}) _{2}(1) 0

e che il polinomio definito in (8)

2( ) 2 _{0}( ) 11 _{1}( ) 17 _{2}( )

è tale che

( _{0}) 2 _{0}( _{0}) 11 _{1}( _{0}) 17 _{2}( _{0}) 2 1 11 0 17 0 2
( _{1}) 2 _{0}( _{1}) 11 _{1}( _{1}) 17 _{2}( _{1}) 2 0 11 1 17 0 11
( _{2}) 2 _{0}( _{2}) 11 _{1}( _{2}) 17 _{2}( _{2}) 2 0 11 0 17 1 17

e quindi è proprio il polinomio interpolante cercato.

Dal punto di vista pratico, (8) ha alcuni problemi. Se dopo aver calcolato il polinomio _{}interpolante in 1 punti ( _{0} _{0}), , ( _{}_{}), desideriamo ottenere il polinomio _{1} interpolante in 2 punti ( _{0} _{0}), , ( _{}_{}), la formula (8) è inefficiente poichè bisogna ricalcolare tutti i polinomi di Lagrange. Fortunatamente, esistono altre maniere di esprimere il polinomio interpolatore, come quella di Newton (cf. [1, p.138], [8, p.294], [13], [11]), che non soffrono di questo problema.

Per quanto concerne l’errore, se la funzione è di classe ^{1} , si puó vedere dal teorema del resto che

( ) _{}( ) ^{( 1)} ( )\frac{_{0}^{}( _{})}{( 1)!}  (9)
dove  con il più piccolo intervallo aperto contente _{0}  _{}.

Ma come scegliere (avendone la possibilità ) i punti in cui interpolare la funzione : [ ] �? Mostriamo due casi notevoli:

1. nodi equispaziati: fissato , i punti sono

\frac{( )}{}0  ;  (10)

2. nodi di Chebyshev-Gauss (scalati) [8, p.294]: fissato , i punti sono

\frac{( )}{2} \frac{( )}{2} _{}0  (11)

con

cos \frac{2 1}{2 2}  0  ;  (12)

3. nodi di Chebyshev-Gauss-Lobatto (scalati) [8, p.295]: fissato , i punti sono

\frac{( )}{2} \frac{( )}{2} _{}0  (13)

con

cos \frac{}{}0  ;  (14)

Figura 2: Pafnuty Lvovich Chebyshev (1821-1894) e Rehuel Lobatto (1797-1866).

In seguito useremo i nodi di Chebyshev-Gauss-Lobatto, che a differenza dei nodi di Chebyshev (cf. [14], [4]) includono gli estremi dell’intervallo. Si consiglia di effettuare gli esercizi successivi sostituendo ai nodi di Chebyshev-Gauss-Lobatto (scalati), quelli di Chebyshev-Gauss.

Viene da credere che aumentando , nel caso dei nodi equispaziati, si abbia che ( ) ( ) per [ ]. Come provato da Runge ciò non è vero. Infatti, per

( ) \frac{1}{1 ^{2}}  [ 5 5]  (15)

si ha che _{}non converge (puntualmente) a . Fortunatamente ciò non succede per i nodi di Chebyshev-Gauss, per cui comunque (per un teorema dovuto a Faber) esistono funzioni continue ma non ^{1} tali che l’interpolante _{}non converge puntualmente a .

2 Esercizi sull’interpolazione polinomiale

1. Si implementi una function che calcola i nodi di Chebyshev-Gauss per un intervallo qualsiasi.

2. Si implementi poi una function che interpoli polinomialmente una funzione che assume valori su un vettore di nodi che vengono forniti in input. Tale funzione Matlab deve inoltre calcolare i valori che l’interpolante polinomiale _{}assume nei nodi test . Si testi il codice, producendo dei grafici, sulle funzioni

• ( ) \frac{1}{1 ^{2}} in [ 5 5].

usando nodi equispaziati (con alcuni valori di tra 6 e 11) e nodi di Chebyshev-Gauss-Lobatto (con i medesimi valori di ), calcolando gli errori in norma infinito tra la funzione e le interpolanti. In particolare, si mostri (mediante grafici e/o valori dell’errore) che per il controesempio di Runge (cf. [9], [16]) l’aumento del numero di nodi equispaziati non migliora la ricostruzione della funzione da parte dell’interpolante.

Figura 3: Grafico che illustra l’interpolazione su nodi equispaziati (la funzione ha la linea continua, il polinomio interpolatore è tratteggiato, i nodi di interpolazioni coi cerchietti).

Figura 4: Grafico della funzione di Runge 1/(1 ^{2}) nell’intervallo [ 5 5].

2.1 Traccia di risoluzione

Si scrivano le functions cheb.m, interpol.m, runge.m definite rispettivamente come

function xc=cheb(a,b,n)
for m=1:1:n
    xc(m)=(a+b)/2-((b-a)/2)*cos(pi*(m-1)/(n-1));
end

function [t]=interpol(x,y,s)
%----------------------------------------------
% Interpolazione
%
% In input:
% x: nodi.
% y: valori nei nodi.
% s: nodi su cui calcolare l'interpolante.
%
% In output:
% t: valori dell'interpolante o approssimante.
%----------------------------------------------
m=length(x)-1;
coeff=polyfit(x,y,m);
t=polyval(coeff,s);

function [fx]=runge(x)
fx=1./(x.^2+1);

Eseguire poi da Matlab/Octave il programma principale esperimento.m

n=11;               % GRADO.
x=-5:10/n:5;        % NODI INTERP. (EQUISPAZ.).
y=runge(x);         % FUNZIONE NEI NODI EQUISP.
s=-5:10/(10*n):5;    % NODI TEST.
t=interpol(x,y,s);  % INTERPOLANTE NEI NODI TEST.
xcheb=cheb(-5,5,n); % NODI CHEB.
ycheb=runge(xcheb); % FUNZIONE NEI NODI CHEB.
tcheb=interpol(xcheb,ycheb,s);  % INTP. CHEB.
plot(s,runge(s),s,t,s,tcheb);   % PLOT INTP. VS RUNGE.
err_eqs=norm(runge(s)-t,inf);     % ABS. ERR. EQUISPAZ.
err_cheb=norm(runge(s)-tcheb,inf); % ABS. ERR. CHEB.
fprintf('\n \t [ABS.ERR.][EQS]: 2.2e [CHEB]: %2.2e',err_eqs,err_cheb);

Provando per diversi valori di , e tralasciando i warnings di Matlab relativi a polyfit, otteniamo:

[N]:  2 [ABS.ERR.][EQS]: 6.46e-001 [CHEB]: 9.62e-001
[N]:  3 [ABS.ERR.][EQS]: 7.07e-001 [CHEB]: 6.46e-001
[N]:  4 [ABS.ERR.][EQS]: 4.38e-001 [CHEB]: 8.29e-001
[N]:  5 [ABS.ERR.][EQS]: 4.33e-001 [CHEB]: 4.58e-001
[N]:  6 [ABS.ERR.][EQS]: 6.09e-001 [CHEB]: 6.39e-001
[N]:  7 [ABS.ERR.][EQS]: 2.47e-001 [CHEB]: 3.11e-001
[N]:  8 [ABS.ERR.][EQS]: 1.04e+000 [CHEB]: 4.60e-001
[N]:  9 [ABS.ERR.][EQS]: 2.99e-001 [CHEB]: 2.04e-001
[N]: 10 [ABS.ERR.][EQS]: 1.92e+000 [CHEB]: 3.19e-001
[N]: 11 [ABS.ERR.][EQS]: 5.57e-001 [CHEB]: 1.32e-001
[N]: 12 [ABS.ERR.][EQS]: 3.66e+000 [CHEB]: 2.18e-001
[N]: 13 [ABS.ERR.][EQS]: 1.07e+000 [CHEB]: 8.41e-002
[N]: 14 [ABS.ERR.][EQS]: 7.15e+000 [CHEB]: 1.47e-001
[N]: 15 [ABS.ERR.][EQS]: 2.10e+000 [CHEB]: 5.33e-002

Si nota subito che

• l’istruzione m=length(x)-1; calcola il grado dell’interpolante che è uguale al numero di punti meno uno;

• i nodi di cheb, sono quelli di Chebyshev-Gauss-Lobatto. Infatti per 1

(1)  \frac{( )}{2} \frac{( )}{2} cos( \frac{1 1}{1} )
\frac{( )}{2} \frac{( )}{2} cos(0)
\frac{( )}{2} \frac{( )}{2}
(16)

mentre porge

( )  \frac{( )}{2} \frac{( )}{2} cos( \frac{1}{1} )
\frac{( )}{2} \frac{( )}{2} cos( )
\frac{( )}{2} \frac{( )}{2} ( 1)
(17)

Osserviamo che a prima vista la routine cheb sembra produrre un set di nodi diversi da quelli di Chebyshev-Gauss-Lobatto (scalati) introdotti in (14). In realtà sono gli stessi come si può facilmente verificare matematicamente o in Matlab

>> N=10;
>> k=0:N;
>> t=-cos(k*pi/N);
>> t=sort(t);
>> t
t =
  Columns 1 through 7
   -1.0000   -0.9511   -0.8090   -0.5878   -0.3090   -0.0000    0.3090
  Columns 8 through 11
    0.5878    0.8090    0.9511    1.0000
>> xc=cheb(-1,1,N+1)
xc =
  Columns 1 through 7
   -1.0000   -0.9511   -0.8090   -0.5878   -0.3090   -0.0000    0.3090
  Columns 8 through 11
    0.5878    0.8090    0.9511    1.0000
>>

• il caso dei nodi equispaziati è semplice da programmare. Infatti la chiamata

>> a=-5;
>> b=5;
>> N=10;
>> h=(b-a)/N;
>> x=a:h:b
x =
  -5  -4  -3  -2  -1   0   1   2   3   4   5
>> length(x)
ans = 11
>>

suggerisce come fissato si possa generare un set di 1 punti equispaziati in [ ] (comprendente gli estremi).

• aiutandosi con l’help è possibile avere un’idea di cosa faccia il comando polyfit

>> help polyfit
 POLYFIT Fit polynomial to data.
    POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
    degree N that fits the data, P(X(I))~=Y(I), in a least-squares sense.
    [P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and a
    structure S for use with POLYVAL to obtain error estimates on
    predictions.  If the errors in the data, Y, are independent normal
    with constant variance, POLYVAL will produce error bounds which
    contain at least 50% of the predictions.

    The structure S contains the Cholesky factor of the Vandermonde
    matrix (R), the degrees of freedom (df), and the norm of the
    residuals (normr) as fields.
    [P,S,MU] = POLYFIT(X,Y,N) finds the coefficients of a polynomial
    in XHAT = (X-MU(1))/MU(2) where MU(1) = mean(X) and MU(2) = std(X).
    This centering and scaling transformation improves the numerical
    properties of both the polynomial and the fitting algorithm.
    Warning messages result if N is >= length(X), if X has repeated, or
    nearly repeated, points, or if X might need centering and scaling.
    See also POLY, POLYVAL, ROOTS.
>>

Quello che non dice la routine è come venga descritto il polinomio interpolatore. Aiutiamoci con il calcolo simbolico eseguito da Matlab, per calcolare il polinomio interpolatore in forma usuale.

>> syms x
>> simplify((-2*(x-1)*(x-3)/15)+(11*(x+2)*(x-3)/(-6))+(17*(x+2)*(x-1)/10))
ans =
-4/15*x^2+61/15*x+36/5
>> format long
>> -4/15
ans =
  -0.26666666666667
>> 61/15
ans =
   4.06666666666667
>> 36/5
ans =
   7.20000000000000

Testiamo l’uso di polyfit per quest’esempio

>> x=[-2 1 3];
>> y=[-2 11 17];
>> a=polyfit(x,y,2);
>> a
a =
  -0.26666666666667   4.06666666666667   7.19999999999999
>>

Quindi, se ( _{})_{1 3}, abbiamo

2( ) _{1} ^{2} _{2} _{3}

Più in generale, se _{}è il polinomio interpolatore di grado , e ( _{}) è il vettore ottenuto utilizzando polyfit, allora

( ) _{1} ^{}_{2} ^{1} _{1}

• aiutandosi con l’help è possibile avere un’idea di cosa faccia il comando polyval

>> help polyval
To get started, select "MATLAB Help" from the Help menu.

POLYVAL Evaluate polynomial.
Y = POLYVAL(P,X), when P is a vector of length N+1 whose elements
are the coefficients of a polynomial, is the value of the
polynomial evaluated at X.
    Y = P(1)*X^N + P(2)*X\^(N-1) + ... + P(N)*X + P(N+1)
If X is a matrix or vector, the polynomial is evaluated at all
points in X.  See also POLYVALM for evaluation in a matrix sense.
Y = POLYVAL(P,X,[],MU) uses XHAT = (X-MU(1))/MU(2) in place of X.
The centering and scaling parameters MU are optional output
computed by POLYFIT.
[Y,DELTA] = POLYVAL(P,X,S) or [Y,DELTA] = POLYVAL(P,X,S,MU) uses
the optional output structure S provided by POLYFIT to generate
error estimates, Y +/- delta.  If the errors in the data input to
POLYFIT are independent normal with constant variance, Y +/- DELTA
contains at least 50% of the predictions.
See also POLYFIT, POLYVALM.
>>

L’help ci suggerisce che quanto fornito da polyfit è effettivamente un polinomio scritto nella forma

(1) ^{}(2) ^{(} 1)  ( )  ( 1)

come notato precedentemente.

Se ( _{})_{1 1} e ( _{})_{1} , la chiamata

t=polyval(a,s);

valuta il polinomio

( ) _{1} ^{}_{1}

nei nodi ( _{}) e pone ( _{}) con _{}_{}( _{}).

• l’errore dell’interpolante della funzione di Runge nei nodi equispaziati non converge a 0. Osserviamo la parziale descrescita degli errori dall’iterata 2 alla 4 per poi aumentare peggiorando da iterazione ad iterazione. Similmente osserviamo la parziale descrescita degli errori dall’iterata 3 alla 5, dall’iterata 5 alla 7, per poi peggiorare sempre più . Il grafico in figura fa capire come l’errore sia particolarmente significativo in prossimità degli estremi 5 e 5.

• l’errore dell’interpolante della funzione di Runge nei nodi di Chebyshev-Gauss-Lobatto (scalati) converge a 0 (si osservi la particolare decrescita dei gradi pari e dispari).

• Ci si può chiedere perchè i nodi di Chebyshev-Gauss-Lobatto (scalati) vadano meglio. Ricordiamo che la funzione di Runge è ^{}([ 5 5]) e sussiste il seguente teorema dovuto a Bernstein (cf. ??p.136]COM)

Teorema 2.1 Se ^{1}([ ]) con [ ] intervallo limitato e chiuso della retta reale, il polinomio _{}di grado di interpolazione della funzione nei nodi di Chebyshev-Gauss di grado 1 converge uniformemente a su [ ], per . Se inoltre ^{2}([ ]) si ha la seguente stima dell’errore

( ^{1/2})

Il comportamento dell’interpolante nei nodi di Chebyshev-Gauss-Lobatto è analogo.

3 Facoltativo e tecnico: sul calcolo dell’interpolante polinomiale

Lanciamo in Matlab 6.1 il codice esperimento, che fa un confronto tra nodi equispaziati e di Chebyshev per 12 (quindi ci sono 13 punti di sampling). Il risultato è il seguente

>> esperimento
     [ABS.ERR.][EQS]: 3.66e+000 [CHEB]: 2.18e-001

Variamo in esperimento la variabile , ponendola uguale a 13. Matlab da’ il seguente messaggio di warning:

>> esperimento
Warning: Polynomial is badly conditioned. Remove repeated data points
    or try centering and scaling as described in HELP POLYFIT.
> In C:\MATLAB6p1\toolbox\matlab\polyfun\polyfit.m at line 74
In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\interpol.m
at line 17
In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\esperimento.m

  at line 5
          [ABS.ERR.][EQS]: 1.07e+000 [CHEB]: 8.41e-002
>>

Poniamo ora 15 è il risultato è anche peggio

>> esperimento
Warning: Polynomial is badly conditioned. Remove repeated data points
         or try centering and scaling as described in HELP POLYFIT.
> In C:\MATLAB6p1\toolbox\matlab\polyfun\polyfit.m at line 74
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\interpol.m
  at line 17
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\esperimento.m
  at line 5
Warning: Polynomial is badly conditioned. Remove repeated data points
         or try centering and scaling as described in HELP POLYFIT.
> In C:\MATLAB6p1\toolbox\matlab\polyfun\polyfit.m at line 74
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\interpol.m
   at line 17
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\esperimento.m
   at line 9
          [ABS.ERR.][EQS]: 2.10e+000 [CHEB]: 5.33e-002
>>

Diamo una spiegazione. Il codice di polyfit risolve il sistema lineare prodotto dalla matrice di Vandermonde (cf. [18]), con una fattorizzazione (cf. [17]). In altre parole se ( _{}) con _{}_{}^{1} ( _{}sono le ascisse dei nodi di sampling prescelti), si scrive con ortonormale e triangolare superiore e quindi si risolve ( sono i coefficenti del polinomio interpolante ( ) _{}_{}^{}) come soluzione del sistema ^{}. Se è malcondizionata in norma 1, più precisamente _{1} ^{1} _{1} 10^{10}, Matlab/Octave segnala un warning. La cosa è più elaborata di così in quanto polyfit usa la cosidetta economy size version della fattorizzazione QR, ma tralasciamo i dettagli per chiarezza. Per sopperire a questi problemi una delle possibili soluzioni consiste nell’usare la cosidetta formulazione baricentrica dei polinomi di Lagrange (cf. [2], [19]). Questa è implementata in Matlab nel codice barylag [19]. Una volta scaricato e scompattato il file zip, e digitato il codice esperimento_bar

n=15;                    % GRADO.
x=(-5:10/n:5)';       % NODI INTERP. (EQUISPAZ.).
y=runge(x);        % FUNZIONE NEI NODI EQUISP.
s=(-5:10/(10*n):5)';  % NODI TEST.
data_eqs=[x y];
t=barylag(data_eqs,s);

xcheb=(cheb(-5,5,n))';   % NODI CHEB.
ycheb=runge(xcheb);   % FUNZIONE NEI NODI CHEB.
data_cheb=[xcheb ycheb];
tcheb=barylag(data_cheb,s);
                      % INTP. CHEB.
plot(s,runge(s),s,t,s,tcheb);   % PLOT INTP. VS RUNGE.
err_eqs=norm(runge(s)-t,inf);     % ABS. ERR. EQUISPAZ.
err_cheb=norm(runge(s)-tcheb,inf); % ABS. ERR. CHEB.
fprintf('\n \t [ABS.ERR.][EQS]: %2.2e [CHEB]: %2.2e',err_eqs,err_cheb);

abbiamo

>> esperimento_bar
          [ABS.ERR.][EQS]: 2.10e+000 [CHEB]: 5.33e-002

Per curiosità paragoniamo i due codici esperimento ed esperimento_bar per 100, ottenendo rispettivamente

>> esperimento
Warning: Polynomial is badly conditioned. Remove repeated data points
         or try centering and scaling as described in HELP POLYFIT.
> In C:\MATLAB6p1\toolbox\matlab\polyfun\polyfit.m at line 74
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\interpol.m
   at line 17
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\esperimento.m
   at line 5
Warning: Polynomial is badly conditioned. Remove repeated data points
         or try centering and scaling as described in HELP POLYFIT.
> In C:\MATLAB6p1\toolbox\matlab\polyfun\polyfit.m at line 74
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\interpol.m
   at line 17
  In D:\DVD_27_OTTOBRE_2007\CS_2008\INTERPOLAZIONE_2008\MFILES\esperimento.m
   at line 9
          [ABS.ERR.][EQS]: 5.02e+019 [CHEB]: 2.14e+007
>>

e

>> esperimento_bar
          [ABS.ERR.][EQS]: 4.57e+005 [CHEB]: 5.62e-009
>>

segno che i risultati ottenuti con polyfit sono inaccettabili, mentre quelli di barylag sono perfettamente in parallelo con quanto detto nelle sezioni precedenti relativamente al fenomeno di Runge (si osservi che nel primo caso l’interpolazione nei nodi di Chebyshev sembra divergere, in disaccordo con la teoria!).

4 Sulla formula del resto

Si consideri la funzione ( ) : log( )

e si supponga di conoscere i valori per 1, 1 1 ed 1 2. Si calcoli il polinomio di secondo grado _{2} che interpola tali punti, e lo si valuti in 1 09. Quindi si utilizzi la formula dell’errore per funzioni di classe ^{1} (teorema del resto)

( ) _{}( ) ^{( 1)} ( )\frac{_{0}^{}( _{})}{( 1)!}  ( _{0}  _{})  (18)

per valutare

log( ) _{2}( )

Risoluzione. Dalla formula del resto (per ( ) log( ), 1 09, _{0} 1, _{1} 1 1, 2 1 2 ed 2) otteniamo:

log(1 09) _{2}(1 09) \frac{2}{^{3}} \frac{(1 09 1)(1 09 1 1)(1 09 1 2)}{6}

in quanto la derivata terza di log( )è \frac{2}{^{3}} . Siccome (1 1 2), dalla decrescenza di \frac{2}{^{3}} deduciamo che

1 1574 \frac{2}{1 2^{3}} \frac{2}{^{3}} \frac{2}{1^{3}} 2

Essendo

(1 09 1)(1 09 1 1)(1 09 1 2) 9 9 10^{5}

ricaviamo quindi che l’errore dell’interpolante è

1 9097 10^{5} 1 1574 9 9 10^{5} /6 log( ) _{2}( ) 2 9 9 10^{5} /6 3 3 10^{5}

Utilizzando la funzione interpol.m precedentemente introdotta, si verifica che

2(1 09) 0 08615260795055

mentre log(1 09) 0 08617769624105

e quindi

log(1 09) _{2}(1 09) 2 51 10^{5}

perfettamente in linea con le stime fornite.

Esercizio. Effettuare un programma Matlab che calcoli il valore assunto in 1 09 dal polinomio interpolatore della funzione log relativamente ai punti 1, 1 1, 1 2. E’ buona l’approssimazione fornita dal polinomio interpolatore? (Suggerimento: si noti che 1. il programma è una variante di quanto visto in esperimento.m, per una adeguata

scelta dei nodi x, dei valori y, del punto s;

2. eliminare la parte relativa ai nodi di Chebyshev-Gauss-Lobatto;

3. il plot non è molto indicativo;

4. siccome il comando norm(v,inf) calcola per _{}la quantità

max_{}_{}

visto che si deve valutare l’errore su un solo punto, il comando norm(...,inf) in esperimento.m può essere sostituito da abs(...) (rifletterci sopra);

5. la funzione log è predefinita in Matlab/Octave e quindi non serve ridefinirla come funzione.

Esercizio. Aiutandosi con Matlab/Octave eseguire un programma che calcoli l’interpolante della funzione

( ) : exp( )

relativamente alle ascisse 1, 1 1 ed 1 2. Si calcoli il polinomio di secondo grado _{2} che interpola tali punti, e lo si valuti in 1 09. Fornire una stima dell’errore e verificarne la bontà rispetto al risultato esatto.

Esercizio. Fissati 1, 1, 15, si plottino i nodi di Chebyshev-Gauss. Si suggerisce di utilizzare la function cheb.m e per il plottaggio un comando del tipo

plot(x,y,'r-o')

Ricordiamo che

• l’opzione r-o

1. disegna un cerchietto rosso per ognuno dei punti ( _{1} _{1}), , ( _{}_{}), , ( _{}_{});

2. per 1 1, unisce i punti ( _{}_{}), ( _{1} _{1}) con un segmento rosso.

• essendo importante disegnare solo le ascisse (e non le ordinate), quali ordinate si può porre

     y=zeros(size(x));

Una volta completato l’esercizio si osservi la particolare disposizione dei nodi di Chebyshev-Gauss. Si accumulano verso gli estremi dell’intervallo?

5 Online

Si possono reperire online varie bibliografie dei matematici sopra citati:

http://en.wikipedia.org/wiki/Chebyshev
http://www-history.mcs.st-andrews.ac.uk/Biographies/Chebyshev.html.
http://it.wikipedia.org/wiki/Lagrange.
http://www-history.mcs.st-andrews.ac.uk/Biographies/Lagrange.html.
http://en.wikipedia.org/wiki/Isaac_Newton

Per quanto riguarda l’interpolazione:

http://it.wikipedia.org/wiki/Interpolazione_\%28matematica%29
http://it.wikipedia.org/wiki/Interpolazione_polinomiale
http://en.wikipedia.org/wiki/Interpolation
http://en.wikipedia.org/wiki/Polynomial_interpolation
http://it.wikipedia.org/wiki/Fenomeno_di_Runge
http://en.wikipedia.org/wiki/Runge\%27s_phenomenon
http://en.wikipedia.org/wiki/Newton_polynomial
http://it.wikipedia.org/wiki/Nodi_di_Chebyshev
http://en.wikipedia.org/wiki/QR_decomposition
http://en.wikipedia.org/wiki/Vandermonde_matrix

6 Frasi celebri

1. To isolate mathematics from the practical demands of the sciences is to invite the sterility of a cow shut away from the bulls. (Chebyshev)

2. God made the integers; all else is the work of man. (Kronecker).

3. Before we take to sea we walk on land, Before we create we must understand. (Lagrange)

4. When we ask advice, we are usually looking for an accomplice. (Lagrange)

5. It took the mob only a moment to remove his head; a century will not suffice to reproduce it. [said about the chemist Lavoisier] (Lagrange)

6. I do not know. [summarising his life’s work] (Lagrange)

References

[1] K. Atkinson, An Introduction to Numerical Analysis, Wiley, (1989).

[2] J.P. Berrut and L.N.Trefethen, Barycentric Lagrange interpolation, (SIAM Review, 2004),

  http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/barycentric.pdf.

[3] V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

[4] Mac Tutor, Chebyshev biography,

http://www-history.mcs.st-andrews.ac.uk/Biographies/Chebyshev.html.

[5] Mac Tutor, Lagrange biography,

http://www-history.mcs.st-andrews.ac.uk/Biographies/Lagrange.html.

[6] The MathWorks Inc., Numerical Computing with Matlab,

http://www.mathworks.com/moler.

[7] A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006.

[8] A. Quarteroni, R. Sacco e F. Saleri, Matematica numerica, Springer Verlag, 1998.

[9] Wikipedia, Fenomeno di Runge,

http://it.wikipedia.org/wiki/Fenomeno_di_Runge.

[10] Wikipedia, Interpolazione polinomial,

http://it.wikipedia.org/wiki/Interpolazione_polinomiale.

[11] Wikipedia, Isaac Newton,

http://en.wikipedia.org/wiki/Isaac_Newton.

[12] Wikipedia, Lagrange,

http://it.wikipedia.org/wiki/Lagrange.

[13] Wikipedia, Newton polynomial,

http://en.wikipedia.org/wiki/Newton_polynomial.

[14] Wikipedia, Nodi di Chebyshev,

http://it.wikipedia.org/wiki/Nodi_di_Chebyshev.

[15] Wikipedia, Polynomial interpolation,

http://en.wikipedia.org/wiki/Polynomial_interpolation.

[16] Wikipedia, Runge’s phenomenon,

http://en.wikipedia.org/wiki/Runge%27s_phenomenon.

[17] Wikipedia, QR decomposition,

http://en.wikipedia.org/wiki/QR_decomposition.

[18] Wikipedia, Vandermonde matrix,

http://en.wikipedia.org/wiki/Vandermonde_matrix.

[19] Greg von Winckel, barylag.m,

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=4478&objectType=FILE.