Interpolacja funkcjami sklejanymi
Kazimierz Jakubczyk
19 lutego 2008
Przykład Rungego
Jedyną możliwością uzyskania lepszego przybliżenia w interpolacji wielomiano-
wej jest zwiększenie stopnia wielomianu interpolacyjnego, czyli liczby węzłów.
Nie jest to jednak możliwe, gdy funkcja interpolowana jest nieznana. Nawet gdy
funkcja interpolowana jest znana i bardzo regularna, ciąg wielomianów interpo-
lacyjnych może nie być do niej zbieżny. Przykładem takiej funkcji jest
1
f(x) =
1 + x2
i ciąg wielomianów interpolacyjnych {Pn} opartych na węzłach równoodległych
z przedziału [-5, 5]:
10
xi = -5 + i (i = 0, 1, . . . , n)
n
(Runge) [2]. Pomimo że f jest nieskończenie różniczkowalna, to ciąg {Pn(x)}
jest zbieżny do f(x) tylko dla |x| < 3, 63 . . . i rozbieżny dla |x| > 3, 63 . . .
Funkcje sklejane
Klasyczne funkcje sklejane wywodzą się z praktyki inżynierskiej. Przez wiele lat
do kreślenia elementów konstrukcyjnych w przemyśle okrętowym i samochodo-
wym używano elastycznej listewki drewnianej nazywanej giętką (ang. spline).
Liniał taki, prowadzony przez zadane punkty za pomocą stosownych ciężarków,
ugina się wzdłuż grzywej  najgładszej . Inna odmiana tej techniki rysowania
krzywych gładkich polega na użyciu tzw. krzywików.
Liniał używany przez kreślarza można traktować jak cienką belkę ulegającą
odkształceniu pod wpływem działania sił zewnętrznych w danych punktach.
Siły sprężyste belki równoważą siły zginające, w wyniku czego belka przyjmuje
kształt krzywej  najgładszej  charakteryzującej się własnością minimalnej
krzywizny [1].
Pojęcie  gładkości krzywej uściślimy dla krzywych o równaniu y = f(x),
gdzie f " C2[a, b]. Krzywizna takiej krzywej w punkcie x wynosi
f (x)
º(x) =

3
1 + f (x)2
1
Niech będzie dany układ punktów x0, x1, . . . , xn dzielących przedział [a, b]
na n części, tzn.:
a = x0 < x1 < . . . < xn = b
Niech również będą dane wartości y0, y1, . . . , yn w punktach tego podziału.
Funkcja (krzywa) f " C2[a, b] interpolujÄ…ca punkty (xi, yi) (i = 0, 1, . . . , n)
jest  najgładsza w tym sensie, że spośród wszystkich funkcji (krzywych) kla-
sy C2[a, b] interpolujących te punkty daje najmniejszą wartość tzw. krzywizny
całkowitej [1]:

b
º(x)2dx
a
W praktyce trudno jest szukać rozwiązania zagadnienia interpolacji w posta-
ci funkcji  najgładszej z powodu skomplikowanego wyrażenia opisującego krzy-
wiznę. Nietrudno zauważyć, że dla małych wartości |f (x)| krzywizna w punkcie
x jest w przybliżeniu równa f (x). Uproszczenie to dostarcza interesującego roz-
wiązania w postaci funkcji sklejanych trzeciego stopnia, które wystarcza w wielu
zastosowaniach, nawet gdy wartości |f (x)| nie są małe. Kryterium minimalnej
krzywizny sprowadza się wówczas do minimalizacji całki [1, 2, 4]:

b
f (x)2dx
a
Definicja 1. FunkcjÄ™ rzeczywistÄ… s " Cm-1[a, b] nazywamy funkcjÄ… sklejanÄ…
lub splajnem (ang. spline) stopnia m z węzłami
a = x0 < x1 < . . . < xn = b
jeżeli w każdym z przedziałów [xi, xi+1] jest wielomianem stopnia co najwyżej m.
Funkcja sklejana stopnia 0 jest funkcją schodkową (przedziałami stałą),
a funkcja sklejana stopnia 1 jest łamaną. W praktyce najczęściej są używane
funkcje sklejane stopnia 3, nazywane też sześciennymi lub kubicznymi. W pew-
nych zadaniach stosuje się funkcje sklejane wyższych stopni, a także inne ich
odmiany, np. hiperboliczne funkcje sklejane i funkcje B-sklejane. Dla stopni nie-
parzystych m = 2k + 1 rozpatruje się często naturalne funkcje sklejane, które
poza przedziałem [a, b] są wielomianami stopnia co najwyżej k.
Obliczanie funkcji sklejanych
Szukamy takiej funkcji sklejanej s trzeciego stopnia, która w danych węzłach xi
ma dane wartości yi i która w każdym z przedziałów [xi, xi+1] jest wielomianem
stopnia nie wyższego niż 3. Niech
Mi = s (xi) (i = 0, 1, . . . , n)
Ponieważ s jest w przedziale [xi, xi+1] wielomianem stopnia co najwyżej trze-
ciego, funkcja s jest tam liniowa [1]:
Mi+1 - Mi
s (x) = (x - xi) + Mi
xi+1 - xi
2
Całkując dwukrotnie obie strony tej równości, otrzymujemy najpierw
Mi+1 - Mi
s (x) = (x - xi)2 + Mi(x - xi) + Ä…i
2(xi+1 - xi)
a potem
Mi+1 - Mi Mi
s(x) = (x - xi)3 + (x - xi)2 + Ä…i(x - xi) + yi
6(xi+1 - xi) 2
Stała całkowania yi wynika z warunku interpolacji s(xi) = yi. Z kolei, korzysta-
jąc z powyższego równania i warunku s(xi+1) = yi+1, otrzymujemy
Mi+1 - Mi Mi
yi+1 = (xi+1 - xi)3 + (xi+1 - xi)2 + Ä…i(xi+1 - xi) + yi
6(xi+1 - xi) 2
a stÄ…d
yi+1 - yi Mi+1 + 2Mi
Ä…i = - (xi+1 - xi)
xi+1 - xi 6
Funkcje s i s dają się więc opisać w przedziale [xi, xi+1] wzorami:
Mi+1 - Mi yi+1 - yi
s (x) = (x - xi)2 + Mi(x - xi) + +
2(xi+1 - xi) xi+1 - xi
Mi+1 + 2Mi
- (xi+1 - xi)
6
Mi+1 - Mi Mi
s(x) = (x - xi)3 + (x - xi)2 +
6(xi+1 - xi) 2

yi+1 - yi Mi+1 + 2Mi
+ - (xi+1 - xi) (x - xi) + yi
xi+1 - xi 6
w których należy wyznaczyć wielkości Mi (i = 0, 1, . . . , n). Funkcja s jest ciągła.
Porównując zatem dwa wyrażenia reprezentujące wartość s (xi-) = s (xi+)
odpowiadające przejściu od przedziału [xi-1, xi] do [xi, xi+1], dostajemy
µiMi-1 + 2Mi + iMi+1 = di (i = 1, 2, . . . , n - 1)
gdzie
xi+1 - xi
i = , µi = 1 - i
xi+1 - xi-1

6 yi+1 - yi yi - yi-1
di = -
xi+1 - xi-1 xi+1 - xi xi - xi-1
Otrzymaliśmy układ n - 1 równań liniowych z niewiadomymi M0, M1, . . . , Mn.
Ponieważ niewiadomych jest n + 1, rozwiązanie tego układu wymaga określenia
dodatkowych warunków, które zapiszemy w postaci dwóch dodatkowych równań
liniowych  jednego na początku, a drugiego na końcu układu:
2M0 + M10 = d0
µnMn-1 + 2Mn = dn
3
Wtedy pełny układ wygląda następująco:
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
2 0 0 M0 d0
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
µ1 2 1 M1 d1
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
µ2 2 2 M2 d2
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
× =
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
. .
. . .
. . . . .
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
. . .
. .
ïÅ‚ śł ïÅ‚ śł ïÅ‚ śł
ðÅ‚
µn-1 2 n-1 ûÅ‚ ðÅ‚ Mn-1 ûÅ‚ ðÅ‚ dn-1 ûÅ‚
0 µn 2 Mn dn
Powróćmy do określenia współczynników dwóch dodatkowych równań. Naj-
prościej przyjąć M0 = Mn = 0, tj.
0 = 0, d0 = 0
µn = 0, dn = 0
Daje to normalną funkcję sklejaną trzeciego stopnia, która na lewo i na prawo od
przedziału [a, b] jest wielomianem pierwszego stopnia. Odpowiada to liniałowi
elastycznemu, którego końce są swobodne.
Inny sensowny wybór współczynników dodatkowych równań odpowiada li-
niałowi o usztywnionych końcach, np. belce o końcach wmurowanych w rurach,
w których mogą się one przesuwać. Inaczej mówiąc, warunki dające naturalną
funkcję sklejaną zastępujemy warunkami brzegowymi postaci

