Kubiczne funkcje sklejane

Jeśli zdecydowaliśmy się na użycie kubicznych funkcji sklejanych to powstaje problem wyznaczenia s_f ∈ S₂ interpolującej daną funkcję f , tzn. takiej, że s_f(x_i) = f (x_i) dla 0 ¬ i ¬ n. W tym celu, na każdym przedziale [x_i, x_i+1] przedstawimy s_f w postaci jej rozwinięcia w szereg Taylora w punkcie x_i,

sf(x) = wi(x) = ai+ bi(x − xi) + ci

(x − xi)² 2 + di

(x − xi)³

6 ,

i podamy algorytm obliczania a_i, b_i, c_i, d_i dla 0 ¬ i ¬ n − 1.

Warunki brzegowe i warunki ciągłości dla s⁰⁰_f w wewnętrznych punktach x_i przedziału [a, b] dają nam w₀⁰⁰(x₀) = 0 = w⁰⁰_n−1(x_n) oraz w⁰⁰_i(x_i+1) = w_i+1⁰⁰ (x_i+1), czyli

c₀ = 0,

c_i+ d_ih_i = c_i+1, 0 ¬ i ¬ n − 2, c_n−1+ d_n−1h_n−1 = 0,

gdzie h_i = x_i+1− x_i. Stąd, przyjmując dodatkowo c_n= 0, otrzymujemy di = c_i+1− c_i

h_i , 0 ¬ i ¬ n − 1. (8.2)

Z warunków ciągłości dla s⁰_f w punktach x_i dostajemy z kolei

b_i+ c_ih_i+ ¹₂d_ih²_i = b_i+1, 0 ¬ i ¬ n − 2,

ROZDZIAŁ 8. INTERPOLACJA FUNKCJAMI SKLEJANYMI 70 i uwzględniając (8.2)

b_i+1 = b_i +¹₂h_i(c_i+1+ c_i), 0 ¬ i ¬ n − 2. (8.3) Warunki ciągłości s_f dają w końcu

a_i+ b_ih_i+ ¹₂c_ih²_i +¹₆d_ih³_i = a_i+1, 0 ¬ i ¬ n − 2. (8.4) Równania (8.2), (8.3) i (8.4) definiują nam na odcinku [a, b] naturalną kubiczną funkcję sklejaną.

Ponieważ poszukiwana funkcja sklejana s_f ma interpolować f , mamy dotatkowych n + 1 warunków interpolacyjnych w_i(x_i) = f (x_i) dla 0 ¬ i ¬ n − 1, oraz w_n−1(x_n) = f (x_n), czyli

a_i = f (x_i), 0 ¬ i ¬ n − 1, (8.5)

a_n−1+ b_n−1h_n−1+¹₂c_n−1h²_n−1+¹₆d_n−1h³_n−1 = f (x_n). (8.6) Wykorzystując (8.5) i (8.6) możemy teraz (8.4) przepisać w postaci

f (x_i+1) = f (x_i) + b_ih_i+¹₂c_ih²_i + ¹₆d_ih³_i,

przy czym wzór ten zachodzi również dla i = n − 1. Po wyrugowaniu b_i i podstawieniu d_i z (8.2), mamy

b_i = f [x_i, x_i+1] − ¹₆h_i(c_i+1+ 2c_i), 0 ¬ i ¬ n − 1, (8.7) gdzie f [xi, xi+1] jest odpowiednią różnicą dzieloną. Możemy teraz powyższe wyrażenie na bi podstawić w równaniach (8.3), aby otrzymać

1

6c_ih_i+ ¹₃c_i+1(h_i+ h_i+1) + ¹₆c_i+2h_i+1 = f [x_i+1, x_i+2] − f [x_i, x_i+1], 0 ¬ i ¬ n − 2.

Wprowadzając oznaczenie

c^∗_i = ¹₆ c_i, (8.8)

możemy to równanie przepisać w postaci h_i

Ostatecznie, aby obliczyć współczynniki a_i, b_i, c_i, d_i należy najpierw rozwiązać układ równań (8.9), a potem zastosować wzory (8.8), (8.2), (8.7) i (8.5).

ROZDZIAŁ 8. INTERPOLACJA FUNKCJAMI SKLEJANYMI 71 Zauważmy, że macierz układu (8.9) jest trójdiagonalna i ma dominującą przekątną. Układ ten można więc rozwiązać kosztem proporcjonalnym do n używając algorytmu przeganiania z U. 3.6.

