Interpolacja wielomianowa

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Zadanie interpolacji, czyli poprowadzenia krzywej zadanego rodzaju przez zestaw danych punktów, jest jednym z podstawowych zadań obliczeniowych. Stosuje się je nagminnie w najróżniejszych dziedzinach życia, np. wtedy, gdy trzeba

na podstawie próbki sygnału dźwiękowego (to znaczy: ciągu wartości amplitud sygnału zmierzonych w kolejnych odstępach czasu), odtworzyć jego przebieg;
przybliżyć wykres skomplikowanej (lub wręcz nieznanej) funkcji na podstawie jej wartości uprzednio stablicowanych w wybranych punktach;

Interpolację stosuje się szczególnie chętnie w samej numeryce. Na przykład idea metody siecznych polega na tym, by funkcję, której miejsca zerowego szukamy, przybliżyć prostą interpolującą tę funkcję w dwóch punktach. Metody numerycznego całkowania oraz rozwiązywania równań różniczkowych także korzystają z interpolacji.

Wielomian $w$ (czerwony) stopnia 6, interpolujący 7 zadanych wartości (zaznaczone na zielono) danej funkcji $f$

Niech Parser nie mógł rozpoznać (nieznana funkcja „\subsetR”): {\displaystyle \displaystyle D\subsetR} i niech $F$ będzie pewnym zbiorem funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:D\toR} . Niech $x_{0}, x_{1}, \dots, x_{n}$ będzie ustalonym zbiorem parami różnych punktów z $D$ , zwanych później węzłami.

Powiemy, że wielomian $w$ interpoluje funkcję $f \in F$ w węzłach $x_{j}$ , gdy

w (x_{j}) = f (x_{j}), 0 \leq j \leq n .

Oznaczmy przez $Π_{n}$ przestrzeń liniową wielomianów stopnia co najwyżej $n$ o współczynnikach rzeczywistych,

Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle \Pi_n\,=\,\{\,w(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0:\; a_j\inR, 0\le j\le n\,\}. }

Zadanie znalezienia wielomianu interpolującego zadane wartości nazywamy zadaniem interpolacji Lagrange'a.

Twierdzenie o istnieniu i jednoznaczności zadania interpolacji Lagrange'a

Dla dowolnej funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:D\toR} istnieje dokładnie jeden wielomian $w_{f} \in Π_{n}$ interpolujący $f$ w węzłach $x_{j}$ , $0 \leq j \leq n .$

Dowód

Wybierzmy w $Π_{n}$ dowolną bazę wielomianów $φ_{j}$ , $0 \leq j \leq n$ ,

Π_{n} = span {φ_{0}, φ_{1}, \dots, φ_{n}} .

Wtedy każdy wielomian z $Π_{n}$ można jednoznacznie przedstawić w postaci rozwinięcia względem wybranej bazy. Warunkiem koniecznym i dostatecznym na to, aby wielomian $w_{f} (\cdot) = \sum_{j = 0}^{n} c_{j} φ_{j} (\cdot)$ interpolował $f$ jest spełnienie układu $n + 1$ równań liniowych

\sum_{j = 0}^{n} c_{j} φ_{j} (x_{i}) = f (x_{i}), 0 \leq i \leq n,

z $n + 1$ niewiadomymi $c_{j}$ , który w postaci macierzowej wygląda następująco:

(\begin{array}{cccc} φ_{0} (x_{0}) & φ_{1} (x_{0}) & \dots & φ_{n} (x_{0}) \\ φ_{0} (x_{1}) & φ_{1} (x_{1}) & \dots & φ_{n} (x_{1}) \\ ⋮ \\ φ_{0} (x_{n}) & φ_{1} (x_{n}) & \dots & φ_{n} (x_{n}) \end{array}) (\begin{array}{c} c_{0} \\ c_{1} \\ ⋮ \\ c_{n} \end{array}) = (\begin{array}{c} f (x_{0}) \\ f (x_{1}) \\ ⋮ \\ f (x_{n}) \end{array}) .

Aby wykazać, że układ ten ma jednoznaczne rozwiązanie [[Algebra liniowa z geometrią analityczną/Wykład 8: Zastosowania wyznacznika. Układy równań liniowych|wystarczy, aby wektor zerowy był jedynym rozwiązaniem układu jednorodnego]]. Rzeczywiście, układ jednorodny odpowiada interpolacji danych zerowych, $f (x_{i}) = 0$ , $\forall i$ . Istnienie niezerowego rozwiązania byłoby więc równoważne istnieniu niezerowego wielomianu stopnia nie większego od $n$ , który miałby $n + 1$ różnych zer $x_{i}$ , co jest niemożliwe.

Zadanie znalezienia dla danej funkcji $f$ jej wielomianu interpolacyjnego stopnia co najwyżej $n$ jest więc dobrze zdefiniowane, tzn. rozwiązanie istnieje i jest wyznaczone jednoznacznie. Zauważmy, że wielomian interpolacyjny $w_{f}$ jako taki nie może być wynikiem obliczeń w naszym modelu obliczeniowym. Możemy natomiast wyznaczyć jego współczynniki $c_{j}$ w wybranej bazie.

Definicja

Niech $(φ_{j})_{j = 0}^{n}$ będzie bazą w przestrzeni $Π_{n}$ wielomianów stopnia co najwyżej $n$ . Zadanie interpolacji wielomianowej polega na obliczeniu dla danej funkcji $f$ współczynników $c_{j}$ takich, że wielomian

w_{f} (\cdot) = \sum_{j = 0}^{n} c_{j} φ_{j} (\cdot)

interpoluje $f$ w punktach $x_{j}$ , $0 \leq j \leq n$ .

Wybór bazy wielomianowej

Jak już wiemy, zadanie interpolacji Lagrange'a sprowadza się do rozwiązania układu równań liniowych. Okazuje się, że w zależności od wyboru sposobu reprezentacji naszego wielomianu (czyli od wyboru bazy wielomianowej $(φ_{j})_{j = 0}^{n}$ ), układ ten może być albo bardzo łatwy do rozwiązania, albo bardzo trudny. Co więcej, jego rozwiązanie w arytmetyce $f l_{ν}$ może napotykać na większe bądź mniejsze trudności (w zależności np. od uwarunkowania macierzy układu, który musimy rozwiązać).

