Interpolazione di Newton
Appunti di calcolo numerico
Il problema con l’approccio di Lagrange è che una volta calcolato il polinomio p_n(x) che interpola i dati
\{(x_i,y_i)\}_{i=0}^n
se vogliamo calcolare p_{n+1}(x) che approssima un punto in più, dobbiamo ripetere tutto il procedimento, e non si può calcolare p_{n+1} a partire da p_n. L’idea di Newton è proprio quella di rendere possibile questa costruzione 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{1}
Esempio 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-21) \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 la si può vedere ancora come un sistema lineare, ma non più rettangolare come quello visto all’inizio del capitolo, ma di forma triangolare (inferiore), infatti si può scrivere:
\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{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{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.
Proprietà delle differenze divise
Dalla relazione Equazione 3 e dalla definizione Equazione 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{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:
\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] & & f[x_0,x_1,x_2,x_3]\\ x_2 & f[x_2] & & f[x_1,x_2,x_3] & & f[x_0,x_1,x_2,x_3,x_4]\\ & & f[x_2,x_3] & & f[x_1,x_2,x_3,x_4]\\ x_3 & f[x_3] & & f[x_2,x_3,x_4]\\ & & f[x_3,x_4]\\ x_4 & f[x_4] \end{matrix}
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 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 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}}\\ & & \frac{13-25}{0-(-1)}={\color{red}\bm{-12}}\\ -1 & f[x_1]=25 & & \frac{-12-26}{0-(-2)}={\color{red}\bm{7}}\\ & & \frac{25-51}{-1-(-2)}=-26 & & \frac{7-5}{0-1}={\color{red}\bm{-2}}\\ -2 & f[x_2]=51 & & \frac{-26-(-16)}{-1-1}=5 & & \frac{-2-1}{0-3}={\color{red}\bm{1}}\\ & & \frac{51-3}{-2-1}=-16 & & \frac{5-9}{-1-3}=1\\ 1 & f[x_3]=3 & & \frac{-16-29}{-2-3}=9\\ & & \frac{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).
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)}].
Formula ricorrente
In questo paragrafo dimostriamo la formula ricorrente Equazione 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{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 Equazione 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}] \\ &\qquad= f[x_0,\ldots,x_{k-1},x_k](x-x_k-(x-x_{k-1}))\\ &\qquad=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{7}
Questa è già una formula ricorrente valida, tuttavia non è scritta nel modo standard, per ricondurla alla Equazione 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 Equazione 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 Equazione 5.
Riassunto sulle differenze divise
Il sistema lineare che definisce il problema di interpolazione secondo il metodo di Newton
\bgroup \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.\egroup
è 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}.
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 & \textrm{[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 & \textrm{[Lagrange]}\\ p_n(x)&= c_0+c_1\omega_1(x)+\ldots+c_n\omega_n(x) & \qquad & \textrm{[Newton]} \end{aligned}
Essi rappresentano tutti una combinazione lineare di polinomi, cioé una somma pesata di coefficienti. Il primo modo dà un polinomio facile da scrivere ma costa la risoluzione di un sistema lineare e non permette una costruzione efficiente, inoltre non è efficiente se bisogna aggiungere un nodo. Il suo costo computazionale è quello della soluzione di un sistema lineare, e quindi di circa \mathcal{O}(n^3).
Il metodo di Lagrange usa polinomi complicati (gli \mathcal{L}_i(x)), tuttavia la soluzione è immediata anche se non è incrementale. Il calcolo del polinomio interpolatore di Lagrange costa n+1 polinomi \mathcal{L}_i(x), ciascuno richiede n+1 somme e n prodotti, sia al numeratore che al denominatore, per una complessità di \mathcal{O}(4n^2).
Infine, il metodo di Newton offre una costruzione incrementale e economica dal punto di vista computazionale. Essa si può ulteriormente migliorare applicando il metodo delle differenze divise. Un problema che si presenta è però quello di perdita di cifre significative per sottrazione, che si verifica quando si sottraggono valori simili all’infittirsi dei nodi di interpolazione.
Il calcolo con il metodo delle differenze divise ha una complessità di \mathcal{O}(\frac{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.
L’algoritmo di Aitken-Neville
Per presentare l’algoritmo di Aitken-Neville 1 introduciamo l’insieme formato da n+1 punti di \Bbb{R}\times\Bbb{R}.
I = \{ ( x_i,f( x_i))\ | \ i=0,1,2,\ldots,n,\, x_i\neq x_j \},
e tutti i suoi sottoinsiemi
I_{i_0,i_1,i_2,\ldots,i_k} = \{ (x_{i_s},f(x_{i_{s}}))\ | \ s=0,1,2,\ldots,k \},
con k=0,1,\ldots,n.
Ad ogni sottoinsieme \boldsymbol{I}[i_0,i_1,i_2,\ldots,i_k] associamo il polinomio p_{i_0,i_1,i_2,\ldots,i_k} di grado k che interpola i punti dell’insieme
p_{i_0,i_1,i_2,\ldots,i_k}(z_{i_j}) = f(x_{i_j}) \qquad j=0,1,\ldots,k.
Per s=1,2,\ldots,k-1 vale ovviamente \begin{aligned} p_{i_0,i_1,\ldots,i_k}(x_{i_s}) & = \dfrac{(x_{i_k}-x_{i_s})p_{i_0,i_1,\ldots,i_{k-1}}(x_{i_s})+ (x_{i_s}-x_{i_0})p_{i_1,i_2,\ldots,i_{k}}(x_{i_s})} {x_{i_k}-x_{i_0}}, \\ &= \dfrac{(x_{i_k}-x_{i_s}) f(x_{i_s})+(x_{i_s}-x_{i_0}) f(x_{i_s})} {x_{i_k}-x_{i_0}}, \\ & = f(x_{i_s}), \end{aligned}
\boxed{ p_{i_0,i_1,\ldots,i_k}(x) = \dfrac{(x_{i_k}-x)p_{i_0,i_1,\ldots,i_{k-1}}(x )+ (x-x_{i_0})p_{i_1,i_2,\ldots,i_{k}}(x)} {x_{i_k}-x_{i_0}} }
Nel caso s=0 abbiamo \begin{aligned} p_{i_0,i_1,\ldots,i_k}(x_{i_0}) & = \dfrac{(x_{i_k}-x_{i_0})p_{i_0,i_1,\ldots,i_{k-1}}(x_{i_0})+ (x_{i_0}-x_{i_0})p_{i_1,i_2,\ldots,i_{k}}(x_{i_0})} {x_{i_k}-x_{i_0}}, \\ & = \dfrac{(x_{i_k}-x_{i_0}) f(x_{i_0})}{x_{i_k}-x_{i_0}} = f(x_{i_0}). \end{aligned}
\boxed{ p_{i_0,i_1,\ldots,i_k}(x) = \dfrac{(x_{i_k}-x)p_{i_0,i_1,\ldots,i_{k-1}}(x)+ (x-x_{i_0})p_{i_1,i_2,\ldots,i_{k}}(x)} {x_{i_k}-x_{i_0}} }
Nel caso s=k abbiamo
\begin{aligned} p[i_0,i_1,\ldots,i_k](x_{i_k}) & = \dfrac{(x_{i_k}-x_{i_k}) p_{i_0,i_1,\ldots,i_{k-1}}(x_{i_k})+ (x_{i_k}-x_{i_0}) p_{i_1,i_2,\ldots,i_{k}}(x_{i_k})} {x_{i_k}-x_{i_0}}, \\ & = \dfrac{(x_{i_k}-x_{i_0}) f(x_{i_k})}{x_{i_k}-x_{i_0}} = f(x_{i_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 di interpolazione dell’algoritmo di Aitken-Neville può essere utilizzata in modo efficiente per calcolare il valore del polinomio interpolatore in un punto generico \bar{x} come segue:
Per i=0,1,\ldots,n si pone p_i(\bar{x})=f( x_i).
Per k=1,2,\ldots,n si calcola
Per i=0,1,\ldots,n-k p_{i,\ldots,i+k}(\bar{x}) = \dfrac{(x_{i+k}-\bar{x})p_{i,\ldots,i+k-1}(\bar{x}) + (\bar{x}-x_i)p_{i+1,\ldots,i+k}(\bar{x})} {x_{i+k}-x_i}.
p_{0,1,\ldots,n}(\bar{x}) restituisce il valore del polinomio interpolatore calcolato in \bar{x}.