Splines

Ángeles Martínez Calomardo e Alvise Sommariva

Università degli Studi di Padova

27 novembre 2013

Interpolazione nodi equispaziati e problemi

E’ noto che nel caso dell’interpolazione polinomiale, dati N + 1 punti a = x_{0} < … < x_{N} = b, e i valori y_{0}, …, y_{N} ivi assunti dalla funzione y = f (x), esiste unico il polinomio p_{N} di grado N tale che

p_{N}(x_{i}) = f_{i}, i = 0, …, N.  (1)

Nel caso di nodi equispaziati

x_{k} = a + k \frac{(b − a)}{N}, k = 0, …, N;  (2)

al crescere di N, non si puó garantire che f − p_{n} tenda a 0 puntualmente (si ricordi il fenomeno di Runge!).

Problema: E’ possibile calcolare una interpolante di tipo polinomiale s_{N} per cui al crescere di N si abbia s_{N} → f puntualmente?

Ángeles Martínez Calomardo e Alvise Sommariva  Splines  2/ 23

Splines

Sia [a, b] ⊂ ℝ chiuso e limitato, e sia a = x_{0} < x_{1} < … < x_{n} = b. Una spline di grado m (o ordine m + 1) è una funzione in C^{m−1}([a, b]) che in ogni intervallo [x_{i}, x_{i+1}], con i = 0, …, n − 1, è un polinomio di grado m.

Alcuni esempi:

spline di grado 1 (lineari): una funzione in C([a, b]) che in ogni intervallo [x_{i}, x_{i+1}], con i = 0, …, n − 1, è un polinomio di grado 1. spline di grado 3 (cubiche): una funzione in C^{2}([a, b]) che in ogni intervallo [x_{i}, x_{i+1}], con i = 0, …, n − 1, è un polinomio di grado 3.

Splines di grado 1 e lineari a tratti

Come anticipato sono funzioni C([a, b]) che in ogni intervallo [x_{i}, x_{i+1}], con i = 0, …, n − 1, è un polinomio di grado 1. Di conseguenza sono funzioni continue lineari a tratti, essendo un segmento il grafico di un polinomio di grado 1 in [x_{i}, x_{i+1}].

Figura: Splines grado 1 (n = 10 intervalli di uguale ampiezza).

Ángeles Martínez Calomardo e Alvise Sommariva  Splines  4/ 23

Interpolazione con splines: alcuni fatti.

Dati n + 1 punti in [a, b], diciamo

a = x_{0} < x_{1} < … < x_{n−1} < x_{n} = b,

supponiamo di dover interpolare le coppie (x_{k}, y_{k}) con k = 0, …, n.

esiste unica la spline intp. s_{1} di grado 1 , cioè tale che s_{1}(x_{k}) = y_{k} per k = 0, …, n. non esiste unica la spline intp. s_{3} di grado 3, cioè tale che s_{3}(x_{k}) = y_{k} per k = 0, …, n. Servono 2 ulteriori condizioni per determinare univocamente la spline s_{3} interpolante. Alcune classiche richieste aggiuntive.

Spline naturale: s_{3}^{(2)}(a) = s_{3}^{(2)}(b) = 0.

Spline periodica: s_{3}^{(1)}(a) = s_{3}^{(1)}(b), s_{3}^{(2)}(a) = s_{3}^{(2)}(b).

Spline vincolata: s_{3}^{(1)}(a) = y′_{a}, s_{3}^{(1)}(b) = y′_{b} (con y′_{a}, y′_{b} assegnati).

Ángeles Martínez Calomardo e Alvise Sommariva  Splines  5/ 23

Splines cubiche di tipo not-a-knot

La spline cubica s_{3}^{NAK} con vincolo not-a-knot è definita come segue:

Interpola le coppie (x_{0}, y_{0}), …, (x_{n}, y_{n});

E’ un polinomio di grado 3 nell’intervallo [x_{0}, x_{2}]; E’ un polinomio di grado 3 nell’intervallo [x_{n−2}, x_{n}]

Si osservi che per le spline cubiche generiche s_{3} non è detto che la restrizione di s_{3} a [x_{0}, x_{1}] e [x_{1}, x_{2}] siano lo stesso polinomio, e similmente che la restrizione a [x_{n−2}, x_{n−1}] e [x_{n−1}, x_{n}] siano lo stesso polinomio. Vediamo il confronto dei gradi di libertá.

