Interpolazione polinomiale.

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

21 novembre 2011

Interpolazione polinomiale

Sia ℙ_{n} lo sp. vett. dei polinomi di grado n in ℝ. Date n + 1 coppie (x_{0},y_{0}), …, (x_{n},y_{n}) con x_{j} ≠ x_{k} se j ≠ k, si calcoli p_{n} ∈ ℙ_{n} t.c.

p_{n}(x_{k}) = y_{k}, k = 0, …, n

Tali polinomi p_{n} si dicono interpolare {(x_{k}, y_{k})}_{k=0,...,n} oppure interpolare i valori y_{k} nei nodi x_{k}.

Figura: Grafico che illustra l’interpolazione di cos (x) in [−5, 5] su nodi equispaziati x_{k} = −5 + 2 · k (k = 0, …, 5).

Alvise Sommariva  Interpolazione polinomiale.  2/ 27

Questioni legate all’interpolazione polinomiale

Alcune questioni risultano di importanza fondamentale: ▶ Esiste tale polinomio? ▶ Se esiste è unico? ▶ Se esiste ed è unico, è possibile calcolarlo? ▶ Se f ∈ C([a, b]), esistono {x_{k}}_{k=0,...,n} ∈ [a, b] che

rappresentano una scelta migliore di altre, nel senso che il polinomio interpolante (x_{k},f (x_{k})) (k = 0, …, n) approssima meglio la funzione?

Qualche esempio

▶ n=1: dati due punti (x_{0},y_{0}), (x_{1},y_{1}) trovare un polinomio di grado 1 che passi per i due punti. Equivale a dire: calcolare retta che passa per due punti assegnati.

▶ n=2: dati tre punti (x_{0},y_{0}), (x_{1},y_{1}), (x_{2},y_{2}) trovare un polinomio di grado 2 che passi per i tre punti. Equivale a dire: calcolare parabola

y = a_{1}x^{2} + a_{2}x + a_{3}

che passa per tre punti assegnati.

Polinomi di Lagrange ed esistenza polinomio interpolatore

Consideriamo i polinomi di Lagrange L_{k} ∈ ℙ_{n} (rel. {x_{k}}_{k=0,...,n})

L_{k}(x) = \frac{(x − x_{0}) · . . . · (x − x_{k−1}) · (x − x_{k+1}) · . . . · (x − x_{n})}{(x_{k} − x_{0}) · . . . · (x_{k} − x_{k−1}) · (x_{k} − x_{k+1}) · . . . · (x_{k} − x_{n})}

Detto δ_{k,j} l’operatore di Kronecker, L_{k}(x_{s}) = δ_{k,s} cioè

L_{k}(x_{k}) = 1
L_{k}(x_{j}) = 0, se k ≠ j.

Allora ha le proprietà desiderate il polinomio p_{n} di grado n p_{n}(x) = ∑_{k=0}^{n}y_{k}L_{k}(x)

poichè p_{n}(x_{s}) = ∑_{k=0}^{n} y_{k}L_{k}(x_{s}) = ∑_{k=0}^{n} y_{k}δ_{k,s} = y_{s} per s = 0, …, n. Si noti che {L_{k}}_{k=0,...,n} è base di ℙ_{n}.

Alvise Sommariva  Interpolazione polinomiale.  5/ 27

Unicità polinomio interpolatore

Supponiamo p_{n}, q_{n} ∈ ℙ_{n} siano tali che

p_{n}(x_{k}) = q_{n}(x_{k}) = y_{k}, k = 0, …, n.

Allora s_{n} = p_{n} − q_{n} ∈ ℙ_{n} ed è

s_{n}(x_{k}) = p_{n}(x_{k}) − q_{n}(x_{k}) = 0, k = 0, …, n

per cui

s_{n}(x) = α(x − x_{0}) · . . . · (x − x_{n}).

Se α ≠ 0 allora s_{n} ∈ ℙ_{n+1}, cosa assurda. Quindi α = 0, il che implica p_{n} − q_{n} = s_{n} = 0 cioè p_{n} = q_{n}. Deduciamo così che il polinomio p_{n} ∈ ℙ_{n} interpolante {(x_{k}, y_{k})}_{k=0,...,n} (con x_{k} ≠ x_{s} se k ≠ s) non solo esiste ma è pure unico ed è

