Przypadek węzłów wielokrotnych

odpo-wiednio na węzłach x1, . . . , xn i x0, . . . , xn−1. Stąd

bn = f(x1, . . . , xn) − f(x⁰, . . . , xn−1) xn− x0

= f(x0, x1, . . . , xn), co kończy dowód. 2

Różnicę dzieloną f(x0, x1, . . . , xn) można łatwo obliczyć na podsta-wie wartości f(xj), 0 ≤ j ≤ n, budując następującą tabelkę:

x0 f(x0)

x1 f(x1) f(x0, x1)

x2 f(x2) f(x1, x2) f(x0, x1, x2)

... ... ... ... ...

xn f(xn) f(xn−1, xn) f(xn−2, xn−1, xn) · · · f(x⁰, x1, . . . , xn).

Zauważmy przy tym, że “po drodze” obliczamy f(xi, xi+1, . . . , xj) dla wszystkich 0 ≤ i < j ≤ n, a więc w szczególności również interesu-jące nas różnice dzielone f(x0, x₁, . . . , x_j). Stąd i z Twierdzenia 6.1 na-tychmiast wynika algorytm obliczania współczynników bj wielomianu interpolacyjnego Newtona. Po wykonaniu następującego programu,

for j := 0 to n do b[j] := f (x[j]);

for j := 1 to n do for k := n downto j do

b[k] := (b[k] − b[k − 1])/(x[k] − x[k − j]), współczynniki bj zostaną umieszczone w tablicy b[j].

6.4 Przypadek węzłów wielokrotnych

Uogólnieniem rozpatrzonego zadania interpolacji jest zadanie interpo-lacji Hermite’a. Zakładamy, że oprócz (różnych) węzłów xj dane są również ich krotności nj, 0 ≤ j ≤ k, przy czym^Pkj=0nj = n + 1. Należy skonstruować wielomian w_f ∈ Πn taki, że

w_f⁽ⁱ⁾(xj) = f⁽ⁱ⁾(xj) dla 0 ≤ i ≤ n^j− 1, 0 ≤ j ≤ k.

Oczywiście zakładamy przy tym, że odpowiednie pochodne funkcji f istnieją.

72 ROZDZIAŁ 6. INTERPOLACJA WIELOMIANOWA Lemat 6.2 Zadanie interpolacji Hermite’a ma jednoznaczne rozwiąza-nie.

Dowód Istnienie i jednoznaczność rozwiązania można uzasadnić tak samo jak w przypadku węzłów jednokrotnych. Przedstawiając wielo-mian w dowolnej bazie otrzymujemy układ n + 1 równań z n + 1 nie-wiadomymi, który dla zerowej prawej strony ma jedynie rozwiązanie zerowe. Inaczej bowiem istniałby wielomian niezerowy stopnia nie więk-szego niż n, który miałby zera o łącznej krotności większej niż n. 2

Nas oczywiście interesuje konstrukcja wielomianu w_f. W tym celu ustawimy węzły xj w ciąg i zdefiniujemy uogólnioną bazę Newtona w Πn jako

p₀(x) = 1,

pj(x) = (x − ¯x⁰)(x − ¯x¹) · · · (x − ¯x^j−1), 1 ≤ j ≤ n.

Uogólnimy również pojęcie różnicy dzielonej na węzły powtarzające się kładąc mo-żemy łatwo obliczyć stosując schemat podobny do tego z przypadku węzłów jednokrotnych.

Twierdzenie 6.2 Współczynniki bj wielomianu interpolacyjnego Her-mite’a w bazie Newtona,

wf(·) = ^Xn

j=0

bjpj(·), dane są przez odpowiednie różnice dzielone, tzn.

bj = f(¯x0,¯x1, . . . ,¯xj), 0 ≤ j ≤ n.

6.4. PRZYPADEK WĘZŁÓW WIELOKROTNYCH 73 Dowód Dowód przeprowadzimy podobnie jak dla węzłów jednokrot-nych. Niech wi,j ∈ Πj−i oznacza wielomian interpolacyjny Hermite’a oparty na (być może powtarzających się) węzłach ¯xi,¯xi+1, . . . ,¯xj. To znaczy, wi,j interpoluje f w węzłach xs takich, że xs występuje w ciągu