s (x0+) = y0, s (xn) = yn

gdzie y0 i yn są danymi nachyleniami końców liniału. Korzystając z równania
przedstawiającego postać funkcji s , łatwo uzyskać współczynniki obu równań:

6 y1 - y0
0 = 1, d0 = - y0
x1 - x0 x1 - x0

6 yn-yn-1

µn = 1, dn = yn-
xn-xn-1 xn-xn-1
Oczywiście, można określać warunki mieszane, przyjmując jeden koniec swo-
bodny, a drugi usztywniony. Rozpatruje się również okresowe funkcje sklejane,
które w końcach przedziału [a, b] przyjmują takie same wartości wraz z pierw-
szymi pochodnymi.
Macierz układu równań wiążącego wielkości M0, M1, . . . , Mn jest trójdiago-
nalna, a ponieważ
µi + i = 1 (i = 1, 2, . . . , n - 1)
oraz 0 i µn sÄ… równe 0 lub 1, główna przekÄ…tna tej macierzy jest dominujÄ…ca.
Dowodzi się, że macierz o tej własności jest nieosobliwa, układ ma więc do-
kładnie jedno rozwiązanie. Można go wyznaczyć, stosując uproszczony wariant
metody Gaussa. Mianowicie, eliminując jedynie elementy leżące tuż pod główną
przekÄ…tnÄ…, otrzymujemy
Mi = qiMi+1 + ui (i = 0, 1, . . . , n - 1)
Mn = un
4
gdzie
0 d0
q0 = - , u0 =
2 2
i di - µiui-1
qi = - , ui = (i = 1, 2, . . . , n - 1)
µiqi-1 + 2 µiqi-1 + 2
dn - µnun-1
un =
µnqn-1 + 2
Ostatnie równanie daje bezpośrednio wartość współczynnika Mn. Wartości po-
zostałych obliczamy w kolejności Mn-1, Mn-2, . . . , M0. Algorytm wyznaczania
wartości M0, M1, . . . , Mn można zaprogramować następująco:
if lewy_koniec_swobodny then
begin
q[0] := 0;
u[0] := 0;
end else
begin
q[0] := -0.5;
u[0] := 3*((y[1]-y[0])/(x[1]-x[0])-yp0)/(x[1]-x[0]);
end;
for i := 1 to n-1 do
begin
lambda := (x[i+1]-x[i])/(x[i+1]-x[i-1]);
mi := 1-lambda;
d := 6*((y[i+1]-y[i])/(x[i+1]-x[i])-
(y[i]-y[i-1])/(x[i]-x[i-1]))/(x[i+1]-x[i-1]));
q[i] := -lambda/(mi*q[i-1]+2);
u[i] := (d-mi*u[i-1])/(mi*q[i-1]+2);
end;
if prawy_koniec_swobodny then M[n] := 0 else
begin
d := 6*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]))/(x[n]-x[n-1]);
M[n] := (d-u[n-1])/(q[n-1]+2);
end;
for i := n-1 downto 0 do
M[i] := q[i]*M[i+1]+u[i];
W algorytmie przyjęto, że dwie zmienne logiczne lewy koniec swobodny i pra-
wy koniec swobodny zawierają informację o warunkach brzegowych: wartość true
oznacza, że stosowny koniec liniału elastycznego jest swobodny, a false, że jest

usztywniony. Ponadto zmienne yp0 i ypn mają  odpowiednio  wartości y0