Richieste: Abbiamo n − 2 polinomi da determinare e quindi 4(n − 2) condizioni. Fornite: Interpolazione x_{0}, …, x_{n}: n + 1 condizioni. Regolarità x_{2}, …, x_{n−2}: 3(n − 3) condizioni.

Quindi il numero di gradi di libertá richiesto e fornito è lo stesso.

Ángeles Martínez Calomardo e Alvise Sommariva  Splines  6/ 23

Splines lineari: stima errore

Sia s_{1,Δ} : [a, b] → ℝ una spline di grado 1, che interpola le coppie (x_{i}, f (x_{i})) dove x_{0} = a < x_{1} < … < x_{i} < x_{i+1} < … < x_{N} = b. Δ = {x_{i}}_{i=0,...,N} indica la suddivisione [a, b] = ∪_{i=0}^{n−1}[x_{i}, x_{i+1}]. Si può mostrare che

|f (x) − s_{1,Δ}(x)| ≤ h_{i,Δ}^{2} \frac{M_{i,Δ}}{8} , x ∈ [x_{i}, x_{i+1}],  (3)

con

h_{i,Δ} = x_{i+1} − x_{i},
M_{i,Δ} := max_{x∈[x_{i},x_{i+1}]} |f^{(2)}(x)|.

Di conseguenza, se f ∈ C^{2}([a, b]), per M = ∥f^{(2)}∥_{∞} = max_{i} M_{i}

max_{x∈[a,b]}|f (x) − s_{1,Δ}(x)| = max_{i=0,...,N−1} max_{x∈[x_{i},x_{i}+1]}|f (x) − s_{1,Δ}(x)|
≤ max_{i=0,...,N−1}h_{i,Δ}^{2} \frac{M_{i,Δ}}{8} ≤ \frac{M}{8} (max_{i}h_{i,Δ})^{2}.

Splines lineari: convergenza uniforme

Sia \overline{h}_{Δ} = max_{i} h_{i,Δ}. Sia fissata la famiglia di suddivisioni \{Δ_{n}\} con la proprietà che \overline{h}_{Δ_{n}} → 0 (suddivisioni sempre più fini). Allora la convergenza della successioni di interpolanti s_{1,Δ_{n}} ad una funzione f ∈ C^{2}([a, b]) è uniforme. Infatti:

∥f − s_{1,Δ_{n}}∥ ≤ \frac{M}{8} (max_{i}h_{i,Δ_{n}})^{2} = \frac{M}{8} \overline{h}_{Δ_{n}}^{2} →^{n}0.

Splines cubiche: stima errore

Consideriamo di seguito spline cubiche s_{3,Δ_{n}} su una suddivisione Δ_{n} = {x_{i}} di [a, b]. Poniamo h_{i,Δ} = x_{i+1} − x_{i}, f^{(k)} derivata k-sima di f .

Teorema. Supponiamo f ∈ C^{4}([a, b]). Posto \overline{h}_{Δ} = max_{i} h_{i,Δ}, si consideri una sudd. Δ = \{x_{i}\} di [a, b] con K_{Δ} = max_{i} \overline{h}_{Δ}/h_{i,Δ}.

Esiste una costante c_{s} indipendente da \overline{h}_{Δ} tale che

∥f^{(s)} − s_{3,Δ}^{(s)}∥_{∞} ≤ c_{s}K_{Δ}\overline{h}_{Δ}^{(4−s)}||f^{(4)}||_{∞}, s = 0, 1, 2, 3 (4)

dove ∥f ∥_{∞} = max_{x∈[a,b]} |f (x)|.

PS. Notare le componenti dell’errore dovute alla suddivisione (cioè K_{Δ} e \overline{h}_{Δ}) e alla funzione f (cioè ∥f^{(4)}∥_{∞}).

Splines cubiche: stima errore

Corollario. Supponiamo f ∈ C^{4}([a, b]). Si consideri la suddivisione Δ = {x_{i}} di [a, b] con \overline{h}_{Δ} = h_{i,Δ} = h_{Δ} (cioè \{x_{i}\} nodi equisp.). Allora esiste una costante c_{0} indipendente da h tale che

