Transcript Introduzione alle splines

Introduzione alle splines
Alessandra Sestini
28 dicembre 2016
1
Splines polinomiali
Sia [a , b] un intervallo della retta reale e sia
a = τ0 < · · · < τL = b
una partizione assegnata su tale intervallo. Chiameremo il seguente insieme
T := {τ0 , . . . , τL }
vettore dei nodi. Posto Ii := [τi , τi+1 ), i = 0, . . . , L − 2 e IL−1 := [τL−1 , τL ],
indichiamo con Pm,T lo spazio delle polinomiali a tratti di grado m e nodi in
T,
Pm,T = {f : [a , b] → R | f| Ii ∈ Πm , i = 0, . . . , L − 1} ,
(1)
dove Πm denota lo spazio dei polinomi di grado ≤ m . Come prima osservazione, notiamo che una generica funzione di Pm,T nei nodi interni τi , i =
1, . . . , L − 1, può essere discontinua, presentando in particolare un salto finito. Osserviamo inoltre che Pm,T forma uno spazio vettoriale di funzioni sul
campo R la cui dimensione è,
dim(Pm,T ) = (m + 1) L ,
essendo in generale i tratti polinomiali che compongono una generica funzione
di Pm,T indipendenti l’uno dall’altro e essendo (m + 1) la dimensione di Πm .
Ovviamente una rappresentazione di funzioni in PmT si ottiene scegliendo una
base di Πm (per esempio la base delle potenze ma anche la base di Bernstein
che è spesso preferita in ambiente grafico) da utilizzare per rappresentare
1
ciascuno dei suoi tratti polinomiali. Osserviamo infine che, comunque sia
definito l’insieme dei nodi T , risulta
Πm ⊂ Pm,T .
Poiché le polinomiali a tratti non sono neanche continue in [a , b], siamo
interessati a definire degli spazi di funzioni che siano intermedi fra Πm e Pm,T
che chiameremo appunto funzioni splines, distinguendole in splines classiche
e generalizzate. Le splines classiche formano un sottospazio vettoriale di Pm,T
che indicheremo con Sm,T e che possiamo definire come segue,
Sm,T := Pm,T ∩ C m−1 [a , b] .
(2)
Si tratta quindi di polinomiali a tratti di grado m (o, come anche si
dice, di ordine k := m + 1) e nodi in T che nei nodi interni hanno regolarità
C m−1 , ossia la massima regolarità con cui due polinomi di grado m si possono
raccordare in un punto senza essere necessariamente coincidenti. Le più
semplici sono naturalmente le splines (classiche) lineari il cui grafico è una
spezzata (con angoli in corrispondenza dei nodi interni), le più note sono le
splines cubiche che sono globalmente C 2 e vengono spesso introdotte come
alternativa ai polinomi per risolvere problemi di interpolazione. Osserviamo
che in effetti risulta
Πm ⊂ Sm,T ⊂ Pm,T .
Le splines generalizzate sono state introdotte in ambiente CAGD per avere a disposizione spazi di funzioni più duttili delle splines classiche delle
quali costituiscono appunto una generalizzazione, permettendo di imporre
regolarità variabile nei nodi interni. Per definirle, oltre a specificarne il grado m (o l’ordine k) e il vettore dei nodi T , occorre specificare un ulteriore vettore M = (m1 , . . . , mL−1 ) detto vettore delle molteplicità e tale che
mi ∈ N , 1 ≤ mi ≤ k, i = 1, . . . , L − 1 . Con tale ausilio, lo spazio Sm,T,M delle
splines generalizzate di grado m, nodi in T e molteplicità in M è definito
come segue,
Sm,T ,M := {f ∈ Pm,T | f (j) (τi− ) = f (j) (τi+ ), 0 ≤ j ≤ m − mi , 1 ≤ i ≤ L − 1}
(3)
Osserviamo che, se mi = k, nel nodo τi non è richiesta alcuna regolarità
e, se tutti gli mi in M sono pari a 1, lo spazio SmT ,M viene a coincidere con
Sm,T . Valgono inoltre le seguenti inclusioni
Πm ⊂ Sm,T ⊂ Sm,T ,M ⊂ Pm,T .
2
Dopo aver definito gli spazi delle splines, siamo interessati a individuarne la dimensione e a capire come sia possibile rappresentare le funzioni che
appartengono loro. A quest’ultimo riguardo, osserviamo che, trattandosi comunque di polinomiali a tratti di grado m con nodi in T , è sempre possibile
utilizzare la forma generale di rappresentazione di funzioni di PmT , anche
se in questo caso sarà necessario inserire dei vincoli sui parametri di rappresentazione in modo da garantire la regolarità richiesta nei nodi interni 1 .
Riferendoci dapprima a Sm,T , osserviamo infatti che il requisito di regolarità
C m−1 nel nodo τi corrisponde a m condizioni lineari linearmente indipendenti,
il che implica che
dim(Sm,T ) = (m + 1) L − m(L − 1) = L + m = (L − 1) + k ,
ossia la dimensione di SmT è data dalla somma del numero di nodi interni
con l’ordine. Con un ragionamento analogo si evince poi che
L−1
X
dim(Sm,T ,M ) = (m + 1) L −
(m − mi + 1) = µ + k ,
i=1
dove
µ :=
L−1
X
mi ,
i=1
ossia la dimensione di Sm,T ,M è data della somma delle molteplicità dei nodi
interni con l’ordine2 . Allo scopo di evitare inutili ridondanze, siamo quindi
interessati ad individuare delle basi di Sm,T e di Sm,T ,M , in particolare basi
che abbiano proprietà utili in ambiente grafico.
La prima base che possiamo definire per SmT viene detta base delle
potenze troncate ed è definita come segue,
m
{ 1, x , . . . , xm , (x − τ1 )m
+ , . . . , (x − τL−1 )+ } ,
(4)
dove per un qualsiasi punto fissato τ e per ogni intero positivo `, la funzione
(x − τ )`+ è definita come segue,
(x − τ )` , se x ≥ τ ,
`
(x − τ )+ :=
0,
altrimenti.
1
Questa è per esempio la forma utilizzata più comunemente per definire le splines
cubiche interpolanti nei nodi.
2
Si osservi che quando tutte le molteplicità sono pari a 1 si riottiene la dimensione di
Sm,T , come deve essere.
3
Osserviamo che le prime k funzioni che formano in (4) la base delle potnze
troncate altro non sono che la base delle potenze per Πm che è un sottospazio
di SmT . Per le altre L − 1 funzioni, le potenze troncate appunto, è immediato verificare che si tratta di funzioni appartenenti a Pm,T con due tratti
polinomiali (dei quali uno identicamente nullo). Verificando che le derivate
sinistre e destre fino all’ordine m − 1 di tali funzioni nel nodo di raccordo fra
i due tratti polinomiali sono coincidenti, si deduce quindi che esse sono globalmente di classe C m−1 e quindi che sono comunque funzioni appartenenti a
Sm,T . Poiché le funzioni in (4) sono in tutto proprio k + (L − 1), per verificare
che esse sono una base di Sm,T occorre verificare che esse sono linearmente
indipendenti, come viene fatto nella seguente proposizione,
Proposizione 1. Le funzioni di Sm,T elencate in (4) sono linearmente indipendenti.
Dimostrazione : Per dimostrare l’asserto facciamo vedere che se una combinazione lineare delle funzioni in (4) è identicamente nulla, allora i suoi
coefficienti sono tutti nulli. Sia quindi
m
X
i=0
ai x i +
L−1
X
bj (x − τj )m
+ ≡ 0.
j=1
In particolare quindi la combinazione lineare si annullerà per ogni x ∈ I0 .
Ma se x appartine a tale sottointervallo la combinazione lineare si riduce
alla sola prima sommatoria per definizione di potenze troncate. Poiché un
polinomio di grado m non nullo può avere al più m radici distinte, segue
che a0 = · · · = am = 0. Quindi avendo dimostrato questo, la combinazione
lineare identicamente nulla si riduce alla seconda sommatoria. Prendendo
però x ∈ I1 essa si riduce solo a b1 (x − τ1 )m che si deve annullare ∀x ∈ I1 . Ne
segue che quindi b1 = 0. Ragionando iterativamente allo stesso modo sui vari
sottointervalli, se ne deduce che anche b1 = · · · = bL−1 = 0, il che conclude
la dimostrazione.
La base delle potenze troncate può essere generalizzate come segue per
definire una base di Sm,T ,M ,
m−mi +1
{ 1, x , . . . , xm , (x − τi )m
, i = 1, . . . , L − 1 } .
+ , .., (x − τi )+
(5)
Esercizio 1. Verificare che le funzioni introdotte in (5) formano effettivamente una base di SmT ,M .
4
Sebbene sia stata introdotta per prima e sia di facile definizione, la base
delle potenze troncate presenta degli aspetti non attraenti per le applicazioni. Infatti innanzi tutto si può osservare che i suoi elementi non hanno
supporto locale, cosa di estrema utilità in ambiente grafico. Inoltre le matrici di collocazione ad essa abbinate con le quali si ha a che fare quando si
utilizzano le splines per risolvere problemi di interpolazione risultano spesso
malcondizionate e quindi il loro utilizzo non è raccomandabile in tale contesto
3
.
Tenute presenti le suddette osservazioni, siamo quindi interessati a introdurre una base diversa per le splines che nel prossimo paragrafo introdurremo in modo costruttivo mediante una formula di recursione. In particolare ci
riferiremo dapprima a Sm,T e dopo generalizzeremo la definizione per Sm,T,M .
2
B–splines (a nodi semplici)
Per poter introdurre recursivamente le B–splines di ordine k e nodi in T
definite in [a , b] = [τ0 , τL ], abbiamo bisogno di associare a T il cosiddetto
vettore esteso dei nodi T
T := {ti , i = 0, . . . , n + k} ,
dove n + 1 = k + (L − 1) = dim(Sm,T ) , ti ≤ ti+1 e
{tk−1 , . . . , tn+1 } = {τ0 , . . . , τL } .
(6)
I primi k − 1 nodi t0 , . . . , tk−2 si dicono nodi ausiliari sinistri cosı̀ come gli
ultimi k − 1 si dicono ausiliari destri e possono essere scelti arbitrariamente
purché non decrescenti (anche tutti coincidenti). Essi costituiscono un ausilio
necessario per la definizione delle B-splines che definiremo recursivamente. A
tale scopo indicheremo con Ni,r (t) la i–esima B–spline di ordine r , r ≥ 1 e
siamo interessati al terminea definire le Ni,k (t), i = 0, . . . , n. Definiamo quindi
innanzitutto le B–splines Ni,1 (t), i = 0, . . . , n + k − 1, t ∈ IR di ordine 1 che
sono le seguenti funzioni non negative “a scalino” ,
1 se t ∈ [ti , ti+1 ) ,
Ni,1 (t) =
(7)
0 altrimenti .
3
Questo è il motivo per esempio per cui non si usano per calcoalare le spline cubiche
interpolanti nei nodi.
5
Osserviamo che formalmente (7) implica che se ti = ti+1 , la Ni,1 (t) è identicamente nulla.
A partire dalle B-splines di ordine 1 e riferendoci sempre al vettore esteso dei
nodi T, è possibile recursivamente definire come segue le B–splines Ni,r (t), i =
0, . . . , n + k − r di ordine r ≤ k,
Ni,r (t) , = ωi,r (t) Ni,r−1 (t) + (1 − ωi+1,r (t)) Ni+1,r−1 (t) , i = 0, ..., n + k − r (8)
dove ωi,r (t) è il seguente polinomio lineare (o, eventualmente, identicamente
nullo),
t−ti
se ti < ti+r−1 ,
ti+r−1 −ti
(9)
ωi,r (t) :=
0
altrimenti.
In questo caso possiamo osservare che dalla definizione discende che la Ni,r (t)
è identicamente nulla se ti = · · · = ti+r e anche che per r > 1 risulta Ni,r (ti ) =
Ni,r (ti+r ) = 0.
Esercizio 2. Scrivere l’espressione analitica e fare il grafico delle B–splines
lineari e quadratiche con vettore esteso dei nodi T = {0 , 1 , 2 , . . . , 4 , 5 , 6} .
Esercizio 3. Verificare che se T = {0, . . . , 0, 1, . . . , 1} con 0 e 1 entrambi
di molteplicità k in T, alllora le B–splines Ni,k , i = 0, . . . , k coincidono con i
polinomi di Bernstein Bik−1 (t), i = 0, . . . , k in [0 , 1].
Inoltre per induzione su r, si dimostrano le seguenti proprietà delle B–
splines,
• P1) Ni,r (t) = 0 , se t ∈
/ [ti , ti+r ] , (supporto locale)
• P2) Ni,r (t) ≥ 0 , ∀t ∈ IR (non negatività)
• P3)
n+k−r
X
Ni,r (t) = 1 , ∀t ∈ [tr−1 , tn+1+k−r ] , r = 1, . . . , k.
i=0
La proprietà P1) giustifica perché il sottoinsieme
Nir := {ti , · · · , ti+r }
(10)
di T viene detto insieme dei nodi attivi per la B–spline Ni,r (t). Poiché le
proprietà P1) e P2) sono immediate da verificare, dimostriamo qui solo la
P3).
6
Dimostrazione : Sia quindi r = 1 e t ∈ [t0 , tn+k ]. Se t = tn+k , tutte le B–
splines di ordine 1 si annullano in t ad eccezione dell’ultima che vale proprio
1. Supponiamo quind t 6= tn+k . Allora esisterà ti ∈ T con i ≤ n + k − 1 e
ti < ti+1 tale che t ∈ [ti , ti+1 ). Quindi dalla definizione segue che Ni,1 (t) = 1 e
Nj,1 (t) = 0 se j 6= i. Quindi segue l’asserto per r = 1. Supponiamo quindi che
esso valga per r − 1 e facciamo vedere che vale anche per r. Dalla definizione
recursiva si ha
n+k−r
X
Ni,r (t) =
n+k−r
X
i=0
[ωi,r (t) Ni,r−1 (t) + (1 − ωi+1,r (t)) Ni+1,r−1 (t)] .
i=0
Ma, se t ∈ [tr−1 , tn+1+k−r ), allora sarà necessariamente N0,r−1 (t) = 0 e
Nn+k−r+1,r−1 (t) = 0. Quindi, raccogliendo a secondo membro rispetto alle
B-splines, si ottiene
n+k−r
X
Ni,r (t) =
n+k−r
X
i=0
Ni,r−1 (t) =
i=1
n+k−r+1
X
Ni,r−1 (t) = 1 .
(11)
i=1
Si osservi che se t = tn+1+k−r e r > 2 la catena di uguaglianze in (11) vale
ancora. Inoltre, se r = 2 e t = tn+k−1 , si ha invece che Nn+k−r+1,r−1 (t) =
Nn+k−1,1 (t) = 1 e ωn+k−1,2 (t) = 0. Quindi in questo caso si ottiene che in
(11) vale l’uguaglianza fra il primo e il terzo termine che in effetti si riduce
all’ultimo addendo pari a 1.
Definiamo allora il seguente spazio lineare Sk,T di funzioni definite in [a , b]
e generato dalle B–splines,
Sk,T := hN0,k , · · · , Nn,k i .
(12)
Siamo interessati a far vedere, dimostrando il Teorema di Curry–Shoenberg,
che esso risulta coincidente con Sm,T , m = k − 1 che ha dimensione proprio
n + 1. Da ciò potremo quindi dedurre che le B–splines Ni,k , i = 0, . . . , n ne
costituiscono proprio una base. Organizziamo la dimostrazione in due parti.
Parte 1: Πk−1 ⊆ Sk,T
Dimostrazione : Dimostriamo l’asserto provando la seguente identitità,
detta identitità di Mardsen,
k−1
(x − θ)
=
n
X
ψj,k (θ) Nj,k (x) ,
j=0
7
∀x ∈ [a , b] ,
∀θ ∈ IR ,
dove ψj,k (θ) è il seguente polinomio di grado k − 1,
ψj,k (θ) :=
k−1
Y
(tj+r − θ) .
r=1
La dimostrazione è ancora per induzione su k. Se k = 1, si ha che ∀j, ψj,1 (θ) =
1, ∀θ ∈ IR. Quindi, avendo dimostrato che le B–splines di ordine k formano in
[a , b] = [tk−1 , tn+1 ] una partizione dell’unità, segue la validità dell’asserto.
Supponiamo ora che esso sia valido per k − 1 e facciamo vedere che allora
vale anche per k. Usiamo quindi a tale scopo la definizione recursiva delle
B–splines, ottenendo
n
X
ψj,k (θ) Nj,k (x) =
j=0
n
X
ψj,k (θ) [ωj,k (x) Nj,k−1 (x)+(1−ωj+1,k (x)) Nj+1,k−1 (x)] .
j=0
Tenendo conto che per x ∈ [a , b) risulta N0,k−1 (x) = Nn+1,k−1 (x) = 0 (ciò
è vero anche in b se k > 2), raccogliendo a secondo membro rispetto alle
B–splines di ordine k − 1 e indicando con M s il membro sinistro, si ottiene,
Ms =
n
X
[ψj,k (θ) ωj,k (x) + ψj−1,k (θ) (1 − ωj,k (x)) ] Nj,k−1 (x) .
j=1
Raccogliendo quindi in ciascun addendo a secondo membro il massimo comun
divisore fra ψj,k e ψj−1,k che è ψj,k−1 e ricordando la definzione di ωj,k (x), si
ha
Ms =
n
X
j=1
n
X
j=1
n
X
ψj,k−1 (θ) [(tj+k−1 − θ) ωj,k (x) + (tj − θ) (1 − ωj,k (x)) ] Nj,k−1 (x) =
ψj,k−1 (θ) [(tj+k−1 − tj ) ωj,k (x) + (tj − θ) ] Nj,k−1 (x) =
n
X
ψj,k−1 (θ) (x − θ) Nj,k−1 (x) = (x − θ)
ψj,k−1 (θ) Nj,k−1 (x) .
j=1
j=1
Per induzione si completa quindi la dimostrazione (analizzare separatamente
il caso k = 2, x = b).
Si osservi che per, l’arbitrarietà di θ, l’identità di Mardsen implica che Πk−1 ⊆
Sk,T .
8
Parte 2: Le potenze troncate (x − ti )k−1
+ , i = k, . . . , n appartengono a Sk,T ,
essendo
(x −
ti )k−1
+
=
n
X
ψj,k (ti ) Nj,k (x) ,
∀x ∈ [a , b] ,
i = k, . . . , n .
(13)
j=i
Dimostrazione : Consideriamo una potenza troncata (x−ti )k−1
+ , i = k, . . . , n.
Da Mardsen sappiamo che
k−1
(x − ti )
=
n
X
ψj,k (ti ) Nj,k (x) ,
∀x ∈ [a , b] .
j=0
Osserviamo però che, se i − k < j < i, risulta per definizione ψj,k (ti ) = 0.
Quindi
(x − ti )
k−1
=
i−k
X
ψj,k (ti ) Nj,k (x) +
j=0
n
X
ψj,k (ti ) Nj,k (x) ∀x ∈ [a , b] .
j=i
Ricordando la definizione di potenza troncata e ricordando che la B–spline
j–esima di ordine k ha supporto compatto pari a [tj , tj+k ], da ciò segue
facilmente la validità della relazione in (13).
Osserviamo che da quanto dimostrato segue necessariamente che Sk−1,T =
Sk,T , che le B–splines Ni,k , i = 0, . . . , n sono polinomiali a tratti di grado k −1
con nodi in T e regolarità C k−2 in [a , b] e che esse costituiscono una base di
Sk−1,T . Inoltre, sempre per induzione su k, si può dimostrare che, per k > 2,
la derivata di Ni,k si può scrivere come segue,
Ni,k−1 (t)
Ni+1,k−1 (t)
dNi,k
(t) = (k − 1) [
−
].
dt
ti+k−1 − ti
ti+k − ti+1
3
B–splines (a nodi multipli)
Le B–splines possono essere generalizzate in modo da definire una base per lo
spazio Sm,T ,M . In particolare occorre generalizzare la definizione del vettore
esteso dei nodi introdotto in (6) come segue,
mL−1
m1
,
z }| {
z
}|
{
T := {tk−1 , . . . , tn+1 } = {τ0 , τ1 , . . . , τ1 , . . . , τL−1 , . . . , τL−1 , τL }
9
dove si assume 1 ≤ mi ≤ k, i = 1, . . . , L − 1 e mi ∈ IN dicesi molteplicità del
nodo τi nel vettore esteso dei nodi. T è quindi un vettore in cui ogni nodo
interno τi di T è ripetuto mi volte e che al solito ingloba i nodi ausiliari.
Osserviamo poi che risulta n − k + 3 = µ + 2, cioè n + 1 = µ + k che è
la dimensione di Sm,T ,M e quindi siamo interessati sempre a definire n +
1 B–splines Ni,k , i = 0, . . . , n, di ordine k. A tale scopo si usa sempre la
definizione recursiva introdotta in (8) nella sezione precedente, dove le Ni,1 (t)
sono sempre definite come in (7). Si osservi che la presenza in T di nodi ti
tali che sia almeno ti = ti+1 anche per k < i < n (corrispondenti a nodi
multipli in T ), fa sı̀ che per qualche r il polinomio ωi,r (t) usato in (8) possa
essere identicamente nullo e ciò determina una riduzione della regolarità di
tutte le B–splines Nj,k il cui insieme di nodi attivi contiene il nodo ti con
molteplicità maggiore di 1. Più precisamente si dimostra che, se in nodo ti
ha molteplicità ` in Njr , allora la B–spline Nj,r in ti è C r−`−1 .
Esercizio 4. Ricavare l’espressione analitica delle B–splines N2,2 , N3,2 e di
N1,3 , N2,3 , N3,3 e farne il grafico in IR con T = {0 , 1 , 2 , 3 , 3 , 4 , 5 , 6} .
Esercizio 5. Ricavare T e M a partire dal vettore esteso dei nodi definito
nell’esercizio precedente per k = 2 e per k = 3 e dire quindi su quale intervallo
[a , b] della retta reale le B–splines di ordine k definibili mediante T formano
una base per Sk−1,T,M .
Per quanto riguarda le B–splines a nodi multipli, ossia definite relativamente al vettore esteso dei nodi T come detto sopra, è possibile dimostrare
che valgono ancora le proprietà P1), P2) e P3). Inoltre si dimostra con
tecniche simili a quelle utilizzate per le B–splines classiche la validità del teorema di Curry–Shoenberg che asserisce in questo caso che lo spazio Sk−1,T ,M
coincide con Sk,T e che le B–splines ne costituiscono ancora una base.
Osservazione 1. Data l’arbitrarietà con cui si possono definire i nodi ausiliari sinistri e destri in T, in genere si tende a fare la seguente scelta di nodi
ausiliari multipli,
t0 = · · · = tk−2 = tk−1 = a ,
b = tn+1 = · · · = tn+2 = tn+k .
(14)
Esercizio 6. Verificare che la scelta (14) garantiscel’annullamento di tutte
le B–splines di ordine k in a e in b, rispettivamente ad eccezione di N0,k che
in a vale 1 e di Nn,k che in b vale 1.
10
Fissato quindi l’ordine k delle splines con le quali si vuole lavorare (ossia il
grado m = k −1) e il vettore dei nodi T e il relativo vettore delle molteplicità
M (e di conseguenza fissato il vettore esteso dei nodi T ), una qualunque spline
s appartenente a Sk−1,T ,M si può univocamente rappresentare mediante la
base delle B–spline,
s(t) =
n
X
ci Ni,k (t),
t ∈ [a , b] .
i=0
Tale rappresentazione delle funzioni di Sk−1,T ,M viene usualmente detta B–
form.
Esercizio 7. Utilizzare il comando spmak del toolbox splines del Matlab per
definire una spline in un fissato spazio Sk−1,T ,M in B–form. Farne quindi il
grafico dopo averla tabulata.
Nelle sezioni seguenti per semplicità indicheremo sempre con Sk,T uno
spazio di splines, facendo quindi riferimento alla loro rappresentazione in
B–form.
4
Curve B–spline
Una curva B–spline X : [a , b] → E d è una curva parametrica le cui componenti sono funzioni spline di un certo grado assegnato m (ordine k = m + 1)
espresse nella base delle B–splines di un certo spazio di splines generalizzate Sm,T ,M , con T = {τi , i = 0, . . . , L} vettore dei break–points assegnato
(a = τ0 < · · · < τL = b) e M = {mi , i = 1, . . . , L − 1} vettore delle molteplicità pure assegnato (con mi ∈ IN , 1 ≤ mi ≤ k). Quindi una curva di tale
tipo può essere scritta come segue
X(t) :=
n
X
di Ni,k (t) ,
t ∈ [a , b] ,
i=0
dove i di ∈ E d sono punti assegnati in modo ordinato dove la spezzata che
li congiunge dicesi poligono di controllo (di de Boor) relativo alla curva X.
Fissato lo spazio Sm,T ,M , abbiamo visto che per definire la corrispondente
base delle B–splines Ni,k , i = 0, . . . , n basta definire il corrispondente vettore
esteso dei nodi T = {t0 , . . . , tn+k } a partire dal grado m, da T e da M e in
11
generale la funzione vettoriale X sarà C k−mi −1 nel break–point τi , , sebbene
naturalmente possano esistere particolari poligoni di controllo tali che una o
tutte le componenti della X risultino in τi più regolari (si pensi che i polinomi di grado m possono essere visti come particolari funzioni appartenenti a
Sm,T ,M ).
Grazie alle proprietà delle B–splines, le curve B–splines godono delle
seguenti proprietà:
• invarianza per trasformazioni affini (conseguenza della proprietà di
partizione dell’unità)
• località (conseguenza della proprietà di supporto compatto)
• strong convex–hull (conseguenza di partizione unità, non negatività e
supporto compatto)
• end–point interpolation se il vettore esteso dei nodi ha nodi ausiliari
multipli
• possibilità di definire curve chiuse (se gli ultimi k − 1 punti di controllo
coincidono con i primi k − 1 e il vettore esteso dei nodi è ciclico (ossia
hi := hn−k+2+i , hn+1+i := hk−1+i , i = 0, . . . , k − 2).
Si noti che le curve di Bezier possono essere ottenute come un caso particolare
di curve B–spline usando T = {a , b} e nodi ausiliari multipli. Le curve B–
spline hanno però il vantaggio di svincolare il grado dal numero di punti di
controllo. Inoltre godono della proprietà di località che risulta davvero utile
per il disegno computerizzato e che permette di dimostrare quanto asserito
nel seguente esercizio,
Esercizio 8. Verificare che se k punti di controllo consecutivi sono allineati
allora la curva X per un tratto (corrispondente a un sottointervallo) si schiaccia sul segmento che li congiunge. Se sono solo k − 1 i punti di controllo
consecutivi allineati allora X tocca il segmento che li congiunge. Infine se il
poligono è degenere con k − 1 punti di controllo consecutivi coincidenti con
lo stesso punto, allora la curva passa per tale punto.
Per quanto riguarda la derivazione, se X è derivabile in [a , b] (quindi ogni
mi < k − 1), si dimostra che la sua derivata può essere rappresentata come
12
segue,
n
X
dX
(1)
= (k − 1)
di Ni,k−1 (t) ,
dt
i=1
t ∈ [a , b] ,
dove
di − di−1
, i = 1, . . . , n .
ti+k−1 − ti
Si osservi che, nel caso ci siano k − 1 punti di controllo consecutivi coincidenti, questa formula implica che nel break–point τi in cui tale punto viene
(τi ) = 0.
interpolato si ha che risulta dX
dt
Le curve B–splines nel piano godono anche della proprietà di variation–
diminishing che, come abbiamo visto parlando di curve di Bézier, significa
che il numero di intersezioni di una curva B–spline piana con una qualunque
retta r è minore o uguale del numero di intersezioni fra r e il suo poligono di
controllo. Tale proprietà deriva dal fatto che le B–splines formano un sistema
di funzioni totalmente positivo in [a , b].
(1)
di
4.1
:=
Algoritmi per curve B–spline
Per valutare in un punto t o in una serie di punti in [a , b] una curva B–spline
X si può utilizzare l’algoritmo di de Boor che ha diverse somiglianze con
l’algoritmo di de Casteljau per curve di Bézier. Una differenza significativa
consiste nel fatto che preliminarmente occorre individuare il sottointervallo
[tj , tj+1 ) cui appartiene il punto t nel quale si vuole calcolare X(t), per cui
si può scrivere
j
X
X(t) =
di Ni,k (t) .
i=j−k+1
L’idea è allora quella di eseguire k − 1 passi in ciascuno dei quali si riscrive
X(t) in termini di B–splines diminuite di un ordine. Più precisamente, posto
(0)
di := di , i = j − k + 1, . . . , j , utilizzando la formula recursiva per le B–
splines, al primo passo si scrive
X(t) =
j
X
(1)
di (t)Ni,k−1 (t) ,
i=j−k+2
dove
(1)
(0)
(0)
di (t) := (1 − ωi,k (t)) di−1 + ωi,k (t) di , i = j − k + 2, . . . , j ,
13
e dove si ricorda che
ωi,s (t) :=
t−ti
ti+s−1 −ti
se ti < ti+s−1 ,
altrimenti.
0
In generale, al passo r–esimo scriveremo
j
X
X(t) =
(r)
di (t)Ni,k−r (t) ,
i=j−k+r+1
con
(r)
(r−1)
(r−1)
di (t) := (1 − ωi,k−r+1 (t)) di−1 + ωi,k−r+1 (t) di
, i = j − k + r + 1, . . . , j .
(k−1)
(t) . In tutto occorre calcolare
Al passo (k −1)–esimo otteniamo X(t) = dj
quindi k(k − 1)/2 nuovi punti di controllo ciascuno definito come combinazione convessa di due punti di controllo consecutivi del passo precedente.
Vediamo ora un altro algoritmo molto potente, detto algoritmo di knot–
insertion che è definibile per curve B–spline. Osserviamo che, se T e T̂ sono
due vettori estesi dei nodi entrambi ammissibili per definire splines di ordine
k su [a , b] con T ⊂ T̂ , si deduce dalla definizione di splines che se una
funzione è combinazione lineare di B–splines associate a T essa deve anche
essere esprimibile come combinazione lineare di B–splines associate a T̂ . Tale
osservazione ci porta a dire che in particolare, se
X(t) =
n
X
di Ni,k (t) , t ∈ [a , b],
i=0
dove le B–splines sono definite utilizzando T, allora deve essere possibile
determinare d̂i , i = 0, . . . , n + 1 tali che
X(t) =
n+1
X
d̂i N̂i,k (t) , t ∈ [a , b],
i=0
dove con N̂i,k (t), i = 0, . . . , n+1 si sono indicate le B–splines di ordine sempre
k costruite utilizzando il vettore esteso dei nodi T̂ = T ∪ {t̂}, con t̂ ∈ (a , b)
(che quindi può corrispondere a un nuovo break–point o può essere utilizzato
per aumentare la molteplicità di un vecchio break–point). Omettendo la
14
dimostrazione, riportiamo qui soltanto l’espressione dei d̂i , i = 0, . . . , n + 1,
supponendo che t̂ ∈ [tr , tr+1 ),

se 0 ≤ i ≤ r − k + 1 ,
 di ,
d̂i :=
(1 − ωi,k (t̂)) di−1 + ωi,k (t̂)) di , se r − k + 2 ≤ i ≤ r ,

