Calcolo numerico/Altre scritture dei polinomi interpolanti

Scrittura dei polinomi interpolanti in funzione dei polinomi nodali

Cerchiamo una scrittura equivalente dei polinomi interpolanti:

$\sum _{i=0}^n{y}_i{l}_i=$

$=\sum _{i=0}^n{y}_i\frac{\prod _{i\ne j,j=0}^n(x-x_j)}{\prod _{i\ne j,j=0}^n(x_i-x_j)}$

e tenendo conto dell'espressione di $w_{n+1}(x)$:

$=\sum _{i=0}^n{y}_i\frac{w_{n+1}(x)}{x-x_i}\ast \frac{1}{\prod _{j\ne i,j=1}^n(x_i-x_j)}$

Osserviamo che:

$w_{n+1}(x{)}^{\prime }=(\prod _{j=0}^n(x-x_j){)}^{\prime }=\sum _{i=0}^n\prod _{j\ne i,j=0}^n(x-x_j)$

e valutando in $x_k$:

$w_{n+1}(x_k{)}^{\prime }=\prod _{j\ne k,j=0}^n(x_k-x_j)$

quindi, tornando all'espressione del polinomio interpolante:

$\sum _i{y}_i{l}_i(x)=$

$=W_{n+1}(x)\ast \sum _{i=0}^n{y}_i\frac{1}{x-x_i}\ast \frac{1}{w_{n^{\prime }+1}(x_i)}$

Considerazioni sui costi computazionali

Poniamo

$p_n(x)=w_{n+1}(x)\ast \sum _{i=0}^n\frac{z_i}{x-x_i}$

dove $z_i=\frac{y_i}{\prod _{i\ne j,j=0}^n(x_i-x_j)}$.

Ci si chiede qual è il costo computazionale della valutazione di $p_n(x)$ in un $x$ assegnato. Sono necessarie:

1. $n$ sottrazioni per calcolare $x-x_i$ (al denominatore);
2. $n^2/2$ sottrazioni per calcolare gli $x_i-x_j$ nella produttoria, e questo equivale a calcolare le entrate della matrice $L=(l_{ij})$ con 0 sulla diagonale principale, tale che $l_{ij}=-l_{ji}$;
3. $n^2$ operazioni per calcolare gli $z_i$, in particolare per calcolare ogni $z_i$, una volta che gli $x_i-x_j$ sono noti, sono necessari solo $n$ prodotti. Siccome gli $z_i$ sono $n$ abbiamo in totale $n^2$ operazioni;
4. $n$ addizioni (defivano dalla sommatoria) e $n$ moltiplicazioni per calcolare $\frac{z_i}{x-x_i}$, divisione tra due numeri già noti;
5. $n$ operazioni per la valutazione del polinomio nodale di cui già conosco i fattori.

$tot=n^2/2+n^2$

più una quantità trascurabile dell'ordine di $n$ operazioni.

Ci si chiede quanto costa calcolare il polinomio interpolante in un nuovo punto $\overline{x}$, dopo aver effettuato il calcolo descritto sopra. Le operazioni riferite ai punti 2 e 3 (calcolo completo degli $z_i$) sono indipendenti dal punto in cui si valuta il polinomio, e non vanno ripetute a ogni passaggio. Si hanno quindi solo $2n$ operazioni additive e $2n$ moltiplicative. Rispetto alla risoluzione del sistema delle $n+1$ equazioni nelle $n+1$ incognite (come nella dimostrazione 1), dove la fattorizzazione LU richiede $n^3$ operazioni, si risparmiano operazioni.

Supponiamo di aggiungere un nuovo punto $(x_{n+1},y_{n+1})$ ai dati: si vuole calcolare un nuovo polinomio interpolante di grado $n+1$.

Gli $z_i$ sono già stati calcolati, si ha che il nuovo polinomio è:

$p_{n+1}(x)=w_{n+2}(x)\ast [\sum _{i=0}^n(\frac{z_i^{(n+1)}}{x-x_i})+\frac{z_{n+1}^{(n+1)}}{x-x_{n+1}}]$

$z_i^{(n+1)}=\frac{y_i}{\prod _{j=0}^n(x_i-x_j)}\ast \frac{1}{x_i-x_{n+1}}$, quindi

$z_i^{(n+1)}=\frac{z_i^{(n)}}{x_i-x_{n+1}}$

Inoltre

$w_{n+2}(x)=w_{n+1}(x)\ast (x-x_{n+1})$

e si ricava che

$p_{n+1}(x)=w_{n+1}(x)\ast (x-x_{n+1})\ast [\sum _{i=0}^n(\frac{z_i^{(n)}}{(x-x_i)(x_i-x_{n+1})})+\frac{z_{n+1}^{(n+1)}}{x-x_{n+1}}]$

quindi le operazioni totali in più, avendo già calcolato $p_n(x)$, sono una sottrazione per $x-x_{n+1}$, $n$ sottrazioni per $x_i-x_{n+1}$, $n$ moltiplicazioni per $z_i^{(n+1)}$, $n$ addizioni per la sommatoria, $n$ moltiplicazioni per $z_i^n/(x-x_i)$, 2 operazioni per l'ultimo termine.