p_{n}(x) = ∑_{k=0}^{n}y_{k}L_{k}(x).

Esempio

Calcoliamo il polinomio di grado 2 che assume nei nodi x_{0} = −2, x_{1} = 1, x_{2} = 3 rispettivamente i valori f_{0} = −2, f_{1} = 11, f_{2} = 17.

L_{0}(x) = \frac{(x − x_{1})(x − x_{2})}{(x_{0} − x_{1})(x_{0} − x_{2})} = \frac{(x − 1)(x − 3)}{(−2 − 1)(−2 − 3)} = \frac{(x − 1)(x − 3)}{15}
L_{1}(x) = \frac{(x − x_{0})(x − x_{2})}{(x_{1} − x_{0})(x_{1} − x_{2})} = \frac{(x − (−2))(x − 3)}{(1 − (−2))(1 − 3)} = \frac{(x + 2)(x − 3)}{−6}
L_{2}(x) = \frac{(x − x_{0})(x − x_{1})}{(x_{2} − x_{0})(x_{2} − x_{1})} = \frac{(x − (−2))(x − 1)}{(3 − (−2))(3 − 1)} = \frac{(x + 2)(x − 1)}{10}

Si vede subito che L_{0}, L_{1}, L_{2} ∈ ℙ_{2} che L_{0}(x_{0}) = L_{0}(−2) = 1, L_{0}(x_{1}) = L_{0}(1) = 0, L_{0}(x_{2}) = L_{0}(3) = 0 L_{1}(x_{1}) = L_{1}(1) = 1, L_{1}(x_{0}) = L_{1}(−2) = 0, L_{1}(x_{2}) = L_{1}(3) = 0 L_{2}(x_{2}) = L_{2}(3) = 1, L_{2}(x_{0}) = L_{2}(−2) = 0, L_{2}(x_{1}) = L_{2}(1) = 0

Alvise Sommariva  Interpolazione polinomiale.  7/ 27

Esempio

Posto p_{2}(x) = −2 L_{0}(x) + 11 L_{1}(x) + 17 L_{2}(x) abbiamo

p(x_{0}) = −2 L_{0}(x_{0})+11 L_{1}(x_{0})+17 L_{2}(x_{0}) = −2·1+11·0+17·0 = −2 = y_{0},

p(x_{1}) = −2 L_{0}(x_{1})+11 L_{1}(x_{1})+17 L_{2}(x_{1}) = −2·0+11·1+17·0 = 11 = y_{1},

p(x_{2}) = −2 L_{0}(x_{2})+11 L_{1}(x_{2})+17 L_{2}(x_{2}) = −2·0+11·0+17·1 = 17 = y_{2}.

cioè p_{2} interpola le coppie {(x_{k}, y_{k})}_{k=0,1,2}.

Teorema del resto

Usualmente i dati y_{k} corrispondono alla valutazione di una funzione f ∈ C([a, b]) in nodi x_{k} ∈ [a, b]. Si desidera che il polinomio interpolatore p_{n} approssimi la funzione f . Risulta quindi importante stimare l’errore compiuto. A tal proposito vale il seguente teorema detto del resto.

Teorema. Sia f ∈ C^{(n+1)}(a, b) e sia p_{n} il polinomio che interpola le coppie {(x_{k}, y_{k})}_{k=0,...,n} con x_{k} ≠ x_{s} se k ≠ s. Allora

f (x) − p_{n}(x) = f^{(n+1)}(ξ)\frac{∏_{i=0}^{n}(x − x_{i})}{(n + 1)!}  (1)

dove ξ ∈ I con I il più piccolo intervallo aperto contenente x_{0}, …, x_{n}.

Teorema del resto: esempio

Supponiamo f (x) = log(x) in [1, 2]. Si vede subito che |f^{(n+1)}(x)| = 1/x^{n+1} e quindi ∥f^{(n+1)}∥_{∞} = 1.

Alvise Sommariva  Interpolazione polinomiale.  10/ 27