W matematyce, jeden byt może być opisany na wiele równoważnych sposobów. W numeryce, każdy z nich może mieć diametralnie różne własności numeryczne: od odporności na błędy zaokrągleń, po koszt rozwiązania.
Dlatego, optymalizacja algorytmów numerycznych zaczyna się często od wyrażenia tego samego --- inaczej.

W naturalny sposób powstaje więc problem wyboru "wygodnej" bazy w $Π_{n}$ . Rozpatrzymy trzy bazy: Lagrange'a, potęgową i Newtona.

Baza Lagrange'a (kanoniczna)

Zdefiniujmy dla $0 \leq j \leq n$ wielomiany

l_{j} (x) = \frac{(x - x_{0}) (x - x_{1}) \dots (x - x_{j - 1}) (x - x_{j + 1}) \dots (x - x_{n})}{(x_{j} - x_{0}) (x_{j} - x_{1}) \dots (x_{j} - x_{j - 1}) (x_{j} - x_{j + 1}) \dots (x_{j} - x_{n})} .

Zauważmy, że każdy z $l_{j}$ jest stopnia dokładnie $n$ oraz

l_{j} (x_{i}) = {\begin{cases} 0 & i \neq j, \\ 1 & i = j . \end{cases}

Teraz widać, że wielomiany te stanowią bazę w $Π_{n}$ , którą nazywamy bazą Lagrange'a. Macierz układu zadania interpolacji jest w takim wypadku identycznością i w konsekwencji $c_{j} = f (x_{j})$ , $\forall j$ . Wielomian interpolacyjny dla funkcji $f$ można więc zapisać jako

w_{f} (x) = \sum_{j = 0}^{n} f (x_{j}) l_{j} (x) .

Koszt kombinatoryczny rozwiązania zadania interpolacji jest przy tym zerowy.

Wzory barycentryczne

Przypuśćmy, że chcielibyśmy obliczyć wartość wielomianu interpolacyjnego $w_{f}$ w punkcie $x$ różnym od $x_{j}$ , $0 \leq j \leq n$ . Podstawiając

w_{j} = \frac{1}{(x_{j} - x_{0}) (x_{j} - x_{1}) \dots (x_{j} - x_{j - 1}) (x_{j} - x_{j + 1}) \dots (x_{j} - x_{n})}

oraz $p_{n} (x) = (x - x_{0}) \dots (x - x_{n})$ mamy pierwszy wzór barycentryczny

w_{f} (x) = p_{n} (x) \sum_{j = 0}^{n} \frac{w_{j} f (x_{j})}{x - x_{j}},

i ostatecznie dostajemy tzw. drugi wzór barycentryczny na wielomian interpolacyjny,

w_{f} (x) = \frac{\sum_{j = 0}^{n} q_{j} (x) f (x_{j})}{\sum_{j = 0}^{n} q_{j} (x)},

gdzie $q_{j} (x) = w_{j} / (x - x_{j})$ . W ostatniej równości wykorzystaliśmy fakt, że $p_{n} (x) \equiv (\sum_{j = 0}^{n} q_{j} (x))^{- 1}$ , co łatwo widzieć, rozpatrując zadanie interpolacji funkcji $f \equiv 1$ . Drugi wzór barycentryczny jest korzystniejszy w implementacji.

Dla wielu układów węzłów, wagi $w_{j}$ są zadane jawnymi wzorami, np. dla węzłów równoodległych (niezależnie od tego, na jakim odcinku!) wagi w drugim wzorze barycentrycznym wynoszą po prostu

w_{j} = (- 1)^{j} (\begin{matrix} n \\ j \end{matrix}) .

Również dla węzłów Czebyszewa istnieją eleganckie wzory na takie współczynnki.

Można pokazać, że wartość $\tilde{w_{f} (x)}$ wielomianu iterpolacyjnego obliczona w arytmetyce $f l_{ν}$ według pierwszego wzoru barycentrycznego spełnia

\tilde{w_{f} (x)} = p_{n} (x) \sum_{j = 0}^{n} \frac{w_{j}}{x - x_{j}} f (x_{j}) (1 + ϵ_{j}),

gdzie $| ϵ_{j} | \leq 5 (n + 1)$ , a więc jest to algorytm numerycznie poprawny. Zachowanie drugiej postaci wzoru barycentrycznego w arytmetyce $f l_{ν}$ jest nieco bardziej skomplikowane.

Baza potęgowa (naturalna)

Znacznie prościej można obliczyć wartość wielomianu interpolacyjnego, (a także jego pochodnych), gdy jest on dany w najczęściej używanej bazie potęgowej, $φ_{j} (x) = x^{j}$ , $\forall j$ . Jeśli bowiem

w_{f} (x) = a_{0} + a_{1} x + \dots + a_{n} x^{n},

to również

w_{f} (x) = (\dots (a_{n} x + a_{n - 1}) x + a_{n - 2}) x + \dots + a_{1}) x + a_{0},

co sugeruje zastosowanie następującego schematu Hornera do obliczenia $w_{f} (x)$ :

Algorytm Algorytm Hornera

<math>\displaystyle v_n = a_n;</math>
for (j=n-1; j >= 0 ; j--)
	<math>\displaystyle v_j\, = \,v_{j+1}\cdot x\,+\,a_j</math>;

Po wykonaniu tego algorytmu $w_{f} (x) = v_{0}$ . Schemat Hornera wymaga wykonania tylko $n$ mnożeń i $n$ dodawań. Ma on również głębszy sens, bo jego produktem ubocznym mogą być także wartości pochodnych naszego wielomianu w $x$ . Algorytm Hornera okazuje się optymalny. Każdy inny algorytm obliczający dokładną wartość wielomianu, gdy danymi są współczynniki wielomianu, wymaga wykonania co najmniej $n$ mnożeń i $n$ dodawań. Algorytm Hornera jest też numerycznie poprawny.

Zauważmy jednak, że w przypadku bazy potęgowej macierz $(x_{i}^{j})_{i, j = 0}^{n}$ układu zadania interpolacji jest pełna. Jest to tzw. macierz Vandermonde'a. Obliczenie współczynników wielomianu interpolacyjnego w bazie potęgowej bezpośrednio z tego układu, stosując jedną ze znanych nam już metod, kosztowałoby rzędu $n^{3}$ operacji arytmetycznych. Co gorsza, w często spotykanym przypadku, gdy węzły interpolacji są równoodległe, ta macierz jest bardzo źle uwarunkowana!

