Interpolazione polinomiale - forma di Lagrange, differenze divise ed errore

Prerequisiti: Teoremi di Rolle e di LagrangeRolle: f continua su [a,b], derivabile in ]a,b[, f(a) = f(b) ⇒ esiste c in ]a,b[ con f'(c) = 0. Lagrange (valor medio): stesse ipotesi senza f(a) = f(b) ⇒ esiste c con f'(c) = (f(b) − f(a))/(b − a): la tangente in c è parallela alla corda. Conseguenza: su un intervallo f è costante se e solo se f' = 0; su un insieme che non è un intervallo è solo costante a tratti (arctan x + arctan(1/x)).Teoremi di Rolle e di Lagrange → (per l'errore), Formula di Taylor con resto di LagrangeSe f ha n+1 derivate continue su [a,b], l'errore tra f e il suo polinomio di Taylor si scrive esattamente come f^(n+1)(c)(x-x0)^(n+1)/(n+1)! per un c opportuno: serve a stimare l'errore su tutto un intervallo e a scegliere il grado del polinomio.Formula di Taylor con resto di Lagrange → (caso a nodi coincidenti), Eliminazione di GaussCon tre operazioni elementari sulle righe (scambio, moltiplicazione per uno scalare non nullo, somma di un multiplo di un'altra riga) ogni matrice si riduce a scala senza cambiare il rango; serve a calcolare ranghi, risolvere sistemi, trovare relazioni di dipendenza e matrici che riducono a scala.Eliminazione di Gauss → (sistema di Vandermonde).

Problema. Di una funzione $f$ si conoscono i valori $y_i = f(x_i)$ in $n+1$ punti distinti $x_0 < x_1 < \dots < x_n$ (nodi), ad esempio perché sono misure o una tabella. Si cerca una funzione semplice, un polinomio, che passi per questi punti, per stimare $f$ in altri punti, e poi derivarla o integrarla (Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di RichardsonUna formula di quadratura approssima $\int_a^b f,dx$ con $\sum w_if(x_i)$ ottenuta integrando il polinomio che interpola $f$ nei nodi. Trapezio $\frac{b-a}2[f(a)+f(b)]$ (grado di esattezza 1, errore $-\frac{(b-a)^3}{12}f''(\xi)$), punto medio $(b-a)f(m)$ (grado 1, errore $+\frac{(b-a)^3}{24}f''(\xi)$), Cavalieri-Simpson $\frac{b-a}6[f(a)+4f(m)+f(b)]$ (grado 3, errore $-\frac{(b-a)^5}{2880}f^{(4)}(\xi)$). Composite con passo $h=\frac{b-a}n$: trapezi errore $-\frac{(b-a)h^2}{12}f''(\xi)$ (ordine 2, rapporto tra errori 4 se $h$ si dimezza), Simpson errore $-\frac{(b-a)h^4}{180}f^{(4)}(\xi)$ (ordine 4, rapporto 16, $n$ pari). Estrapolazione di Richardson: se $Q(h)=I+Ch^p+\dots$ allora $I\approx\frac{2^pQ(h/2)-Q(h)}{2^p-1}$; per i trapezi $p=2$ e $\frac{4T_{h/2}-T_h}3=S_{h/2}$; ripetuta dà la tabella di Romberg. Stima dell'errore a posteriori: $I-Q(h/2)\approx\frac{Q(h/2)-Q(h)}{2^p-1}$.Formule di quadratura - trapezi, Cavalieri-Simpson ed estrapolazione di Richardson →). Un polinomio è facile da valutare; un polinomio di grado $n$ è determinato da $n+1$ numeri, proprio come i dati.

Teorema (esistenza e unicità). Per $n+1$ coppie $(x_i, y_i)$ con $x_i$ a due a due distinti esiste un unico polinomio $p_n$ di grado $\le n$ con $p_n(x_i) = y_i$ per ogni $i$.

Dimostrazione. Scrivendo $p_n(x) = a_0 + a_1x + \dots + a_nx^n$, le condizioni $p_n(x_i) = y_i$ formano un sistema lineare $n+1$ per $n+1$ nelle incognite $a_j$, con matrice di Vandermonde $V_{ij} = x_i^{\,j}$, il cui determinante è $\prod_{i<j}(x_j - x_i) \ne 0$ perché i nodi sono distinti. Quindi il sistema ha una e una sola soluzione. (Unicità anche direttamente: la differenza di due soluzioni è un polinomio di grado $\le n$ con $n+1$ zeri, quindi nulla.) ∎

Risolvere il sistema di Vandermonde è però mal condizionato (la matrice ha numero di condizionamento enorme: circa $7\cdot10^2$ per 5 nodi equispaziati in $[0,1]$, $1{,}5\cdot10^7$ per 10 nodi, $4\cdot10^{11}$ per 15, $10^{16}$ per 20; vedi Condizionamento di un problema e stabilità di un algoritmoL'errore finale ha tre origini (dati, troncamento, arrotondamento). Un problema è ben posto se la soluzione esiste, è unica e dipende con continuità dai dati; il numero di condizionamento K = |x f'(x)/f(x)| misura di quanto un errore relativo sui dati si amplifica nel risultato (si perdono circa log10 K cifre). La stabilità riguarda invece l'algoritmo: un algoritmo instabile amplifica gli errori anche su un problema ben condizionato (ricorrenza I_n = 1 − n I_{n−1}); un problema mal condizionato non si salva con nessun algoritmo. Uno schema numerico converge se è consistente e stabile (derivata numerica: troncamento contro arrotondamento).Condizionamento di un problema e stabilità di un algoritmo →). Per questo si usano due altre scritture dello stesso polinomio.

Forma di Lagrange

Si costruiscono i polinomi di grado $n$ $$L_i(x) = \prod_{\substack{j=0\\ j\ne i}}^{n}\frac{x - x_j}{x_i - x_j}, \qquad L_i(x_j) = \delta_{ij} = \begin{cases}1 & j = i\\ 0 & j \ne i\end{cases}$$ (il numeratore si annulla in tutti i nodi tranne $x_i$; il denominatore normalizza a $1$ in $x_i$). Allora $$p_n(x) = \sum_{i=0}^{n} y_i\,L_i(x).$$ Infatti $p_n(x_j) = \sum_i y_i\delta_{ij} = y_j$, e il grado è $\le n$. Dato che $\sum_i L_i(x) \equiv 1$ (interpola la funzione costante $1$), i $L_i$ sono "pesi" che sommano a $1$.

Esempio svolto. $f(x) = e^x$ con nodi $x_0 = 0$, $x_1 = 0{,}5$, $x_2 = 1$, valori $y = 1;\ 1{,}648721;\ 2{,}718282$. Si stima $f(0{,}25)$. Con $x = 0{,}25$:

• $L_0 = \dfrac{(x - 0{,}5)(x-1)}{(0 - 0{,}5)(0-1)} = \dfrac{(-0{,}25)(-0{,}75)}{0{,}5} = 0{,}375$;
• $L_1 = \dfrac{(x - 0)(x-1)}{(0{,}5 - 0)(0{,}5-1)} = \dfrac{(0{,}25)(-0{,}75)}{-0{,}25} = 0{,}75$;
• $L_2 = \dfrac{(x - 0)(x-0{,}5)}{(1 - 0)(1 - 0{,}5)} = \dfrac{(0{,}25)(-0{,}25)}{0{,}5} = -0{,}125$.

Somma dei pesi: $0{,}375 + 0{,}75 - 0{,}125 = 1$ ✓. Quindi $$p_2(0{,}25) = 0{,}375\cdot 1 + 0{,}75\cdot 1{,}648721 - 0{,}125\cdot 2{,}718282 = 0{,}375 + 1{,}236541 - 0{,}339785 = 1{,}271756.$$ Il valore esatto è $e^{0{,}25} = 1{,}284025$: errore $0{,}01227$.

La forma di Lagrange è comoda per la teoria, ma se si aggiunge un nodo si deve ricalcolare tutto, e valutarla costa $O(n^2)$.

Forma di Newton e differenze divise

Definizione (differenze divise). $f[x_i] = f(x_i)$ e ricorsivamente $$f[x_i, \dots, x_{i+k}] = \frac{f[x_{i+1}, \dots, x_{i+k}] - f[x_i, \dots, x_{i+k-1}]}{x_{i+k} - x_i}.$$

Teorema (forma di Newton). Il polinomio di interpolazione è $$p_n(x) = f[x_0] + f[x_0,x_1](x - x_0) + f[x_0,x_1,x_2](x-x_0)(x-x_1) + \dots + f[x_0,\dots,x_n]\prod_{j=0}^{n-1}(x - x_j).$$

Perché funziona: $p_k = p_{k-1} + c_k\prod_{j<k}(x - x_j)$ aggiunge ogni volta un termine che si annulla nei nodi precedenti, quindi non rovina l'interpolazione già fatta; il coefficiente $c_k$ è scelto per far passare il nuovo nodo. Proprietà utili:

• Le differenze divise sono simmetriche nei nodi (non dipendono dall'ordine).
• $f[x_0,\dots,x_k] = \dfrac{f^{(k)}(\xi)}{k!}$ per un certo $\xi$ nel più piccolo intervallo contenente i nodi (si ricava dall'errore sotto). In particolare, se $f$ è un polinomio di grado $d$, tutte le differenze divise di ordine $> d$ sono nulle, e quelle di ordine $d$ valgono il coefficiente direttivo.
• Aggiungere un nodo costa solo una nuova riga della tabella: si aggiunge un termine alla forma di Newton senza toccare i precedenti.

Esempio svolto: tabella delle differenze divise di $e^x$. Per i nodi $0$, $0{,}5$, $1$: $$f[0,0{,}5] = \frac{1{,}648721 - 1}{0{,}5} = 1{,}297443,\qquad f[0{,}5,1] = \frac{2{,}718282 - 1{,}648721}{0{,}5} = 2{,}139121,$$ $$f[0,0{,}5,1] = \frac{2{,}139121 - 1{,}297443}{1 - 0} = 0{,}841679.$$ Allora $p_2(x) = 1 + 1{,}297443\,x + 0{,}841679\,x(x - 0{,}5)$ e $$p_2(0{,}25) = 1 + 0{,}324361 + 0{,}841679\cdot 0{,}25\cdot(-0{,}25) = 1{,}324361 - 0{,}052605 = 1{,}271756,$$ lo stesso valore di Lagrange (stesso polinomio, scritto in un'altra base).

Caso limite (nodi coincidenti). Se tutti i nodi tendono a $x_0$, $f[x_0,\dots,x_k] \to \frac{f^{(k)}(x_0)}{k!}$ e $p_n$ diventa il polinomio di Taylor; la formula dell'errore diventa il resto di Lagrange. Quindi Taylor è il caso limite dell'interpolazione.

Errore di interpolazione

Teorema (errore). Sia $f \in C^{n+1}([a,b])$, con nodi $x_0,\dots,x_n \in [a,b]$. Per ogni $x \in [a,b]$ esiste $\xi = \xi(x) \in\,]a,b[$ tale che $$f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,\omega_{n+1}(x), \qquad \omega_{n+1}(x) = \prod_{i=0}^{n}(x - x_i).$$

Dimostrazione. Se $x$ è un nodo entrambi i membri sono nulli. Sia $x$ non nodo e si pone $$g(t) = f(t) - p_n(t) - K\,\omega_{n+1}(t), \qquad K = \frac{f(x) - p_n(x)}{\omega_{n+1}(x)}.$$ $g$ si annulla negli $n+1$ nodi (perché $f = p_n$ e $\omega = 0$ lì) e anche in $t = x$ (per la scelta di $K$): ha $n+2$ zeri distinti. Applicando il teorema di Rolle (Teoremi di Rolle e di LagrangeRolle: f continua su [a,b], derivabile in ]a,b[, f(a) = f(b) ⇒ esiste c in ]a,b[ con f'(c) = 0. Lagrange (valor medio): stesse ipotesi senza f(a) = f(b) ⇒ esiste c con f'(c) = (f(b) − f(a))/(b − a): la tangente in c è parallela alla corda. Conseguenza: su un intervallo f è costante se e solo se f' = 0; su un insieme che non è un intervallo è solo costante a tratti (arctan x + arctan(1/x)).Teoremi di Rolle e di Lagrange →) ripetutamente, $g'$ ha almeno $n+1$ zeri, $g''$ ne ha $n$, …, $g^{(n+1)}$ ne ha almeno $1$, in un punto $\xi$. Ora $p_n^{(n+1)} \equiv 0$ (grado $\le n$) e $\omega_{n+1}^{(n+1)} \equiv (n+1)!$, quindi $$0 = g^{(n+1)}(\xi) = f^{(n+1)}(\xi) - K\,(n+1)! \quad\Longrightarrow\quad K = \frac{f^{(n+1)}(\xi)}{(n+1)!}$$ ∎.

Maggiorazione. Se $|f^{(n+1)}| \le M_{n+1}$ su $[a,b]$, $$|f(x) - p_n(x)| \le \frac{M_{n+1}}{(n+1)!}\,|\omega_{n+1}(x)| \le \frac{M_{n+1}}{(n+1)!}\max_{[a,b]}|\omega_{n+1}|.$$ L'errore ha due fattori: uno dipende da $f$ (la sua derivata $(n+1)$-esima), l'altro solo dai nodi (il polinomio nodale $\omega$). Il primo non si controlla; il secondo sì, scegliendo bene i nodi.

Esempio (continuazione). Per $e^x$, $n=2$, $M_3 = e^1 = 2{,}718$ su $[0,1]$. Per $x = 0{,}25$: $|\omega_3(0{,}25)| = |0{,}25\cdot(-0{,}25)\cdot(-0{,}75)| = 0{,}046875$. Maggiorazione: $\frac{2{,}718}{6}\cdot 0{,}046875 = 0{,}02124$. L'errore effettivo $0{,}01227$ è inferiore ✓.

Esempio: $f(x) = \frac1x$, nodi $1, 2, 4$, valore in $x=3$. $f[1,2] = \frac{0{,}5 - 1}{1} = -0{,}5$, $f[2,4] = \frac{0{,}25 - 0{,}5}{2} = -0{,}125$, $f[1,2,4] = \frac{-0{,}125 + 0{,}5}{3} = 0{,}125$. Il polinomio è $p_2(x) = 1 - 0{,}5(x-1) + 0{,}125(x-1)(x-2)$ e $p_2(3) = 1 - 1 + 0{,}125\cdot 2\cdot 1 = 0{,}25$, contro $\frac13 = 0{,}3333$: errore $0{,}0833$. Con la formula: $f''' = -\frac{6}{x^4}$ e $\omega_3(3) = 2\cdot1\cdot(-1) = -2$, quindi $f - p_2 = \frac{-6/\xi^4}{6}(-2) = \frac{2}{\xi^4}$; ponendolo uguale a $0{,}0833$ si trova $\xi = 2{,}21 \in\,]1,4[$ ✓. Aggiungendo il nodo $x_3 = 5$ ($y = 0{,}2$) basta una differenza divisa in più: $f[1,2,4,5] = -0{,}025$, e $p_3(3) = 0{,}25 + (-0{,}025)(3-1)(3-2)(3-4) = 0{,}25 + 0{,}05 = 0{,}30$ (più vicino a $0{,}3333$).

La scelta dei nodi: fenomeno di Runge e nodi di Chebyshev

Più nodi non significano più accuratezza. Per $f(x) = \dfrac{1}{1+25x^2}$ su $[-1,1]$ con nodi equispaziati l'errore massimo cresce con $n$:

(errore massimo su $[-1,1]$.) Con $n = 8$ equispaziati si ha $p_8(0{,}95) = -0{,}892$ mentre $f(0{,}95) = 0{,}0424$: le oscillazioni agli estremi dell'intervallo sono enormi. È il fenomeno di Runge: $\max|\omega_{n+1}|$ per nodi equispaziati è piccolo al centro e grande ai bordi, e le derivate di $f$ crescono fattorialmente.

Per ridurre il problema si usano i nodi di Chebyshev sull'intervallo $[-1,1]$: $$x_k = \cos\frac{(2k+1)\pi}{2n+2},\quad k = 0,\dots,n,$$ (le proiezioni sull'asse $x$ di punti equispaziati sulla semicirconferenza: più fitti ai bordi). Per questa scelta $\max_{[-1,1]}|\omega_{n+1}| = 2^{-n}$, che è il minimo possibile su tutte le scelte di $n+1$ nodi. Confronto di $\max|\omega|$:

Su un intervallo generico $[a,b]$ i nodi si ottengono con $\frac{a+b}{2} + \frac{b-a}{2}x_k$ e vale $$\max_{[a,b]}|\omega_{n+1}| = 2\left(\frac{b-a}{4}\right)^{n+1}, \qquad |f - p_n| \le \frac{2M_{n+1}}{(n+1)!}\left(\frac{b-a}{4}\right)^{n+1}.$$

Esempio: quanti nodi per $\sin x$ su $[0,\frac\pi2]$ con errore $< 10^{-6}$? Qui $M_{n+1} = 1$. Con nodi di Chebyshev: $\frac{2}{(n+1)!}\left(\frac\pi8\right)^{n+1}$ vale $1{,}0\cdot10^{-5}$ per $n=5$ e $5{,}7\cdot10^{-7}$ per $n=6$: servono $n = 6$ (7 nodi). Con nodi equispaziati ($|\omega| \le \frac{n!}{4}h^{n+1}$ con $h = \frac\pi{2n}$, quindi maggiorante $\frac{h^{n+1}}{4(n+1)}$): $3{,}0\cdot10^{-6}$ per $n=6$ e $2{,}0\cdot10^{-7}$ per $n=7$: servono $n = 7$.

In pratica. Si evitano i polinomi di grado alto su nodi equispaziati; per molti dati si usano polinomi a tratti (Interpolazione a tratti, spline cubiche e minimi quadratiInvece di un polinomio di grado alto si usano polinomi di grado basso su ogni intervallo. Lineare a tratti: errore ≤ h²/8·max|f''|. Spline cubica: polinomio di grado 3 su ogni intervallo, di classe C² nei nodi; i momenti M_i = S''(x_i) risolvono un sistema tridiagonale (h_{i−1}/6)M_{i−1} + ((h_{i−1}+h_i)/3)M_i + (h_i/6)M_{i+1} = differenza tra le pendenze dei due intervalli adiacenti, con due condizioni agli estremi (naturale M_0 = M_n = 0, vincolata, not-a-knot); errore O(h⁴), nessun fenomeno di Runge. Minimi quadrati: con più dati che parametri si minimizza ‖Ac − y‖², cioè si risolvono le equazioni normali AᵀA c = Aᵀy (proiezione ortogonale di y sulle colonne di A); la retta di regressione minimizza la somma dei quadrati dei residui e i residui sono ortogonali alle colonne.Interpolazione a tratti, spline cubiche e minimi quadrati →) o l'approssimazione ai minimi quadrati.

Codice

import numpy as np

def diffdiv(x, y):                       # coefficienti di Newton f[x0], f[x0,x1], ...
    x = np.asarray(x, float); c = np.array(y, float)
    n = len(x)
    for j in range(1, n):
        c[j:] = (c[j:] - c[j-1:-1]) / (x[j:] - x[:n-j])
    return c

def newton_val(x, c, t):                 # schema di Horner generalizzato
    p = c[-1]
    for k in range(len(c) - 2, -1, -1):
        p = p * (t - x[k]) + c[k]
    return p

x = [0, 0.5, 1]; y = np.exp(x)
c = diffdiv(x, y)
print(c)                                 # differenze divise
print(newton_val(x, c, 0.25), np.exp(0.25))

x = [1, 2, 4]; y = [1, 0.5, 0.25]        # 1/x
c = diffdiv(x, y)
print(c, newton_val(x, c, 3))
print(np.polyfit(x, y, 2))               # stessi nodi: coefficienti in x^2, x, 1

Output:

[1.         1.29744254 0.84167857]
1.2717557244677156 1.2840254166877414
[ 1.    -0.5    0.125] 0.25
[ 0.125 -0.875  1.75 ]

In MATLAB: p = polyfit(x, y, n) con $n$ pari al numero di nodi meno uno restituisce il polinomio di interpolazione (coefficienti in potenze di $x$), polyval(p, t) lo valuta. Per molti nodi polyfit avvisa del cattivo condizionamento (matrice di Vandermonde).

Errori tipici

• Confondere grado e numero di nodi: con $n+1$ nodi il grado è $\le n$, e l'errore contiene $f^{(n+1)}$.
• Dimenticare che $\xi$ dipende da $x$ e che non si conosce: nell'esercizio si maggiora $|f^{(n+1)}|$ sull'intervallo.
• Usare molti nodi equispaziati sperando in un errore minore (Runge).
• Sbagliare gli indici delle differenze divise: il denominatore è $x_{i+k} - x_i$ (nodi estremi), non la differenza di due nodi adiacenti.

Domande d'esame

1. Perché il polinomio di interpolazione di grado $\le n$ su $n+1$ nodi distinti è unico? Risposta: matrice di Vandermonde con $\det = \prod_{i<j}(x_j - x_i) \ne 0$; oppure la differenza di due soluzioni ha $n+1$ zeri e grado $\le n$.
2. Enunciare la formula dell'errore e indicare da cosa dipende. Risposta: $f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod(x - x_i)$; dipende da $f$ (derivata $(n+1)$-esima) e dai nodi.
3. Calcolare la differenza divisa del secondo ordine di $f(x)=x^2$ sui nodi $0,1,3$. Risposta: $f[0,1] = 1$, $f[1,3] = 4$, $f[0,1,3] = \frac{4 - 1}{3} = 1$ (il coefficiente di $x^2$).
4. Cos'è il fenomeno di Runge e come si limita? Risposta: crescita dell'errore con $n$ per nodi equispaziati; si usano nodi di Chebyshev o polinomi a tratti.