In totale abbiamo $2n$ addizioni e $2n$ moltiplicazioni.

Forma di Newton

Il polinomio interpolante è unico ma ci sono diverse riscritture. La forma di Newton permette di calcolare i polinomi interpolanti aumentando le operazioni additive e diminuendo quelle moltiplicative.

Il polinomio interpolante

$p_n(x)=\sum _{k=0}^n{a}_k\ast {\varphi }_k(x)$

viene riscritto in funzione della base dei polinomi nodali, ponendo

${\varphi }_k=\prod _{j=1}^k(x-x_j)$

Gli $a_k$ vengono posti uguali alla differenza divisa di ordine $k$.

Vogliamo dimostrare che

$p_n(x)=f[x_0]+f[x_0,x_1]\ast (x-x_0)+f[x_0,x_1,x_2]\ast (x-x_0)\ast (x-x_1)+\dots +f[x_0,x_1,\dots ,x_n)\ast (x-x_0)\ast (x-x_1)\ast (x-x_{n-1})$

è un polinomio interpolante, dimostriamo prima il lemma di Neville.

Lemma 5.1 (formula di Neville)

Dato $I=\{i_1,\dots ,i_k\}$ un insieme di indici, e definendo $p(i_0,\dots ,i_k)$ il polinomio valutato nei nodi $(x_{i_j},y_{i_j})$ per $j=0,\dots ,n$, si ha che:

$p_k(x)=\frac{(x-x_{i_0})\ast p_{k-1}(x)(i_1,\dots ,i_k)-(x-x_{i_k})\ast P_{k-1}(x)(i_0,\dots ,i_{k-1})}{x_{i_k}-x_{i_0}}$

è un polinomio interpolante.

Dim. Il lemma dà una rappresentazione ricorsiva del polinomio interpolante.

Per calcolare l'espressione di $p(i_0)$ si tiene conto che il termine $(x-x_{i_0})\ast p_{k-1}(i_1,\dots ,i_k)$ è nullo, quindi rimane

$p(i_0)=\frac{-(x-x_{i_k})\ast P_{k-1}(x_{i_0})(i_0,\dots ,i_{k-1})}{x_{i_k}-x_{i_0}}=p_{k-1}(x_{i_0})(i_0,\dots ,i_{k-1})=y_{i_0}$

Valutando invece $p(i_k)$ si annulla il secondo termine, e rimane:

$p(i_k)=\frac{(x_{i_k}-x_{i_0})\ast p_{k-1}(i_1,\dots ,i_k)(x_{i_k})}{x_{i_k}-x_{i_0}}=$

$p_{k-1}(i_1,\dots ,i_k)(x_{i_k})=y_{i_k}$

per un generico $j=1,2,\dots ,k-1$, si ha:

$p(i_j)=p(x_{i_j})=\frac{(x_{i_j}-x_{i_0})\ast p_{k-1}(i_1,\dots ,i_k)(x_{i_j})-(x_{i_j}-x_{i_k})\ast p_{k-1}(i_0,\dots ,i_{k-1})(x_{i_j})}{x_{i_k}-x_{i_0}}$

$=\frac{(x_{i_j}-x_{i_0})\ast y_{i_j}-(x_{i_j}-x_{i_k})\ast y_{i_j}}{x_{i_k}-x_{i_0}}$

$=y_{i_j}$

Abbiamo quindi verificato tutte le condizioni di interpolazione.

Teorema 5.3 (dimostrazione della forma di Newton)

$p_n(x)=f[x_0]+f[x_0,x_1]\ast (x-x_0)+f[x_0,x_1,x_2]\ast (x-x_0)\ast (x-x_1)+\dots +f[x_0,x_1,\dots ,x_n]\ast (x-x_0)\ast \dots \ast (x-x_{n-1})$

è un polinomio interpolante.

Dimostrazione

Dimostriamo l'asserto per induzione.

1. Se $n=0$, otteniamo che $p_0(x)=f[x_0]=y_0$, per definizione.
2. Per ipotesi induttiva, assumo che sia vero che il polinomio di Newton è il polinomio interpolante per il grado $n-1$ rispetto alla base nodale, e lo dimostro per $n$.
3. Supponiamo che $p_n(x)$ sia il polinomio interpolante di grado $n$ rispetto alla base delle ${\varphi }_k(x)$ date dai polinomi nodali di indice $k$, allora possiamo scriverlo come:

$p_n(x)=p_{n-1}(x)+a_n\ast (x-x_0)(x-x_1)\ast \dots \ast (x-x_n)$

vogliamo dimostrare che $a_n=f[x_0,\dots ,x_n]$. Consideriamo la formula di interpolazione di Neville. Allora

$p_n(x)=\frac{(x-x_0)\ast p_{n-1}(1,\dots ,n)-(x-x_n)\ast p_{n-1}(x)(0,\dots ,n-1)}{x_n-x_0}$

