interpolacji

Wst ˛ep do metod numerycznych

7. Interpolacja

P. F. Góra

http://th-www.if.uj.edu.pl/zfs/gora/

Interpolacja

Dana jest funkcja w postaci stabelaryzowanej

x_i x₁ x₂ x₃ . . . x_n

f_i = f (x_i) f₁ f₂ f₃ . . . f_n

Punkty x_i nazywamy w ˛ezłami interpolacji. Problem: chcemy znale´z´c łatwy

sposób na wyliczanie warto´sci funkcji pomi ˛edzy w ˛ezłami (interpolacja) lub poza obszarem obejmuj ˛acym w ˛ezły (ekstrapolacja).

Interpolacja odcinkami liniowa -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4

Interpolacja odcinkami liniowa na 11 w ˛ezłach. “Prawdziwa” funkcja interpolowana zaznaczona jest lini ˛a

kropkowan ˛a. Istniejeniesko ´nczenie wiele funkcji ci ˛agłych, które s ˛a sobie równe w sko ´nczonej liczbie

w ˛ezłów!

Mo˙zna to zrobi´c prowadz ˛ac łaman ˛a pomi ˛edzy poszczególnymi punktami w ˛ezłowymi (xi, fi),

ale jest to nieeleganckie i w pewnych przypadkach mo˙ze powodowa´c problemy, gdy˙z

Interpolacja wielomianowa

Tabela taka, jak na pocz ˛atku wykładu, wyznacza jednoznacznie wielomian stopnia n−1. Istotnie, rozpatrzmy wielomian

a_n−1xn−1 + a_n−2xn−2 + · · · + a₁x + a₀ . (1) Je´sli do wielomianu (1) podstawimy za x kolejno x₁, x₂,. . . ,xn, przyjmuj ˛ac,

˙ze warto´s´c wielomianu w tych punktach wynosi odpowiednio f₁, f₂,. . . ,fn,

otrzymamy a_n−1xn−1₁ + a_n−2xn−2₁ + · · · + a₁x₁ + a₀ = f₁ a_n−1xn−1₂ + a_n−2xn−2₂ + · · · + a₁x₂ + a₀ = f₂ · · · · a_n−1xn−1_n + a_n−2xn−2_n + · · · + a₁x_n + a₀ = f_n (2)

Układ (2) zapisany w postaci macierzowej ma posta´c       xn−1₁ xn−2₁ · · · x₁ 1 xn−1₂ xn−2₂ · · · x₂ 1 · · · · xn−1_n xn−2_n · · · x_n 1            a_n−1 a_n−2 ... a₀      =      f₁ f₂ ... f_n      (3)

Rozwi ˛azaniem układu równa ´n (3) s ˛a współczynniki wielomianu (1). Ma-cierz układu (3) nosi nazw ˛e maMa-cierzy Vandermonde’a.

Mo˙zna pokaza´c∗, ˙ze wyznacznik macierzy Vandermonde’a ma posta´c

Y

16i<j6n

(x_i − x_j) , (4)

∗_{Dzi ˛ekuj ˛e internetowemu komentatorowi juventusik plus za wskazanie bł ˛edu w}

a zatem wyznacznik macierzy Vandermonde’a jest ró˙zny od zera, je˙zeli ˙zadne punkty x₁, x₂, . . . , x_n nie pokrywaj ˛a si ˛e. Tym samym problem zna-lezienia współczynników wielomianu (1) ma jednoznaczne rozwi ˛azanie. Gdyby´smy chcieli rozwi ˛azywa´c układ równa ´n (3) za pomoc ˛a ju˙z pozna-nych metod, mogłoby si ˛e wydawa´c, ˙ze koszt interpolacji wielomianowej wynosi O(n3). W rzeczywisto´sci układ ten mo˙zna, korzystaj ˛ac z syme-trii macierzy Vandermode’a, rozwi ˛aza´c w czasie O(n2)_{. Przekonany si ˛e} o tym konstruuj ˛ac interpolacj ˛e wielomianow ˛a w całkiem inny sposób.

