Personal tools

MN09

From Studia Informatyczne


Spis treści

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  (czerwony) stopnia 6, interpolujący 7 zadanych wartości (zaznaczone na zielono) danej funkcji
Enlarge
Wielomian \displaystyle w (czerwony) stopnia 6, interpolujący 7 zadanych wartości (zaznaczone na zielono) danej funkcji \displaystyle f

Niech \displaystyle D\subsetR i niech \displaystyle F będzie pewnym zbiorem funkcji \displaystyle f:D\toR. Niech \displaystyle x_0,x_1,\ldots,x_n będzie ustalonym zbiorem parami różnych punktów z \displaystyle D, zwanych później węzłami.

Powiemy, że wielomian \displaystyle w interpoluje funkcję \displaystyle f\in F w węzłach \displaystyle x_j, gdy

\displaystyle w(x_j)\,=\,f(x_j),\qquad 0\le j\le n.

Oznaczmy przez \displaystyle \Pi_n przestrzeń liniową wielomianów stopnia co najwyżej \displaystyle n o współczynnikach rzeczywistych,

\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 \displaystyle f:D\to R istnieje dokładnie jeden wielomian \displaystyle w_f\in\Pi_n interpolujący \displaystyle f w węzłach \displaystyle x_j, \displaystyle 0\le j\le n.

Dowód

Wybierzmy w \displaystyle \Pi_n dowolną bazę wielomianów \displaystyle \varphi_j, \displaystyle 0\le j\le n,

\displaystyle \Pi_n\,=\, \mbox{span} \{\,\varphi_0,\varphi_1,\ldots,\varphi_n\,\}.

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

\displaystyle \sum_{j=0}^n c_j\varphi_j(x_i)\,=\,f(x_i),\qquad 0\le i\le n,

z \displaystyle n+1 niewiadomymi \displaystyle c_j, który w postaci macierzowej wygląda następująco:

\displaystyle    \left(\begin{array} {cccc}     \varphi_0(x_0) & \varphi_1(x_0) & \cdots & \varphi_n(x_0) \\     \varphi_0(x_1) & \varphi_1(x_1) & \cdots & \varphi_n(x_1) \\     \vdots \\     \varphi_0(x_n) & \varphi_1(x_n) & \cdots & \varphi_n(x_n)      \end{array} \right)\left(\begin{array} {c}      c_0 \\ c_1 \\ \vdots \\ c_n \end{array} \right)\,=\,     \left(\begin{array} {c}      f(x_0) \\ f(x_1) \\ \vdots \\ f(x_n) \end{array} \right).

Aby wykazać, że układ ten ma jednoznaczne rozwiązanie wystarczy, aby wektor zerowy był jedynym rozwiązaniem układu jednorodnego. Rzeczywiście, układ jednorodny odpowiada interpolacji danych zerowych, \displaystyle f(x_i)=0, \displaystyle \forall i. Istnienie niezerowego rozwiązania byłoby więc równoważne istnieniu niezerowego wielomianu stopnia nie większego od \displaystyle n, który miałby \displaystyle n+1 różnych zer \displaystyle x_i, co jest niemożliwe.

image:End_of_proof.gif

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

Definicja

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

\displaystyle    w_f(\cdot)\,=\,\sum_{j=0}^n c_j\varphi_j(\cdot)

interpoluje \displaystyle f w punktach \displaystyle x_j, \displaystyle 0\le j\le 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 \displaystyle (\varphi_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 \displaystyle fl_\nu 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 \displaystyle \Pi_n. Rozpatrzymy trzy bazy: Lagrange'a, potęgową i Newtona.

Baza Lagrange'a (kanoniczna)

Zdefiniujmy dla \displaystyle 0\le j\le n wielomiany

\displaystyle    l_j(x)\,=\,\frac      {(x  -x_0)(x  -x_1)\cdots(x  -x_{j-1})(x  -x_{j+1})\cdots(x  -x_n)}      {(x_j-x_0)(x_j-x_1)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_n)}.

Zauważmy, że każdy z \displaystyle l_j jest stopnia dokładnie \displaystyle n oraz