Baza Newtona

Rozwiązaniem pośrednim, które łączy prostotę obliczenia współczynników z prostotą obliczenia wartości $w_{f} (x)$ i ewentualnie jego pochodnych, jest wybór bazy Newtona,

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned p_0(x) &= 1, \\ p_j(x) &= (x-x_0)(x-x_1)\cdots(x-x_{j-1}),\qquad 1\le j\le n. \endaligned}

W tym przypadku współczynniki rozwinięcia wielomianu interpolacyjnego będziemy oznaczać przez $b_{j}$ ,

w_{f} = \sum_{j = 0}^{n} b_{j} p_{j} .

Zwróćmy od razu uwagę na ważną własność bazy Newtona. Jeśli $w_{f, j} \in Π_{j}$ jest wielomianem interpolacyjnym dla funkcji $f$ opartym na węzłach $x_{0}, x_{1}, \dots, x_{j}$ , $0 \leq j \leq n$ , to $w_{f, 0} = b_{0}$ oraz

w_{f, j} = w_{f, j - 1} + b_{j} p_{j}, 1 \leq j \leq n .

Wartość $w_{f} (x)$ możemy obliczyć, stosując prostą modyfikację algorytmu Hornera:

Algorytm Algorytm Hornera dla bazy Newtona

<math>\displaystyle v_n = b_n;</math>
for (j=n-1; j >= 0 ; j--)
	<math>\displaystyle v_j\, = \,v_{j+1}\cdot (x-x_j)\,+\,b_j</math>;

Ponadto układ równań zadania interpolacji jest trójkątny dolny, o specyficznej strukturze, dzięki czemu można stworzyć elegancki algorytm, który teraz przedstawimy.

Algorytm różnic dzielonych

Różnicę dzieloną funkcji $f$ opartą na różnych węzłach $t_{0}, t_{1}, \dots, t_{s}$ , gdzie $s \geq 1$ , definiuje się indukcyjnie jako

f (t_{0}, t_{1}, \dots, t_{s}) = \frac{f (t_{1}, t_{2}, \dots, t_{s}) - f (t_{0}, t_{1}, \dots, t_{s - 1})}{t_{s} - t_{0}} .

Zachodzi następujące ważne twierdzenie.

Twierdzenie O różnicach dzielonych

Współczynniki $b_{j}$ wielomianu interpolacyjnego Newtona dla danej funkcji $f$ dane są przez różnice dzielone $f$ w węzłach $x_{0}, x_{1}, \dots, x_{j}$ , tzn.

b_{j} = f (x_{0}, x_{1}, \dots, x_{j}), 0 \leq j \leq n .

Dowód

Dla $0 \leq i \leq j \leq n$ , oznaczmy przez $w_{i, j}$ wielomian z $Π_{j - i}$ interpolujący $f$ w węzłach $x_{i}, x_{i + 1}, \dots, x_{j}$ . Wtedy ma miejsce następująca równość ( $i < j$ ):

w_{i, j} (x) = \frac{(x - x_{i}) w_{i + 1, j} (x) - (x - x_{j}) w_{i, j - 1} (x)}{x_{j} - x_{i}}, \forall x .

Aby ją pokazać wystarczy, że prawa strona tej równości, którą oznaczymy przez $v (x)$ , przyjmuje wartości $f (x_{s})$ dla $x = x_{s}$ , $i \leq s \leq j$ . Rzeczywiście, jeśli $i + 1 \leq s \leq j - 1$ to

v (x_{s}) = \frac{(x_{s} - x_{i}) f (x_{s}) - (x_{s} - x_{j}) f (x_{s})}{x_{j} - x_{i}} = f (x_{s}) .

Ponadto

v (x_{i}) = \frac{- (x_{i} - x_{j})}{x_{j} - x_{i}} f (x_{i}) = f (x_{i}),

oraz podobnie $v (x_{j}) = f (x_{j})$ . Stąd $v$ jest wielominem z $Π_{j - i}$ interpolującym $f$ w węzłach $x_{s}$ , $i \leq s \leq j$ , czyli $w_{i, j} = v$ .

Dalej postępujemy indukcyjnie ze względu na stopień $n$ wielomianu interpolacyjnego. Dla $n = 0$ mamy oczywiście $b_{0} = f (x_{0})$ . Niech $n \geq 1$ . Ponieważ, jak łatwo zauważyć,

w_{0, n} (x) = w_{0, n - 1} (x) + b_{n} p_{n} (x),

z założenia indukcyjnego mamy $b_{j} = f (x_{0}, \dots, x_{j})$ dla $0 \leq j \leq n - 1$ . Aby pokazać podobną równość dla $b_{n}$ , zauważmy, że

w_{0, n} (x) = \frac{(x - x_{0}) w_{1, n} (x) - (x - x_{n}) w_{0, n - 1} (x)}{x_{n} - x_{0}} .

Zauważmy teraz, że $b_{n}$ jest współczynnikiem przy $x^{n}$ w wielomianie $w_{0, n}$ . Z założenia indukcyjnego wynika, że współczynniki przy $x^{n - 1}$ w wielomianach $w_{1, n}$ i $w_{0, n - 1}$ są ilorazami różnicowymi opartymi odpowiednio na węzłach $x_{1}, \dots, x_{n}$ i $x_{0}, \dots, x_{n - 1}$ . Stąd

b_{n} = \frac{f (x_{1}, \dots, x_{n}) - f (x_{0}, \dots, x_{n - 1})}{x_{n} - x_{0}} = f (x_{0}, x_{1}, \dots, x_{n}),

co kończy dowód.

Różnicę dzieloną $f (x_{0}, x_{1}, \dots, x_{n})$ można łatwo obliczyć na podstawie wartości $f (x_{j})$ , $0 \leq j \leq n$ , budując następującą tabelkę:

\begin{array}{llllll} x_{0} & f (x_{0}) \\ x_{1} & f (x_{1}) & f (x_{0}, x_{1}) \\ x_{2} & f (x_{2}) & f (x_{1}, x_{2}) & f (x_{0}, x_{1}, x_{2}) \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \\ x_{n} & f (x_{n}) & f (x_{n - 1}, x_{n}) & f (x_{n - 2}, x_{n - 1}, x_{n}) & \dots & f (x_{0}, x_{1}, \dots, x_{n}) . \end{array}

