Calcolo numerico/Interpolazione polinomiale a tratti (spline)

Definizione di spline

Definizione

Consideriamo una funzione $s(x)\colon [a,b]\to \mathbb {R}$, diciamo che $s(x)$ è una spline di grado $k$ se, dato un intervallo $[a,b]$ su cui consideriamo il partizionamento

$\{a=x_{0}<x_{1}<\dots <x_{n}=b\}$

$s$ è definita come

$s(x)={\begin{cases}s_{0}(x)&\iff x\in [x_{0},x_{1}]\\s_{1}(x)&\iff x\in [x_{1},x_{2}]\\s_{n-1}(x)&\iff x\in [x_{n-1},x_{n}]\end{cases}}$

dove le $s_{i}$ sono polinomi di grado $k$ per $i=0,\dots ,n-1$, localmente $C^{\infty }$, e tali che globalmente $s_{x}\in C_{k-1}$.

Nel caso $k=1$, consideriamo l'intervallo $[a,b]$ e una partizione, in ciascun sottointervallo si ha un polinomio di grado 1, e le $s_{i}$ sono localmente polinomi di grado 1 e $s$ è di classe $C^{0}$ sull'intervallo; si ha l'interpolazione lineare a tratti.

Nel caso $k=3$ si ha l'interpolazione spline cubica, si considera una partizione di $[a,b]$, in ciascun sottointervallo definisco un polinomio di grado 3. $s$ è di classe $C^{2}$ sull'intervallo.

Sistema per la determinazione della spline

$s_{i}(x)$ sono polinomi di grado 3, con $i=0,\dots ,n-1$. Abbiamo $4n$ coefficienti da determinare.

Dobbiamo imporre le seguenti condizioni:

• $n+1$ condizioni di interpolazione, perché la spline deve valere $y_{i}$ in $x_{i}$;
• $n-1$ condizioni di raccordo con continuità, cioè $s_{i-1}(x_{i})=s_{i}(x_{i})$ per $i=1,\dots n-1$;
• $2n-2$ condizioni sulle derivate, cioè $s_{i-1}'(x_{i})=s_{i}'(x_{i})$ e $s_{i-1}''(x_{i})=s_{i}''(x_{i})$ perché $s$ dev'essere di classe $C^{2}$.

Quindi mancano due condizioni per avere un sistema determinato di $4n$ equazioni in $4n$ incognite.

In base alle due condizioni aggiuntive che si impongono si hanno diversi tipi di spline:

Spline naturale

imponiamo che $s''(a)=s''(b)=0$.

Spline completa e vincolata

$s_{0}'(a)=f'(a)$ e $s_{n-1}'(b)=f'(b)$.

Spline periodica

supponiamo $T=b-a$, allora $s_{0}'(a)=s_{n-1}'(b)$ e analogamente $s_{0}''(a)=s_{n-1}''(b)$.

Abbiamo un sistema di $4n$ equazioni in $4n$ incognite, che però ha grandi dimensioni ed è spesso mal condizionato.

Metodo dei momenti

Definizione della funzione: sia $M_{i}$ la derivata seconda di $s$ in $x_{i}$, per $i=1,\dots ,n$ (incognite da calcolare), dove $M_{i}$ sono detti momenti. Siccome $s_{i}$ ha grado 3, la sua derivata è un polinomio di grado 1, inoltre, imponendo le condizioni:

${\begin{cases}s_{i}''(x_{i})=M_{i}\\s_{i+1}''(x_{i+1})=M_{i+1}\end{cases}}$

si ottiene

$s_{i}''(x)=M_{i+1}(x-x_{i})/h_{i}-M_{i}(x-x_{i+1})/h_{i}$

dove $h_{i}=x_{i+1}-x_{i}$.

Integrando due volte otteniamo:

$s_{i}'(x)={\frac {M_{i+1}(x-x_{i})^{2}}{2h_{i}}}-M_{i}{\frac {(x-x_{i+1})^{2}}{2H_{i}}}+\alpha _{i}$

Integrando nuovamente

$s_{i}(x)=M_{i+1}{\frac {(x-x_{i})^{3}}{6H_{i}}}-M_{i}{\frac {(x-x_{i+1})^{3}}{6H_{i}}}+\alpha *(x-x_{i})+\beta _{i}\quad {\hbox{equazione }}\ast$