Alcune scelte dei nodi

Consideriamo l’intervallo chiuso e limitato [a, b]. Vediamo alcuni sets di nodi.

▶ Equispaziati: x_{k} = a + k \frac{(b−a)}{n}, k = 0, …, n. ▶ Gauss-Chebyshev (scalati): fissato n, i punti sono

x_{k} = \frac{(a + b)}{2} + \frac{(b − a)}{2} t_{k}, k = 0, …, n  (2)

con

t_{k} = cos (\frac{2k + 1}{2n + 2} π) , k = 0, …, n;  (3)

▶ Gauss-Chebyshev-Lobatto (scalati): fissato n, i punti sono

x_{k} = \frac{(a + b)}{2} + \frac{(b − a)}{2} t_{k}, k = 0, …, n  (4)

con

t_{k} = − cos (\frac{kπ}{n}) , k = 0, …, n.  (5)
Alvise Sommariva  Interpolazione polinomiale.  11/ 27

Interpretazione geometrica Gauss-Legendre-Lobatto

I punti di Gauss-Chebyshev-Lobatto si possono interpretare geometricamente.

▶ Si considerano gli n + 1 punti della circonferenza

{(cos(θ_{k}), sin(θ_{k}))}_{k=0,...,n} con θ_{k} = (k · π)/n ▶ Si proiettano {(cos(θ_{k}), sin(θ_{k}))}_{k=0,...,n} sull’asse delle ascisse.

Le coordinate x di tali punti corrispondono ai punti di Gauss-Chebyshev-Lobatto.

Interpretazione geometrica Gauss-Chebyshev-Lobatto

Figura: Interpretazione geometrica dei punti di Gauss-Chebyshev-Lobatto.

Alvise Sommariva  Interpolazione polinomiale.  13/ 27

Convergenza ed esempio di Runge

L’interpolante polinomiale in un set di nodi prefissati non converge sempre puntualmente alla funzione da approssimare. Infatti, per la funzione di Runge

f (x) = \frac{1}{1 + x^{2}}, x ∈ [−5, 5]  (6)

si ha che il polinomio interpolatore p_{n} in nodi equispaziati non converge (puntualmente) a f . Fortunatamente ciò non succede per i nodi di Gauss-Chebyshev(-Lobatto).

Purtroppo, per un teorema dovuto a Faber, esistono comunque funzioni continue f (ma non C^{1}!!) tali che l’interpolante p_{n} nei nodi di Gauss-Chebyshev(-Lobatto) non converge puntualmente a f . Questo teorema è generalizzabile ad una famiglia arbitraria di nodi.

Esempio di Runge

Figura: Grafico della funzione di Runge 1/(1 + x^{2}) nell’intervallo [−5, 5] e delle sue interpolanti di grado 12 nei nodi equispaziati e di Gauss-Chebyshev-Lobatto.

Figura: Grafico errore della funzione di Runge 1/(1 + x^{2}) nell’intervallo [−5, 5] con le sue interpolanti di grado 12 nei nodi equispaziati e di Gauss-Chebyshev-Lobatto.

Esempio di Runge

Se p_{n} ∈ ℙ_{n} è polinomio intp. nei nodi eqsp. o di G.-C.-L., vediamo al variare del grado quali sono gli errori ∥1/(1 + x^{2}) − p_{n}(x)∥_{∞}:

Deg Err. Eqs. Err. GCL
2 6.46e − 001 9.62e − 001
3 7.07e − 001 6.46e − 001
4 4.38e − 001 8.29e − 001
5 4.33e − 001 4.58e − 001
6 6.09e − 001 6.39e − 001
7 2.47e − 001 3.11e − 001
8 1.04e + 000 4.60e − 001
9 2.99e − 001 2.04e − 001
10 1.92e + 000 3.19e − 001
11 5.57e − 001 1.32e − 001
12 3.66e + 000 2.18e − 001
13 1.07e + 000 8.41e − 002
14 7.15e + 000 1.47e − 001

Alvise Sommariva  Interpolazione polinomiale.  17/ 27

Importante: l’esempio di Runge mostra che esiste una funzione C^{∞}([−5, 5]) tale che al crescere del numero di nodi equispaziati n non sia garantita nemmeno la convergenza puntuale!