\displaystyle l_j(x_i)\,=\,\left\{\,\begin{array} {ll}                         0 & \quad i\ne j, \\ 1 & \quad i=j. \end{array} \right.

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

\displaystyle 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 \displaystyle w_f w punkcie \displaystyle x różnym od \displaystyle x_j, \displaystyle 0\le j\le n. Podstawiając

\displaystyle w_j\,=\,\frac 1 {(x_j-x_0)(x_j-x_1)\cdots(x_j-x_{j-1})       (x_j-x_{j+1})\cdots(x_j-x_n)}

oraz \displaystyle p_n(x)=(x-x_0)\cdots(x-x_n) mamy pierwszy wzór barycentryczny

\displaystyle    w_f(x)\,=\,p_n(x)\sum_{j=0}^n\frac{w_jf(x_j)}{x-x_j},

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

\displaystyle w_f(x)\,=\,\frac{\sum_{j=0}^n q_j(x)f(x_j)}{\sum_{j=0}^n q_j(x)},

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

Dla wielu układów węzłów, wagi \displaystyle 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

\displaystyle w_j = (-1)^j \begin{pmatrix}  n \\ j \end{pmatrix} .

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

Można pokazać, że wartość \displaystyle \widetilde{w_f(x)} wielomianu iterpolacyjnego obliczona w arytmetyce \displaystyle fl_\nu według pierwszego wzoru barycentrycznego spełnia

\displaystyle  \widetilde{w_f(x)} = p_n(x) \sum_{j=0}^n\frac{w_j}{x-x_j}f(x_j)(1+\epsilon_j),

gdzie \displaystyle |\epsilon_j| \leq 5(n+1), a więc jest to algorytm numerycznie poprawny. Zachowanie drugiej postaci wzoru barycentrycznego w arytmetyce \displaystyle fl_\nu 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, \displaystyle \varphi_j(x)=x^j, \displaystyle \forall j. Jeśli bowiem

\displaystyle w_f(x)\,=\,a_0+a_1x+\cdots+ a_nx^n,

to również

\displaystyle w_f(x)\,=\,(\cdots(a_nx+a_{n-1})x+a_{n-2})x+\cdots+a_1)x+a_0,

co sugeruje zastosowanie następującego schematu Hornera do obliczenia \displaystyle 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 \displaystyle w_f(x)=v_0. Schemat Hornera wymaga wykonania tylko \displaystyle n mnożeń i \displaystyle n dodawań. Ma on również głębszy sens, bo jego produktem ubocznym mogą być także wartości pochodnych naszego wielomianu w \displaystyle 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 \displaystyle n mnożeń i \displaystyle n dodawań. Algorytm Hornera jest też numerycznie poprawny.

Zauważmy jednak, że w przypadku bazy potęgowej macierz \displaystyle (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 \displaystyle 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 \displaystyle w_f(x) i ewentualnie jego pochodnych, jest wybór bazy Newtona,

\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 \displaystyle b_j,

\displaystyle w_f\,=\,\sum_{j=0}^n b_jp_j.

Zwróćmy od razu uwagę na ważną własność bazy Newtona. Jeśli \displaystyle w_{f,j}\in\Pi_j jest wielomianem interpolacyjnym dla funkcji \displaystyle f opartym na węzłach \displaystyle x_0,x_1,\ldots,x_j, \displaystyle 0\le j\le n, to \displaystyle w_{f,0}=b_0 oraz

\displaystyle    w_{f,j}\,=\,w_{f,j-1}\,+\,b_jp_j,\qquad 1\le j\le n.

Wartość \displaystyle 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 \displaystyle f opartą na różnych węzłach \displaystyle t_0,t_1,\ldots,t_s, gdzie \displaystyle s\ge 1, definiuje się indukcyjnie jako

\displaystyle    f(t_0,t_1,\ldots,t_s)\,=\,\frac      {f(t_1,t_2,\ldots,t_s)\,-\,f(t_0,t_1,\ldots,t_{s-1})}{t_s\,-\,t_0}.

Zachodzi następujące ważne twierdzenie.

Twierdzenie O różnicach dzielonych

Współczynniki \displaystyle b_j wielomianu interpolacyjnego Newtona dla danej funkcji \displaystyle f dane są przez różnice dzielone \displaystyle f w węzłach \displaystyle x_0,x_1,\ldots,x_j, tzn.

\displaystyle b_j\,=\,f(x_0,x_1,\ldots,x_j),\qquad 0\le j\le n.

Dowód

Dla \displaystyle 0\le i\le j\le n, oznaczmy przez \displaystyle w_{i,j} wielomian z \displaystyle \Pi_{j-i} interpolujący \displaystyle f w węzłach \displaystyle x_i,x_{i+1},\ldots,x_j. Wtedy ma miejsce następująca równość (\displaystyle i<j):

\displaystyle    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}, \qquad\forall x.

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