i yn, gdy końce liniału są usztywnione.
Aby znalez
ć wartość w = s(t) funkcji sklejanej s w punkcie t " [a, b], trze-
ba najpierw ustalić, do którego z przedziałów [xi, xi+1] punkt t należy. W po-
niższym algorytmie przeglądanie to odbywa się od strony prawej do lewej. Po
znalezieniu odpowiedniego przedziału obliczana jest według schematu Hornera
wartość właściwego wielomianu trzeciego stopnia:
5
i := n;
repeat
i := i-1;
dt := t-x[i];
until (i = 0) or (dt >= 0);
dx := x[i+1]-x[i];
w := (((M[i+1]-M[i])/(6*dx)*dt+
M[i]/2)*dt+
((y[i+1]-y[i])/dx-(M[i+1]+2*M[i])/6*dx))*dt+
y[i];
Algorytm można łatwo rozbudować tak, by uwzględnione zostały również
przedziały (-", x0) i (xn, "). Należy wówczas przyjąć, że s jest w tych prze-
działach wielomianem stopnia co najwyżej pierwszego, czyli że jest naturalną
funkcjÄ… sklejanÄ….
Jakość interpolacji
Istnienie i jednoznaczność funkcji sklejanej trzeciego stopnia interpolującej da-
ne punkty wynika z faktu, że macierz rozpatrywanego wyżej układu równań
liniowych jest nieosobliwa. InterpolujÄ…ca funkcja sklejana jest w pewnym sensie
najgładszą funkcją interpolującą [3].
Twierdzenie 1. Jeżeli f " C2[a, b], a = x0 < x1 < . . . < xn = b i s jest
funkcją sklejaną trzeciego stopnia interpolującą f w węzłach xi (i = 0, 1, . . . , n),
to

b b
[s (x)]2 dx [f (x)]2 dx
a a
Dowód. Niech q a" f - s. Mamy:

b b b b
[f (x)]2 dx = [s (x)]2 dx + [g (x)]2 dx + 2 s (x) g (x) dx
a a a a
Wystarczy pokazać, że ostatnia całka po prawej stronie znika. Całkując przez
części, otrzymujemy

n
b xi

s (x)g (x) dx = s (x)g (x) dx =
a xi-1
i=1


n
xi

= s (xi)g (xi) -s (xi-1)g (xi-1) - s (x)g (x) dx =
xi-1
i=1

n
xi

= s (b)g (b) - s (a)g (a) - ci g (x) dx =
xi-1
i=1
n

= s (b)g (b) - s (a)g (a) - ci [g(xi) - g(xi-1)]
i=1
gdzie ci są stałymi takimi, że ci = s (x) dla x " [xi-1, xi] (funkcja s jest
przedziałami stała). Z definicji funkcji g wynika, że g(xi) = 0. Zatem

b
s (x)g (x) dx = s (b)g (b) - s (a)g (a)
a
6
Wyrażenie po prawej stronie równości jest równe zero, gdy warunki brzegowe
określają naturalną funkcję sklejaną (swobodne końce):
s (a) = s (b) = 0
Jeżeli warunki określające naturalną funkcję sklejaną zastąpimy warunkami
alternatywnymi (usztywnione końce):

s (a) = f (a) = y0, s (b) = f (b) = yn
to wyrażenie to będzie również miało wartość zero:
s (b)g (b) - s (a)g (a) = s (b) [f (b) - s (b)] - s (a) [f (a) - s (a)] = 0
Wbrew zachowaniu się wielomianów interpolacyjnych (por. przykład Run-
gego [2]) funkcje sklejane są zbieżne do funkcji, którą interpolują, jeżeli tylko
wybieramy coraz drobniejsze podziały rozpatrywanego przedziału [a, b]. Wyni-
ka to z następującego twierdzenia dotyczącego oszacowania błędu interpolacji
funkcjami sklejanymi trzeciego stopnia [2] (bez dowodu).
Twierdzenie 2. Jeżeli f " C2[a, b], a = x0 < x1 < . . . < xn = b i s jest
funkcją sklejaną trzeciego stopnia interpolującą f w węzłach xi (i = 0, 1, . . . , n),
to dla każdego x " [a, b] zachodzi nierówność
|f(x) - s(x)| 5M max (xi - xi-1)2
1 i n
gdzie
M = max |f (¾)|
a ¾ b
Literatura
[1] Ahlberg J.H., Nilson E.N., Walsh J.L.: The Theory of Splines and Their
Applications, Academic Press, New York and London 1967.
[2] Jankowscy J. i M.: Przegląd metod i algorytmów numerycznych, część 1,
WNT, Warszawa 1981.
[3] Kincaid D., Cheney W.: Analiza numeryczna, WNT, Warszawa 2006.
[4] Stoer J.: Wstęp do metod numerycznych, tom 1, PWN, Warszawa 1979.
7