8 Interpolazione di Newton
Un limite dell’approccio di Lagrange1 è che, una volta calcolato il polinomio p_n(x) che interpola i dati,
\{(x_i,y_i)\}_{i=0}^n
per aggiungere un ulteriore punto e calcolare p_{n+1}(x) occorre ripetere l’intero procedimento: in generale, non è possibile ottenere direttamente p_{n+1} a partire da p_n. L’idea di Newton2 consente invece di costruire il polinomio in modo incrementale, scrivendo
p_{n+1}(x)=p_n(x)+g(x)
per un polinomio g(x) da determinare. Se valutiamo p_{n+1} sui nodi x_0,x_1,\ldots,x_{n+1} avremo che:
\begin{aligned} y_0 & = p_{n+1}(x_0)= p_n(x_0) + g(x_0) & \implies & g(x_0)=0\\ y_1 & = p_{n+1}(x_1)= p_n(x_1) + g(x_1) & \implies & g(x_1)=0\\ \vdots & & \vdots\\ y_n & = p_{n+1}(x_n)= p_n(x_n) + g(x_n) & \implies & g(x_n)=0\\ y_{n+1} & = p_{n+1}(x_{n+1})= p_n(x_{n+1}) + g(x_{n+1}) & \implies & g(x_{n+1})=y_{n+1}-p_n(x_{n+1}) \end{aligned}
Dunque g(x) è un polinomio che deve valere zero sui nodi x_0,\ldots,x_{n} e quindi si scriverà come
g(x)=c\omega_{n+1}(x)= c\cdot (x-x_0)(x-x_1)\cdot\ldots\cdot(x-x_n)
Allora
p_{n+1}(x)=p_n(x)+c\omega_{n+1}(x)
con \omega_{n+1}(x_i)=0 per i=0,\ldots,n. Il grado complessivo di p_{n+1}(x) è proprio n+1 per il contributo di \omega_{n+1}(x). Valutando p_{n+1}(x) sul nodo x_{n+1} si determina facilmente la costante c:
\begin{aligned} p_{n+1}(x_{n+1}) &=p_n(x_{n+1})+c\omega_{n+1}(x_{n+1})=y_{n+1} \\ a_{n+1}= c &= \dfrac{y_{n+1}-p_n(x_{n+1})}{\omega_{n+1}(x_{n+1})}. \end{aligned} \tag{8.1}
Esempio 8.1 Vediamo un esempio di costruzione incrementale col metodo di Newton. Consideriamo i dati \{(0,1),(2,3),(1,0)\}. I passi da fare sono:
- n=0. Il primo polinomio interpolante è p_0(x) = y_0 = 1.
- n=1. Il polinomio p_1(x) = p_0(x)+ c\omega_1(x). Si ha immediatamente che \omega_1(x)=x-x_0=x e quindi
\begin{aligned} p_1(x) &= p_0(x)+ c\omega_1(x) = 1 + cx = 3 \quad\implies \quad c=1\\ &= 1+x \end{aligned}
n=2. Il polinomio p_2(x)=p_1(x)+c\omega_2(x) con
\omega_2(x)=(x-x_0)(x-x_1)=x(x-2)=x^2-2x.
Allora si calcola c a partire da p_2 mediante
\begin{aligned} p_2(x) &=p_1(x)+c\omega_2(x)=(1+x) + c(x^2-2x)\\ p_2(x_2) &=p_1(1)+c\omega_2(1)=(1+1) + c(1^2-2\cdot 1) \quad\implies c=2\\ p_2(x) &= 1+x+2(x^2-2x)=1-3x+2x^2 \end{aligned}
Il polinomio interpolante è quindi 1-3x+2x^2.
Riprendiamo in modo più formale e schematico quanto appena visto. Per prima cosa definiamo il polinomio \omega_{k}(x) di grado k come
\omega_k(x):=\prod_{j=0}^{k-1}(x-x_j).
Si nota immediatamente che \omega_{k+1}(x) = (x-x_k)\omega_k(x) e la costruzione incrementale si può scrivere come
\begin{aligned} p_0(x) &= a_0 =f(x_0)=y_0\\ p_1(x) &= p_0(x)+a_1\omega_1(x) = a_0+a_1\omega_1(x)\\ p_2(x) &= p_1(x)+a_2\omega_2(x) = a_0+a_1\omega_1(x)+a_2\omega_2(x)\\ & \vdots \\ p_n(x) &= p_{n-1}(x)+a_n\omega_n(x) = a_0+a_1\omega_1(x)+\ldots+a_n\omega_n(x) \end{aligned}
La costruzione di Newton può essere formulata nuovamente come un sistema lineare, questa volta triangolare inferiore anziché rettangolare, come si vede dalla seguente espressione:
\begin{aligned} f(x_0) &= a_0\\ f(x_1) &= a_0+a_1\omega_1(x_1)\\ f(x_2) &= a_0+a_1\omega_1(x_2)+a_2\omega_2(x_2)\\ f(x_3) &= a_0+a_1\omega_1(x_3)+a_2\omega_2(x_3)+a_3\omega_3(x_3)\\ & \vdots & \\ f(x_n) &= a_0+a_1\omega_1(x_n)+a_2\omega_2(x_n)+\ldots+a_n\omega_n(x_n). \end{aligned}
Equivalentemente, in forma matriciale si ha \bm{M}\bm{a}=\bm{b}:
\begin{bmatrix} 1 & 0 & 0 & \ldots & 0 \\ 1 & \omega_1(x_1) & 0 & \ldots & 0 \\ \vdots & & & \vdots & \\ 1 & \omega_1(x_n) & \omega_2(x_n) & \ldots & \omega_n(x_n) \end{bmatrix} \begin{bmatrix} a_0\\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} f(x_0)\\ f(x_1) \\ \vdots \\ f(x_n) \end{bmatrix}. \tag{8.2}
Il determinante di \bm{M} è dato dal prodotto dei polinomi \omega_i e vale
\det \bm{M} =\omega_0(x_0)\omega_1(x_1)\cdots\omega_n(x_n)= \prod_{k=0}^n\omega_k(x_k) =\prod_{k=0}^n \prod_{j=0}^{k-1}(x_k-x_j).
Di conseguenza, nell’ipotesi fatta di avere nodi distinti, cioè x_i\neq x_j se i\neq j, il determinante è non nullo ed il problema ammette un’unica soluzione.
Seguendo la tradizione, si usa indicare il termine a_k come la differenza divisa di ordine k e si scrive con un simbolo speciale
a_k = f[x_0,x_1,\ldots,x_k] \tag{8.4}
e come conseguenza, il polinomio interpolatore di Newton si scrive come
p_n(x)=\sum_{k=0}^n f[x_0,x_1,\ldots,x_k]\omega_k(x).
La differenza divisa f[x_0,x_1,\ldots,x_k] è il coefficiente del monomio di grado massimo (k) del polinomio interpolatore dei primi k+1 punti, e non dipende dall’ordine con cui sono considerati i nodi. Come conseguenza, il polinomio di interpolazione non cambia se si cambia la sequenza dei nodi. Il modo efficiente di usare le differenze divise è quello di costruirne la tavola in modo ricorsivo. Siano x_i i nodi, la prima differenza divisa è f[x_i]=f(x_i), e i passi successivi dipendono solo dai precedenti.
8.1 Proprietà delle differenze divise
Dalla relazione 8.3 e dalla definizione 8.4 si intuisce che c’è una ricorsione che permette di calcolare un coefficiente a partire dai precedenti:
\begin{aligned} a_0 &= f[x_0]:=f(x_0)\\ a_1 &= f[x_0,x_1]:=\dfrac{f[x_1]-f[x_0]}{x_1-x_0}\\ a_2 &= f[x_0,x_1,x_2]:=\dfrac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_1}\\ & \vdots \\ a_k &= f[x_0,x_1,\ldots,x_k]:=\dfrac{f[x_1,\ldots,x_k]-f[x_0,x_1,\ldots,x_{k-1}]}{x_k-x_0}. \end{aligned}
Il coefficiente di grado massimo del polinomio interpolatore p_n(x) è proprio l’ultima differenza divisa f[x_0,\ldots,x_n]. Infatti il polinomio \omega_n è monico in quanto prodotto di fattori della forma (x-x_i). In generale si possono definire le differenze divise per una sequenza di punti qualsiasi:
f[x_i,x_{i+1},\ldots,x_{i+k}]:=\dfrac{f[x_{i+1},\ldots,x_{i+k}]-f[x_i,\ldots,x_{i+k-1}]}{x_{i+k}-x_i} \tag{8.5}
e da questa definizione si vede che per calcolare una differenza di ordine k servono le differenze divise di ordine j=0,\ldots,k-1. Questa dipendenza si esprime bene con la seguente tavola delle differenze divise:
La tavola si calcola una colonna alla volta da sinistra verso destra, utilizzando i valori precedenti per calcolare quelli nuovi. Osservando con attenzione la tavola, si nota che i coefficienti del polinomio interpolatore si trovano sul lato superiore del triangolo delle differenze divise. Notiamo anche che è necessario calcolare tutta la tavola, perché, ad esempio, per calcolare f[x_0,x_1,x_2] servono prima f[x_0,x_1], f[x_1,x_2] ed in seconda istanza f[x_0], f[x_1] ed f[x_2].
Esempio 8.2 Interpoliamo i punti \{(-1,1),(2,7),(0,1)\}. La costruzione della tavola delle differenze divise è la seguente:
\begin{matrix} x_0 & f[x_0]\\ & & f[x_0,x_1]\\ x_1 & f[x_1] & & f[x_0,x_1,x_2]\\ & & f[x_1,x_2] \\ x_2 & f[x_2] \\ \end{matrix}
che calcolata sui nodi specifici dell’esempio dà
\begin{matrix} -1\; & {\color{red}\bm{1}} \\ & & f[x_0,x_1]=\dfrac{1-7}{-1-2}={\color{red}\bm{2}}\\ 2\; & 7 & & f[x_0,x_1,x_2]=\dfrac{2-3}{-1-0}={\color{red}\bm{1}}\\ & & f[x_1,x_2]= \dfrac{7-1}{2-0}=3 \\ 0 \; & 1 \\ \end{matrix}
I coefficienti in rosso sul triangolo della tavola sono i coefficienti del polinomio interpolatore cercato, che si scriverà come combinazione lineare dei polinomi della base \omega_k:
\begin{aligned} \omega_0(x) &= 1\\ \omega_1(x) &= (x+1)\\ \omega_2(x) &= \omega_1(x)(x-2)=(x+1)(x-2)=x^2-x-2. \end{aligned}
Il polinomio interpolatore è quindi
\begin{aligned} p_2(x) &= {\color{red}\bm{1}}\omega_0(x)+{\color{red}\bm{2}}\omega_1(x)+{\color{red}\bm{1}}\omega_2(x)\\ &= 1+2(x+1)+1(x^2-x-2)= 1+x+x^2. \end{aligned}
Esempio 8.3 Vediamo un altro esempio con più punti da interpolare: \{(0,13),(-1,25),(-2,51),(1,3),(3,61)\}. La tavola delle differenze divise è:
\begin{matrix} 0 & f[x_0]={\color{red}\bm{13}}\\ & & \dfrac{13-25}{0-(-1)}={\color{red}\bm{-12}}\\ -1 & f[x_1]=25 & & \dfrac{-12-26}{0-(-2)}={\color{red}\bm{7}}\\ & & \dfrac{25-51}{-1-(-2)}=-26 & & \dfrac{7-5}{0-1}={\color{red}\bm{-2}}\\ -2 & f[x_2]=51 & & \dfrac{-26-(-16)}{-1-1}=5 & & \dfrac{-2-1}{0-3}={\color{red}\bm{1}}\\ & & \dfrac{51-3}{-2-1}=-16 & & \dfrac{5-9}{-1-3}=1\\ 1 & f[x_3]=3 & & \dfrac{-16-29}{-2-3}=9\\ & & \dfrac{3-61}{1-3}=29\\ 3 & f[x_4]=61 \end{matrix}
Il polinomio di interpolazione è allora
p_4(x) = 13-12\omega_1(x)+7\omega_2(x)-2\omega_3(x)+1\omega_4(x).
8.2 Invarianza per permutazioni dei nodi
La differenza divisa f[x_0,\ldots,x_k] non dipende dall’ordine in cui si considerano i nodi di interpolazione. Sia dunque \sigma\in S_{k+1} una permutazione degli interi 0,1,2,\ldots,k, dimostriamo che
f[x_0,\ldots,x_k]=f[x_{\sigma(0)},\ldots,x_{\sigma(k)}]
Si procede per induzione: la base dell’induzione è per k=1 e dalla definizione si ha che
\begin{aligned} f[x_0,x_1] &= \dfrac{f(x_1)-f(x_0)}{x_1-x_0}=\dfrac{-(f(x_0)-f(x_1))}{-(x_0-x_1)}\\ &= \dfrac{f(x_0)-f(x_1)}{x_0-x_1}=f[x_1,x_0]. \end{aligned}
Per l’ipotesi induttiva, si può assumere che l’argomento è valido per k>1 e resta da dimostrare il passo induttivo che porta da k-1 a k. Possiamo scrivere il polinomio interpolatore p_k(x) che deve coincidere sulla sequenza di punti x_0,x_1,\ldots, x_{k} sia sulla sequenza permutata x_{\sigma(0)},\ldots,x_{\sigma(k)}:
\begin{aligned} p_k(x) &= f[x_0]\omega_0(x)+\ldots+f[x_0,\ldots,x_k]\omega_k(x) \\ &= f[x_{\sigma(0)}]\hat{\omega}_0(x)+\ldots+f[x_{\sigma(0)},\ldots,x_{\sigma(k)}]\hat{\omega}_k(x). \end{aligned}
Il termine di grado massimo (ovvero k) di p_k(x) è dato dal termine di grado massimo di \omega_k(x) e \hat{\omega}_k(x) ed è esattamente x^k per entrambi, in altre parole i polinomi \omega_j (e \hat{\omega}_j) sono polinomi monici. Allora per il principio di sovrapposizione dei polinomi, il coefficiente che moltiplica x^k di p_k(x) deve essere lo stesso, e quindi f[x_0,\ldots,x_k]=f[x_{\sigma(0)},\ldots,x_{\sigma(k)}].
8.3 Formula ricorrente
In questo paragrafo dimostriamo la formula ricorrente 8.5 per le differenze divise che permette di costruire la tavola delle differenze divise. Questa dimostrazione si basa sull’invarianza per permutazioni mostrata in precedenza. Scriviamo il polinomio p(x) che interpola i nodi
x_0,\ldots,x_{k-2},x_{k-1},x_k
e il polinomio \hat{p}(x) che interpola gli stessi nodi ma con gli ultimi due scambiati, ovvero
x_0,\ldots,x_{k-2},x_k,x_{k-1}
In questa situazione abbiamo
\begin{aligned} \omega_0(x) &= 1 \\ \omega_1(x) &= \omega_0(x)(x-x_0) \\ & \vdots \\ \omega_{k-2}(x)&= \omega_{k-3}(x)(x-x_{k-3}) \\ \omega_{k-1}(x)&=\omega_{k-2}(x)(x-x_{k-2}) \\ \color{blue}\omega_k(x) &\color{blue}=\omega_{k-1}(x)(x-{\color{red}x_{k-1}}) \end{aligned}
\begin{aligned} \hat\omega_0(x) &= 1 \\ \hat\omega_1(x) &= \omega_0(x)(x-x_0) \\ & \vdots \\ \hat\omega_{k-2}(x) &= \hat\omega_{k-3}(x)(x-x_{n-3}) \\ \hat\omega_{k-1}(x) &=\hat\omega_{k-2}(x)(x-x_{k-2}) \\ \color{blue}\hat\omega_k(x) &\color{blue} = \hat\omega_{k-1}(x)(x-{\color{red}x_{k-2}}) \end{aligned}
e quindi
\omega_j(x)=\hat\omega_j(x),\qquad j=0,1,\ldots,k-1
Ovviamente abbiamo già visto che il polinomio interpolante è lo stesso, cioè p(x)=\hat p(x) scriviamolo per esteso:
\begin{aligned} \hat p(x) &= f[x_0]\omega_0(x) \\ &+f[x_0,x_1]\omega_1(x) \\ & \vdots \\ & +f[x_0,\ldots,x_{k-2}] \omega_{k-2}(x)\\ & +f[x_0,\ldots,x_{k-2},x_{k-1}] \omega_{k-1}(x)\\ & +f[x_0,\ldots,x_{k-1},x_k]\omega_k(x) \end{aligned}
\begin{aligned} p_k(x) &= f[x_0]\omega_0(x) \\ &+ f[x_0,x_1]\omega_1(x) \\ & \vdots \\ &+f[x_0,\ldots,x_{k-2}]\hat{\omega}_{k-2}(x)\\ &+f[x_0,\ldots,x_{k-2},x_{k}]\hat{\omega}_{k-1}(x)\\ &+f[x_0,\ldots,x_k,x_{k-1}]\hat{\omega}_k(x). \end{aligned}
Sottraendo membro a membro, i primi k-2 addendi sono uguali in entrambi i casi, rimane solamente
\begin{aligned} 0 &= f[x_0,\ldots,x_{k-2},x_{k-1}]\omega_{k-1}(x)+ f[x_0,\ldots,x_{k-1},x_k]\omega_k(x) \\ & -f[x_0,\ldots,x_{k-2},\,x_{k}]\hat{\omega}_{k-1}(x)- f[x_0,\ldots,x_k,x_{k-1}]\hat{\omega}_k(x). \end{aligned} \tag{8.6}
L’espressione precedente si può semplificare osservando che
\begin{aligned} \omega_{k-1}(x)&=\hat{\omega}_{k-1}(x),\\ \omega_k(x)&=\omega_{k-1}(x)(x-x_{k-1}),\\ \hat\omega_k(x)&=\omega_{k-1}(x)(x-x_{k-2}), \end{aligned}
e la 8.6 diventa, dopo aver diviso per il fattore comune \omega_{k-1},
\begin{aligned} 0 &= f[x_0,\ldots,x_{k-2},x_{k-1}] +f[x_0,\ldots,x_{k-1},x_k](x-x_{k-1}) \\ & -f[x_0,\ldots,x_{k-2},\,x_{k}] -f[x_0,\ldots,x_k,x_{k-1}](x-x_k). \end{aligned}
Per la proprietà di invarianza per permutazione, l’ultimo addendo dell’espressione precedente diventa f[x_0,\ldots,x_{k-1},x_k](x-x_k) e si può mettere a fattore comune con il secondo addendo, così da ottenere:
\begin{aligned} f[x_0,\ldots,x_{k-2},x_{k-1}] -f[x_0,\ldots,x_{k-2},\,x_{k}] &= f[x_0,\ldots,x_{k-1},x_k](x-x_k-(x-x_{k-1}))\\ &= f[x_0,\ldots,x_k](x_{k-1}-x_k). \end{aligned}
Ora possiamo raccogliere il termine f[x_0,\ldots,x_k] e scrivere
f[x_0,\ldots,x_k] = \dfrac{f[x_0,\ldots,x_{k-2},x_{k-1}] -f[x_0,\ldots,x_{k-2},\,x_{k}]}{x_{k-1}-x_k}. \tag{8.7}
Questa è già una formula ricorrente valida, tuttavia non è scritta nel modo standard, per ricondurla alla 8.5 basta operare una nuova permutazione dei nodi. Chiamiamo per comodità y la sequenza di nodi x_0,\ldots,x_{k-2}, a=x_{k-1} e b=x_k, allora l’espressione 8.7 si riscrive come
f[y,a,b]=\dfrac{f[y,a]-f[y,b]}{a-b},
e permutando l’ordine possiamo scrivere
f[a,y,b]=\dfrac{f[a,y]-f[y,b]}{a-b}
che è proprio la ricorrenza 8.5.
Riassunto sulle differenze divise
Il sistema lineare che definisce il problema di interpolazione secondo il metodo di Newton
\left\{\begin{aligned} f(x_0) &= a_0 \omega_{k+1} \\ f(x_1) &= a_0 + a_1 \omega_1 (x_1) \omega{k+1} \\ f( x_2 ) &= a_0 + a_1 \omega_1 ( x_2 ) +a_2 \omega_2 ( x_2 ) \omega_{k+1} \\ &\vdots \\ f(x_n) &= a_0 + a_1 \omega_1 (x_n) +\cdots + a_n \omega_n (x_n) \end{aligned}\right.
è triangolare inferiore e può essere risolto esplicitamente con la tecnica delle sostituzioni in avanti.
La soluzione del sistema lineare può essere scritta come
\begin{aligned} a_0 &= f(x_0), \\ a_1 &= \dfrac{f(x_1)-f(x_0)}{\omega_1 (x_1)}, \\ a_2 &= \dfrac{f( x_2 )-f(x_0)-\omega_1 (x_1) \dfrac{f(x_1)-f(x_0)}{\omega_1 (x_1)}} {\omega_2 ( x_2 )}, \\ &\vdots \end{aligned}
Si noti che ogni a_i dipende dai nodi di interpolazione ( x_j,f( x_j)) per j=0,1,\ldots i e non dai successivi.
Indichiamo questa dipendenza con il simbolo speciale (notare le parentesi quadre)
a_j = f[x_0,x_1,\ldots, x_j],
che si chiama differenza divisa di ordine j.
Tramite le differenze divise possiamo scrivere il polinomio interpolante come segue
\boxed{p(x) = \sum_{k=0}^n f[x_0,x_1,\ldots,x_k ] \omega_k (x )}
Le differenze divise godono di una serie notevole di proprietà interessanti.
Proprietà 1: La differenza divisa f[x_0,x_1,\ldots,x_k ] è il coefficiente del termine di grado massimo k del polinomio interpolatore che interpola i k+1 nodi di ascisse x_0, x_1, \ldots, x_k.
Proprietà 2: La differenza divisa f[x_0,x_1,\ldots,x_k ] non dipende dall’ordine con cui si prendono i nodi di interpolazione, cioè f[x_0,x_1,\ldots,x_k ] = f[x_{i_0},x_{i_1},\ldots,x_{i_k}],
dove i_0,i_1,\ldots,i_k è una qualunque permutazione dei numeri 0,1,\ldots,k.
Proprietà 3 Definizione ricorsiva delle differenze divise di ogni ordine
f[ x_i] = f( x_i),\qquad f[\ldots, x_i, x_j] = \dfrac{f[\ldots, x_i]-f[\ldots, x_j]}{ x_i- x_j}.
Il procedimento ricorsivo è quello illustrato in Figura 8.1.
Costi computazionali
In conclusione, abbiamo visto tre modi di calcolare il polinomio interpolatore dei dati, quello naïve, quello di Lagrange e quello di Newton, rispettivamente:
\begin{aligned} p_n(x)&= a_0+a_1x+\ldots+a_nx^n & \qquad & \text{[sistema lineare]}\\ p_n(x)&= b_0\mathcal{L}_0(x)+b_1\mathcal{L}_1(x)+\ldots+b_n\mathcal{L}_n(x) & \qquad & \text{[Lagrange]}\\ p_n(x)&= c_0+c_1\omega_1(x)+\ldots+c_n\omega_n(x) & \qquad & \text{[Newton]} \end{aligned}
In tutti e tre i casi, il polinomio è espresso come combinazione lineare di una base. Il primo approccio porta a un polinomio semplice da scrivere, ma richiede la soluzione di un sistema lineare e non consente di aggiornare efficientemente il risultato quando si aggiunge un nodo. Il costo computazionale è quindi quello della soluzione del sistema, pari a circa \mathcal{O}(n^3).
Il metodo di Lagrange impiega polinomi base \mathcal{L}_i(x) la cui espressione è articolata; consente tuttavia di scrivere direttamente la soluzione, pur non offrendo una costruzione incrementale. Per calcolare il polinomio interpolatore occorre formare n+1 polinomi \mathcal{L}_i(x); ciascuno richiede n+1 somme e n prodotti sia al numeratore sia al denominatore. Il costo complessivo è dell’ordine di \mathcal{O}(4n^2).
Infine, il metodo di Newton consente una costruzione incrementale con un costo computazionale contenuto. L’uso delle differenze divise rende il procedimento ancora più efficiente. Può però verificarsi una perdita di cifre significative per cancellazione, quando, con l’infittirsi dei nodi di interpolazione, si sottraggono valori molto vicini.
Il calcolo con il metodo delle differenze divise ha una complessità di \mathcal{O}(\dfrac{3}{2}n^2), ma ha complessità lineare se da p_n si vuol passare a p_{n+1} e non serve ripercorrere tutto il calcolo precedentemente fatto per p_n.
In generale una funzione f si può approssimare come una combinazione lineare del tipo
f(x)=a_0b_0(x)+a_1b_1(x)+\ldots+a_nb_n(x)+E_n(x),
che è una combinazione lineare degli elementi della base \{b_i\}_{i=0}^n più un errore E_n. Abbiamo visto il caso in cui la base è una base polinomiale canonica 1,x,x^2,\ldots, il caso della base di Lagrange \mathcal{L}_i(x) e la base di Newton \omega_i(x). Ci sono anche altre basi polinomiali, ma si possono scegliere anche altre forme funzionali, molto usate sono basi di funzioni esponenziali, trigonometriche o funzioni razionali.
8.4 L’algoritmo di Aitken-Neville
Per presentare l’algoritmo di Aitken-Neville 3
introduciamo l’insieme formato da n+1 punti di \mathbb{R}\times\mathbb{R}
\bm{I} = \{ (x_i,f(x_i))\ | \ i=0,1,2,\ldots,n,\, x_i\neq x_j \},
e tutti i suoi sottoinsiemi
I_{i,i+1,i+2,\ldots,j} = \{ (x_{k},f(x_{k}))\ | \ k=i,i+1,i+2,\ldots,j \},
Ad ogni sottoinsieme I_{0,1,2,\ldots,k} associamo il polinomio p_{0,1,2,\ldots,k} di grado k che interpola i punti dell’insieme
p_{0,1,2,\ldots,k}(x_j) = f(x_j) \qquad j=0,1,\ldots,k.
Osservazione 8.1. Conoscendo i due polinomi interpolatori
- p_{0,1,2,\ldots,k-1}(x) e
- p_{1,2,\ldots,k-1,k}(x)
possiamo facilmente costruire il polinomio interpolatore
- p_{0,1,2,\ldots,k-1,k}(x),
che è dato dalla formula
p_{0,1,\ldots,k}(x) = \dfrac{(x_k-x) p_{0,1,\ldots,k-1}(x)+ (x-x_0) p_{1,2,\ldots,k}(x)} {x_k-x_0}
Per s=1,2,\ldots,k-1 vale ovviamente
\begin{aligned} p_{0,1,\ldots,k}(x_s) & = \dfrac{(x_k-x_s) p_{0,1,\ldots,k-1}(x_s)+ (x_s-x_0) p_{1,2,\ldots,k}(x_s)} {x_k-x_0}, \\ &= \dfrac{(x_k-x_s) f(x_s)+(x_s-x_0) f(x_s)} {x_k-x_0}, \\ & = f(x_s), \end{aligned}
Nel caso s=0 abbiamo
\begin{aligned} p_{0,1,\ldots,k}(x_0) & = \dfrac{(x_k-x_0) p_{0,1,\ldots,k-1}(x_0)+ (x_0-x_0) p_{1,2,\ldots,k}(x_0)} {x_k-x_0}, \\ & = \dfrac{(x_k-x_0) f(x_0)}{x_k-x_0} = f(x_0). \end{aligned}
Nel caso s=k abbiamo
\begin{aligned} p_{0,1,\ldots,k}(x_k) & = \dfrac{(x_k-x_k) p_{0,1,\ldots,k-1}(x_k)+ (x_k-x_0) p_{1,2,\ldots,k}(x_k)} {x_k-x_0}, \\ & = \dfrac{(x_k-x_0) f(x_k)}{x_k-x_0} = f(x_k). \end{aligned}
Osservazioni finali sull’algoritmo di Aitken-Neville
I polinomi p_i(x) con i=0,1,\ldots,n sono semplicemente le costanti f(x_0), f(x_1),\ldots f(x_n); cioè:
p_i(x) = f(x_i),\qquad i=0,1,\ldots,n
Analogamente
p_{i,i+1}(x) = \dfrac{(x_{i+1}-x)p_i(x) + (x-x_i)p_{i+1}(x)} {x_{i+1}-x_{i}},
per i=0,1,\ldots,n-1.
La formula ricorrente si usa in modo efficiente per calcolare il valore del polinomio interpolatore in un punto generico \bar{x} senza mai costruirne i coefficienti: si riempie una tavola triangolare colonna per colonna, partendo dai valori nei nodi.
Lo schema di calcolo è lo stesso triangolo delle differenze divise (Figura 8.1), con i polinomi al posto delle differenze:
In sintesi
- La forma di Newton costruisce il polinomio incrementalmente: aggiungere un nodo costa \mathcal{O}(n) e non richiede di rifare i conti.
- Le differenze divise sono invarianti per permutazione dei nodi e si calcolano con una ricorrenza a triangolo, in \mathcal{O}(n^2).
- Vale f[t_0,\ldots,t_n]=\dfrac{f^{(n)}(\zeta)}{n!}: le differenze divise sono derivate discrete, e questo spiega la forma dell’errore.
- L’algoritmo di Aitken-Neville valuta il polinomio in un punto senza costruirne i coefficienti: preferibile quando serve un solo valore.
Alexander Craig Aitken 1895–1967↩︎