D’altra parte sussiste il seguente teorema dovuto a Bernstein:

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

∥f − p_{n}∥_{∞} = O(n^{−1/2}).

Interpolazione di Newton

Supponiamo note le coppie (x_{0},f (x_{0})), …, (x_{n},f (x_{n})) con x_{k} ≠ x_{j} se k ≠ j. Allora posto w_{k}(x) = (x − x_{0}) · . . . · (x − x_{k−1}),

f [x_{j}, x_{j+1}, …, x_{s}] = \frac{f [x_{j+1}, …, x_{s}] − f [x_{j}, …, x_{s−1}]}{x_{s} − x_{j}}, f [x_{j}] = f (x_{j})

si dimostra che il polinomio p_{n} ∈ ℙ_{n} interpolante le coppie citate è

p_{n}(x) = ∑_{k=0}^{n}w_{k}(x)f [x_{0}, …, x_{k}].

Il vantaggio di questa scrittura è che se aggiungiamo un nuovo punto da interpolare, diciamo (x_{n+1},f (x_{n+1})):

p_{n+1}(x) = ∑_{k=0}^{n+1}w_{k}(x)f [x_{0}, …, x_{k}] = ∑_{k=0}^{n}w_{k}(x)f [x_{0}, …, x_{k}] +
+ w_{n+1}f [x_{0}, …, x_{n+1}] = p_{n}(x) + w_{n+1}f [x_{0}, …, x_{n+1}]

Differenze divise

Le quantità f [x_{j}, …, x_{s}] si chiamano differenze divise. Di solito le quantità utili f [x_{0}, …, x_{k}] per il calcolo del polinomio interpolatore, si scrivono per via tabulare. Esempio:

x_{0} f [x_{0}]
x_{1} f [x_{1}] f [x_{0}, x_{1}]
x_{2} f [x_{2}] f [x_{1}, x_{2}]  f [x_{0}, x_{1}, x_{2}]
⋮ ⋮  ⋮  ...
x_{n} f [x_{n}] f [x_{n−1}, x_{n}] f [x_{n−2}, x_{n−1}, x_{n}] … f [x_{0}, …, x_{n}]

Differenze divise

Se vogliamo intp. le coppie dell’esempio precedente, cioè (x,y) ove

x = [−2; 1; 3], y = [−2; 11; 17]

otteniamo la tabella

−2 −2
1 11 4.\overline{3}
3 17 3.0 −0.2\overline{6}

e il polinomio intp. p_{2} risulta

p(x) = −2 + 4.\overline{3} · (x + 2) − 0.2\overline{6} · (x + 2)(x − 1) In effetti p_{2}(−2) = −2, p_{2}(1) = −2 + 4.\overline{3} · 3 − 0.2\overline{6} · 0 = 11 e p_{2}(3) = −2 + 4.\overline{3} · 5 − 0.2\overline{6} · 10 = −2 + 21.\overline{6} − 2.\overline{6} = 17.

Alvise Sommariva  Interpolazione polinomiale.  21/ 27

Interpolazione in Matlab/Octave

Supponiamo di dover interpolare le coppie {(x_{k}, y_{k})}_{k=0,...,n} e supponiamo sia x = [x_{0}, …, x_{n}], y = [y_{0}, …, y_{n}]. I coefficienti del polinomio interpolatore sono ottenibili dal comando polyfit. A tal proposito l’help di Matlab suggerisce:

>> h e l p p o l y f i t

 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 .

…

 See also POLY , POLYVAL , ROOTS .

>>

Per capire qualcosa in più eseguiamo il seguente codice

>> x=[−2 1 3 ] ;
>> y=[−2 11 1 7 ] ;
>> a=p o l y f i t ( x , y , 2 )

a =

−0.2667  4.0667  7.2000

>>

In effetti, calcolando manualmente il polinomio interpolatore si ha, semplificando quando ottenuto coi polinomi di Lagrange che è

p_{2}(x) = (−4/15) · x^{2} + (61/15) · x + (36/5) ≈ 0.2\overline{6}x^{2} + 4.0\overline{6}x + 7.2.