Koszt znalezienia wszystkich współczynników naturalnej kubicznej funkcji sklejanej interpolującej f jest więc też proporcjonalny do n.

Na końcu oszacujemy jeszcze błąd interpolacji naturalnymi kubicznymi funkcjami sklejanymi na przedziale [a, b]. Będziemy zakładać, że f jest dwa razy różniczkowalna w sposób ciągły.

Twierdzenie 8.3 Jeśli f ∈ C_M² ([a, b]) to

kf − s_fk_C([a,b]) ¬ 5 M max

1¬i¬n(xi− x_i−1)².

W szczególności, dla podziału równomiernego, x_i = a + _nⁱ(b − a) dla 0 ¬ i ¬ n, mamy kf − s_fk_C([a,b]) ¬ 5 M

b − a n

2

.

Dowód. Wykorzystamy obliczoną wcześniej postać interpolującej funkcji sklejanej s_f. Niech x będzie punktem w przedziale [x_i, x_i+1]. Wtedy

s_f(x) = f (x_i) + f (x_i+1) − f (x_i) hi

− h_i(c^∗_i+1+ 2c^∗_i)

!

(x − x_i) + 3c^∗_i(x − x_i)²+c^∗_i+1− c^∗_i hi

(x − x_i)³. Z rozwinięcia f w szereg Taylora w punkcie x_idostajemy f (x) = f (x_i)+f⁰(x_i)(x−x_i)+¹₂f⁰⁰(ξ₁)(x−x_i)² oraz _h¹

i(f (x_i+1) − f (x_i)) = f⁰(x_i) + ¹₂h_if⁰⁰(ξ₂) dla pewnych ξ₁, ξ₂. Stąd f (x) − s_f(x) = f⁰⁰(ξ₁)

2 (x − x_i)²− f⁰⁰(ξ₂)

2 − (c^∗_i+1+ 2c^∗_i)

!

h_i(x − x_i)

− 3 c^∗_i(x − xi)²− c^∗_i+1− c^∗_i

h_i (x − xi)³, czyli

|f (x) − s_f(x)| ¬ ^M + 2|c^∗_i+1| + 6|c^∗_i|^h²_i = ^M + 8 max

1¬i¬n−1|c^∗_i|^h²_i. (8.10) Niech teraz max_1¬i¬n−1|c^∗_i| = |c^∗_s|. Przyjmując u₁ = 0 = w_n−1, z postaci układu (8.9) otrzymujemy

|c^∗_s| ¬ 2|c^∗_s| − |c^∗_s|(us+ ws) ¬ |usc^∗_s−1+ 2c^∗_s+ wsc^∗_s+1| = |f [xs−1, xs, xs+1]| ¬

   

f⁰⁰(ξ3) 2

    ¬ 1

2M, a stąd i z (8.10)

|f (x) − sf(x)| ¬ 5M h²_i, co kończy dowód. 

ROZDZIAŁ 8. INTERPOLACJA FUNKCJAMI SKLEJANYMI 72

Uwagi i uzupełnienia

U. 8.1 Zamiast terminu funkcje sklejane używa się też często terminów splajny (ang. spline-sklejać), albo funkcje gięte.

U. 8.2 W Rozdziale 7 pokazaliśmy, że aproksymacja kawałkami wielomianowa stopnia r jest optymalna (co do rzędu zbieżności) w klasie C_M^r+1([a, b]), wśród wszystkich aproksymacji korzystających jedynie z informacji o wartościach funkcji w danej liczbie n punktów. Okazuje się, że podobne optymalne własności wykazują funkcje sklejane, ale w innych klasach funkcji.

Niech

W_M^r (a, b) = ⁿf ∈ W^r(a, b) : Z b

a

f^(r)(x)^2 dx ¬ M^o.

Ustalmy węzły a = x₀ < · · · < x_n = b. Dla f ∈ W_M^r (a, b), niech s_f będzie naturalną funkcją sklejaną interpolującą f w x_j dla 0 ¬ j ¬ n, a a_f dowolną inną aproksymacją korzystającą jedynie z informacji o wartościach f w tych węzłach , tzn.

a_f = φ(f (x₀), . . . , f (x_n)).

Załóżmy, że błąd aproksymacji mierzymy nie w normie Czebyszewa, ale w normie średniokwadratowej, zde-finiowanej jako