Imposizione delle condizioni di interpolazione: per determinare i nuovi parametri $\alpha _{i},\beta _{i}$, imponiamo le condizioni di interpolazione e otteniamo:

$s_{i}(x_{i})=-M_{i}{\frac {(x_{i}-x_{i-1})^{3}}{6H_{i}}}+\beta _{i}=f(x_{i})$

Inoltre

$s_{i}(x_{i+1})=M_{i+1}*(x_{i+1}-x_{i})^{3}/(6H_{i})+\alpha _{i}(x_{i+1}-x_{i})+\beta _{i}=f(x_{i+1})$

dove l'ultimo passaggio vale per la condizione di raccordo con continuità $s_{i}(x_{i+1}=s_{i+1}(x_{i+1})$. Quindi, per $i=0,\dots ,n-1$ otteniamo:

${\begin{cases}M_{i}{\frac {H_{i}^{2}}{6}}+\beta _{i}=f(x_{i})\\M_{i+1}{\frac {H_{i}^{2}}{6}}+\alpha _{i}H_{i}+\beta _{i}=f(x_{i+1})\end{cases}}$

Risolviamo il sistema:

${\begin{cases}M_{i}{\frac {H_{i}^{2}}{6}}-f(x_{i})=-\beta _{i}\\M_{i+1}{\frac {H_{i}^{2}}{6}}+\alpha _{i}H_{i}-f(x_{i+1})=-\beta _{i}\end{cases}}$

ed eguagliando le due equazioni otteniamo:

$M_{i}{\frac {H_{i}^{2}}{6}}-f(x_{i})=M_{i+1}{\frac {H_{i}^{2}}{6}}+\alpha _{i}H_{i}-f(x_{i+1})$

$f(x_{i+1})-f(x_{i})+(M_{i}-M_{i+1}){\frac {H_{i}^{2}}{6}}=\alpha _{i}H_{i}$

$\alpha _{i}={\frac {f(x_{i+1})-f(x_{i})}{H_{i}}}+(M_{i}-M_{i+1}){\frac {H_{i}}{6}}$

invece ricaviamo $\beta _{i}$ dalla prima equazione:

$\beta _{i}=-M_{i}{\frac {H_{i}^{2}}{6}}+f(x_{i})$

e come soluzione del sistema otteniamo

${\begin{cases}\beta _{i}=f(x_{i})-M_{i}H_{i}^{2}/6\\\alpha _{i}={\frac {f(x_{i+1})-f(x_{i})}{H_{i}}}-H_{i}/6(M_{i+1}-M_{i})\end{cases}}$

Imposizione delle condizioni di raccordo con continuità delle derivate:

$s_{i-1}'(x)={\frac {M_{i}(x-x_{i-1})^{2}}{2h_{i-1}}}-M_{i-1}{\frac {(x-x_{i})^{2}}{2H_{i-1}}}+\alpha _{i}$

$={\frac {M_{i}(x-x_{i-1})^{2}}{2h_{i-1}}}-M_{i-1}{\frac {(x-x_{i})^{2}}{2H_{i-1}}}+{\frac {f(x_{i})-f(x_{i-1})}{H_{i-1}}}-H_{i-1}/6(M_{i}-M_{i-1})$

Quindi, valutando in $x_{i}$ si ottiene:

$s_{i-1}'(x_{i})={\frac {M_{i}h_{i-1}}{2}}+{\frac {f(x_{i})-f(x_{i-1})}{H_{i-1}}}-H_{i-1}/6(M_{i}-M_{i-1})$

analogamente

$s_{i}'(x_{i})={\frac {M_{i+1}h_{i}}{2}}+{\frac {f(x_{i+1})-f(x_{i})}{H_{i}}}-H_{i}/6(M_{i+1}-M_{i})$

ed eguagliando le due espressioni:

$M_{i-1}*H_{i-1}/6+M_{i}(H_{i-1}+H_{i})/3+M_{i+1}H_{i}/6={\frac {f(x_{i+1})-f(x_{i})}{H_{i}}}-{\frac {f(x_{i})-f(x_{i-1})}{H_{i-1}}}$