<flash>file=Interpolacja.swf|width=550|height=300</flash>

Wyznaczenie wielomianu

w

interpolującego zestaw punktów

(0, 2) (1, 5) (- 1, 7)

algorytmem różnic dzielonych

Zauważmy przy tym, że "po drodze" obliczamy $f (x_{i}, x_{i + 1}, \dots, x_{j})$ dla wszystkich $0 \leq i < j \leq n$ , a więc w szczególności również interesujące nas różnice dzielone $f (x_{0}, x_{1}, \dots, x_{j})$ . Stąd i z twierdzenia o różnicach dzielonych wynika algorytm obliczania współczynników $b_{j}$ wielomianu interpolacyjnego w bazie Newtona. Po wykonaniu następującego algorytmu,

Algorytm Metoda różnic dzielonych

for (j = 0; j <= n; j++)
	<math>\displaystyle b_j</math> = <math>\displaystyle f(x_j)</math>; 
for (j = 0; j <= n; j++)
	for (k = n; k >= j; k--)
		<math>\displaystyle b_j</math> = <math>\displaystyle (b_k-b_{k-1})/(x_k - x_{k-j})</math>;

współczynniki $b_{j}$ na końcu algorytmu zawierają wspólczynniki wielomianu interpolacyjnego w bazie Newtona. Czy gdybyś zobaczył ten algorytm na samym początku tego wykładu, zgadłbyś, do czego może służyć?!

<flash>file=Interpolacjainsitu.swf|width=550|height=300</flash>

Wyznaczenie tego samego wielomianu

w

, interpolującego zestaw punktów

(0, 2) (1, 5) (- 1, 7)

algorytmem różnic dzielonych --- wykonanym tym razem in situ.

Okazuje się, że przy realizacji w $f l_{ν}$ algorytmu różnic dzielonych istotną rolę odgrywa porządek węzłów. Można pokazać, że --- o ile węzły są uporządkowane nierosnąco lub niemalejąco --- algorytm liczenia $f (t_{0}, \dots, t_{n})$ jest numerycznie poprawny ze względu na dane interpolacyjne $f (t_{j})$ , a cały algorytm różnic dzielonych daje w arytmetyce $f l_{ν}$ współczynniki wielomianu interpolacyjnego, będące niewiekim zaburzeniem wartości dokładnych.

Uwarunkowanie

Danymi w zadaniu interpolacji są zarówno wartości interpolowanej funkcji, jak i węzły interpolacji. Traktując węzły jako sztywno zadane parametry zadania i dopuszczając jedynie zaburzenia wartości funkcji, można pokazać, że jeśli zamiast $f$ rozpatrzyć jej zaburzenie $f + Δ f$ , gdzie $| Δ f | \leq ϵ$ , to

| w_{f} (x) - w_{f + Δ f} (x) | \leq cond (x, f) | w_{f} (x) | ϵ,

gdzie

cond (x, f) = \frac{\sum_{j = 0}^{n} | l_{j} (x) f (x_{j}) |}{| p_{n} (x) |} \geq 1 .

Znacznie rzadziej rozważa się uwarunkowanie zadania interpolacji ze względu na zaburzenie węzłów. Warto zaznaczyć, że zaburzenie danych interpolacji tylko w jednym punkcie może mieć wpływ na przebieg całego wielomianu interpolacyjnego, co ukazuje poniższy przykład:

Przykład

Pokażemy zmianę kilku bazowych wielomianów Lagrange'a stopnia 10 (dla węzłów równoodległych w $[0, 1]$ ) w sytuacji, gdy trzeci węzeł interpolacji zostanie zaburzony o 0.01.

Wybrane wielomiany bazowe Lagrange'a oparte na węzłach równoodległych (zielone) kontra te same wielomiany, oparte na tych samych węzłach, z jednym wyjątkiem: węzeł $x_{3} = 0.2$ został zmieniony na $x_{3} = 0.21$ (czerwone).

Jak widać, to lokalne zaburzenie danych może powodować wyraźne globalne zaburzenie całego wielomianu interpolacyjnego (zwróć uwagę na prawy koniec przedziału!).

Biblioteki

MATLAB i Octave mają wbudowaną funkcję wyznaczającą wielomian, interpolujący zadane wartości: jeśli x jest wektorem zawierającym $N$ węzłów, a y --- wektorem zawierającym wartości w węzłach, to

c = polyfit(x,y,N-1);