\displaystyle 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

\displaystyle v(x_i)\,=\,\frac{-(x_i-x_j)}{x_j-x_i}f(x_i)\,=\,f(x_i),

oraz podobnie \displaystyle v(x_j)=f(x_j). Stąd \displaystyle v jest wielominem z \displaystyle \Pi_{j-i} interpolującym \displaystyle f w węzłach \displaystyle x_s, \displaystyle i\le s\le j, czyli \displaystyle w_{i,j}=v.

Dalej postępujemy indukcyjnie ze względu na stopień \displaystyle n wielomianu interpolacyjnego. Dla \displaystyle n=0 mamy oczywiście \displaystyle b_0=f(x_0). Niech \displaystyle n\ge 1. Ponieważ, jak łatwo zauważyć,

\displaystyle w_{0,n}(x)\,=\,w_{0,n-1}(x)+b_n p_n(x),

z założenia indukcyjnego mamy \displaystyle b_j=f(x_0,\ldots,x_j) dla \displaystyle 0\le j\le n-1. Aby pokazać podobną równość dla \displaystyle b_n, zauważmy, że

\displaystyle 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 \displaystyle b_n jest współczynnikiem przy \displaystyle x^n w wielomianie \displaystyle w_{0,n}. Z założenia indukcyjnego wynika, że współczynniki przy \displaystyle x^{n-1} w wielomianach \displaystyle w_{1,n} i \displaystyle w_{0,n-1} są ilorazami różnicowymi opartymi odpowiednio na węzłach \displaystyle x_1,\ldots,x_n i \displaystyle x_0,\ldots,x_{n-1}. Stąd

\displaystyle b_n\,=\,\frac{f(x_1,\ldots,x_n)-f(x_0,\ldots,x_{n-1})}{x_n-x_0}     \,=\,f(x_0,x_1,\ldots,x_n),

co kończy dowód.

image:End_of_proof.gif

Różnicę dzieloną \displaystyle f(x_0,x_1,\ldots,x_n) można łatwo obliczyć na podstawie wartości \displaystyle f(x_j), \displaystyle 0\le j\le n, budując następującą tabelkę:

\displaystyle \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) \\   \vdots &\vdots &\vdots &\vdots &\ddots \\   x_n & f(x_n) & f(x_{n-1},x_n) & f(x_{n-2},x_{n-1},x_n) &\cdots     & f(x_0,x_1,\ldots,x_n).\end{array}


Wyznaczenie wielomianu \displaystyle w interpolującego zestaw punktów \displaystyle (0,2)\displaystyle (1,5)\displaystyle (-1,7) algorytmem różnic dzielonych

Zauważmy przy tym, że "po drodze" obliczamy \displaystyle f(x_i,x_{i+1},\ldots,x_j) dla wszystkich \displaystyle 0\le i < j\le n, a więc w szczególności również interesujące nas różnice dzielone \displaystyle f(x_0,x_1,\ldots,x_j). Stąd i z twierdzenia o różnicach dzielonych wynika algorytm obliczania współczynników \displaystyle 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 = 1; j <= n; j++)
	for (k = n; k >= j; k--)
		<math>\displaystyle b_k</math> = <math>\displaystyle (b_k-b_{k-1})/(x_k - x_{k-j})</math>;

współczynniki \displaystyle 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ć?!



Wyznaczenie tego samego wielomianu \displaystyle w, interpolującego zestaw punktów \displaystyle (0,2)\displaystyle (1,5)\displaystyle (-1,7) algorytmem różnic dzielonych --- wykonanym tym razem in situ.