di−1 ,
se r + 1 ≤ i ≤ n + 1 .
Un’applicazione interessante di questo algoritmo consiste nel portare ogni
nodo ad avere molteplicità k ottenendo quindi la rappresentazione Bézier–
spline della curva stessa (ossia polinomiale a tratti).
5
Interpolazione e approssimazione ai minimi
quadrati in Sk,T,M
Le splines costituiscono una generalizzazione estremamente duttile dei polinomi pur mantenendone la facilità di impego. Questo spiega perchè esse
abbiano largo impiego sia nell’ambito dell’approssimazione che in quello delle
applicazioni grafiche.
Per quanto riguarda l’interpolazione, siamo innanzi tutto interessati a
stabilire se esiste ed è unica la spline s appartenente a Sk,T tale che
s(xj ) = fj , j = 0, . . . , n ,
(15)
dove a ≤ x0 < · · · < xn ≤ b (con a = tk−1 e b = tn+1 ), e dove fj denota
il valore noto assunto in xj dalla funzione f (x) da interpolare. Al solito le
ascisse xi , i = 0, . . . , n, si dicono ascisse di interpolazione e si suppongono
distinte. Ovviamente, usando la B–form per rappresentare s, il problema si
riconduce a stabilire se è nonsingolare la seguente matrice A di dimensione
(n + 1) × (n + 1),