¯xi, . . .¯xj, a jego krotność jest liczbą powtórzeń xs w tym ciągu.

Zauważmy najpierw, że dla ¯xi 6= ¯xj zachodzi znany nam już wzór (6.7),

wi,j(x) = (x − ¯xⁱ)wi+1,j(x) − (x − ¯x^j)wi,j−1(x)

¯xj − ¯xⁱ . (6.8)

Rzeczywiście, oznaczmy przez v(x) prawą stronę powyższej równości.

Dla k mniejszego od krotności danego węzła xsw ciągu ¯xi, . . .¯xj, mamy w^(k−1)_i+1,j(xs) = w_i,j−1^(k−1)(xs), a ponieważ

v^(k)(x) = k(w^(k−1)_i+1,j(x) − w^(k−1)i,j−1(x))

¯xj − ¯xⁱ

+(x − ¯xⁱ)w_i+1,j^(k) (x) − (x − ¯x^j)w_i,j−1^(k) (x)

¯xj− ¯xi

, to

v^(k)(xs) = (xs− ¯xⁱ)w_i+1,j^(k) (xs) − (x^s− ¯x^j)w_i,j−1^(k) (xs)

¯xj − ¯xⁱ .

Korzystając z tego wzoru sprawdzamy, że v spełnia odpowiednie wa-runki interpolacyjne, a stąd wi,j = v.

Dalej postępujemy indukcyjnie ze względu na n. Dla n = 0 mamy b₀ = f(x0). Dla n ≥ 1 wystarczy pokazać, że bⁿ = f(¯x0,¯x1, . . . ,¯xn). W tym celu rozpatrzymy dwa przypadki.

Jeśli ¯x0 = ¯xn to mamy jeden węzeł x0 o krotności n + 1. Wielomian interpolacyjny jest wtedy postaci

w_f(x) =

Xn j=0

f^(j)(x0)

j! (x − x⁰)^j, a stąd bn = f⁽ⁿ⁾(x0)/(n!) = f(x0, . . . , x₀

| {z }

n+1

). Jeśli zaś ¯x0 6= ¯xj to równość bn = f(¯x0,¯x1, . . . ,¯xn) wynika z (6.8) po podstawieniu i = 0 i j = n, oraz z założenia indukcyjnego. 2

74 ROZDZIAŁ 6. INTERPOLACJA WIELOMIANOWA

Uwagi i uzupełnienia

U. 6.1 Algorytm Hornera okazuje się optymalny. Każdy inny algorytm ob-liczający dokładną wartość wielomianu znając jego współczynniki wymaga wykonania co najmniej n mnożeń i n dodawań. Algorytm Hornera jest też numerycznie poprawny, zob. Ćw. 6.1.

U. 6.2 Mając obliczone współczynniki bj wielomianu interpolacyjnego w postaci Newtona, można łatwo przejść do jego współczynników a_j w bazie potęgowej. Algorytm wynika bowiem z następującego wzoru rekurencyjnego.

Dla i = n, n − 1, . . . , 0, niech a⁽ⁱ⁾j , i ≤ j ≤ n, będą odpowiednimi współ-czynnikami w rozwinięciu

b_i+b_i+1(x−xi)+· · ·+bn(x−xi) · · · (x−xn−1) = a⁽ⁱ⁾_i +a⁽ⁱ⁾_i+1x+· · ·+a⁽ⁱ⁾n xⁿ⁻ⁱ. Wtedy a⁽ⁱ⁾n = bn oraz

a⁽ⁱ⁾_j = a⁽ⁱ⁺¹⁾_j − xi∗ a⁽ⁱ⁺¹⁾_j+1 , i≤ j ≤ n − 1.

Uwzględniając fakt, że poszukiwanymi współczynnikami są aj = a⁽⁰⁾_j , 0 ≤ j≤ n, algorytm może być zrealizowany następująco:

for j := 0 to n do a[j] := b[j];

for i := n − 1 downto 0 do for j := i to n − 1 do

a[j] := a[j] − x[i] ∗ a[j + 1].