Okazuje się, że przy realizacji w \displaystyle fl_\nu 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 \displaystyle f(t_0,\ldots,t_n) jest numerycznie poprawny ze względu na dane interpolacyjne \displaystyle f(t_j), a cały algorytm różnic dzielonych daje w arytmetyce \displaystyle fl_\nu 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 \displaystyle f rozpatrzyć jej zaburzenie \displaystyle f+\Delta f, gdzie \displaystyle |\Delta f| \leq \epsilon, to

\displaystyle |w_f(x) - w_{f+\Delta f}(x)| \leq  \mbox{cond} (x,f)|w_f(x)|\epsilon,

gdzie

\displaystyle  \mbox{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 \displaystyle [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ł  został zmieniony na  (czerwone).
Enlarge
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ł \displaystyle x_3 = 0.2 został zmieniony na \displaystyle 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!).


MATLAB i Octave mają wbudowaną funkcję wyznaczającą wielomian, interpolujący zadane wartości: jeśli x jest wektorem zawierającym \displaystyle 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 \displaystyle 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 \displaystyle X, także musimy użyć specjalnej funkcji,

Y = polyval(c,X);

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

Przykład

Interpolujemy tabelkę

\displaystyle x 2 1 0
\displaystyle 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 \displaystyle 1\cdot x^2 + 0\cdot x + 1. Wartość tego wielomianu dla \displaystyle 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ę) \displaystyle 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 \displaystyle x_j dane są również ich krotności \displaystyle n_j, \displaystyle 0\le j\le k, przy czym \displaystyle \sum_{j=0}^k n_j=n+1. Należy skonstruować wielomian \displaystyle w_f\in\Pi_n taki, że

\displaystyle w_f^{(i)}(x_j)\,=\,f^{(i)}(x_j)\qquad \mbox{ dla } \quad        0\le i\le n_j-1, 0\le j\le k.

Oczywiście zakładamy przy tym, że odpowiednie pochodne funkcji \displaystyle 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 \displaystyle n+1 równań z \displaystyle 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ż \displaystyle n, który miałby zera o łącznej krotności większej niż \displaystyle n.

image:End_of_proof.gif

Nas oczywiście interesuje konstrukcja wielomianu \displaystyle w_f. W tym celu ustawimy węzły \displaystyle x_j w ciąg

\displaystyle (\bar x_0,\bar x_1,\ldots,\bar x_n)\,=\,   (\underbrace{x_0,\ldots,x_0}_{n_0},\underbrace{x_1,\ldots,x_1}_{n_1},    \ldots,\underbrace{x_k,\ldots,x_k}_{n_k})

i zdefiniujemy uogólnioną bazę Newtona w \displaystyle \Pi_n jako

\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

\displaystyle f(\bar x_i,\bar x_{i+1},\ldots,\bar x_j)\,=\,     \frac{f^{(j-i)}(\bar x_i)}{(j-i)!}

dla \displaystyle \bar x_i=\bar x_{i+1}=\cdots=\bar x_j, oraz

\displaystyle f(\bar x_i,\bar x_{i+1},\ldots,\bar x_j)\,=\,\frac    {f(\bar x_{i+1},\ldots,\bar x_j)-f(\bar x_i,\ldots,x_{j-1})}     {\bar x_j-\bar x_i}

