Interpolazione polinomiale nella base di Lagrange

Un'alternativa alla base di Newton che permette di calcolare il polinomio interpolante senza cadere in un problema mal condizionato consiste nel calcolare i coefficienti del polinomio interpolante nella base di Lagrange.

Gli elementi della base di Lagrange sono così definiti:

$$\mathbf{L}_{in}(x)= \frac {\prod_{\substack{j=0 \\ j \neq i}}^{n}{(x-x_j)} } {\prod_{\substack{j=0 \\ j \neq i}}^{n}{(x_i-x_j)} }$$ (3.3)

• Questo oggetto è ben definito se per $i \neq j$ si ha $x_i \neq x_j$ e questa è la condizione fondamentale che si richiede alle ascisse di interpolazione.

• $L_{in}(x) \in \Pi_n$.

• I polinomi $L_{in}(x)$ con $i=0,\ldots, n$ sono esattamente $n+1$ e si può dimostrare che costituiscono una base per $\Pi_n$.

• Il denominatore è una costante reale

Osserviamo anche che

$$\left \{ \begin{array}{l} L_{in}(x_j)=0 \qquad j \ne i \\ L_{in}(x_i)=1 \end{array} \right.$$

Vogliamo ora trovare i coefficienti del polinomio interpolante nella base di Lagrange

$$p_n(x)=\sum_{i=0}^n \beta_i L_{in}(x)$$

con le condizioni $p_n(x_k)=f_k \quad k=0,\ldots,n$
ma

$$p_n(x_k)=\sum_{i=0}^n \beta_i L_{in}(x_k)=\beta_k L_{kn}(x_k)=\beta_k=f_k$$

quindi i coefficienti nella base di Lagrange coincidono con gli $f_i$.
In conclusione:

$$p_n(x)=\sum_{i=0}^n f_i L_{in}(x)$$

è il polinomio interpolante.

Questo risultato è importante dal punto di vista implementativo: i coefficienti del polinomio interpolante sono noti, essendo gli $f_i$ i dati del problema, l'algoritmo quindi deve solo calcolare il valore dei polinomi di base nelle ascisse di tabulazione.

In particolare la procedura LagrangeInterpola costruisce una matrice $L$ che ha $n$ righe (dove $n$ è il grado del polinomio interpolante) e $steps$ colonne dove con $steps$ indichiamo il numero di punti di tabulazione; questa matrice contiene nella posizione $(i,j)$ il valore di $\mathbf{L}_{in}$ nell'ascissa di tabulazione di posizione $j$. In simboli se indichiamo con $xtab$ il vettore che contiene i punti di tabulazione e con $l_{ij}$ l'elemento di posto $(i,j)$ della matrice $L$ si ha $L \in \mathbf{R}^{(n \; x \; steps)}$ e

$$l_{ij}=\mathbf{L_{in}}(xtab_j)$$

Una volta costruita questa matrice, i valori assunti dal polinomio interpolante nelle ascisse di tabulazione si ottengono facilmente dal prodotto del vettore che contiene gli $f_i$ e la matrice $L$. Infatti

$$p_n(xtab_j)=\sum_{i=0}^n{f_i L_{in}(xtab_j)}$$

e questo è proprio il valore assunto dalla posizione $j$ del vettore $fL$ se $f$ è il vettore degli $f_i$. Da questo prodotto si ottiene un vettore di lunghezza $steps$ se $f$ è pensato come vettore colonna infatti le dimensioni di $f$ e di $L$ sono in questo caso $(1 \; x \; n)$ e $(n \; x \; steps)$ quindi il risultato del prodotto avrà dimensioni $(1 \; x \; steps)$.

La procedura LagrangeInterpola riceve una funzione da interpolare, il vettore delle ascisse di interpolazione e il numero di punti di tabulazione da calcolare. La matrice $L$ è ottenuta richiamando la procedura ausilaria Lag che calcola in parallelo il valore del nominatore dei vari $L_{in}$ nelle ascisse di tabulazione e poi dividendo tutti i suoi elementi per il valore del denominatore che dipende solo dalle ascisse si interpolazione e non dal punto in cui deve essere calcolato $L_{in}$ (questi calcoli sono svolti dalla procedura lag_denom).

I valori del polinomio nei punti di tabulazione sono ottenuti dal prodotto del vettore $ordinate$ (che contiene gli $f_i$) per la matrice $L$.

Dato che il polinomio di grado $n$ interpolante una funzione è unico, anche se i calcoli effettuati sono diversi, la procedura LagrangeInterpola trova gli stessi valori della procedura Interpola.

Valgono perciò le considerazioni esposte nel paragrafo 3.4 e quindi è conveniente passare anche a questa procedura un vettore $x$ contenente le ascisse di Cebyshev.