e inserendo l'espressione della differenza divisa:

$M_{i-1}*H_{i-1}/6+M_{i}(H_{i-1}+H_{i})/3+M_{i+1}H_{i}/6=f[x_{i+1},x_{i},x_{i-1}]*(H_{i}-H_{i-1})$

e dividendo ambo i membri per $H_{i}-H_{i-1}$ otteniamo:

$M_{i-1}*{\frac {H_{i-1}}{6(h_{i}-h_{i-1})}}+M_{i}/3+M_{i+1}{\frac {H_{i}}{6(h_{i}-h_{i-1})}}=f[x_{i+1},x_{i},x_{i-1}]$

e moltiplicando per 6

$M_{i-1}*{\frac {H_{i-1}}{H_{i}-H_{i-1}}}+2M_{i}+M_{i+1}{\frac {H_{i}}{H_{i}-H_{i-1}}}=6f[x_{i+1},x_{i},x_{i-1}]$

Chiamiamo

$\gamma _{i}={\frac {H_{i-1}}{H_{i}-H_{i-1}}}$

e simmetricamente

$\delta _{i}={\frac {H_{i}}{H_{i}-H_{i-1}}}$

e sono due quantità positive, che hanno la proprietà

$\gamma _{i}+\delta _{i}=1,\,\forall i$

per qualsiasi tipo di nodi.

Abbiamo $n-1$ equazioni in $n+1$ incognite, date da:

$M_{i-1}*\gamma _{i}+2M_{i}+M_{i+1}\delta _{i}=6f[x_{i+1},x_{i},x_{i-1}]$

Imposizione delle condizioni iniziali: la situazione dipende dalle condizioni imposte.

Spline naturale: imponiamo che $M_{0}=s''(x_{0})=0$, e $M_{n}=s''(x_{n})=0$. Allora otteniamo un sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$ dove $A$ è una matrice tridiagonale della forma:

$A={\begin{array}{ccccc}2&\delta _{1}&0&0&0\\\gamma _{2}&2&\delta _{2}&0&0\\0&\gamma _{3}&2&\dots &0\\0&0&\dots &2&\delta _{i-1}\\0&0&0&\gamma _{n-1}&2\end{array}}$

${\boldsymbol {x}}={\begin{array}{c}M_{1}\\M_{2}\\\dots \\M_{n-2}M_{n-1}\\\end{array}}$

è il vettore delle incognite, mentre il termine noto è

${\boldsymbol {x}}={\begin{array}{c}f[x_{0},x_{1},x_{2}]\\f[x_{1},x_{2},x_{3}]\\f[x_{2},x_{3},x_{4}]\\\dots \\f[x_{n-2},x_{n-1},x_{n}]\end{array}}$

La matrice è diagonalmente dominante sulle righe.

Se i nodi sono equispaziati, cioè se $H_{i}=H\forall i$, la matrice è simmetrica. Per risolvere il sistema si può usare l'eliminazione gaussiana senza pivot, oppure metodi diretti che convergono per la dominanza diagonale stretta.

Spline completa: imponiamo le condizioni

$s'(x_{0})=f'(a),\quad s'(x_{n})=f'(b)$

Tenendo conto che

$s'(x)={\frac {M_{i}(x-x_{i-1})^{2}}{2h_{i-1}}}-M_{i-1}{\frac {(x-x_{i})^{2}}{2H_{i-1}}}+{\frac {f(x_{i})-f(x_{i-1})}{H_{i-1}}}-H_{i-1}/6(M_{i}-M_{i-1})$

si ha

$s'(x_{0})=-M_{0}H_{0}/2+f[x_{0},x_{1}]+H_{0}/6(M_{0}-M_{1})$

$f'(x_{0})=H_{0}f[x_{0},x_{0}]=s'(x_{0})=-M_{0}H_{0}/3-M_{1}H_{0}/6+H_{0}f[x_{0},x_{1}]$

che si riscrive come

$H_{0}/3M_{0}-H_{0}/6M_{1}=(f[x_{0},x_{1}]-f[x_{0},x_{0}])*H_{0}=f[x_{0},x_{0},x_{1}]*H_{0}$

dividendo per $H_{0}$ e moltiplicando per 6:

$2M_{0}+M_{1}=6f[x_{0},x_{0},x_{1}]$

Analogamente imponendo $f'(x_{n})=s'(x_{n})$ si ottene:

$M_{n-1}+2M_{n}=6f[x_{n-1},x_{n},x_{n}]$

In questo caso abbiamo un sistema di $n+1$ equazioni della forma $A{\boldsymbol {x}}={\boldsymbol {b}}$ con:

$A={\begin{array}{ccccc}2&1&0&0&0\\\gamma _{1}&2&\delta _{1}&0&0\\0&\gamma _{2}&2&\delta _{2}&0\\0&0&\gamma _{n-1}&2&\delta _{n-1}\\0&0&0&1&2\end{array}}$

Il vettore delle incognite è

${\begin{array}{c}M_{0}\\M_{1}\\\dots \\M_{n-1}\\M_{n}\end{array}}$

Il termine noto è:

$6*{\begin{array}{c}f[x_{0},x_{0},x_{1}]\\f[x_{0},x_{1},x_{2}]\\\dots \\f[x_{n-2},x_{n-1},x_{n}]\\f[x_{n-1}x_{n},x_{n}]\end{array}}$

Abbiamo una matrice diagonalmente dominante in senso stretto, non è simmetrica neanche nel caso dei nodi equispaziati.

Spline periodica: $s'(x_{0})=s'(x_{n})$ e analogamente $s''(x_{0})=s''(x_{n})$.

La seconda si traduce nella condizione $M_{0}=M_{n}$. Prendendo l'espressione di $s_{i}'$, si impone l'uguaglianza e si ha $f(x_{n})=f(x_{0})$. Introduciamo un nodo fittizio $x_{n+1}$ esterno all'intervallo, imponiamo che $H_{n}=H_{0}$. Imponendo tutte le condizioni si ottiene

$\delta _{n}M_{1}+\gamma _{n}M_{n-1}+2M_{n}=6f[x_{n-1},x_{n},x_{n+1}]$

$\gamma _{1}M_{n}+2M_{1}+\delta _{1}M_{2}=6f[x_{0},x_{1},x_{2}]$

Quindi il sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$ è tale che:

$A={\begin{array}{cccccc}2&\delta _{1}&0&0&0&\gamma _{1}\\\gamma _{2}&2&\delta _{2}&0&0&0\\0&\gamma _{3}&2&\delta _{3}&0&0\\0&0&0&\gamma _{n-1}&2&\delta _{n-1}\\\gamma _{n}&0&0&0&\delta _{n}&2\end{array}}$

Il vettore delle incognite è

${\boldsymbol {x}}={\begin{array}{c}M_{1}\\M_{2}\\\dots \\M_{n-1}\\M_{n}\end{array}}$

mentre il vettore del termine noto è

${\begin{array}{c}f[x_{0},x_{1},x_{2}]\\f[x_{1},x_{2},x_{3}]\\\dots \\f[x_{n-2},x_{n-1},x_{n}]\\f[x_{n-1},x_{n},x_{n+1}]\end{array}}$

Questa matrice $n\times n$ non è tridiagonale a causa di $\delta _{n}$ in posizione $A_{n1}$ e $\gamma _{1}$ in posizione $A_{1n}$. La matrice in ogni caso è diagonalmente dominante in senso stretto. Se i nodi sono equispaziati si ha la simmetria.

Conclusione

Nei tre casi, la matrice è diagonalmente dominante in senso stretto, e per il primo teorema di Gersch-Gorin è non singolare.

Teorema 5.7

Se i nodi $A=x_{0}<x_{1}<\dots <x_{n-1}<x_{n}=b$ sono distinti, allora esiste ed è unica la spline cubica che interpola $f$ funzione assegnata.

Teorema 5.8

Supponiamo che $f$ sia $C_{4}([a,b])$, consideriamo un partizionamento di $[a,b]$ in sottointervalli $H_{i}$, poniamo $H=\max _{i}H_{i}$ e chiamiamo $\beta ={\frac {H}{\min _{i}H_{i}}}$, e sia $s$ la spline cubica interpolante.

$\vert f_{r}-s_{r}\vert ^{\infty }\leq c_{r}^{4-r}\vert f^{(4)}\vert$