daje współczynniki wielomianu interpolacyjnego (Ostatni argument jest równy $N - 1$ , bo taki powinien być stopień wielomianu interpolacyjnego Lagrange'a!).

Co ciekawe (i budzące trochę zgrozy!) --- wielomian (zarówno w MATLABie, jak w Octave) jest wyznaczany w bazie naturalnej, przez rozwiązanie układu równań z macierzą Vandermonde'a, a więc w sposób najgorszy z możliwych. Nie sądzisz, że czas najwyższy, aby to zmienić? Napisz odpowiedni kod i wyślij do Octave-forge!

Aby teraz wyznaczyć wartości takiego wielomianu w zadanych punktach $X$ , także musimy użyć specjalnej funkcji,

Y = polyval(c,X);

Domyślamy się, że implementuje ona algorytm Hornera.

Przykład

Interpolujemy tabelkę


$x$	2	1	0
$y$	5	2	1

wielomianem stopnia co najwyżej 2.

octave:1> x = [2, 1, 0]
x =
  2  1  0

octave:2> y = [5, 2, 1]
y =
  5  2  1

octave:3> c = polyfit(x,y,2)
c =
  1  0  1
  
octave:4> polyval(c,3)
ans =  10

Zgodnie z przewidywaniami, otrzymaliśmy wielomian $1 \cdot x^{2} + 0 \cdot x + 1$ . Wartość tego wielomianu dla $x = 3$ rzeczywiście jest równa 10.

A co się stanie, gdy będziemy szukać wielomianu stopnia niższego?

octave:6> c1 = polyfit(x,y,1)
c1 =
   2.00000   0.66667

Też "coś" zostało obliczone --- wielomian (jak domyślamy się) $2 \cdot x + \frac{2}{3}$ . Nie dziwi, że ten wielomian nie jest wielomianem interpolacyjnym (dlaczego?) --- więc czym może być? Okazuje się, że to coś to wielomian nalepiej pasujący do danych w sensie aproksymacji średniokwadratowej, o czym będzie mowa w innym wykładzie.

Warto jeszcze może wiedzieć, że polyfit można także wywołać dla jeszcze wyższego stopnia wielomianu, jednak, co niespodziewane, wynikiem nie będzie wielomian stopnia 2, uzyskany poprzednio:

octave:7> c3 = polyfit(x,y,3)
c3 =
   0.21429   0.35714   0.42857   1.00000

Wynika to stąd, że gdy dopuszczalny stopień wielomianu jest wyższy niż wymagany w zadaniu interpolacji Lagrange'a, zadanie interpolacji ma nieskończenie wiele rozwiązań. Funkcja polyfit wybiera z nich to, które spełnia warunek, że norma euklidesowa wektora współczynników wielomianu jest najmniejsza z możliwych.

Pragnąc wykorzystać interpolację we własnym programie w C, najlepiej samemu zaprogramować bądź drugi wzór barycentryczny, bądź algorytm różnic dzielonych --- w zależności od potrzeb.

Przypadek węzłów wielokrotnych

Uogólnieniem rozpatrzonego zadania interpolacji jest zadanie interpolacji Hermite'a. Zakładamy, że oprócz (różnych) węzłów $x_{j}$ dane są również ich krotności $n_{j}$ , $0 \leq j \leq k$ , przy czym $\sum_{j = 0}^{k} n_{j} = n + 1$ . Należy skonstruować wielomian $w_{f} \in Π_{n}$ taki, że

w_{f}^{(i)} (x_{j}) = f^{(i)} (x_{j}) dla 0 \leq i \leq n_{j} - 1, 0 \leq j \leq k .

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

Lemat

Zadanie interpolacji Hermite'a ma jednoznaczne rozwiązanie.

Dowód

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

Nas oczywiście interesuje konstrukcja wielomianu $w_{f}$ . W tym celu ustawimy węzły $x_{j}$ w ciąg

({\bar{x}}_{0}, {\bar{x}}_{1}, \dots, {\bar{x}}_{n}) = (\underset{n_{0}}{\underset{⏟}{x_{0}, \dots, x_{0}}}, \underset{n_{1}}{\underset{⏟}{x_{1}, \dots, x_{1}}}, \dots, \underset{n_{k}}{\underset{⏟}{x_{k}, \dots, x_{k}}})

i zdefiniujemy uogólnioną bazę Newtona w $Π_{n}$ jako

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned p_0(x) &= 1, \\ p_j(x) &= (x-\bar x_0)(x-\bar x_1)\cdots (x-\bar x_{j-1}), \qquad 1\le j\le n. \endaligned}

Uogólnimy również pojęcie różnicy dzielonej na węzły powtarzające się, kładąc

f ({\bar{x}}_{i}, {\bar{x}}_{i + 1}, \dots, {\bar{x}}_{j}) = \frac{f^{(j - i)} ({\bar{x}}_{i})}{(j - i)!}

dla ${\bar{x}}_{i} = {\bar{x}}_{i + 1} = \dots = {\bar{x}}_{j}$ , oraz

f ({\bar{x}}_{i}, {\bar{x}}_{i + 1}, \dots, {\bar{x}}_{j}) = \frac{f ({\bar{x}}_{i + 1}, \dots, {\bar{x}}_{j}) - f ({\bar{x}}_{i}, \dots, x_{j - 1})}{{\bar{x}}_{j} - {\bar{x}}_{i}}

dla ${\bar{x}}_{i} \neq {\bar{x}}_{j}$ . Zauważmy, że przy tej definicji różnice $f ({\bar{x}}_{i}, \dots, {\bar{x}}_{j})$ możemy łatwo obliczyć stosując schemat podobny do tego z przypadku węzłów jednokrotnych.

Twierdzenie O różnicach dzielonych dla interpolacji Hermite'a

Współczynniki $b_{j}$ wielomianu interpolacyjnego Hermite'a w bazie Newtona,

w_{f} (\cdot) = \sum_{j = 0}^{n} b_{j} p_{j} (\cdot),

dane są przez odpowiednie różnice dzielone, tzn.

b_{j} = f ({\bar{x}}_{0}, {\bar{x}}_{1}, \dots, {\bar{x}}_{j}), 0 \leq j \leq n .

Dowód

Dowód przeprowadzimy podobnie jak dla węzłów jednokrotnych. Niech $w_{i, j} \in Π_{j - i}$ oznacza wielomian interpolacyjny Hermite'a oparty na (być może powtarzających się) węzłach ${\bar{x}}_{i}, {\bar{x}}_{i + 1}, \dots, {\bar{x}}_{j}$ . To znaczy, $w_{i, j}$ interpoluje $f$ w węzłach $x_{s}$ takich, że $x_{s}$ występuje w ciągu ${\bar{x}}_{i}, \dots {\bar{x}}_{j}$ , a jego krotność jest liczbą powtórzeń $x_{s}$ w tym ciągu.

Zauważmy najpierw, że dla ${\bar{x}}_{i} \neq {\bar{x}}_{j}$ zachodzi znany nam już wzór,

w_{i, j} (x) = \frac{(x - {\bar{x}}_{i}) w_{i + 1, j} (x) - (x - {\bar{x}}_{j}) w_{i, j - 1} (x)}{{\bar{x}}_{j} - {\bar{x}}_{i}} .

Rzeczywiście, oznaczmy przez $v (x)$ prawą stronę powyższej równości. Dla $k$ mniejszego od krotności danego węzła $x_{s}$ w ciągu ${\bar{x}}_{i}, \dots {\bar{x}}_{j}$ , mamy $w_{i + 1, j}^{(k - 1)} (x_{s}) = w_{i, j - 1}^{(k - 1)} (x_{s})$ , a ponieważ

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned v^{(k)}(x)&=\frac{k\,(w_{i+1,j}^{(k-1)}(x)-w_{i,j-1}^{(k-1)}(x))} {\bar x_j-\bar x_i} \\ && \qquad +\, \frac{(x-\bar x_i)w_{i+1,j}^{(k)}(x)-(x-\bar x_j)w_{i,j-1}^{(k)}(x)} {\bar x_j-\bar x_i}, \endaligned}

to

v^{(k)} (x_{s}) = \frac{(x_{s} - {\bar{x}}_{i}) w_{i + 1, j}^{(k)} (x_{s}) - (x_{s} - {\bar{x}}_{j}) w_{i, j - 1}^{(k)} (x_{s})}{{\bar{x}}_{j} - {\bar{x}}_{i}} .