Wzór interpolacyjny Lagrange’a

Zamiast szuka´c rozwi ˛azania równania (3), postulujemy, ˙ze poszukiwany wzór interpolacyjny ma posta´c f (x) = n X j=1 l_j(x) f_j + E(x) , (5a) gdzie l_j(x) = (x − x1) . . . (x − xj−1)(x − xj+1) . . . (x − xn) (x_j − x₁) . . . (x_j − x_j−1)(x_j − x_j+1) . . . (x_j − xn) . (5b)

E(x) w (5a) jest nazwywane reszt ˛a lub bł ˛edem interpolacji. Zauwa˙zmy, ˙ze l_j(x) jest wielomianem stopnia n−1 oraz

l_j(x_k) = δ_jk . (6)

Je˙zeli f (x) jest wielomianem stopnia co najwy˙zej n−1, E(x) znika to˙zsamo´sciowo. Mó-wimy, ˙ze interpolacja (5a) ma dokładno´s´c n−1.

Przykłady -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4

Interpolacja pewnej funkcji (oznaczonej lini ˛a kropkowan ˛a) oparta na 7 (lewy panel) i 9 (prawy panel) w ˛ezłach. W tym wypadku zwi ˛ekszanie

Przykłady -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4

Interpolacja tej samej funkcji, co poprzednio, oparta na 11 i 13 w ˛ezłach. Dalsze zwi ˛ekszanie liczby w ˛ezłów prowadzi do pogorszenia jako´sci interpolacji, zwłaszcza w pobli˙zu kra ´nców przedziału zawieraj ˛acego w ˛ezły.

Oscylacje Rungego -2 0 2 4 6 8 10 12 14 -4 -3 -2 -1 0 1 2 3 4 oscylacje Rugego

Wielomiany wysokiego stopnia s ˛a “sztywne”. Je˙zeli narzuci´c im warunek, ˙ze maj ˛a

prze-chodzi´c przez ustalone z góry punkty, mog ˛a to kompensowa´c znacznymi wahaniami

po-mi ˛edzy w ˛ezłapo-mi. Zjawisko to nazywa si ˛e oscylacjapo-mi Rungego i oznacza, ˙ze

Dalsze przykłady -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4

Mo˙ze si ˛e zdarzy´c, ˙ze dalsze powi ˛ekszanie liczby w ˛ezłów poprawi jako´s´c interpolacji (interpolacja oparta na 31 w ˛ezłach, lewy panel), jednak w ko ´ncu zbyt du˙za liczba w ˛ezłów doprowadza do katastrofy (interpolacja

oparta na 41 w ˛ezłach, prawy panel — oscylacje Rungego osi ˛agaj ˛a amplitud ˛e ∼ 200).

Bł ˛ad interpolacji

Wielomian (5b) wygodnie jest niekiedy zapisywa´c w postaci l_j(x) = pn(x) (x − x_j)p0_n(x_j) (7a) gdzie p_n(x) = n Y i=1 (x − x_i) , p0_n(x_j) = dpn dx     x=x_j (7b) Oznaczmy y(x) = Pn j=1

l_j(x) f_j (jest to wielomianowa cz ˛e´s´c wzoru inter-polacyjnego Lagrange’a (5a), bez reszty).

Rozpatrzmy funkcj ˛e

F (z) = f (z) − y(z) − [f (x) − y(x)] pn(z)

p_n(x) . (8)

F (z), jako funkcja zmiennej z, ma n+1 miejsc zerowych: x₁, x₂, . . . , x_n oraz x. Zakładaj ˛ac, ˙ze funkcja f (·) jest dostatecznie gładka, stosujemy do funkcji F (z) n-krotnie twierdzenie Rolle’a i stwierdzamy, ˙ze pochodna

F(n)(z) = f(n)(z) − y(n)(z) − [f (x) − y(x)] n!

