W

W

y

y

k

k

ł

ł

a

a

d

d

6

6

:

:

I

I

n

n

t

t

e

e

r

r

p

p

o

o

l

l

a

a

c

c

j

j

a

a

f

f

u

u

n

n

k

k

c

c

j

j

a

a

m

m

i

i

s

s

k

k

l

l

e

e

j

j

a

a

n

n

y

y

m

m

i

i

(

(

s

s

p

p

l

l

a

a

j

j

n

n

y

y

)

)

W tym wykładzie omówimy problem interpolacji przy pomocy tzw. funkcji sklejanych,

zwanych też (żargonowo) splajnami. W przeciwieństwie do metod interpolacyjnych

opisanych w Wykładzie nr 1, gdzie stosowaliśmy jeden globalny wielomian dla całego

przedziału interpolacji, w metodzie splajnów stosowane są funkcje zdefiniowane jako

wielomiany niskiego stopnia osobno dla każdego odcinka pomiędzy sąsiednimi węzłami

interpolacyjnymi. Te lokalne wielomiany są jednak tak dobrane, aby – oprócz warunków

interpolacji – spełniać warunki sklejenia w taki sposób, aby cały splajn był funkcją o

odpowiedniej regularności. Skoncentrujemy się przede wszystkim na zagadnieniu interpolacji

za pomocą splajnu kubicznego, tj. funkcji ciągłej wraz z pochodnymi do rzędu drugiego

włącznie i zbudowanej z wielomianów 3-ego stopnia.

Należy wspomnieć, że funkcje sklejane (i ich daleko idące uogólnienia), maja wiele

zastosowań praktycznych, w szczególności stanowią podstawowe narzędzie

współczesnego projektowanie geometrycznego (CAD).

Rozważmy układ

n

węzłów interpolacyjnych

0

0

1

1

1

1

{( ,

), ( ,

),...,(

,

)}

n

n

x y

x y

x

y





, gdzie

0

1

1

..

n

a

x

x

x

b





  



.

Splajnem kubicznym C = C(x) nazywamy funkcje określoną na przedziale [a,b] i taką, że:

1.

2

( )

([ , ])

C x

C

a b



, tj. jest ona ciągła wraz z pierwszą i drugą pochodną w [a,b].

2.

1

3

2

,3

,2

,1

,0

[

,

]

( ) :

( )

k

k

k

k

k

k

k

x x

C x

C x

a x

a x

a x

a













,

0,..,

2

k

n





, tj. wewnątrz każdego

podprzedziału funkcja ta jest pewnym wielomianem 3-ego stopnia.

3. Dla każdego węzła funkcja C(x) spełnia warunki interpolacji tj.

(

)

,

0,..,

1

k

k

C x

y

k

n







.

Z powyższego wynika, że wielomiany lokalne muszą spełniać warunki interpolacyjne

(

)

,

0,..,

1

k

k

k

C x

y

k

n







,

oraz warunki sklejenia zapewniające założoną regularność funkcji C, a mianowicie

1

1

1

(

)

(

)

(

)

(

)

(

)

(

)

k

k

k

k

k

k

k

k

k

k

k

k

C

x

C x

C

x

C x

C

x

C x











 







 







dla

1,..,

2

k

n





(tj. w węzłach wewnętrznych).

Zauważmy, że liczba postawionych warunków jest równa 4n - 6. Całkowita liczba nieznanych
współczynników lokalnych wielomianów jest natomiast równa 4(n - 1) = 4

n

- 4. Zatem,

problem wyznaczenia splajna kubicznego jest niedookreślony.

Potrzebujemy nałożyć dwa dodatkowe warunki tak, aby zagadnienie wyznaczenia
funkcji C miało jednoznaczne rozwiązanie.

Splajn kubiczny

x

0

x

1

x

2

x

k

x

k-1

x

k+1

x

n-2

x

n-1

y

0

y

1

y

2

y

k-1

y

k

y

k+1

y

n-2

y

n-1

C

1

(x)

C