U. 6.3 Okazuje się, że przy realizacji w fl_ν algorytmu różnic dzielonych istotną rolę odgrywa porządek węzłów. Można pokazać, że algorytm liczenia f (t0, . . . , tn) jest numerycznie poprawny ze względu na dane interpolacyjne f⁽ⁱ⁾(t_j), o ile węzły są uporządkowane nierosnąco lub niemalejąco, zob. Ćw.

6.5.

U. 6.4 Zauważmy, ze pojęcie różnicy dzielonej formalnie zdefiniowaliśmy je-dynie dla ciągu węzłów postaci x0, . . . , x₀, x₁, . . . , x₁, . . . , x_k, . . . , x_k, gdzie xj

są parami różne. Tą definicję można rozszerzyć do dowolnego ciągu węzłów.

Można bowiem powiedzieć, że f(t₀, t₁, . . . , t_n) jest współczynnikiem przy xⁿ wielomianu wt0,...,tn ∈ Πn interpolującego f w węzłach tj (uwzględniając krotności). Równoważnie,

f (t₀, t₁, . . . , t_n) = w⁽ⁿ⁾_t0,...,tn n! .

6.4. PRZYPADEK WĘZŁÓW WIELOKROTNYCH 75

Ćwiczenia

Ćw. 6.1 Pokazać numeryczną poprawność algorytmu Hornera obliczania wartości wielomianu ze względu na dane współczynniki a_j tego wielomianu.

Ćw. 6.2 Pokazać, że algorytm Hornera obliczania wartości w(ξ) wielomianu danego w postaci potęgowej jest jednocześnie algorytmem dzielenia tego wie-lomianu przez jednomian (x − ξ). Dokładniej, jeśli w(x) =^Pnj=0ajx^j to

w(x) = ^ Xn j=1

v_jx^j−1∗ (x − ξ) + v0, gdzie vj są zdefiniowane tak jak w algorytmie Hornera.

Ćw. 6.3 Zaproponować możliwie tani algorytm obliczający współczynniki wielomianu w bazie Newtona, dysponując współczynnikami tego wielomianu w bazie potęgowej.

Ćw. 6.4 Pokazać, że dla różnych węzłów t0, . . . , t_n mamy f (t₀, t₁, . . . , t_n) =

Xn j=0

f (t_j)

(tj− t0) · · · (tj− tj−1)(tj− tj+1) · · · (tj− tn). Ćw. 6.5 Niech t₀ < t₁ < . . . < t_n. Pokazać, że realizując w fl_ν algo-rytm różnic dzielonych, zamiast dokładnej wartości f(t₀, t₁, . . . , t_n) otrzymu-jemy fl_ν^f (t₀, t₁, . . . , t_n)^, która jest dokładną różnicą dzieloną dla danych f (t_j)(1 + εj), gdzie |εj| ≤ Kν, 0 ≤ j ≤ n.

Wskazówka Rozpatrzyć najpierw proste przypadki n = 1, 2. Dla dowolnego n skorzystać z Ćw. 6.4.

Ćw. 6.6 Uzasadnić, że iloraz różnicowy jest funkcją symetryczną, tzn. dla dowolnej permutacji ti0, . . . , t_is węzłów t0, . . . , t_s mamy

f (t_i0, . . . , t_is) = f(t0, . . . , t_s).

Ćw. 6.7 Niech w będzie wielomianem stopnia dokładnie n. Pokazać, że dla ustalonych t1, t2, . . . , ts funkcja

w_t1,...,ts(t) = w(t, t₁, t₂, . . . , t_s), t∈ R,

jest wielomianem stopnia dokładnie n − s dla s ≤ n, oraz jest zerem dla s≥ n + 1.

76 ROZDZIAŁ 6. INTERPOLACJA WIELOMIANOWA Ćw. 6.8 Pokazać, że jeśli funkcja f : D → R jest n-krotnie różniczkowalna na D w sposób ciągły, f ∈ C⁽ⁿ⁾(D), to

f (t0, t1, . . . , tn−1, t0) = lim

tn→t0

f (t0, t1, . . . , tn−1, tn).

Ćw. 6.9 Stosując algorytm różnic dzielonym znależć wielomian w postaci Hermite’a interpolujący funkcję f(x) = x⁴ w trzech dwukrotnych węzłach równych 0, 1 i 2.