Quindi, se a = (a_{k})_{k=1,...,3}, abbiamo p_{2}(x) = a_{1} x^{2} + a_{2} x + a_{3} e più in generale, se p_{n} è il polinomio interpolatore di grado n, e a = (a_{k}) è il vettore ottenuto utilizzando polyfit, allora

p_{n}(x) = a_{1} x^{n} + a_{2} x^{n−1} + … + a_{n+1}.
Alvise Sommariva  Interpolazione polinomiale.  23/ 27

Interpolazione in Matlab/Octave

Per valutare in un vettore di ascisse X = [X_{k}]_{k=1,...,m} un polinomio

p_{n}(x) = a_{1} x^{n} + a_{2} x^{n−1} + … + a_{n+1}

i cui coefficienti sono i coeff. di in un vettore P = [a_{k}]_{k=1,...,n+1} usiamo il comando polyval. Dall’help:

>> h e l p p o l y v a l

To get started , select ”MATLAB Help” from the

Help menu .

POLYVAL Evaluate polynomial .
Y = POLYVAL ( P , X ) , when P is a vector of l e n g t h
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)

>>

Dati i vettori x = [x_{k}]_{k=1,...,n}, y = [y_{k}]_{k=1,...,n} sia p_{n−1}(x_{k}) = y_{k}, per k = 1, …, n. Sia s = [s_{k}]_{k=1,...,m} e desideriamo calcolare t = [t_{k}]_{k=1,...,m} per cui t_{k} = p_{n−1}(s_{k}), per ogni k. A tal proposito introduciamo la funzione:

f u n c t i o n t=interpol ( x , y , s )
m=l e n g t h ( x ) −1;
coeff=p o l y f i t ( x , y , m ) ;
t=p o l y v a l ( coeff , s ) ;

Interpoliamo la funzione di Runge 1/(1 + x^{2}) in [−5, 5] sia su nodi equispaziati che di tipo Gauss-Chebyshev-Lobatto. A tal proposito definiamo la funzione di Runge

f u n c t i o n [ fx]=runge ( x )
fx =1./( x .ˆ 2+1) ;

e una funzione che genera n nodi di Gauss-Chebyshev-Lobatto nell’intervallo [a, b]:

f u n c t i o n xc=cheb ( a , b , n )
f o r m =1:1: n
xc ( m )=(a+b ) /2 −((b−a ) /2) ∗cos ( p i ∗(m−1) /( n−1) ) ;

end

Esempio di Runge

Quindi scriviamo il file esperimento.m

n=11; % GRADO.
% NODI TEST .
s= −5:10/(10∗n ) : 5 ;
% NODI EQSP . : ASCISSE/ORDINATE + INTP . TEST .
x=−5:10/n : 5 ; y=runge ( x ) ;

t=interpol ( x , y , s ) ;

% NODI GCL . : ASCISSE/ORDINATE+INTP . TEST .

xgcl=cheb ( −5 ,5 , n+1) ; ygcl=runge ( xgcl ) ;
tt =interpol ( xgcl , ygcl , s ) ;

% PLOT INTP . VS RUNGE.
p l o t ( s , runge ( s ) , s , t , s , tt ) ;
% ERRORI ASSOLUTI .

ee=norm ( runge ( s )−t , inf ) ; ec=norm ( runge ( s )−tt , inf ) ;

f p r i n t f ( ’ \n\ t [ ERR . ] [ EQS] : % 2 . 2 e [ GCL] : % 2 . 2 e ’ , ee , ec ) ;

Esempio di Runge: risultati

Al variare di n: ▶ Otteniamo la tabella degli errori vista in precedenza. ▶ Notiamo che la scelta di n non può essere eccessiva. Provare

n = 30!! ▶ Risulta evidente che non sussiste la convergenza puntuale al

crescere di n, di p_{n} a 1/(1 + x^{2}), qualora si utilizzino nodi equispaziati.

▶ Risulta evidente che sussiste la convergenza puntuale al crescere di n, di p_{n} a 1/(1 + x^{2}), qualora si utilizzino nodi di Gauss-Chebyshev-Lobatto.