0

(x)

C

k-1

(x)

C

k

(x)

C

n-2

(x)

y

x

Zwykle (ale nie zawsze) warunki dodatkowe mają formę warunków „brzegowych”
nałożonych na pierwszą lub drugą pochodną funkcji C, a mianowicie:

(

0

(

)

C x







lub

0

(

)

C x







) i (

1

(

)

n

C x









lub

1

(

)

n

C x









)

W powyższych warunkach liczby α, β, γ, i δ są oczywiście zadane.

Ważnym przypadkiem szczególnym jest tzw. splajn naturalny. Jest to taki splajn kubicznym
który spełnia warunki postaci

0

(

)

0

C x





,

1

(

)

0

n

C x







Splajn naturalny posiada interesującą własność. Okazuje się, że spośród wszystkich funkcji
ciągłych wraz z dwiema pierwszymi pochodnymi i interpolujących zadany układ węzłów
splajn naturalny ma najmniejszą wartość całki z kwadratu drugiej pochodnej na
przedziale interpolacji [a,b], tj.

2

( )

min

b

a

C x

dx















Dokładniej, mamy następujące

TWIERDZENIE: Niech

2

0 ,

1

( [

] )

n

f

C

x

x





będzie dowolna. Załóżmy, że

( )

0

C a





i

( )

0

C b





lub

( )

( )

C a

f a







i

( )

( )

C b

f b







. Wówczas

2

2

[

( )]

[

( )]

b

b

a

a

C x

dx

f

x

dx











Dowód:

0

0

( )[

( )

( )]

( )[

( )

( )]

( )[

( )

( )]

( )[

( )

( )]

( )[

( )

( )]

( )[

( )

( )]

przez

czesci

z zal

b

b

x b

x a

ozenia

z zalozenia

a

a

b

a

C x

f

x

C x dx

C x

f x

C x

C

x

f x

C x dx

C b

f b

C b

C a

f a

C a

C

x

f x

C x dx














































































1

1

,3

1

1

1

,3

0

0

6

0

(

)

(

),

0

1

,3

,1 ..

0

, ,

( ) [

( )

( )]

6

[

( )

( )]

6

[ ( )

( )]

0

k

k

k

k

k

k

k

k

k

a

con

n

n

x

x

k

k

x

x

k

k

n

x

st

bo C x

f x

x

k

n

x

k

x

k

C x

f x

C x dx

a

f x

C x dx

a

f x

C x












































 





 















Otrzymaliśmy równość

2

( )

( )

[

( )]

b

b

a

a

C x f

x dx

C x

dx













.

Dalej mamy

2

2

2

2

2

2

[

( )]

0

[

( )

( )]

[

( )]

2

( )

( )

[

( )]

[

( )]

[

( )]

b

a

b

b

b

b

a

a

a

a

b

b

a

a

C x

dx

f

x

C x

dx

f

x

dx

f

x C x dx

C x

dx

f

x

dx

C x

dx



















































czyli

2

2

[

( )]

[

( )]

b

b

a

a

f

x

dx

C x

dx











. Koniec dowodu.

Pozostaje kwestia jak w efektywny sposób wyznaczyć (skonstruować) splajn kubiczny dla
zadanego układu węzłów. W teorii, moglibyśmy zbudować i rozwiązać układ równań
(liniowych) dla nieznanych współczynników lokalnych wielomianów C

k

(k = 0,1,..,n-2).

Układ taki zawierałby 4n-4 równań, a macierz współczynników miałaby dość „paskudną”
strukturę.

Okazuje się (jak zwykle?), że istnieje alternatywna metoda inteligentna!

Zacznijmy od spostrzeżenia, że druga pochodna poszukiwanej funkcji sklejanej C jest funkcją
kawałkami liniową (mówimy też – jest splajnem liniowym). Można ją zapisać następująco:

,

1

[

]

1

1

1

1

( )

( )

(

)

(

)

k

k

k

x x

k

k

x

k

x

k

k

k

k

k

C x

C x