æ

Rozdział 7

Interpolacja a aproksymacja funkcji

Gdy mamy do czynienia z funkcją, która jest “skomplikowana” to często dobrze jest zastąpić ją funkcją “prostszą”. Mówimy wtedy o aproksyma-cji (przybliżaniu) funkaproksyma-cji. Funkcję musimy również aproksymać wtedy, gdy nie jesteśmy w stanie uzyskać pełnej o niej informacji. Na przykład, gdy funkcja reprezentuje pewien proces fizyczny to często zdarza się, że dysponujemy jedynie ciągiem próbek, czyli wartościami tej funkcji w pewnych punktach. Jasne jest, że chcielibyśmy przy tym, aby błąd aproksymacji był możliwie mały.

Z tego punktu widzenia, intepolacja wielomianowa może być trak-towana jako jeden ze sposobów aproksymacji funkcji, opartym na prób-kowaniu. Naturalnym staje się więc pytanie o błąd takiej aproksymacji.

7.1 Błąd interpolacji wielomianowej

Niech x0, x₁, . . . , x_n będą (niekoniecznie różnymi) węzłami należącymi do pewnego (być może nieskończonego) przedziału D ⊂ R. Dla da-nej funkcji f : D → R, przez w^f oznaczymy, tak jak w poprzednim rozdziale, wielomian interpolacyjny stopnia co najwyżej n interpolu-jący f w zadanych węzłach. W przypadku węzłów wielokrotnych jest to oczywiście wielomian interpolacyjny Hermite’a.

77

78ROZDZIAŁ 7. INTERPOLACJA A APROKSYMACJA FUNKCJI Lemat 7.1 Dla dowolnego punktu ¯x∈ D błąd interpolacji w ¯x wyraża się wzorem

f(¯x) − w^f(¯x) = (¯x − x⁰)(¯x − x¹) · · · (¯x − xⁿ)f(x0, x₁, . . . , x_n,¯x). (7.1) Jeśli ponadto f ∈ C⁽ⁿ⁺¹⁾(D), czyli pochodna f⁽ⁿ⁺¹⁾ w D istnieje i jest ciągła, to

f(¯x) − w^f(¯x) = (¯x − x⁰)(¯x − x¹) · · · (¯x − xⁿ)f⁽ⁿ⁺¹⁾(ξ) (n + 1)! ,

gdzie ξ = ξ(¯x) jest pewnym punktem należącym do najmniejszego prze-działu zawierającego punkty x0, x1, . . . , xn,¯x.

Dowód Możemy założyć, że ¯x nie jest żadnym z węzłów x_j, 0 ≤ j ≤ n. Niech ¯wf ∈ Πⁿ⁺¹ będzie wielomianem interpolacyjnym funkcji f opartym na węzłach x0, . . . , xn i dodatkowo na węźle ¯x. Mamy wtedy

¯

w_f(x) = wf(x) + (x − x⁰)(x − x¹) · · · (x − xⁿ)f(x0, x₁, . . . , x_n,¯x), a ponieważ z warunku interpolacyjnego f(¯x) = ¯wf(¯x), to mamy też (7.1).

Aby pokazać drugą część lematu, rozpatrzmy funkcję ψ : D → R, ψ(x) = f(x) − ¯wf(x)

= f(x) − w^f(x) − (x − x⁰)(x − x¹) · · · (x − xⁿ)f(x0, . . . , xn,¯x).

Z warunków interpolacyjnych na ¯wf ∈ Πn+1 wynika, że funkcja ψ ma punkty zerowe o łącznej krotności co najmniej n + 2. Wykorzystując twierdzenie Rolle’a wnioskujemy stąd, że ψ^′ ma zera o łącznej krotności co najmniej n + 1, ψ^′′ma zera o łącznej krotności co najmniej n, itd. W końcu funkcja ψ⁽ⁿ⁺¹⁾ zeruje się w co najmniej jednym punkcie ξ = ξ(¯x) należącym do najmniejszego przedziału zawierającego x0, x1, . . . , xn,¯x.