Korzystając z tego wzoru sprawdzamy, że $v$ spełnia odpowiednie warunki interpolacyjne, a stąd $w_{i, j} = v$ .

Dalej postępujemy indukcyjnie ze względu na $n$ . Dla $n = 0$ mamy $b_{0} = f (x_{0})$ . Dla $n \geq 1$ wystarczy pokazać, że $b_{n} = f ({\bar{x}}_{0}, {\bar{x}}_{1}, \dots, {\bar{x}}_{n})$ . W tym celu rozpatrzymy dwa przypadki.

Jeśli ${\bar{x}}_{0} = {\bar{x}}_{n}$ , to mamy jeden węzeł $x_{0}$ o krotności $n + 1$ . Wielomian interpolacyjny jest wtedy postaci

w_{f} (x) = \sum_{j = 0}^{n} \frac{f^{(j)} (x_{0})}{j!} (x - x_{0})^{j},

a stąd $b_{n} = f^{(n)} (x_{0}) / (n!) = f (\underset{n + 1}{\underset{⏟}{x_{0}, \dots, x_{0}}})$ . Jeśli zaś ${\bar{x}}_{0} \neq {\bar{x}}_{j}$ , to równość $b_{n} = f ({\bar{x}}_{0}, {\bar{x}}_{1}, \dots, {\bar{x}}_{n})$ wynika z wcześniej wyprowadzonych wzorów oraz z założenia indukcyjnego.

Uwaga

Zauważmy, ze pojęcie różnicy dzielonej formalnie zdefiniowaliśmy jedynie dla ciągu węzłów postaci $x_{0}, \dots, x_{0}, x_{1}, \dots, x_{1}, \dots, x_{k}, \dots, x_{k}$ , gdzie $x_{j}$ są parami różne. Tą definicję można rozszerzyć do dowolnego ciągu węzłów. Można bowiem powiedzieć, że $f (t_{0}, t_{1}, \dots, t_{n})$ jest współczynnikiem przy $x^{n}$ wielomianu $w_{t_{0}, \dots, t_{n}} \in Π_{n}$ interpolującego $f$ w węzłach $t_{j}$ (uwzględniając krotności). Równoważnie,

f (t_{0}, t_{1}, \dots, t_{n}) = \frac{w_{t_{0}, \dots, t_{n}}^{(n)}}{n!} .

Błąd interpolacji

Gdy mamy do czynienia z funkcją, która jest "skomplikowana", często dobrze jest zastąpić ją funkcją "prostszą". Mówimy wtedy o aproksymacji funkcji. 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, 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.

Podobnie ma się sprawa w przypadku implementacji funkcji elementarnych ( $\sin, \exp, . . .$ ) w bibliotece funkcji matematycznych, czy wręcz w procesorze. Tam również najchętniej poszukiwalibyśmy sposobu taniego przybliżenia wartości dokładnej funkcji. I rzeczywiście, często w tym celu stosuje się m.in. specjalnie konstruowaną aproksymację wielomianową.

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

Niech $x_{0}, x_{1}, \dots, x_{n}$ będą (niekoniecznie różnymi) węzłami należącymi do pewnego (być może nieskończonego) przedziału Parser nie mógł rozpoznać (nieznana funkcja „\subsetR”): {\displaystyle \displaystyle D\subsetR} . Dla danej funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:D\toR} , przez $w_{f}$ rozważamy, tak jak w całym wykładzie, wielomian interpolacyjny stopnia co najwyżej $n$ interpolujący $f$ w zadanych węzłach. W przypadku węzłów wielokrotnych jest to oczywiście wielomian interpolacyjny Hermite'a; gdy węzły są jednokrotne, mamy do czynienia z interpolacją Lagrange'a.

Lemat Postać błędu interpolacji

Dla dowolnego punktu $\bar{x} \in D$ błąd interpolacji w $\bar{x}$ wyraża się wzorem

f (\bar{x}) - w_{f} (\bar{x}) = (\bar{x} - x_{0}) (\bar{x} - x_{1}) \dots (\bar{x} - x_{n}) f (x_{0}, x_{1}, \dots, x_{n}, \bar{x}) .

Jeśli ponadto $f \in C^{(n + 1)} (D)$ , czyli pochodna $f^{(n + 1)}$ w $D$ istnieje i jest ciągła, to

f (\bar{x}) - w_{f} (\bar{x}) = (\bar{x} - x_{0}) (\bar{x} - x_{1}) \dots (\bar{x} - x_{n}) \frac{f^{(n + 1)} (ξ)}{(n + 1)!},

gdzie $ξ = ξ (\bar{x})$ jest pewnym punktem należącym do najmniejszego przedziału zawierającego punkty $x_{0}, x_{1}, \dots, x_{n}, \bar{x}$ .

Dowód

Możemy założyć, że $\bar{x}$ nie jest żadnym z węzłów $x_{j}$ , $0 \leq j \leq n$ . Niech ${\bar{w}}_{f} \in Π_{n + 1}$ będzie wielomianem interpolacyjnym funkcji $f$ opartym na węzłach $x_{0}, \dots, x_{n}$ i dodatkowo na węźle $\bar{x}$ . Mamy wtedy

{\bar{w}}_{f} (x) = w_{f} (x) + (x - x_{0}) (x - x_{1}) \dots (x - x_{n}) f (x_{0}, x_{1}, \dots, x_{n}, \bar{x}),

a ponieważ z warunku interpolacyjnego $f (\bar{x}) = {\bar{w}}_{f} (\bar{x})$ , to mamy też pierwszą równość w lemacie.

Aby pokazać drugą część lematu, rozpatrzmy funkcję Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle \psi:D\toR} ,

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \lefteqn{\psi(x) \;=\; f(x)-\bar w_f(x)} \\ &= \, f(x)-w_f(x)-(x-x_0)(x-x_1)\cdots(x-x_n) f(x_0,\ldots,x_n,\bar x). \endaligned}