x

x

x

x

C x

C x

x

x

x

x



































lub

1

1

( )

k

k

k

k

k

k

k

x

x

x

x

C x

m

m

h

h















gdzie

1

(

) ,

k

k

k

k

k

m

C x

h

x

x











Całkując dwukrotnie powyższą formułę otrzymamy postać lokalnego wielomianu

k

C

3

3

1

1

1

( )

(

)

(

)

(

)

(

)

6

6

k

k

k

k

k

k

k

k

k

k

k

m

m

C x

x

x

x

x

p x

x

q x

x

h

h























przy czym symbole

p

k

i

q

k

oznaczają stałe całkowania. Czytelnik będzie uprzejmy upewnić

się, że druga pochodna funkcji

k

C ma istotnie odpowiednią postać.

x

0

x

1

x

2

x

k

x

k-1

x

k+1

x

n-2

x

n-1

y

0

y

1

y

2

y

k-1

y

k

y

k+1

y

n-2

y

n-1

y

x

L

0

(x)

L

1

(x)

L

k-1

(x)

L

k

(x)

L

n-2

(x)

Póki co, stałe całkowania były dowolne. Teraz dobierzemy je tak, aby spełnić warunki
interpolacji

1

1

(

)

,

(

)

k

k

k

k

k

k

C x

y

C x

y









Otrzymujemy wartości stałych

p

k

i

q

k

2

1

1

6

6

k

k

k

k

k

k

k

k

k

k

y

y

m h

p h

p

m h

h











2

1

1

1

6

6

1

1

1

k

k

k

k

k

k

k

k

k

k

y

y

m h

q h

q

m h

h



















i w konsekwencji ostateczna postać lokalnego wielomianu wyraża się formułą

3

3

1

1

1

1

6

6

1

1

1

( )

(

)

(

)

(

)

(

)

6

6

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

m

m

y

y

C x

x

x

x

x

m h

x

x

m h

x

x

h

h

h

h































 





















przy czym indeks k przyjmuje wartości od 0 do n-2.

Pozostało obliczyć wartości drugiej pochodnej splajnu w węzłach, czyli wielkości

0

1

1

,

,...,

n

m m

m



.

Zauważmy, że nie wykorzystaliśmy jeszcze warunku „dopasowania” pierwszej pochodnej
sąsiadujących wielomianów lokalnych. Różniczkując otrzymaną wyżej formułę dla C

k

mamy

2

2

1

1

1
6

1

1

( )

(

)

(

)

(

)

2

2

k

k

k

k

k

k

k

k

k

k

k

k

k

m

m

y

y

C x

x

x

x

x

m

m h

h

h

h













 













Rozważmy węzeł

x

k

. Z powyższego wzoru wynika dla

k

x

x



, że

1

1

1

1

1

3

6

3

6

1

1

(

)

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

d

y

y

C x

m h

m h

m h

m h

d

h











 





 





Wartość pierwszej pochodnej wielomianu C

k-1

w węźle

x

k

otrzymamy następująco: najpierw

w ogólnej formule dla 1-szej pochodnej wielomianu C

k

podmienimy formalnie k na k-1, a

następnie podstawimy

k

x

x



. Oto rezultat (sprawdzić!)

1

1

1

1

1

1

3

6

3

6

1

1

1

1

1

1

1

1

1

(

)

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

d

y

y

C

x

m h

m h

m h

m h

d

h







































Z warunków ciągłości 1-szej pochodnej w węzłach mamy

1

(

)

(

)

k

k

k

k

C

x

C x









, co po

podstawieniu otrzymanych wzorów i prostych przekształceniach prowadzi do następującego
układu równań dla wielkości

0

1

1

,

,...,

n

m m

m



1

1

1

1

2(

)

,

1,2,..,

2

k

k

k

k

k

k

k

k

h m

h

h m

h m

u

k

n





















gdzie oznaczyliśmy

1

1

1

1

6(

)

6

k

k

k

k

k

k

k

k

k

y

y

y

y

u