Wobec tego, że w⁽ⁿ⁺¹⁾_f ≡ 0, a (n + 1)-sza pochodna wielomianu (x − x₀) · · · (x − xⁿ) wynosi (n + 1)!, mamy

0 = ψ⁽ⁿ⁺¹⁾(ξ) = f⁽ⁿ⁺¹⁾(ξ) − (n + 1)!f(x⁰, . . . , xn,¯x).

Stąd

f(x₀, x₁, . . . , x_n,¯x) = f⁽ⁿ⁺¹⁾(ξ) (n + 1)! ,

7.1. BŁĄD INTERPOLACJI WIELOMIANOWEJ 79 co kończy dowód. 2

Zwykle interesuje nas nie tyle błąd w ustalonym punkcie ¯x ∈ D, ale na całym przedziale D. Zakładając teraz, że przedział D jest domknięty, czyli

D = [a, b]

dla pewnych −∞ < a < b < +∞, błąd ten będziemy mierzyć w normie jednostajnej (Czebyszewa). Dla funkcji ciągłej g : [a, b]→ R, norma ta jest zdefiniowana jako

kgk^C([a,b]) = max

x∈D |g(x)|.

Niech F_M^r ([a, b]), gdzie r ≥ 0, będzie klasą funkcji

F_M^r ([a, b]) = { f ∈ C^(r+1)([a, b]) : kf^(r+1)k^C([a,b]) ≤ M }, gdzie 0 < M < ∞. Mamy następujące twiedzenie.

Twierdzenie 7.1 Załóżmy, że każdą funkcję f ∈ FM^r ([a, b]) aproksy-mujemy jej wielomianem interpolacyjnym wf ∈ Π^r opartym na r + 1 węzłach x0, . . . , xr ∈ [a, b]. Wtedy maksymalny błąd takiej aproksymacji wynosi

e(F_M^r ([a, b]); x0, x1, . . . , xr) = max

f ∈F_M^r([a,b])kf − wfk^C([a,b])

= M

(r + 1)! · max

a≤x≤b|(x − x0) · · · (x − x^r)|.

Dowód Oszacowanie górne wynika bezpośrednio z Lematu 7.1, bo-wiem dla f ∈ FM^r ([a, b]) mamy

kf − wfk^C([a,b]) = max

a≤x≤b|f(x) − wf(x)|

= max

a≤x≤b|(x − x⁰) · · · (x − x^r)||f^(r+1)(ξ(x))|

(r + 1)!

≤ M

(r + 1)!max

x∈D |(x − x0) · · · (x − xr)|.

80ROZDZIAŁ 7. INTERPOLACJA A APROKSYMACJA FUNKCJI Z drugiej strony zauważmy, że dla wielomianu v(x) = Mx^r+1/(r+1)!

mamy v ∈ FM^r([a, b]) oraz kv − w^vk^C([a,b]) = M

(r + 1)! · max

a≤x≤b|(x − x⁰) · · · (x − x^r)|, co kończy dowód. 2

Zauważmy, że błąd aproksymacji e(F_M^r([a, b]); x0, . . . , xr) w istotny sposób zależy od wyboru węzłów xj. Naturalne jest więc teraz następu-jące pytanie. W których punktach xj przedziału [a, b] należy obliczać wartości funkcji, aby błąd był minimalny? Problem ten sprowadza się oczywiście do minimalizacji wielkości maxa≤x≤b|(x − x0) · · · (x − x^r)|

względem węzłów xj.

Twierdzenie 7.2 Błąd aproksymacji F_M^r ([a, b])(x0,· · · , x^r) jest mini-malny gdy węzły

x^∗_j = b− a

2 ·cos^2j + 1 2r + 2π



+ a+ b

2 , 0 ≤ j ≤ r.

Ponadto, dla węzłów optymalnych x^∗_j mamy

e(F_M^r ([a, b]); x^∗₀, . . . , x^∗_r) = 2M (r + 1)!

b− a 4

!r+1

.

Dowód tego twierdzenia opiera się na własnościach pewnego waż-nego ciągu wielomianów, który teraz przedstawimy.

W dokumencie Trzynaście Wykładów z Matematyki Obliczeniowej (Stron 76-85)