Z warunków interpolacyjnych na ${\bar{w}}_{f} \in Π_{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 $ψ^{(n + 1)}$ zeruje się w co najmniej jednym punkcie $ξ = ξ (\bar{x})$ należącym do najmniejszego przedziału zawierającego $x_{0}, x_{1}, \dots, x_{n}, \bar{x}$ . Wobec tego, że $w_{f}^{(n + 1)} \equiv 0$ , a $(n + 1)$ -sza pochodna wielomianu $(x - x_{0}) \dots (x - x_{n})$ wynosi $(n + 1)!$ , mamy

0 = ψ^{(n + 1)} (ξ) = f^{(n + 1)} (ξ) - (n + 1)! f (x_{0}, \dots, x_{n}, \bar{x}) .

Stąd

f (x_{0}, x_{1}, \dots, x_{n}, \bar{x}) = \frac{f^{(n + 1)} (ξ)}{(n + 1)!},

co kończy dowód.

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

D = [a, b]

dla pewnych $- \infty < a < b < + \infty$ , błąd ten będziemy mierzyć w normie jednostajnej (Czebyszewa). Dla funkcji ciągłej Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle g:[a,b]\toR} , norma ta jest zdefiniowana jako

‖ g ‖_{C ([a, b])} = \max_{x \in D} | g (x) | .

Niech $F_{M}^{r} ([a, b])$ , gdzie $r \geq 0$ , będzie klasą funkcji

F_{M}^{r} ([a, b]) = {f \in C^{(r + 1)} ([a, b]) : ‖ f^{(r + 1)} ‖_{C ([a, b])} \leq M},

gdzie $0 < M < \infty$ . Mamy następujące twiedzenie.

Twierdzenie O najgorszym możliwym błędzie interpolacji w klasie

Załóżmy, że każdą funkcję $f \in F_{M}^{r} ([a, b])$ aproksymujemy jej wielomianem interpolacyjnym $w_{f} \in Π_{r}$ opartym na $r + 1$ węzłach $x_{0}, \dots, x_{r} \in [a, b]$ . Wtedy maksymalny błąd takiej aproksymacji wynosi

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned e(F^r_M([a,b]);x_0,x_1,\ldots,x_r) &= \max_{f\in F^r_M([a,b])} \|f-w_f\|_{ C([a,b])} \\ &= \frac M{(r+1)!}\cdot \max_{a\le x\le b}|(x-x_0)\cdots(x-x_r)|. \endaligned}

Dowód

Oszacowanie górne wynika bezpośrednio z lematu o postaci błędu interpolacji, bowiem dla $f \in F_{M}^{r} ([a, b])$ mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \|f-w_f\|_{ C([a,b])}&=\max_{a\le x\le b}|f(x)-w_f(x)| \\ &= \max_{a\le x\le b}|(x-x_0)\cdots(x-x_r)| \frac{|f^{(r+1)}(\xi(x))|}{(r+1)!} \\ &\le & \frac{M}{(r+1)!}\max_{x\in D}|(x-x_0)\cdots(x-x_r)|. \endaligned}

Z drugiej strony zauważmy, że dla wielomianu $v (x) = M \frac{x^{r + 1}}{(r + 1)!}$ mamy $v \in F_{M}^{r} ([a, b])$ oraz

‖ v - w_{v} ‖_{C ([a, b])} = \frac{M}{(r + 1)!} \cdot \max_{a \leq x \leq b} | (x - x_{0}) \dots (x - x_{r}) |,

co kończy dowód.

Zjawisko Rungego

Rozważmy zadanie interpolacji funkcji

f (x) = \frac{1}{1 + x^{2}}

w $N$ równoodległych węzłach na przedziale $[- 5, 5]$ . Okazuje się, że dla dużych wartości $N$ , wielomian interpolacyjny ma poważne kłopoty z aproksymacją tej funkcji przy krańcach przedziału:

Zjawisko Rungego: interpolacja w $N = 17$ węzłach równoodległych dla $f (x) = \frac{1}{1 + x^{2}}$

Z kolei wielomian oparty na węzłach Czebyszewa znacznie lepiej przybliża tę funkcję.

Zjawisko Rungego: interpolacja w węzłach równoodległych, kontra interpolacja w węzłach Czebyszewa

Rzeczywiście, węzły Czebyszewa zagęszczają się w pobliżu krańców odcinka.

Zjawisko Rungego: interpolacja w węzłach Czebyszewa

Wiąże się to z zachowaniem się samych wielomianów bazowych: wielomiany oparte na węzłach równoodległych właśnie silnie oscylują w pobliżu krańców przedziału (jasne: nasz wielomian jest wysokiego stopnia, musi mieć dużo zer, a z drugiej strony, jako wielomian wysokiego stopnia, chce szybko uciec do nieskończoności, dlatego "wije się" jak może). Natomiast wielomiany bazowe oparte na węzłach Czebyszewa są najspokojniejsze: wiją się, ale z umiarem, bo zagęszczone przy krańcach węzły skutecznie je "duszą".

Zauważmy, że błąd aproksymacji $e (F_{M}^{r} ([a, b]); x_{0}, \dots, x_{r})$ w istotny sposób zależy od wyboru węzłów $x_{j}$ . Naturalne jest więc teraz następujące pytanie: w których punktach $x_{j}$ 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 $\max_{a \leq x \leq b} | (x - x_{0}) \dots (x - x_{r}) |$ względem węzłów $x_{j}$ .

Twierdzenie O optymalnym doborze węzłów

Błąd aproksymacji w klasie funkcji $F_{M}^{r} ([a, b]) (x_{0}, \dots, x_{r})$ jest minimalny gdy węzły interpolacji są zadane jako węzły Czebyszewa na $(a, b)$ , tzn.

x_{j}^{*} = \frac{b - a}{2} \cdot \cos (\frac{2 j + 1}{2 r + 2} π) + \frac{a + b}{2}, 0 \leq j \leq r .

Ponadto, dla optymalnych węzłów $x_{j}^{*}$ mamy