p_n(x) (9)

ma co najmniej jedno miejsce zerowe w najmniejszym przedziale domkni ˛e-tym zawieraj ˛acym punkty x₁, x₂, . . . , x_n oraz x. Oznaczmy to miejsce ze-rowe przez ξ. Zauwa˙zmy, ˙ze y(n)(z) = 0, gdy˙z y(z) jest wielomianem stopnia n−1. Ostatecznie otrzymujemy

0 = F(n)(ξ) = f(n)(ξ) − [f (x) − y(x)] | {z } E(x) n! p_n(x) (10) czyli E(x) = pn(x) n! f (n)_{(ξ) . (11)}

ξ jest pewnym punktem wewn ˛etrzmym przedziału zawieraj ˛acego w ˛ezły i x (to ostatnie jest wa˙zne w wypadku ekstrapolacji). Nie wiemy, którym punk-tem, zatem dla bezpiecze ´nstwa nale˙załoby bra´c najwi ˛eksz ˛a (co do mo-dułu) warto´s´c f(n)(ξ). Trudno´s´c w szacowaniu bł ˛edu interpolacji polega na trudno´sci w szacowaniu wysokich pochodnych interpolowanych funk-cji. Praktyka pokazuje, ˙ze wysokie pochodne nawet “porz ˛adnych” funkcji (niewielomianowych) mog ˛a przybiera´c znaczne warto´sci.

Interpolacja Hermite’a

Je˙zeli znamy nie tylko warto´sci funkcji interpolowanej w w ˛ezłach, ale tak˙ze

warto´sci pochodnej,

x_i x₁ x₂ x₃ . . . x_n

f_i = f (x_i) f₁ f₂ f₃ . . . f_n f_i0 = f0(x_i) f₁0 f₂0 f₃0 . . . f_n0

narzuca to 2n warunków na wielomian interpolacyjny†. Mo˙zna skonstru-owa´c wówczas interpolacj ˛e wielomianow ˛a rz ˛edu 2n − 1, postaci

†_{Jest to zagadnienie o stosunkowo niewielkich zastosowaniach praktycznych, ale za to}

y(x) = n X i=1 h_i(x)f_i + n X i=1 ¯_hi(x)f0 i + E(x) , (12a) gdzie h_i(x) = 1 − 2(x − x_i)l_i0(x_i) l_i2(x) , (12b) ¯_h i(x) = (x − xi)l2i (x) . (12c)

l_j(x) oznacza to samo, co w interpolacji Lagrange’a, natomiast

E(x) = p 2 n(x) (2n)!f (2n)_{(ξ) ,} (12d) gdzie ξ jest pewnym punktem wewn ˛etrzym przedziału rozpi ˛etego na w ˛e-złach i warto´sci x. Wielomian interpolacyjny (12a) zgadza si ˛e z interpolo-wan ˛a funkcj ˛a oraz jej pochodn ˛a w w ˛ezłach.

Interpolacja za pomoc ˛a funkcji sklejanych

Interpolacja wielomianowa jest koncepcyjnie najprostszym sposobem in-terpolacji, prowadzi´c jednak mo˙ze, jak to pokazali´smy, do niepo˙z ˛adanych zachowa ´n. Najbardziej popularnym sposobem unikni ˛ecia oscylacji Run-gego, zwi ˛azanych z trudno´sci ˛a szacowania bł ˛edu interpolacji wielomiano-wej, jest interpolacja za pomoc ˛a funkcji sklejanych.

Funkcj ˛a sklejan ˛a rz ˛edu k, czyli “splajnem” (ang. spline), nazywam funkcj ˛e, która

1. lokalnie jest wielomianem rz ˛edu k,

2. jest (k−1)-krotnie ró˙zniczkowalna w w ˛ezłach (z czego wynika, ˙ze jej pochodne rz ˛edu k−2 i ni˙zszych s ˛a ci ˛agłe).