∥f − s_{3,Δ}∥_{∞} ≤ c_{0} h_{Δ}^{4}||f^{(4)}||_{∞}  (5)

Inoltre esistono delle costanti c_{1}, c_{2}, c_{3} indipendenti da h tali che

∥f^{(1)} − s_{3,Δ}^{(1)}∥_{∞} ≤ c_{1} h_{Δ}^{3}||f^{(4)}||_{∞},  (6)
∥f^{(2)} − s_{3,Δ}^{(2)}∥_{∞} ≤ c_{2} h_{Δ}^{2}||f^{(4)}||_{∞},  (7)
∥f^{(3)} − s_{3,Δ}^{(3)}∥_{∞} ≤ c_{3} h_{Δ}^{1}||f^{(4)}||_{∞},  (8)

dove ∥f ∥_{∞} = max_{x∈[a,b]} |f (x)|.

Nel caso di spline interpolanti e vincolate si ha c_{0} = 5/384, c_{1} = 1/24, c_{2} = 3/8.

Ángeles Martínez Calomardo e Alvise Sommariva  Splines  10/ 23

Splines cubiche: convergenza uniforme

Sia \overline{h}_{Δ} = max_{i} h_{i,Δ}. Sia fissata la famiglia di suddivisioni \{Δ_{n}\} con la proprietà che K = max_{n}K_{Δ_{n}} < +∞ e \overline{h}_{Δ_{n}} → 0 (suddivisioni sempre più fini). Allora la convergenza della successioni di interpolanti s_{3,Δ_{n}} ad una funzione f ∈ C^{2}([a, b]) è uniforme. Infatti da (10):

∥f^{(s)} − s_{3,Δ_{n}}^{(s)}∥_{∞} ≤ c_{s}K\overline{h}_{Δ_{n}}^{(4−s)}||f^{(4)}||_{∞}→^{n}0.

NB: Per s = 3, significa che la successione di funzioni a gradino s_{3,Δ_{n}}^{(3)} converge a f^{(3)} uniformemente.

Splines lineari: esempio di Runge.

La funzione di Runge f (x) = 1/(1 + x^{2}) appartiene a C^{∞}([−5, 5]) e di conseguenza, nonostante la convergenza puntuale dell’interpolazione polinomiale in nodi equispaziati non sia garantita, ciò non si può dire per l’interpolazione mediante splines lineari o cubiche (sappiamo sussistere la conv. unif.!). Vediamo una stima dell’errore più precisa. Nel caso lineare, se Δ = {x_{i}} è la sudd. equisp. di ampiezza h_{Δ}, se x ∈ [x_{i}, x_{i+1}]

|f (x) − s_{1,Δ}(x)| ≤ h_{i,Δ}^{2} \frac{M_{i,Δ}}{8} = h_{Δ}^{2} \frac{M_{i,Δ}}{8}  (9)

dove

M_{i,Δ} = max_{x∈(x_{i},x_{i+1})}|f^{(2)}(x)|.

con

f^{(2)}(x) = (8 x^{2})/(x^{2} + 1)^{3} − 2/(x^{2} + 1)^{2}.

Splines lineari: esempio di Runge.

Dal grafico di |f^{(2)}| deduciamo che

M_{i} ≤ M = ∥f^{(2)}∥_{∞} = |f^{(2)}(0)| = 2

e quindi se x ∈ [x_{i}, x_{i+1}]

|f (x) − s_{1,Δ}(x)| ≤ h_{Δ}^{2} \frac{M_{i}}{8} ≤ h_{Δ}^{2} \frac{M}{8} ≤ \frac{h_{Δ}^{2}}{4} ⇒ ∥f − s_{1,Δ}∥_{∞} ≤ \frac{h_{Δ}^{2}}{4} .

Splines cubiche: esempio di Runge.

Usiamo le notazioni precedenti. Sappiamo che ∥f − s_{3,Δ}∥_{∞} ≤ c_{0} h_{Δ}^{4}||f^{(4)}||_{∞} e si può vedere che

f^{(4)}(x) = 24/(x^{2} + 1)^{3} − 288 · x^{2}/(x^{2} + 1)^{4} + 384 · x^{4}/(x^{2} + 1)^{5}

