Calcolo numerico/Altre scritture dei polinomi interpolanti

Scrittura dei polinomi interpolanti in funzione dei polinomi nodali

Definizione

Definiamo polinomio nodale

$W_{n+1}(x)=\prod _{i=0}^{n}(x-x_{i}).$

Cerchiamo una scrittura equivalente dei polinomi interpolanti:

$\sum _{i=0}^{n}y_{i}l_{i}=$

$=\sum _{i=0}^{n}y_{i}{\frac {\prod _{i\neq j,j=0}^{n}(x-x_{j})}{\prod _{i\neq 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}}}*{\frac {1}{\prod _{j\neq i,j=1}^{n}(x_{i}-x_{j})}}$

Osserviamo che:

$w_{n+1}(x)'=(\prod _{j=0}^{n}(x-x_{j}))'=\sum _{i=0}^{n}\prod _{j\neq i,j=0}^{n}(x-x_{j})$

e valutando in $x_{k}$:

$w_{n+1}(x_{k})'=\prod _{j\neq 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)*\sum _{i=0}^{n}y_{i}{\frac {1}{x-x_{i}}}*{\frac {1}{w_{n+1}'(x_{i})}}$

Considerazioni sui costi computazionali

Poniamo

$p_{n}(x)=w_{n+1}(x)*\sum _{i=0}^{n}{\frac {z_{i}}{x-x_{i}}}$

dove $z_{i}={\frac {y_{i}}{\prod _{i\neq 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 ${\bar {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)*[\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})}}*{\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)*(x-x_{n+1})$

e si ricava che

$p_{n+1}(x)=w_{n+1}(x)*(x-x_{n+1})*[\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}*\phi _{k}(x)$

viene riscritto in funzione della base dei polinomi nodali, ponendo

$\phi _{k}=\prod _{j=1}^{k}(x-x_{j})$

Gli $a_{k}$ vengono posti uguali alla differenza divisa di ordine $k$.

Definizione

Si indica con $f[x_{0},\dots ,x_{k}]$ la differenza divisa di ordine $k$ ($k+1$ nodi). La differenza divisa si definisce induttivamente, ponendo:

$f[x_{0}]=f(x_{0}){\hbox{per definizione}}$

$f[x_{0},x_{1}]={\frac {f(x_{1})-f(x_{0})}{x_{1}-x_{0}}}\quad {\hbox{rapporto incrementale}}$

$f[x_{0},\dots ,x_{k}]={\frac {f[x_{1},\dots ,x_{k}]-f[x_{0},\dots ,x_{k-1}]}{x_{k}-x_{0}}}$

Vogliamo dimostrare che

$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},x_{1},\dots ,x_{n})*(x-x_{0})*(x-x_{1})*(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}})*p_{k-1}(x)(i_{1},\dots ,i_{k})-(x-x_{i_{k}})*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}})*p_{k-1}(i_{1},\dots ,i_{k})$ è nullo, quindi rimane

$p(i_{0})={\frac {-(x-x_{i_{k}})*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}})*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}})*p_{k-1}(i_{1},\dots ,i_{k})(x_{i_{j}})-(x_{i_{j}}-x_{i_{k}})*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}})*y_{i_{j}}-(x_{i_{j}}-x_{i_{k}})*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}]*(x-x_{0})+f[x_{0},x_{1},x_{2}]*(x-x_{0})*(x-x_{1})+\dots +f[x_{0},x_{1},\dots ,x_{n}]*(x-x_{0})*\dots *(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 $\phi _{k}(x)$ date dai polinomi nodali di indice $k$, allora possiamo scriverlo come:

$p_{n}(x)=p_{n-1}(x)+a_{n}*(x-x_{0})(x-x_{1})*\dots *(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})*p_{n-1}(1,\dots ,n)-(x-x_{n})*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*K_{1}-1*K_{2}}{x_{k}-x_{0}}}$

$={\frac {1*f[x_{1},\dots ,x_{n}]-1*f[x_{0},\dots ,x_{n-1}]}{x_{n}-x_{0}}}=f[x_{0},\dots ,x_{n}]\quad {\hbox{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+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:

${\rm {{for}(i=1:n-1)}}$

${\rm {{for}(j=n:i)}}$

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

${\rm {end}}$

${\rm {end}}$

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*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})*[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})*[a_{1}+(x-x_{1})*[a_{2}+(x-x_{2})*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.