Najcz ˛e´sciej u˙zywa sie funkcji sklejanych rz ˛edu 3, czyli “splajnów kubicz-nych” (ang. cubic splines).

Cubic splines

Załó˙zmy‡, ˙ze oprócz warto´sci funkcji w w ˛ezłach, znamy tak˙ze warto´sci

drugiej pochodnej interpolowanej funkcji w w ˛ezłach. Mamy wi ˛ec tabelk˛e postaci:

x_i x₁ x₂ x₃ . . . x_n

f_i = f (x_i) f₁ f₂ f₃ . . . f_n f_i00 = f00(x_i) f₁00 f₂00 f₃00 . . . f_n00

W ka˙zdym przedziale [x_j, x_j+1], j = 1, 2, . . . , n − 1, konstruujemy wielo-mian trzeciego stopnia

y_j(x) = A f_j + B f_j+1 + C f_j00 + D f_j+100 , (13a) gdzie A = xj+1 − x x_j+1 − x_j, B = x − x_j x_j+1 − x_j, (13b) C = 1 6(A 3_−A)(x j+1−xj)2, D = 1 6(B 3_−B)(x j+1−xj)2.(13c)

Łatwo sprawdzi´c, ˙ze y_j(x_j) = f_j, y_j(x_j+1) = f_j+1, a poniewa˙z, co mo˙zna wykaza´c prostym rachunkiem,

d2y_j

dx2 = A f

00

j + B fj+100 , (14)

tak˙ze warto´sci drugiej pochodnej (13a) zgadzaj ˛a si ˛e z zadanymi warto-´sciami w w ˛ezłach.

Jest jednak pewien problem: w rzeczywisto´sci nie znamy warto´sci dru-gich pochodnych f_j00. Nie skorzystali´smy jednak jeszcze z wymogu ci ˛ a-gło´sci pierwszej pochodnej (13a) w w ˛ezłach. W tym celu ˙z ˛adamy, aby pochodna y_j(x) obliczana w prawym kra ´ncu przedziału równała si ˛e po-chodnej y_j+1(x) obliczanej w lewym kra ´ncu odpowiedniego przedziału. Gdy to zrobimy, otrzymamy równanie

x_j − x_j−1 6 f 00 j−1 + x_j+1 − x_j−1 3 f 00 j + x_j+1 − x_j 6 f 00 j+1 = fj+1 − fj x_j+1 − x_j − f_j − f_j−1 x_j − x_j−1 . (15)

Naturalny splajn kubiczny

Je˙zeli mamy n w ˛ezłów interpolacji, mamy n−2 wewn ˛etrznych punktów zszycia, w których mo˙zemy ˙z ˛ada´c ci ˛agło´sci pochodnej. W takim wypadku (15) stanowi układ n−2 równa ´n z n niewiadomymi. Trzeba poda´c jakie´s dodatkowe warunki.

Najcz ˛e´sciej przyjmuje si ˛e, ˙ze f₁00 = f_n00 = 0. Jest to wówczas tak zwany

naturalny splajn kubiczny. Je˙zeli sytuacja tego wymaga — lub je˙zeli mamy jakie´s przesłanki, aby tak zrobi´c — mo˙zemy narzuci´c inne warunki na dru-gie pochodne na brzegach lub na kombinacje liniowe drugich pochodnych.

Równoodległe w ˛ezły

Je˙zeli w ˛ezły interpolacji s ˛a równoodległe, x_j+1 − x_j = h, równanie (15) przybiera szczególnie prost ˛a posta´c:

       4 1 1 4 1 1 4 1 . . . . 1 4 1 1 4               f₂00 f₃00 f₄00 ... f_n−200 f_n−100        = 6 h2        f1 − 2f2 + f3 f2 − 2f3 + f4 f3 − 2f4 + f5 ... fn−3 − 2fn−2 + fn−1 fn−2 − 2fn−1 + fn        (16)

