Splines

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

21 novembre 2011

Interpolazione nodi equispaziati e problemi

Si è visto 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 da una funzione y = f (x), esiste uno ed un solo 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 (x) − p_{n}(x) tenda a 0 (si ricordi il fenomeno di Runge!). Sorge spontaneo chiedersi se qualora si possegga un gran numero di punti, anche equispaziati, sia possibile calcolare un’approssimante di tipo polinomiale per cui al crescere di N si abbia p_{N} → f .

Alvise Sommariva  Splines  2/ 25

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: Esempio di splines grado 1 (n = 10 intervalli di uguale ampiezza).

Alvise Sommariva  Splines  4/ 25

Interpolazione con splines: grado 1. Esistenza e unicità

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, con una spline s_{1} di grado 1 cioè determinare s_{1} tale che s_{1}(x_{k}) = y_{k} per k = 0, …, n. Poichè per due punti assegnati (x_{j}, y_{j}), (x_{j+1}, y_{j+1}) passa un unico segmento è evidente che tale funzione s_{1} è unica. Alternativamente, basta osservare che in ogni subintervallo esiste ed è unico il polinomio di grado 1 che interpola 2 punti assegnati (x_{j}, y_{j}), (x_{j+1}, y_{j+1}).

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, con una spline s_{3} di grado 3 cioè determinare s_{3} tale che s_{3}(x_{k}) = y_{k} per k = 0, …, n. Questo caso è più complesso di quello delle splines di grado 1. Mostriamo che esistono splines s_{3} interpolanti ma non sono uniche.

In ogni sottointervallo [x_{j}, x_{j+1}] una spline s_{3} è un polinomio p_{3,j}(x) = ∑_{i=0}^{3} a_{j}x^{j} di grado 3 e quindi definita da 4 coefficienti a_{k}. In totale i gradi di libertà sono 4n.

Le condizioni di interpolazione nei punti x_{0}, …, x_{n} sono n + 1. Le condizioni di regolarità fino alla derivata seconda nei punti x_{1}, …, x_{n−1} sono 3n − 3 ed quindi le condizioni imposte sono in totale 4n − 2 = (n + 1) + (3n − 3). Dal confronto coi gradi di libertà , è quindi evidente che servano 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).

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.

Alvise Sommariva  Splines  9/ 25

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}]. Sia p_{i,Δ} = s_{1,Δ}|_{[x_{i},x_{i+1}]} ∈ ℙ_{1} dove ℙ_{1} è l’insieme dei polinomi di primo grado. Dal teorema dell’errore dell’interpolazione polinomiale nel caso N = 1 abbiamo che puntualmente

f (x) − p_{i,Δ}(x) = f^{(2)}(ξ)\frac{(x − x_{i})(x − x_{i+1})}{2}, ξ ∈ (x_{i}, x_{i} + 1) (3)

da cui, posto h_{i,Δ} = (x_{i+1} − x_{i}), M_{i,Δ} = max_{x∈(x_{i},x_{i+1})} |f^{(2)}(x)| è immediata la stima d’errore per f ∈ C^{2}([a, b])

|f (x) − s_{1,Δ}(x)| = |f (x) − p_{i,Δ}(x)| = |f^{(2)}(ξ)|\frac{|x − x_{i}| · |x − x_{i+1}|}{2}
≤ M_{i} \frac{h_{i,Δ}^{2}}{2} , x ∈ [x_{i}, x_{i+1}]
Alvise Sommariva  Splines  10/ 25

Splines lineari: stima errore

In realtà si può mostrare, con tecniche più precise, che

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

Quindi, essendo f ∈ C^{2}([a, b]), per il teorema di Weierstrass esiste ed è finito M = max_{x∈[a,b]} |f^{(2)}(x)|, ed evidentemente

M_{i,Δ} := max_{x∈[x_{i},x_{i+1}]}|f^{(2)}(x)| ≤ max_{x∈[a,b]}|f^{(2)}(x)| := M. Di conseguenza

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}.
Alvise Sommariva  Splines  11/ 25

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_{k} 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  (5)

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

Teorema. 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)}||_{∞}  (6)

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)}||_{∞},  (7)
∥f^{(2)} − s_{3,Δ}^{(2)}∥_{∞} ≤ c_{2} h_{Δ}^{2}||f^{(4)}||_{∞},  (8)
∥f^{(3)} − s_{3,Δ}^{(3)}∥_{∞} ≤ c_{3} h_{Δ}^{1}||f^{(4)}||_{∞},  (9)

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.

Alvise Sommariva  Splines  14/ 25

Splines cubiche: 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_{3,Δ_{n}} ad una funzione f ∈ C^{2}([a, b]) è uniforme. Infatti da (5):

∥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}  (10)

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 ’ ) ;
Alvise Sommariva  Splines  19/ 25

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.12 2.36 2.36 1.46 . 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.50 e−02

>>

Alvise Sommariva  Splines  24/ 25

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.91 e−03

>>

Alvise Sommariva  Splines  25/ 25