e (F_{M}^{r} ([a, b]); x_{0}^{*}, \dots, x_{r}^{*}) = \frac{2 M}{(r + 1)!} {(\frac{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.

Wielomiany Czebyszewa

Ciąg ${T_{k}}_{k \geq 0}$ wielomianów Czebyszewa (pierwszego rodzaju) zdefiniowany jest indukcyjnie jako

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned T_0(x) &= 1, \\ T_1(x) &= x, \\ T_{k+1}(x) &= 2xT_k(x)-T_{k-1}(x),\qquad \mbox{ dla } \quad k\ge 1. \endaligned}

Pafnutij Lwowicz Czebyszew
Zobacz biografię

Zauważmy, że $T_{k}$ jest wielomianem stopnia dokładnie $k$ o współczynniku przy $x^{k}$ równym $2^{k - 1}$ ( $k \geq 1$ ). Ponadto wielomian $T_{k}$ można dla $| x | \leq 1$ przedstawić w postaci

T_{k} (x) = \cos (k \arccos x) .

Rzeczywiście, łatwo sprawdzić, że jest to prawdą dla $k = 0, 1$ . Stosując podstawienie $\cos t = x$ , $0 \leq t \leq π$ , oraz wzór na sumę cosinusów otrzymujemy dla $k \geq 1$

\cos ((k + 1) t) = 2 \cdot \cos t \cos (k t) - \cos ((k - 1) t),

co jest równoważne formule rekurencyjnej dla $T_{k + 1}$ .

Kilka pierwszych wielomianów Czebyszewa na odcinku $[- 1, 1]$

Ze wzoru $T_{k} (x) = \cos (k \arccos x)$ wynikają również inne ważne własności wielomianów Czebyszewa. Norma wielomianu Czebyszewa na $[- 1, 1]$ wynosi

‖ T_{k} ‖_{C ([- 1, 1])} = \max_{- 1 \leq x \leq 1} | T_{k} (x) | = 1

i jest osiągana w $k + 1$ punktach tego przedziału równych

y_{j} = \cos (\frac{j}{k} π), 0 \leq j \leq k,

przy czym $T_{k} (y_{j}) = (- 1)^{j}$ .

W końcu, $k$ -ty wielomian Czebyszewa $T_{k}$ ma dokładnie $k$ pojedynczych zer w $[- 1, 1]$ równych

z_{j} = \cos (\frac{2 j + 1}{2 r} π), 0 \leq j \leq k - 1 .

Miejsca zerowe wielomianu Czebyszewa będziemy nazywać węzłami Czebyszewa. Konsekwencją wymienionych własności jest następująca własność ekstremalna wielomianów Czebyszewa.

Przez ${\overline{Π}}_{k}$ oznaczymy klasę wielomianów stopnia $k$ o współczynniku wiodącym równym $1$ , tzn.

{\overline{Π}}_{k} = {w \in Π_{k} : w (x) = x^{k} + \dots} .

Twierdzenie O minimaksie

Niech $k \geq 1$ . W klasie ${\overline{Π}}_{k}$ minimalną normę jednostajną na przedziale $[- 1, 1]$ ma wielomian $w^{*} = 2^{1 - k} T_{k}$ , tzn.

\min_{w \in {\overline{Π}}_{k}} ‖ w ‖_{C ([- 1, 1])} = ‖ w^{*} ‖_{C ([- 1, 1])} = \frac{1}{2^{k - 1}} .

Wielomian stopnia 9 oparty na węzłach Czebyszewa kontra oparty na węzłach równoodległych. Zwróć uwagę na wielkie oscylacje tego drugiego pry końcach odcinka.

Możemy teraz przeprowadzić dowód twierdzenia o optymalnym doborze węzłów:

Dowód

Dowód wynika teraz bezpośrednio z twierdzenia o minimaksie. Zauważmy bowiem, że wielomian $(x - x_{0}) (x - x_{1}) \dots (x - x_{r})$ jest w klasie ${\overline{Π}}_{r + 1}$ . Stąd dla $[a, b] = [- 1, 1]$ optymalnymi węzłami są zera $z_{j}$ wielomianu Czebyszewa, przy których

(x - z_{0}) (x - z_{1}) \dots (x - z_{r}) = \frac{T_{r + 1} (x)}{2^{r}} .

Jeśli przedział $[a, b]$ jest inny niż $[- 1, 1]$ , należy dokonać liniowej zamiany zmiennych tak, aby przeszedł on na $[- 1, 1]$ . Bezpośrednie sprawdzenie pokazuje, że w klasie ${\overline{Π}}_{r + 1}$ minimalną normę Czebyszewa na przedziale $[a, b]$ ma wielomian

w_{a, b}^{*} (x) = (\frac{b - a}{2})^{r + 1} w^{*} (\frac{2 x - (a + b)}{b - a}) .

Stąd

‖ w_{a, b}^{*} ‖_{C ([a, b])} = (\frac{b - a}{2})^{r + 1} \frac{1}{2^{r}} = 2 (\frac{b - a}{4})^{r + 1}

i węzły

x_{j}^{*}

są optymalne.

Wielomiany Czebyszewa znajdują bardzo wiele, czasem zaskakujących, zastosowań w różnych działach numeryki, m.in. w konstrukcji metod iteracyjnych rozwiązywania równań liniowych.

Równie interesujący jest fakt, że wielomian interpolacyjny oparty na węzłach Czebyszewa jest prawie optymalnym przybliżeniem wielomianowym zadanej funkcji:

Twierdzenie Jacksona, o prawie optymalnej interpolacji w węzłach Czebyszewa

Dla $f \in C [- 1, 1]$ , wielomian interpolacyjny $w_{f}$ stopnia co najwyżej $n$ , oparty na węzłach Czebyszewa, spełnia

| | f - w_{f} | |_{C [- 1, 1]} \leq (2 + \frac{2}{π} \log (n + 1)) | | f - w_{f}^{*} | |_{C [- 1, 1]}

gdzie $w_{f}^{*}$ jest wielomianem stopnia co najwyżej $n$ , najlepiej aproksymującym $f$ w sensie normy jednostajnej.

Jeśli więc $n \leq 5$ , to wielomian oparty na węzłach Czebyszewa jest co najwyżej 3.02 razy, a gdy $n \leq 20$ --- maksymalnie 4 razy gorszy od optymalnego. Można więc powiedzieć, że jest prawie optymalny.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 6.1--6.3 w

D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

MN09

Spis treści

Interpolacja wielomianowa

Wybór bazy wielomianowej

Baza Lagrange'a (kanoniczna)

Wzory barycentryczne

Baza potęgowa (naturalna)

Baza Newtona

Algorytm różnic dzielonych

Uwarunkowanie

Biblioteki

Przypadek węzłów wielokrotnych

Błąd interpolacji

Zjawisko Rungego

Wielomiany Czebyszewa

Literatura

Menu nawigacyjne

Działania na stronie

Opcje strony

Narzędzia osobiste

Nawigacja

Szukaj

Narzędzia