Macierz po lewej stronie tego równania posiada łatwy do znalezienia roz-kład Cholesky’ego. Z doroz-kładno´sci ˛a do czynników “6”, prawa strona zawiera drugie ilorazy ró˙znicowe interpolowanej funkcji.

Praktyczne zastosowanie naturalnych splajnów kubicznych przebiega dwu-etapowo:

1. Rozwi ˛azujemy układ równa ´n (15) (lub, je˙zeli mo˙zna, (16)) na n−2 nie-znanych wielko´sci nf_j00on−1

j=2. Poniewa˙z układ ten jest trójdiagonalny,

koszt obliczeniowy wynosi O(n). Krok ten wykonujemy tylko raz na pocz ˛atku całej procedury.

2. W celu znalezienia warto´sci pomi ˛edzy w ˛ezłami, wykorzystujemy rów-nanie (13a) tyle razy, ile warto´sci chcemy znale´z´c. W ka˙zdym prze-dziale [x_j, x_j+1] u˙zywamy odpowiedniego wielomianu y_j(x)!

Przykład -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4

Interpolacja za pomoc ˛a naturalnych splajnów kubicznych z 21 w ˛ezłami. Wynik

interpola-cji praktycznie pokrywa si ˛e z funkcj ˛a interpolowan ˛a (t ˛a sam ˛a, co w poprzednich

Interpolacja na płaszczy´znie — splajny bikubiczne

Przypu´s´cmy, ˙ze pewn ˛a funkcj ˛e dwu zmiennych, f (x, y), mamy stabelary-zowan ˛a w w ˛ezłach dwuwymiarowej siatki kwadratowej:

f₁₁=f (x₁, y₁) f₂₁=f (x₂, y₁) f₃₁=f (x₃, y₁) · · · f_n1=f (x_n, y₁) f₁₂=f (x₁, y₂) f₂₂=f (x₂, y₂) f₃₂=f (x₃, y₂) · · · f_n2=f (x_n, y₂) f₁₃=f (x₁, y₃) f₂₃=f (x₂, y₃) f₃₁=f (x₃, y₃) · · · f_n1=f (x_n, y₃) · · · · f_1n=f (x₁, y_n) f_2n=f (x₂, y_n) f_3n=f (x₃, y_n) · · · f_n1=f (x_n, y_n) (17) “Wiersze” tej siatki odpowiadaj ˛a ustalonym warto´sciom zmiennej y.

Przypu´s´cmy, ˙ze chcemy znale´z´c warto´s´c funkcji f (x∗, y∗)_{, gdzie (x}∗, y∗) jest wewn ˛etrznym punktem siatki. W tym celu post ˛epujemy jak nast ˛epuje:

1. Przeprowadzamy splajn wzdłu˙z ka˙zdego “wiersza”. W ka˙zdym wier-szu warto´s´c zmiennej y jest ustalona, wi ˛ec jest to za ka˙zdym razem zwykły splajn jednowymiarowy. Ka˙zdy splajn poci ˛aga koszt nume-ryczny rz ˛edu O(n), a zatem obliczenie splajnów wzdłu˙z wszystkich wierszy poci ˛aga koszt rz ˛edu O(n2)_.

2. Obliczamy warto´s´c ka˙zdego z powy˙zszych splajnów w punkcie x = x∗. W ten sposób dostajemy n warto´sci funkcji w punktach (x∗, y₁)_, (x∗, y₂), . . . , (x∗, y_n).

3. Przez powy˙zsze punkty przeprowadzmy splajn w kierunku y (przy usta-lonej warto´sci x = x∗) i wyliczamy warto´s´c tego splajnu w punkcie (x∗, y∗)_{. Wymaga to dodatkowych O(n) operacji, zatem cały koszt} jest zdominowany przez O(n2).

Interpolacja za pomoc ˛a funkcji wymiernych

“Sztywno´sci” interpolacji wielomianowej mo˙zna unikn ˛a´c interpoluj ˛ac za po-moc ˛a funkcji wymiernych, to znaczy ilorazów wielomianów:

r(x) = Pµ(x)

Q_ν(x) . (18)

Funkcje wymierne z łatwo´sci ˛a modeluj ˛a wi ˛eksze bogactwo zachowa ´n, ni˙z wielomiany. Zagadnienie interpolacji wymiernej jest opracowane od strony teoretycznej gorzej ni˙z interpolacji wielomianowej, a poniewa˙z problem in-terpolacji wymiernej nie ma jednoznacznego rozwi ˛azania, istnieje szereg konkurencyjnych podej´s´c.

W tym wykładzie skorzystam z opublikowanego w 2007 algorytmu Floatera i Hormanna (zobacz tak˙ze tutaj).

Algorytm Floatera i Hormanna

Niech x₀, x₁, . . . , x_n b ˛ed ˛a wzajemnie ró˙znymi punktami (w ˛ezłami inter-polacji) i niech f_j = f (x_j) b ˛ed ˛a stabelaryzowanymi warto´sciami pew-nej funkcji w w ˛ezłach. Wybieramy parametr interpolacji d, 0 6 d 6 n. Niech p_i(x) b ˛edzie wielomianem interpoluj ˛acym rozpi ˛etym na punktach x_i, . . . , x_i+d. Wówczas r(x) = n−d P i=0 λ_i(x) p_i(x) n−d P i=0 λ_i(x) , (19a) gdzie λ_i(x) = (−1) i (x − x_i)· · · (x − x_i+d) . (19b)

Mo˙zna pokaza´c, ˙ze r(x) nie ma biegunów na osi rzeczywistej oraz ˙ze mo˙zna go zapisa´c w nast ˛epuj ˛acej postaci barycentrycznej:

r(x) = n X k=0 w_k x − x_kfk n X k=0 w_k x − x_k (20a) w_k = X i∈J_k (−1)i i+d Y j=i,j6=k 1 x_k − x_j (20b) gdzie J_k = _{{i ∈ I : k − d 6 i 6 k}, I = {0, 1, . . . , n − d}.}

Zastanówmy si ˛e, czy (20) istotnie daje interpolacj ˛e, to znaczy czy zga-dza si ˛e z funkcj ˛a interpolowan ˛a w w ˛ezłach. Niech x → x_l, gdzie x_l jest którym´s w ˛ezłem, a wi ˛ec zerem mianownika którego´s z ułamków wyst ˛epu-j ˛acych w liczniku i mianowniku (20a). Wówczas w obu sumach dominowa´c b ˛edzie tylko człon z k = l, a zatem

r(x → x_l) → w_l x − x_lfl w_l x − x_l = f_l , (21)

Praktyczne zastosowanie algorytmu Floatera i Hormanna wygl ˛ada tak: • Sprawdzamy, czy x jest blisko (z dokładno´sci ˛a do bł ˛edu obci ˛ecia) w

˛e-zła x_k; je´sli tak, wynikiem jest stabelaryzowana warto´s´c funkcji f_k; • je˙zeli nie, obliczamy r(x) według wzoru (20a). Wagi w_k obliczamy

tylko raz, na pocz ˛atku całej procedury.

Jak dobra´c parametr d? Praktyka pokazuje, ˙ze w wi ˛ekszo´sci typowych przypadków wystarcza bra´c d = 3, aczkolwiek niekiedy potrzebne jest na-wet d = 8. Je˙zeli interpolowana funkcja jest dostatecznie gładka, bł ˛ad interpolacji nie przekracza O hd+1_max , gdzie hmax jest najwi ˛eksz ˛a

Je˙zeli w ˛ezły interpolacji s ˛a równoodległe, ∀i : x_i − x_i−1 = h, wyra˙zenia na wagi przyjmuj ˛a szczególnie prost ˛a posta´c:

w_k = (−1) k−d hd X i∈J_k 1 (k − i)!(i + d − k)! (22)

Poniewa˙z ostateczny wynik nie zmieni si ˛e, je´sli wszystkie wagi przemno-˙zymy przez t ˛e sam ˛a stał ˛a (obliczamy stosunek dwóch wyra˙ze ´n!), dla rów-noodległych w ˛ezłów interpolacji mo˙zemy wybra´c wagi całkowite postaci

w_k = (−1)k−d X i∈J_k  d k − i  . (23)

Przykład -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 3 4

Interpolacja funkcjami wymiernymi wg algorytmu Floatera i Hormanna z 21 równoodle-głymi w ˛ezłami i parametrem d = 3. Dokładno´s´c interpolacji jest efektywnie taka sama,

Inne rodzaje interpolacji

W specyficznych sytuacjach, głównie w kontekscie analizy sygnałów, sto-suje si ˛e te˙z zupełnie inne rodzje interpolacji: Interpolacj ˛e trygonometrycz-n ˛a i falkow ˛a (waveletow ˛a). Mianowicie, “sygnał” (ci ˛ag zdyskretyzowanych warto´sci w w ˛ezłach) rozkłada si ˛e na ci ˛ag funkcji trygonometrycznych lub falek, a nast ˛epnie warto´sci pomi ˛edzy w ˛ezłami (lub brakuj ˛ace warto´sci w w ˛e-złach!) znajduje si ˛e korzystaj ˛ac ze znalezionego rozkładu. Ten zespół za-gadnie ´n wykracza poza ramy niniejszego wykładu.

Ró˙zniczkowanie numeryczne

Zdarza si ˛e, ˙ze maj ˛ac tylko stabelaryzowane warto´sci funcji z jakim´s sko ´n-czonym krokiem — tak, jak w wypadku interpolacji — chcemy numerycznie wyznaczy´c warto´s´c pochodnej funkcji w w ˛ezłach. Jest to zagadnienie bar-dzo podatne na bł ˛edy i nale˙zy go unika´c, ale czasami nie ma wyj´scia. . . Pochodna zdefiniowana jest jako granica ilorazu ró˙znicowego:

f0(x) = lim

h→0

f (x + h) − f (x)

h . (24)

Na potrzeby numeryczne, maj ˛ac do dyspozycji tylko zdyskretyzowane war-to´sci funkcji, granic ˛e mo˙zemy zast ˛api´c wyra˙zeniem sko ´nczonym na (co najmniej) trzy ró˙zne sposoby (dla uproszczenia zakładamy, ˙ze odległo´s´c mi ˛edzy w ˛ezłami jest stała):

Iloraz ró˙znicowy “do przodu”:

f_j0 = fj+1 − fj

h (25a)

Wsteczny iloraz ró˙znicowy:

f_j0 = fj − fj−1

h (25b)

Symetryczny iloraz ró˙znicowy:

f_j0 = fj+1 − fj−1

Wszystkie te przybli˙zenia mog ˛a da´c ró˙zne wyniki. Co wi ˛ecej, niesyme-ryczne przybli˙zenia (25a),(25b) mog ˛a wprowadza´c pewien bł ˛ad systema-tyczny, zale˙zny od wypukło´sci (od drugiej pochodnej) analizowanej funkcji.

Przybli˙zenie symetryczne (25c) jest pod tym wzgl ˛edem najbezpieczniej-sze.

Jako´s´c uzyskanego przybli˙zenia pochodnej mocno zale˙zy od kroku inter-polacji: im wi ˛ekszy krok, tym przybli˙zenie pochodnej jest gorsze.

Przykład

Niech f (x) = 1₂x2. Wówczas poszczególne ilorazy ró˙znicowe daj ˛a

f_{do przodu}0 = x + 1₂h , (26a)

f_wsteczny0 = x − 1₂h , (26b)

f_symetryczny0 = x . (26c)

Symteryczny iloraz ró˙znicowy daje w tym wypadku warto´s´c dokładn ˛a, ale jest to przypadek, wynikaj ˛acy ze szczególnie prostej postaci funkcji (z tego,

Przykłady -1.5 -1 -0.5 0 0.5 1 1.5 -4 -3 -2 -1 0 1 2 3 4 do przodu wsteczny symetryczny -1.5 -1 -0.5 0 0.5 1 1.5 -4 -3 -2 -1 0 1 2 3 4 do przodu wsteczny symetryczny

Wyniki ró˙zniczkowania numerycznego funkcji w 9 (lewy panel) i 21 w ˛ezłach (prawy

panel) za pomoc ˛a wzorów (25a)-(25c). Linia kropkowana oznacza dokładn ˛a warto´s´c

pochodnej. Widoczne s ˛a systematyczne ró˙znice pomi ˛edzy przybli˙zeniami “do przodu”

Ró˙zniczkowanie splajnów

Najlepszym, a zarazem numerycznie tanim, sposobem ró˙zniczkowania nu-merycznego jest przeprowadzenie splajnu, a nast ˛epnie zró˙zniczkowanie go w w ˛ezłach. Korzystamy przy tym z wszystkich własno´sci splajnów, a wi ˛ec z semi-analitycznych wzorów i z pewno´sci, ˙ze pochodna jest ci ˛ a-gła w w ˛ezłach, co pozwala nam unikn ˛a´c niejednoznaczno´sci zwi ˛azanej ze stosowaniem wzorów (25).

Dla uproszczenia w dalszym ci ˛agu b ˛edziemy zakłada´c, ˙ze mamy do czy-nienia z naturalnymi splajnami kubicznymi z równoodległymi w ˛ezłami.

Ró˙zniczkuj ˛ac wyra˙zenie (13a), wła´sciwe dla przedziału [x_j, x_j+1], i obli-czaj ˛ac pochodn ˛a w prawym kra ´ncu przedziału, otrzymujemy

y_j+10 = 1 h(fj+1 − fj) + 1 6h(2f 00 j+1 + fj00) . (27)

Równie dobrze mogliby´smy jednak wzi ˛a´c wielomian dla przedziału [x_j+1, x_j+2]_{, zró˙zniczkowa´c go i obliczy´c pochodn ˛a w lewym kra ´ncu} prze-działu: y_j+10 = 1 h(fj+2 − fj+1) − 1 6h(2f 00 j+1 + fj+200 ) . (28)

Wiemy, ˙ze wyra˙zenia (27), (28) s ˛a sobie równe. Istotnie, porównuj ˛ac ich prawe strony, dostajemy wyra˙zenie

f_j00 + 4f_j+100 + f_j+200 = 6

h2(fj − 2fj+1 + fj+2) , (29) odpowiadaj ˛ace (j+1) wierszowi równania (16).

Wyra˙zenia (27), (28) s ˛a sobie równe w arytmetyce dokładnej. W prak-tycznych obliczeniach numerycznych mog ˛a wyst ˛api´c jakie´s ró˙znice. Aby je zminimalizowa´c, bierzemy ´sredni ˛a prawych stron(27), (28). Ostatecznie jako przybli˙zenie pochodnej dostajemy

y0(x_j) = fj+1 − fj−1 2h − 1 12h  f_j+100 − f_j−100  , (30) j = 2, 3, . . . , n − 1

Wyra˙zenie (30) ma posta´c symetrycznego ilorazu ró˙znicowego z popraw-kami wynikaj ˛acymi ze znajomo´sci przybli˙zenia drugiej pochodnej. Wiel-ko´sci f_k00 wyliczamy z równania (16).

Przykład -1.5 -1 -0.5 0 0.5 1 1.5 -4 -3 -2 -1 0 1 2 3 4 pochodna splajnu

Ró˙zniczkowanie numeryczne na podstawie naturalnego splajnu kubicznego z 21 równoodległymi w ˛ezłami. Linia kropkowana oznacza prawdziwy przebieg pochodnej