d

d

h

h

































.

Wiemy już, że do wyznaczenia funkcji sklejanej C potrzebujemy dwóch dodatkowych
warunków. Jeśli zdecydujemy się na określenie wartości pierwszej pochodnej w węzłach
skrajnych to warunki te przyjmą postać

1

1

3

6

0

0

0

0

1

0

0

0

0

1

0

(

)

2

6(

)

C x

h m

h m

d

h m

h m

d







 















1

1

3

6

1

2

1

2

2

2

2

2

2

1

2

(

)

2

6(

)

n

n

n

n

n

n

n

n

n

n

n

C x

h

m

h

m

d

h

m

h

m

d













































Jeśli natomiast określimy wartości brzegowe drugiej pochodnej to dodatkowe równania są
bardzo proste, a mianowicie

0

m





i

1

n

m







. W szczególności, jeśli

0

 

 

to

otrzymamy splajn kubiczny naturalny.

Podsumowując, kompletny układ równań liniowych dla wartości drugiej pochodnej splajna
kubicznego w węzłach może być zapisany następująco:

0

0

0

0

1

0

1

1

1

1

1

2

2

2

1

2

2

6(

)

,

0

2(

)

,

1,..,

2

2

6(

) ,

1

k

k

k

k

k

k

k

k

n

n

n

n

n

n

m

h m

h m

d

k

h

m

h

h m

h m

u

k

n

m

h

m

h

m

d

k

n

lub

lub



































































 

W notacji macierzowo-wektorowej:



Tm r

. Zauważmy, że macierz

T

jest trójdiagonalna.

Niezerowe elementy tej macierzy można zapisać jako elementy trzech wektorów a,b and c.

0

0

1

1

2

2

2

1

1

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

k

k

k

n

n

n

n

n

c

b

a

c

b

a

c

b

a

c

b

a

c









































































T

Wartości zapisane w tych wektorach zależą od przyjętego wariantu „warunków
brzegowych” i przedstawiają się następująco:

0

1

1

1

2

0

,

1,..,

2

0

(

)

k

k

n

n

n

nie uzyw

a

a

h

k

n

a

lub

h

y

a

an













  











 



0

0

0

1

0

,

1,..,

2

(

)

0

k

k

n

l

b

b

h

b

h

nie uzuwan

u

b

y

k

n

b





 



  













0

0

0

1

1

1

2

1

2

2(

) ,

1,..,

2

1

2

k

k

k

n

n

n

lub

l

c

c

h

c

h

h

k

n

c

c

ub

h















 

















Do efektywnego rozwiązania układu liniowego z macierzą trójdiagonalną stosujemy
specjalny wariant metody eliminacji Gaussa zwany algorytmem przeganiania (Thomasa).
Przedstawimy ten algorytm zakładając, że rozwiązywany układ równań ma postać

0

0

0

1

0

1

1

1

2

1

1

1

,

1,..,

2

k

k

k

k

k

k

k

n

n

n

n

n

k

n

c m

b m

r

a m

c m

b m

r

a

m

c

m

r















































Metoda przeganiania (algorytm Thomasa)

Rozważmy dwa pierwsze równania układu

0

0

0

1

0

1

0

1

1

1

2

c m

b m

r

a m

c m

b m

r

















i załóżmy, że

0

0

c



. Najpierw „normalizujemy” pierwsze równanie dzieląc je przez c

0

,

następnie eliminujemy m

0

metodą przeciwnych współczynników

0

0

0

0

0

0

0

0

1

0

1

1

0

1

1

2

,

/

,

/

b c

r c

a m

c m

b

m

m

m

r





























i w końcu sprowadzamy zmodyfikowane drugie równanie do postaci „znormalizowanej”

1

0 1

1

1

2

1

1

0

1

0 1

1

1

0

1

1

1

1

0

1

1

2

1

1

0 1

1

(

)

/ : (

)

,

,

c

a m

b m

r

a

c

a

r

a

b

m

a

a

m

c

c



























 

