kgk_L2(a,b) = s

Z b a

(g(x))²dx.

Wtedy

sup

f ∈W_M^r (a,b)

kf − s_fk_L2(a,b) ¬ sup

f ∈W_M^r (a,b)

kf − a_fk_L2(a,b).

Aproksymacja naturalnymi funkcjami sklejanymi jest więc optymalna w klasie W_M^r (a, b).

Można również pokazać, że interpolacja s^∗_f naturalnymi funkcjami sklejanymi na węzłach równoodległych, x_j = a+_nⁱ(b−a) dla 0 ¬ j ¬ n, jest optymalna co do rzędu w klasie W_M^r (a, b), wśród wszystkich aproksymacji korzystających jedynie z informacji o wartościach funkcji w n + 1 dowolnych punktach, oraz

max

f ∈W_M^r (a,b)

kf − s^∗_fk_L2(a,b)  n^−r.

U. 8.3 Tak jak wielomiany, naturalne funkcje sklejane interpolujące dane funkcje można reprezentować przez ich współczynniki w różnych bazach. Do tego celu można na przykład użyć bazy kanonicznej K_j, 0 ¬ j ¬ n, zdefiniowanej równościami

K_j(x_i) =

( 0 i 6= j, 1 i = j,

przy której s_f(x) =^Pn_j=0f (xj)K_j(x). Baza kanoniczna jest jednak niewygodna w użyciu, bo funkcje K_j w ogólności nie zerują się na żadnym podprzedziale, a tym samym manipulowanie nimi jest utrudnione.

Częściej używa się bazy B-sklejanej B_j, 0 ¬ j ¬ n. W przypadku funkcji kubicznych, r = 2, jest ona zdefiniowana przez następujące warunki:

B_j(x_j) = 1, dla 0 ¬ j ¬ n,

Bj(x) = 0, dla x ¬ x_j−2, j ­ 2, oraz dla x ­ xj+2, j ¬ n − 2.

Dla B₀ i B₁ dodatkowo żądamy, aby

B₀⁰⁰(x₀) = 0 = B₁⁰⁰(x₀), B₁(x₀) = 0, a dla B_n−1 i B_n podobnie

B_n−1⁰⁰ (x_n) = 0 = B⁰⁰_n(x_n), B_n−1(x_n−1) = 0.

ROZDZIAŁ 8. INTERPOLACJA FUNKCJAMI SKLEJANYMI 73 Wtedy B_j nie zerują się tylko na przedziale (x_j−2, xj+2), Wyznaczenie współczynników rozwinięcia w bazie {B_i}ⁿ_i=0 funkcji sklejanej interpolującej f wymaga rozwiązania układu liniowego z macierzą trójdiagonalną {B_j(x_i)}ⁿ_i,j=0, a więc koszt obliczenia tych współczynników jest proporcjonalny do n.

U. 8.4 Oprócz naturalnych funkcji sklejanych często rozpatruje się też okresowe funkcje sklejane. Są to funkcje ˜s : R → R spełniające warunki (i), (ii) z Rozdziału 8.1, oraz warunek:

(iii)’ ˜s⁽ⁱ⁾ jest dla 0 ¬ i ¬ r − 1 funkcją okresową o okresie (b − a), tzn. ˜s⁽ⁱ⁾(x) = ˜s⁽ⁱ⁾(x + b − a).

Klasę okresowych funkcji sklejanych rzędu r oznaczymy przez ˜S_r. Funkcje te mają podobne własności jak naturalne funkcje sklejane. Dokładniej, niech

W˜^r(a, b) = ^f ∈ W^r(a, b) : f⁽ⁱ⁾(a) = f⁽ⁱ⁾(b), 0 ¬ i ¬ r − 1 ,

tzn. ˜W^r(a, b) jest klasą funkcji z W^r(a, b), które można przedłużyć do funkcji, krórych wszystkie pochodne do rzędu r − 1 włącznie są (b − a)-okresowe na R. Wtedy dla dowolnej funkcji f ∈ ˜W^r(a, b) zerującej się w węzłach x_j, oraz dla dowolnej ˜s ∈ ˜S_r mamy

Z b a

f^(r)(x) ˜s^(r)(x) dx = 0.

Jest to odpowiednik Lematu 8.1 w przypadku okresowym. W szczególności wynika z niego jednoznaczność rozwiązania zadania interpolacyjnego dla okresowych funkcji f (tzn. takich, że f (a) = f (b)), jak również odpowiednia własność minimalizacyjna okresowych funkcji sklejanych. Dokładniej, jeśli f ∈ ˜W^r(a, b) oraz

˜s_f ∈ ˜S_r interpoluje f w węzłach x_j, 0 ¬ j ¬ n, to Z b

a

f^(r)(x)^2dx ­ Z b

a



s˜^(r)_f (x)^2dx. (8.11)

U. 8.5 Rozpatrzmy przestrzeń W^r(a, b). Utożsamiając funkcje f₁, f₂ ∈ W^r(a, b) takie, że f₁(x) − f₂(x) ∈ Π_r−1, zdefiniujemy w W^r(a, b) normę

kf k = s

Z b a

f^(r)(x)
²dx.

Dla ustalonych węzłów a = x₀ < · · · < x_n = b, niech S_r będzie podprzestrzenią W^r(a, b) składającą się z naturalnych funkcji sklejanych rzędu r opartych na węzłach x_j dla 0 ¬ j ¬ n. Oczywiście dim S_r = n + 1, co wynika z jednoznaczności rozwiązania w S_r zadania interpolacji. Okazuje się, że wtedy optymalną dla f ∈ W^r(a, b) funkcją w podprzestrzeni S jest naturalna funkcja sklejana s_f interpolująca f w węzłach x_j, tzn.

kf − s_fk = min

s∈Sr

kf − sk.

Rzeczywiście, ponieważ norma w przestrzeni W^r(a, b) generowana jest przez iloczyn skalarny hf₁, f₂i =

Z b a

f₁^(r)(x)f₂^(r)(x) dx,

jest to przestrzeń unitarna. Znane twierdzenie mówi, że w przestrzeni unitarnej najbliższą danej f funkcją w dowolnej domkniętej podprzestrzeni V jest rzut prostopadły f na V , albo równoważnie, taka funkcja vf ∈ V_n+1, że iloczyn skalarny

hf − v_f, vi = 0 dla wszystkich v ∈ V, cf. U. 7.7. W naszym przypadku, ostatnia równość jest równoważna

Z b a

(f − v_f)^(r)(x) s^(r)(x) dx = 0.

To zaś jest na mocy Lematu 8.1 prawdą gdy v_f interpoluje f w punktach x_j, czyli v_f = s_f.

ROZDZIAŁ 8. INTERPOLACJA FUNKCJAMI SKLEJANYMI 74

Ćwiczenia

Ćw. 8.1 Dla z ∈ R, niech

z₊ = max(0, z) =

( z, z ­ 0, 0, z < 0.

Pokaż, że każdą funkcję sklejaną s rzędu r można przedstawić w postaci s(x) = w(x) +

n

X

j=0

aj

(x − x_j)₊^2r−1,

gdzie a_j są pewnymi współczynnikami rzeczywistymi, a w ∈ Π_2r−1. Ponadto, jeśli s jest naturalna, to w ∈ Π_r−1, a współczynniki a_j spełniają równania

n

X

j=0

ajxⁱ_j = 0 dla 0 ¬ i ¬ r − 1.

Ćw. 8.2 Niech h > 0 i c ∈ R. Wyznacz współczynniki kubicznej funkcji sklejanej s opartej na pięciu węzłach

−2h, −h, 0, h, 2h i spełniającej dodatkowo następujące warunki interpolacyjne:

s(0) = c, s^(k)(±2h) = 0 dla k = 0, 1, 2.

Ćw. 8.3 Wykorzystując Ćw. 8.2 wykaż, że dla dowolnego układu węzłów a ¬ x₁ < · · · < x_n¬ b mamy max^kf⁰⁰k_C([a,b]): f ∈ C_M² ([a, b]), f (x_i) = 0, 1 ¬ i ¬ n ­ M

24

b − a n + 1

2

.

Stąd i z U. 7.13 wywnioskuj, że dla każdej aproksymacji a_f korzystającej jedynie z wartości f w n punktach, max

f ∈C_M² ([a,b])

kf − a_fk_C([a,b]) ­ M 24

b − a n + 1

2

.

Ćw. 8.4 Pokaż, że funkcje B_j, 0 ¬ j ¬ n, zdefiniowane w U. 8.3 istnieją i tworzą bazę w przestrzeni S_r. Ćw. 8.5 Pokaż jednoznaczność rozwiązania zadania interpolacyjnego w przypadku okresowych funkcji skle-janych zdefiniowanych w U. 8.4, oraz własność minimalizacyjną (8.11).

Wskazówka. Zastosuj technikę dowodową podobną do tej z przypadku naturalnych funkcji sklejanych.

Ćw. 8.6 Pokaż, że ogólnie postawione zadanie aproksymacyjne z U. 8.5 ma rozwiązanie.

Rozdział 9

Całkowanie numeryczne

Zajmiemy się teraz zadaniem całkowania numerycznego. Polega ono na obliczeniu (a raczej przybli-żeniu) całki oznaczonej

S(f ) =

Z b a

f (x) dx,

gdzie −∞ < a < b < +∞, a f należy do pewnej klasy F funkcji rzeczywistych określonych i całkowalnych w sensie Riemanna na całym przedziale [a, b].

Będziemy zakładać, że mamy możliwość obliczania wartości funkcji f , a w niektórych przypadkach również jej pochodnych, o ile istnieją. Dokładna całka S(f ) będzie więc w ogólności przybliżana wartością A(f ), która zależy tylko od wartości f i ew. jej pochodnych w skończonej liczbie punktów.

9.1 Co to są kwadratury?

Kwadraturami nazywamy funkcjonały liniowe Q_n : F → R postaci Q_n(f ) =

n

X

i=0

a_if (x_i), albo ogólniej

Q_n(f ) =

k

X

i=0 ni−1

X

j=0

a_i,jf^(j)(x_i), (9.1)

gdzie x_i są punktami z [a, b], współczynniki a_i i a_i,j są rzeczywiste, oraz n + 1 =^Pk_i=0n_i. Zauważmy, że obliczenia kwadratur są dopuszczalne w naszym modelu obliczeniowym, mogą więc służyć jako sposób przybliżania całki.

Jeden z możliwych sposobów konstrukcji kwadratur jest następujący. Najpierw wybieramy węzły xj (pojedyncze lub wielokrotne), budujemy wielomian interpolacyjny odpowiadający tym węzłom, a następnie całkujemy go. Ponieważ wielomian interpolacyjny zależy tylko od danej informacji o f , otrzymana w ten sposób wartość też będzie zależeć tylko od tej informacji, a w konsekwencji funkcjonał wynikowy będzie postaci (9.1). Są to tzw. kwadratury interpolacyjne.

Definicja 9.1 Kwadraturę Q^I_n opartą na węzłach o łącznej krotności n + 1 nazywamy kwadraturą interpolacyjną, jeśli

Q^I_n(f ) =

Z b a

w_f,n(x) dx,

gdzie w_f,njest wielomianem interpolacyjnym funkcji f stopnia co najwyżej n, opartym na tych węzłach.

75

ROZDZIAŁ 9. CAŁKOWANIE NUMERYCZNE 76 Współczynniki kwadratur interpolacyjnych można łatwo wyliczyć. Rozpatrzmy dla uproszczenia przypadek, gdy węzły są jednokrotne. Zapisując wielomian interpolacyjny w postaci jego rozwinięcia w bazie kanonicznej Lagrange’a l_i (cf. (6.3)), otrzymujemy

Q^I_n(f ) =

Kwadratura prostokątów jest oparta na jednym węźle x₀ = ¹₂(a + b), R(f ) = (b − a)f

a + b 2



.

Kwadratura trapezów jest oparta na jednokrotnych węzłach x₀ = a, x₁ = b i jest równa polu odpo-wiedniego trapezu,

T (f ) = b − a 2

f (a) + f (b)^.

Kwadratura parabol (Simpsona) jest oparta na jednokrotnych węzłach x₀ = a, x₁ = ¹₂(a + b), x₂ = b, i jest równa polu pod parabolą interpolującą f w tych węzłach,

P (f ) = b − a

Zauważmy, że kwadratury trapezów i parabol są oparte na węzłach jednokrotnych i równoodle-głych, przy czym x₀ = a i x_n= b. Ogólnie, kwadratury interpolacyjne oparte na węzłach równoodle-głych, x_i = a + _nⁱ(b − a) dla 0 ¬ i ¬ n, nazywamy kwadraturami Newtona–Cotesa.

W dokumencie Matematyka Obliczeniowa (kurs podstawowy dla studentów UW) Leszek Plaskota Instytut Matematyki Stosowanej i Mechaniki Uniwersytet Warszawski (Stron 72-79)