il cui massimo modulo in [−5, 5] è |f^{(4)}(0)| = 24. Quindi per splines cubiche intp. vincolate, essendo c_{0} = 5/384

∥f − s_{3,Δ}∥_{∞} ≤ c_{0} h_{Δ}^{4}||f^{(4)}||_{∞} ≤ (120/384)h_{Δ}^{4}.

Splines lineari e cubiche: confronto esempio di Runge.

Figura: Grafico stime errori intp. splines lineari e cubiche vincolate (ascissa: x = −log_{10}(h)).

>> n = 0 : 0 . 2 : 1 5 ; h=10.ˆ (−n ) ;
>> e1=(h . ˆ 2 ) / 4 ; e3 =(120/384) ∗(h . ˆ 4 ) ;
>> s e m i l o g y ( n , e1 , ’ r− ’ , n , e3 , ’ k− ’ ) ;
>> l e g e n d ( ’ l i n e a r e ’ , ’ c u b i c a ’ ) ;
Ángeles Martínez Calomardo e Alvise Sommariva  Splines  15/ 23

Splines interpolanti in Matlab/Octave: spline

Alcuni comandi per calcolare interpolazione spline in Matlab/Octave

SPLINE  Cubic s p l i n e data  interpolation .
YY = SPLINE ( X , Y , XX ) uses  cubic s p l i n e interpolation  to f i n d
     YY , the  values  of the  underlying f u n c t i o n Y at the
     points  in the  vector  XX .
The  vector X specifies  the  points  at which the  data Y is
     given . . . .
Ordinarily , the not−a−knot end conditions  are  used . However ,
     i f Y contains  two more values  than X has  entries , then
     the  first  and  last  value  in Y are  used as the  endslopes

f o r the  cubic s p l i n e .

…

Splines interpolanti in Matlab/Octave: spline

Example : This  generates a sine  curve , then  samples  the

s p l i n e over a finer mesh :
x = 0 : 1 0 ; y = s i n ( x ) ;
xx = 0 : . 2 5 : 1 0 ; yy = s p l i n e ( x , y , xx ) ;
p l o t ( x , y , ’ o ’ , xx , yy )

Example : This  illustrates  the  use of  clamped  or  complete
     s p l i n e interpolation  where end slopes  are  prescribed .
     Zero  slopes  at the  ends of an  interpolant  to the  values
     of a certain  distribution  are  enforced :

x = −4:4; y = [ 0 . 1 5 1 . 1 2 2 . 3 6 2 . 3 6 1 . 4 6 . 4 9 . 0 6 0 ] ;
cs = s p l i n e ( x , [ 0 y 0 ] ) ;
xx = l i n s p a c e ( −4 ,4 ,101) ;
p l o t ( x , y , ’ o ’ , xx , ppval ( cs , xx ) , ’− ’ ) ;

See  also  INTERP1 , PPVAL , SPLINES ( The  Spline  Toolbox ) .

Splines interpolanti in Matlab/Octave: interp1

INTERP1 1−D interpolation ( table  lookup ) .
YI = INTERP1 ( X , Y , XI ) interpolates  to f i n d YI , the  values  of
    the  underlying f u n c t i o n Y at the  points  in the  vector  XI
     . The  vector X specifies  the  points  at which the  data Y
    is  given . . . .
YI = INTERP1 ( X , Y , XI , ’ method ’ ) specifies  alternate  methods .
The  default  is  linear  interpolation . Available  methods  are :
’ n e a r e s t ’ − nearest  neighbor  interpolation
’ l i n e a r ’ − linear  interpolation
’ s p l i n e ’ − piecewise  cubic s p l i n e interpolation ( SPLINE )
’ p c h i p ’  − piecewise  cubic  Hermite  interpolation ( PCHIP )

’ c u b i c ’  − same as ’ p c h i p ’

’ v 5 c u b i c ’ − the  cubic  interpolation  from  MATLAB 5 , which
    does  not  extrapolate  and  uses ’ s p l i n e ’ i f Xmis  not
     equally  spaced . . . .

Splines interpolanti in Matlab/Octave: interp1

For  example , generate a coarse  sine  curve  and  interpolate
over a finer  abscissa :