dla \displaystyle \bar x_i\ne\bar x_j. Zauważmy, że przy tej definicji różnice \displaystyle f(\bar x_i,\ldots,\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 \displaystyle b_j wielomianu interpolacyjnego Hermite'a w bazie Newtona,

\displaystyle w_f(\cdot)\,=\,\sum_{j=0}^n b_jp_j(\cdot),

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

\displaystyle b_j\,=\,f(\bar x_0,\bar x_1,\ldots,\bar x_j),\qquad 0\le j\le n.

Dowód

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

Zauważmy najpierw, że dla \displaystyle \bar x_i\ne\bar x_j zachodzi znany nam już wzór,

\displaystyle    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 \displaystyle v(x) prawą stronę powyższej równości. Dla \displaystyle k mniejszego od krotności danego węzła \displaystyle x_s w ciągu \displaystyle \bar x_i,\ldots\bar x_j, mamy \displaystyle w_{i+1,j}^{(k-1)}(x_s)=w_{i,j-1}^{(k-1)}(x_s), a ponieważ

\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

\displaystyle 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 \displaystyle v spełnia odpowiednie warunki interpolacyjne, a stąd \displaystyle w_{i,j}=v.

Dalej postępujemy indukcyjnie ze względu na \displaystyle n. Dla \displaystyle n=0 mamy \displaystyle b_0=f(x_0). Dla \displaystyle n\ge 1 wystarczy pokazać, że \displaystyle b_n\,=\,f(\bar x_0,\bar x_1,\ldots,\bar x_n). W tym celu rozpatrzymy dwa przypadki.

Jeśli \displaystyle \bar x_0=\bar x_n, to mamy jeden węzeł \displaystyle x_0 o krotności \displaystyle n+1. Wielomian interpolacyjny jest wtedy postaci

\displaystyle w_f(x)\,=\,\sum_{j=0}^n \frac{f^{(j)}(x_0)}{j!}(x-x_0)^j,

a stąd \displaystyle b_n=f^{(n)}(x_0)/(n!)=f(\underbrace{x_0,\ldots,x_0}_{n+1}). Jeśli zaś \displaystyle \bar x_0\ne\bar x_j, to równość \displaystyle b_n\,=\,f(\bar x_0,\bar x_1,\ldots,\bar x_n) wynika z wcześniej wyprowadzonych wzorów oraz z założenia indukcyjnego.

image:End_of_proof.gif

Uwaga

Zauważmy, ze pojęcie różnicy dzielonej formalnie zdefiniowaliśmy jedynie dla ciągu węzłów postaci \displaystyle x_0,\ldots,x_0,x_1,\ldots,x_1,\ldots,x_k,\ldots,x_k, gdzie \displaystyle x_j są parami różne. Tą definicję można rozszerzyć do dowolnego ciągu węzłów. Można bowiem powiedzieć, że \displaystyle f(t_0,t_1,\ldots,t_n) jest współczynnikiem przy \displaystyle x^n wielomianu \displaystyle w_{t_0,\ldots,t_n}\in\Pi_n interpolującego \displaystyle f w węzłach \displaystyle t_j (uwzględniając krotności). Równoważnie,

\displaystyle f(t_0,t_1,\ldots,t_n)\,=\,\frac{w^{(n)}_{t_0,\ldots,t_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 (\displaystyle \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 \displaystyle x_0,x_1,\ldots,x_n będą (niekoniecznie różnymi) węzłami należącymi do pewnego (być może nieskończonego) przedziału \displaystyle D\subsetR. Dla danej funkcji \displaystyle f:D\toR, przez \displaystyle w_f rozważamy, tak jak w całym wykładzie, wielomian interpolacyjny stopnia co najwyżej \displaystyle n interpolujący \displaystyle 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 \displaystyle \bar x\in D błąd interpolacji w \displaystyle \bar x wyraża się wzorem

\displaystyle    f(\bar x)-w_f(\bar x)\,=\,(\bar x-x_0)(\bar x-x_1)      \cdots(\bar x-x_n)f(x_0,x_1,\ldots,x_n,\bar x).

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

\displaystyle f(\bar x)-w_f(\bar x)\,=\,(\bar x-x_0)(\bar x-x_1)      \cdots(\bar x-x_n)\frac{f^{(n+1)}(\xi)}{(n+1)!},

gdzie \displaystyle \xi=\xi(\bar x) jest pewnym punktem należącym do najmniejszego przedziału zawierającego punkty \displaystyle x_0,x_1,\ldots,x_n,\bar x.

Dowód

Możemy założyć, że \displaystyle \bar x nie jest żadnym z węzłów \displaystyle x_j, \displaystyle 0\le j\le n. Niech \displaystyle \bar w_f\in\Pi_{n+1} będzie wielomianem interpolacyjnym funkcji \displaystyle f opartym na węzłach \displaystyle x_0,\ldots,x_n i dodatkowo na węźle \displaystyle \bar x. Mamy wtedy

\displaystyle \bar w_f(x)\,=\,w_f(x)\,+\,(x-x_0)(x-x_1)\cdots(x-x_n)      f(x_0,x_1,\ldots,x_n,\bar x),

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

Aby pokazać drugą część lematu, rozpatrzmy funkcję \displaystyle \psi:D\toR,

\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 \displaystyle \bar w_f\in\Pi_{n+1} wynika, że funkcja \displaystyle \psi ma punkty zerowe o łącznej krotności co najmniej \displaystyle n+2. Wykorzystując twierdzenie Rolle'a wnioskujemy stąd, że \displaystyle \psi' ma zera o łącznej krotności co najmniej \displaystyle n+1, \displaystyle \psi'' ma zera o łącznej krotności co najmniej \displaystyle n, itd. W końcu funkcja \displaystyle \psi^{(n+1)} zeruje się w co najmniej jednym punkcie \displaystyle \xi=\xi(\bar x) należącym do najmniejszego przedziału zawierającego \displaystyle x_0,x_1,\ldots,x_n,\bar x. Wobec tego, że \displaystyle w_f^{(n+1)}\equiv 0, a \displaystyle (n+1)-sza pochodna wielomianu \displaystyle (x-x_0)\cdots(x-x_n) wynosi \displaystyle (n+1)!, mamy

\displaystyle 0 \,=\, \psi^{(n+1)}(\xi)\,=\,f^{(n+1)}(\xi)-(n+1)!         f(x_0,\ldots,x_n,\bar x).

Stąd

\displaystyle f(x_0,x_1,\ldots,x_n,\bar x)\,=\,    \frac{f^{(n+1)}(\xi)}{(n+1)!},
co kończy dowód. image:End_of_proof.gif

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

\displaystyle D\,=\,[a,b]

dla pewnych \displaystyle -\infty<a<b<+\infty, błąd ten będziemy mierzyć w normie jednostajnej (Czebyszewa). Dla funkcji ciągłej \displaystyle g:[a,b]\toR, norma ta jest zdefiniowana jako

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

Niech \displaystyle F^r_M([a,b]), gdzie \displaystyle r\ge 0, będzie klasą funkcji

\displaystyle F^r_M([a,b])\,=\,\{\,f\in C^{(r+1)}([a,b]):\,         \|f^{(r+1)}\|_{ C([a,b])}\le M\,\},

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

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

Załóżmy, że każdą funkcję \displaystyle f\in F^r_M([a,b]) aproksymujemy jej wielomianem interpolacyjnym \displaystyle w_f\in\Pi_r opartym na \displaystyle r+1 węzłach \displaystyle x_0,\ldots,x_r\in [a,b]. Wtedy maksymalny błąd takiej aproksymacji wynosi

\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 \displaystyle f\in F^r_M([a,b]) mamy

\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 \displaystyle v(x)=M\frac{x^{r+1}}{(r+1)!} mamy \displaystyle v\in F^r_M([a,b]) oraz

\displaystyle \|v-w_v\|_{ C([a,b])}\,=\,\frac M{(r+1)!}\cdot    \max_{a\le x\le b}|(x-x_0)\cdots(x-x_r)|,
co kończy dowód. image:End_of_proof.gif

Zjawisko Rungego i dobór węzłów interpolacji

Rozważmy zadanie interpolacji funkcji

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

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

Zjawisko Rungego: interpolacja w  węzłach równoodległych dla
Enlarge
Zjawisko Rungego: interpolacja w \displaystyle N=17 węzłach równoodległych dla \displaystyle 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
Enlarge
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
Enlarge
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 \displaystyle e(F^r_M([a,b]);x_0,\ldots,x_r) w istotny sposób zależy od wyboru węzłów \displaystyle x_j. Naturalne jest więc teraz następujące pytanie: w których punktach \displaystyle x_j przedziału \displaystyle [a,b] należy obliczać wartości funkcji, aby błąd był minimalny? Problem ten sprowadza się oczywiście do minimalizacji wielkości \displaystyle \max_{a\le x\le b}|(x-x_0)\cdots(x-x_r)| względem węzłów \displaystyle x_j.

Twierdzenie O optymalnym doborze węzłów

Błąd aproksymacji w klasie funkcji \displaystyle F^r_M([a,b])(x_0,\ldots,x_r) jest minimalny gdy węzły interpolacji są zadane jako węzły Czebyszewa na \displaystyle (a,b), tzn.

\displaystyle x_j^*\,=\,\frac{b-a}2\cdot          \cos\left(\frac{2j+1}{2r+2}\pi\right)\,+\,            \frac{a+b}2,\qquad 0\le j\le r.

Ponadto, dla optymalnych węzłów \displaystyle x_j^* mamy

\displaystyle e(F_M^r([a,b]);x_0^*,\ldots,x_r^*)\,=\,     \frac{2M}{(r+1)!}\left(\frac{b-a}4\right)^{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 \displaystyle \{T_k\}_{k\ge 0} wielomianów Czebyszewa (pierwszego rodzaju) zdefiniowany jest indukcyjnie jako

\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ę
Enlarge
Pafnutij Lwowicz Czebyszew
Zobacz biografię

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

\displaystyle    T_k(x)\,=\,\cos(k\arccos x).

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

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

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

Kilka pierwszych wielomianów Czebyszewa na odcinku
Enlarge
Kilka pierwszych wielomianów Czebyszewa na odcinku \displaystyle [-1,1]

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

\displaystyle \|T_k\|_{ C([-1,1])}\,=\,\max_{-1\le x\le 1} |T_k(x)|      \,=\,1

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

\displaystyle    y_j\,=\,\cos\Big(\frac jk\pi\Big),\qquad 0\le j\le k,

przy czym \displaystyle T_k(y_j)=(-1)^j.

W końcu, \displaystyle k-ty wielomian Czebyszewa \displaystyle T_k ma dokładnie \displaystyle k pojedynczych zer w \displaystyle [-1,1] równych

\displaystyle z_j\,=\,\cos\Big(\frac{2j+1}{2r}\pi\Big),       \qquad 0\le j\le 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 \displaystyle \overline{\Pi}_k oznaczymy klasę wielomianów stopnia \displaystyle k o współczynniku wiodącym równym \displaystyle 1, tzn.

\displaystyle \overline{\Pi}_k\,=\,\{\,w\in\Pi_k:\,      w(x)=x^k+\cdots\,\}.

Twierdzenie O minimaksie

Niech \displaystyle k\ge 1. W klasie \displaystyle \overline{\Pi}_k minimalną normę jednostajną na przedziale \displaystyle [-1,1] ma wielomian \displaystyle w^*=2^{1-k}T_k, tzn.

\displaystyle \min_{w\in\overline{\Pi}_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.
Enlarge
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 \displaystyle (x-x_0)(x-x_1)\cdots(x-x_r) jest w klasie \displaystyle \overline\Pi_{r+1}. Stąd dla \displaystyle [a,b]=[-1,1] optymalnymi węzłami są zera \displaystyle z_j wielomianu Czebyszewa, przy których

\displaystyle (x-z_0)(x-z_1)\cdots(x-z_r)\,=\,\frac{T_{r+1}(x)}{2^r}.

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

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

Stąd

\displaystyle \|w_{a,b}^*\|_{C([a,b])}\,=\,\Big(\frac{b-a}{2}\Big)^{r+1}     \frac 1{2^r}\,=\,2\,\Big(\frac{b-a}{4}\Big)^{r+1}
i węzły \displaystyle x^*_j są optymalne. image:End_of_proof.gif

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 \displaystyle f\in C[-1,1], wielomian interpolacyjny \displaystyle w_f stopnia co najwyżej \displaystyle n, oparty na węzłach Czebyszewa, spełnia

\displaystyle ||f-w_f||_{C[-1,1]}  \leq \left(2+\frac{2}{\pi}\log(n+1)\right) ||f-w_f^*||_{C[-1,1]}

gdzie \displaystyle w_f^* jest wielomianem stopnia co najwyżej \displaystyle n, najlepiej aproksymującym \displaystyle f w sensie normy jednostajnej.


Jeśli więc \displaystyle n \leq 5, to wielomian oparty na węzłach Czebyszewa jest co najwyżej 3.02 razy, a gdy \displaystyle 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.