N0,k (x0 ) · · · Nn,k (x0 )


..
..
..
A := 
(16)
,
.
.
.
N0,k (xn ) · · · Nn,k (xn )
detta matrice di collocazione delle B–splines nelle ascisse xi , i = 0, . . . , n.
Infatti le condizioni di interpolazione in (15) corrispondono ad imporre il
seguente sistema lineare,
Ac = f ,
15
dove c := (c0 , , · · · , cn )T è il vettore delle incognite corrispondenti ai coefficienti di s in B–form e f := (f0 , , · · · , fn )T è il termine noto del sistema.
Ora, a differenza di quanto accade per l’interpolazione polinomiale, è facile
dimostrare che nel caso delle splines la sola ipotesi di ascisse di interpolazione distinte non è sufficiente a garantire la nonsingolarità della matrice A
e quindi l’esistenza e l’unicità della spline interpolante. Se per esempio si
scelgono tutte le ascisse xi < tn , i = 0, . . . , n, l’ultima colonna di A è tutta
nulla e quindi necessariamente A risulterà singolare. E’ quindi evidente che,
per garantire la nonsingolarità di A sarà necessario assumere delle restrizioni
sulle ascisse di interpolazione legate alla distribuzione del vettore dei nodi.
Nel seguito, per k > 1, faremo sempre riferimento al caso di molteplicità
mi < k, i = 1, . . . .L − 1. Il teorema fondamentale a questo riguardo risulta
essere il cosiddetto teorema di Whitney–Schoenberg,
Teorema 1. Sia k > 1 e mi < k, i = 1, . . . , L − 1. La matrice di collocazione
A delle B–splines definita in (16) risulta nonsingolare sse le ascisse di interpolazione xi , i = 0, . . . , n, oltre ad essere distinte (e ordinate in senso crescente), verificano le seguenti disuguaglianze (dette di Whitney–Schoenberg)
,
Ni,k (xi ) > 0 , i = 0, . . . , n .
(17)
Dimostrazione : Facciamo intanto vedere che le condizioni in (17) sono
necessarie, osservando che esse corrispondono a chiedere ti < xi < ti+k , con la
possibilità di x0 = t0 (xn = tn+k ) nel caso di nodi ausiliari multipli.Escludendo
per ora questi ultimi due casi particolari, supponiamo allora per assurdo che
esista xi ∈
/ (ti , ti+k ). Se risulta xi ≤ ti , si ha che Nj,k (xi ) = 0 , ∀j ≥ i. Ma
essendo le ascisse di interpolazione ordinate in senso crescente, ciò implica
che Nj,k (xr ) = 0 , ∀j ≥ i, ∀r ≤ i . La matrice di collocazione A introdotta in
(16) sarebbe quindi una matrice a blocchi del tipo
A1,1 0
A=
,
(18)
A2,1 A2,2
con il blocco nullo di dimensione (i + 1) × (n − i + 1) . Essendo il blocco A2,2
di dimensione (n − i) × (n − i + 1), le ultime (n − i + 1) colonne della matrice
A sarebbero necessariamente linearmente dipendenti e quindi la matrice A
sarebbe singolare. Se invece risulta xi ≥ ti+k , si ha che Nj,k (xi ) = 0, ∀j ≤ i
e quindi Nj,k (xr ) = 0 , ∀j ≤ i, ∀r ≥ i . A allora avrebbe la struttura a blocchi
16
del tipo
A=
A1,1 A1,2
0 A2,2
,
(19)
con il blocco nullo di dimensione (n − i + 1) × (i + 1) . Essendo il blocco A2,2
di dimensione ((n − i + 1) × (n − i), le ultime (n − i + 1) righe della matrice
A sarebbero necessariamente linearmente dipendenti e quindi la matrice A
sarebbe singolare.
Facciamo ora vedere che le condizioni in (17) sono anche sufficienti. Facciamo però innanzi tutto vedere che, senza perdita di generalità, si può assumere ti+1 ≤ xi ≤ ti+k−1 , i = 0, . . . , n. Infatti nel caso esistesse una
ascissa di interpolazione xi tale che per esempio xi < ti+1 , allora la matrice
A definita in (16) sarebbe una matrice triangolare a blocchi con due blocchi diagonali quadrati rispettivamente di dimensione i + 1 e n − i e quindi
la sua nonsingolarità si otterrebbe dimostrando la nonsingolarità di queste
due sottomatrici. Il raginamento potrebbe essere ripetuto iterativamente fino ad arrivare a considerare valida l’ipotesi suddetta. Osserviamo che nel
caso k > 2 si può assumere addirittura, senza perdita di generalità, che
ti+1 < xi < ti+k−1 , i = 0, . . . , n (dato che, se xi = ti+1 o xi = ti+k−1 , è
ancora valido il ragionamento precedente).
Per k = 2 assumiamo quindi xi = ti+1 . Allora risulta Ni,2 (xi ) = 1 e
Ni,2 (xj ) = 0 se j 6= i. Quindi la matrice A coincide con la matrice identità che
è ovviamente nonsingolare. Consideriamo quindi il caso k ≥ 3, supponendo
ti+1 < xi < ti+k−1 , i = 0, . . . , n. Osserviamo che ciò implica anche che
mi ≤ k − 2, i = 1, . . . , L − 1 e quindi che Sk,T,M ⊂ C 1 [a , b]. Ora, se per
assurdo la matrice A fosse singolare, esisterebbe un vettore c ∈ IRn+1 non
identicamente nullo tale che Ac = 0, cioé esisterebbe una spline s ∈ Sk,T,M
non identicamente nulla che si annullerebbe nelle (n + 1) ascisse distinte
xi , i = 0, . . . n tali che ti+1 < xi < ti+k−1 . Dal teorema di Rolle seguirebbe
allora che esisterebbero n ascisse distinte yi , i = 0, . . . , n − 1, tali che ti+1 <
xi < yi < xi+1 < ti+k , ossia tali che ti+1 < yi < ti+k , i = 0, . . . , n − 1 nelle
quali dovrebbe annullarsi s0 ∈ Sk−1,T 0 ,M con T 0 = {t1 , . . . , tn+k−1 }. Quindi,
se k = 3 si arriverebbe ad un assurdo in quanto le ascisse yi , i = 0, . . . , n − 1,
verificano le condizioni di Whitney–Schoenberg relativamente all’ordine k −1
per T 0 . Se k > 3 si può quindi fare la dimostrazione per induzione.
Osservazione 2. Se k > 1 e mi < k, i = 1, . . . , L − 1 e le ascisse di interpolazione sono scelte come le cosiddette ascisse di Greville definite mediante
17
i nodi come segue,
xi := (ti+1 + · · · + ti+k−1 )/(k − 1) ,
i = 0,... ,n,
allora le condizioni di Whitney–Schoenberg sono verificate. Infatti in tal caso,
essedo i nodi in T non decrescenti e al più di molteplicità k − 1 necessariamente ti < xi < ti+k e quindi valgono le disuguaglianze in (17).
Osserviamo inoltre che nel caso lineare (k = 2) a nodi semplici n = L e
quindi i nodi τ0 = t1 , . . . , τL = tn+1 sono esattamente tanti quante sono le condizioni di interpolazione. In tal caso, scegliendo xi = τi , i = 0, . . . , n, le condizioni (17) sono sicuramente verificate e quindi la spline lineare a nodi semplici
interpolante nei nodi risulta univocamente determinata (il suo grafico altro
non è che la spezzata che congiunge i punti di coordinate (xi , fi ), i = 0, . . . , n).
Se però k > 2, risulta L = n − k + 2 < n e quindi se ci si limita a considerare
condizioni di interpolazione nei nodi non si hanno sufficienti condizioni per
caratterizzare la spline. Le splines cubiche interpolanti nei nodi per esempio richiedono 2 condizioni ausiliari opportune che vengono assegnate agli
estremi a = tk−1 e b = tn+1 .
Per quanto riguarda l’approssimazione ai minimi quadrati, date m n
ascisse xi , i = 0, . . . , m siamo interessati a determinare la spline s ∈ Sk,T,M
tale che sia minima la seguente quantità detta residuo,
m
X
(fi − s(xi ))2 .
i=0
Rappresentando s in B–form attraverso il vettore dei coefficienti c = (c0 ,
· · · , cn )T , tale problema può essere formalizzato come segue,
min
kf − Ack2
n+1
c∈IR
dove ora A è una matrice rettangolare di dimensione (m + 1) × (n + 1) che
sempre colloca le B–splines in tutte le ascisse xi , i = 0, . . . , m. Sappiamo
quindi dall’algebra lineare che il problema ammette soluzione unica se la matrice A ha rango massimo (n + 1) (e tale soluzione si calcola utilizzando la
fattorizzazione QR di A). Il rango massimo pari a (n + 1) risulta garantito
se è possibile estrarre (n + 1) righe da A tali che la corrispondente matrice quadrata An+1 risulti non singolare. Ciò significa che deve esistere una
sottosequenza xij , j = 0, . . . , n di xi , i = 0, . . . , m, di ascisse distinte che verifichino le condizioni di Whitney–Schoenberg, ossia tali che tj < xij < tj+k ,
per j = 0, . . . , n.
18
Riferimenti bibliografici
[1] C de Boor (2001), A Practical Guide to Splines, Revised edition,
Applied Mathematical Sciences 27, Springer–Verlag, New York.
19