Dall'ipotesi di induzione segue che: $p_{n-1}(1,\dots ,n)$ è il polinomio interpolante di grado $n-1$ sui nodi $x_1,\dots ,x_n$, e per ipotesi di induzione il coefficiente $K_1$ del termine di grado massimo è $f[x_1,\dots ,x_n]$. $p_{n-1}(0,\dots ,n-1)$ è il polinomio interpolante rispetto ai nodi $(x_0,\dots ,x_{n-1})$ e il coefficiente $K_2$ del suo termine di grado massimo è $f[x_0,\dots ,x_{n-1}]$. Quindi, nel polinomio $p_n$, in base alla formula di Neville e a quanto osservato, il coefficiente del termine di grado massimo è:

$K=\frac{1\ast K_1-1\ast K_2}{x_k-x_0}$

$=\frac{1\ast f[x_1,\dots ,x_n]-1\ast f[x_0,\dots ,x_{n-1}]}{x_n-x_0}=f[x_0,\dots ,x_n]\text{per definizione}$

Osservazione 5.1

La differenza divisa non dipende dall'ordine dei punti. Infatti, confrontando la scrittura del polinomio interpolante in forma di Newton e di Lagrange, per il coefficiente $a_n$ si ha:

$a_n=f[x_0,\dots ,x_n]=\sum _{i=0}^n\frac{f(x_i)}{p_{n^{\prime }+1}(x_i)}$

siccome nella sommatoria non c'è dipendenza dall'ordine dei punti, anche la differenza divisa non ne dipende in aritmetica esatta.

Algoritmo alla Neville

Consideriamo i due vettori colonna $x=(x_0,x_1,\dots ,x_n)$ e $y=(f_0,f_1,\dots ,f_n)$ dove $f_i=f(x_i)$. L'algoritmo alla Neville porta alla costruzione di una matrice triangolare inferiore $A$, dove le entrate sono determinate nel seguente modo:

$\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{(}\mathrm{i}\mathrm{=}\mathrm{1}\mathrm{:}\mathrm{n}\mathrm{-}\mathrm{1}\mathrm{)}$

$\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{(}\mathrm{j}\mathrm{=}\mathrm{n}\mathrm{:}\mathrm{i}\mathrm{)}$

$Y(j)=(y(j)-y(j-1))/(x(j)-x(j-i))\text{differenza divisa}$

$\mathrm{e}\mathrm{n}\mathrm{d}$

$\mathrm{e}\mathrm{n}\mathrm{d}$

Ad esempio, presi 6 punti $(x_0,x_1,x_2,x_3,x_4,x_5$ la matrice che risulta dall'algoritmo di Neville è:

$\begin{array}{ccccccc}{x}_0& f(x_0)& & & & & \\ x_1& f(x_1)& f[x_1,x_0]& 0& 0& 0& 0\\ x_2& f(x_2)& f[x_2,x_1]& f[x_2,x_1,x_0]& 0& 0& 0\\ x_3& f(x_3)& f[x_3,x_2]& f[x_3,x_2,x_1]& f[x_3,\dots ,x_0]& 0& 0\\ x_4& f(x_4)& f[x_4,x_3]& f[x_4,x_3,x_2]& f[x_4,\dots ,x_1]& f[x_4,\dots ,x_0]& 0\\ x_5& f(x_5)& f[x_5,x_4]& f[x_5,x_4,x_3]& f[x_5,\dots ,x_2]& f[x_5,\dots ,x_1]& f[x_5,x_0]\end{array}$

Il ciclo viene fatto a ritroso perché, calcolando la triangolare inferiore con la colonna partendo dal basso, si possono sovrascrivere direttamente i vettori. In questo algoritmo, l'entrata $A_{ij}$ diventa la differenza divisa tra gli elementi $A_{i,j-1}$ e $A_{i-1,j-1}$, per ogni $j$.

Numero di operazioni effettuate: per ogni elemento faccio una divisione e due sottrazioni, e la triangolare inferiore ha $n^2/2$ elementi, quindi vengono effettuate $2\ast n^2/2=n^2$ operazioni additive e $n^2/2$ operazioni moltiplicative.

Algoritmo di Hornern

Per la valutazione del polinomio di Newton in un punto $x$ si usa l'algoritmo di Horner: ad esempio, nel caso di un polinomio di grado 3, il polinomio di Newton si rappresenta come:

$p_n=a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+a_3(x-x_0)(x-x_1)(x-x_2)$

Raccogliamo $x-x_0$ tra gli ultimi tre addendi:

$p_n=a_0+(x-x_0)\ast [a_1+a_2(x-x_1)+a_3(x-x_1)(x-x_2)]$

Poi raccogliamo $x-x_1$ tra gli ultimi due addendi:

$p_n=a_0+(x-x_0)\ast [a_1+(x-x_1)\ast [a_2+(x-x_2)\ast a_3]]$

e per eseguire le operazioni partiamo dall'ultima parentesi. Tuttavia le operazioni necessarie per la valutazione del polinomio in un punto sono di più rispetto a quelle del polinomio di Lagrange.