W trakcie obliczeń pojawiają się dwie pomocnicze wielkości

1



and

1



.

Zauważmy, że nowy zredukowany układ równań z niewiadomymi

1

1

,...,

n

m

m



wygląda „tak

samo” jak oryginalny, tj. ma strukturę 3-diagonalną i pierwsze z równań ma postać 3-
diagonalną i pierwsze równanie ma postać znormalizowaną. Układ ten jest zatem gotowy do
kontunuowania procedury eliminacji kolejnych niewiadomych. Po k krokach procedury
eliminacyjnej dochodzimy do etapu, w którym układ zwiera niewiadome o numerach od k do
n-1. Następny krok polega na eliminacji niewiadomej

m

k

metodą przeciwnych

współczynników. W szczegółach wygląda to następująco

1

1

1

1

1

1

2

1

1

1

1

1

2

1

1

1

1

1

1

1

1

1

1

1

2

1

1

1

/

(

)

,

,

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

k

m

m

a

a

m

c m

b m

r

c

a

m

b m

r

b

r

a

m

m

c

b

a

c

a

























































































































Podczas obliczeń pojawia się kolejna para wielkości pomocniczych

1

k





and

1

k





.

Ostatecznie, po n-1 krokach eliminacji proces osiąga ostatnie równanie układu. Ostatni krok
eliminacji przebiega następująco

2

2

1

2

1

1

2

1

1

1

1

2

1

1

1

1

2

1

1

2

1

1

2

1

/

(

)

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

m

m

a

a

r

b

m

c

m

c

m

r

c

a

m

r

a

b



























































































Zauważmy, że ostatnie równanie układu zawiera jedynie dwie niewiadome (przedostatnią i
ostatnią), przez co wyniku eliminacji m

n-2

otrzymujemy równanie z tylko jedną niewiadomą

m

n-1

. Zauważmy również, że podczas procesu eliminacji otrzymaliśmy rekursywną formułę

wiążącą dwie niewiadome o kolejnych numerach. Ogólna postać tej formuły wynika z
pierwszego równania (znormalizowanego) zredukowanego układu po k krokach eliminacji, a
mianowicie

1

,

2,

3,...,1,0

k

k

k

k

m

m

k

n

n

 







 



W ten sposób wszystkie wyeliminowane wcześniej niewiadome mogą być wyznaczone w
pętli chodzącej wspak.

Metodę przeganiania (algorytm Thomasa) można podsumować następująco:

0

0

0

0

0

0

1

1

1

1

1

1

1

1

1

2

1

1

1

2

1

1

1 (

)

/

;

/

;

0,..,

3

:

;

;

;

2 (

)

;

2,..,

!

:

!

0

k

k

k

k

k

k

k

k

k

k

k

k

n

n

n

n

n

n

n

k

k

k

k

ETAP

sweep

up

b

c

r

c

for k

n

do

b

c

a

r

a

c

a

end

ETAP

sweep

down

r

a

m

c

a

for k

n

ta petla c

do

m

m

hodzi wst

en

ecz



















































































 





;

d

UWAGA:

Ponieważ metoda przeganiania jest pewnym wariantem metody eliminacji Gaussa (bez
wyboru elementu głównego – vide Wykład 7), to na macierz trójdiagonalną należy nałożyć
pewne ograniczenia gwarantujące powodzenie obliczeń. Chodzi przede wszystkim o
gwarancję, że wszystkie operacje dzielenia będą wykonalne (nie wystąpi dzielenie przez
zero).

Można pokazać, że warunki wystarczające dla powodzenia przebiegu obliczeń metodą
przeganiania mają następującą postać:

0

1

0

0

1

1

1)

0,

0 ,

0 ,

0 ,

1,..,

2

2)

,

1,..,

2

,

n

k

k

k

k

k

n

n

c

c

a

b

k

n

c

a

b

i

n

c

b

c

a































Przynajmniej jedna z ty

warunki

ch niero

diagonalnej domina

wnosci musi byc OS

cji

TRA!