x = 0 : 1 0 ; y = s i n ( x ) ; xi = 0 : . 2 5 : 1 0 ;
yi = i n t e r p 1 ( x , y , xi ) ; p l o t ( x , y , ’ o ’ , xi , yi )

See  also  INTERP1Q , INTERPFT , SPLINE , INTERP2 , INTERP3 ,
INTERPN .

Quindi dall’help di spline e interp1 deduciamo che: spline: calcola spline cubica intp. con vincolo not-a-knot e con qualche assegnazione anche vincolate; interp1: calcola spline lineare o spline cubica intp. con vincolo not-a-knot e con qualche assegnazione anche di altro tipo.

Splines interpolanti in Matlab/Octave: esempio di Runge

Vediamo il comportamento della spl. lineare intp. la funzione di Runge in 13 nodi equisp., ricordando che l’interpolante polinomiale p_{12} su 13 nodi, aveva ∥f − p_{12}∥ ≈ 3.66 + 00. Il file è runge lin.m.

n=13;
f=inline ( ’ 1./(1+ x . ˆ 2 ) ’ ) ; % RUNGE.
t=l i n s p a c e ( −5 ,5 ,1000) ; % NODI TEST .
x=l i n s p a c e ( −5 ,5 , n ) ; % NODI SUDDIV .
y=f e v a l ( f , x ) ; % ( x , y ) COPPIE INTP .
st=i n t e r p 1 ( x , y , t , ’ l i n e a r ’ ) ;
ft=f e v a l ( f , t ) ;
err=norm ( ft−st , inf ) ; % ERRORE | f−s3 |
f p r i n t f ( ’ \n \ t e r r : %2.2 e \n ’ , err ) ;

Otteniamo, come supposto, un migliore risultato rispetto a p_{12} (stime: ≤ ((10/12)^{2})/4 ≈ 0.17.)

>> runge_lin

err : 6 . 5 0 e−02

>>

Splines interpolanti in Matlab/Octave: esempio di Runge

Vediamo il comportamento delle spline cubiche not-a-knot interpolanti la funzione di Runge in 13 nodi equispaziati (avevamo ∥f − p_{12}∥ ≈ 3.66 + 00). Il file è runge cub.m.

n=13;
f=inline ( ’ 1./(1+ x . ˆ 2 ) ’ ) ; % RUNGE.
t=l i n s p a c e ( −5 ,5 ,1000) ; % NODI TEST .
x=l i n s p a c e ( −5 ,5 , n ) ; % NODI SUDDIV .
y=f e v a l ( f , x ) ; % ( x , y ) COPPIE INTP .
st=i n t e r p 1 ( x , y , t , ’ s p l i n e ’ ) ;
ft=f e v a l ( f , t ) ;
err=norm ( ft−st , inf ) ; % ERRORE | f−s3 |
f p r i n t f ( ’ \n \ t e r r : %2.2 e \n ’ , err ) ;

Otteniamo, come supposto, un migliore risultato rispetto a p_{12} e alle splines lineari (errore come h^{4} e non h^{2})

>> runge_lin

err : 6 . 9 1 e−03

>>

Esercizio

Un esempio interessante è l’approssimazione del cerchio con splines. Essendo l’equazione del cerchio

f (t) = ( \begin{matrix}cos(t) \\ sin(t) \end{matrix} )  (10)

per t = [0, 2 π]. Se

t_{k} = \frac{2kπ}{N}, k = 0, …, N,

l’interpolante spline puó essere facilmente calcolata via il comando Matlab interp1 applicato alla coppia di vettori t, u e t, v dove al solito t = (t_{k}), u = (u(t_{k})), v = (v(t_{k})).

Esercizio

Per testare il comportamento di tale ricostruzione, fissati dei nodi test

t_{k}^{∗} = \frac{2kπ}{100} , k = 0, …, 100,

si calcolino per M = 10 le splines s_{x}, s_{y} lineari (o cubiche not-a-knot) interpolanti le coppie

(t_{k}, cos(t_{k})), t_{k} = \frac{2kπ}{M}, k = 0, …, M,

(t_{k}, sin(t_{k})), t_{k} = \frac{2kπ}{M}, k = 0, …, M,

si plottino le coppie (s_{x}(t_{k}^{∗}), s_{y}(t_{k}^{∗})). Ripetere l’esercizio per M = 100.