MN09

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

Niech Parser nie mógł rozpoznać (nieznana funkcja „\subsetR”): {\displaystyle \displaystyle D\subsetR} i niech $\displaystyle F$ będzie pewnym zbiorem funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \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\leq j\leq n.$

Oznaczmy przez $\displaystyle \Pi _{n}$ przestrzeń liniową wielomianów stopnia co najwyżej $\displaystyle 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 $\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\leq j\leq n.$

Dowód

Wybierzmy w $\displaystyle \Pi _{n}$ dowolną bazę wielomianów $\displaystyle \varphi _{j}$, $\displaystyle 0\leq j\leq 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\leq i\leq 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.

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\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 $\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\leq j\leq 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\neq 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\leq j\leq 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_{j}f(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_{1}x+\cdots +a_{n}x^{n},$

to również

$\displaystyle w_{f}(x)\,=\,(\cdots (a_{n}x+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,

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 $\displaystyle b_{j}$,

$\displaystyle w_{f}\,=\,\sum _{j=0}^{n}b_{j}p_{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\leq j\leq n$, to $\displaystyle w_{f,0}=b_{0}$ oraz

$\displaystyle w_{f,j}\,=\,w_{f,j-1}\,+\,b_{j}p_{j},\qquad 1\leq j\leq 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\geq 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\leq j\leq n.$

Dowód

Dla $\displaystyle 0\leq i\leq j\leq 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\leq s\leq j$. Rzeczywiście, jeśli $\displaystyle i+1\leq s\leq 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\leq s\leq 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\geq 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\leq j\leq 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.

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\leq j\leq 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}}$

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

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\leq i<j\leq 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ć?!

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

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ł $\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\leq j\leq 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\leq i\leq n_{j}-1,0\leq j\leq 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$.

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

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

$\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}\neq {\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_{j}p_{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\leq j\leq 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}\neq {\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ż

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

$\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\geq 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}\neq {\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.

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_{t_{0},\ldots ,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 ($\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 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 $\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\leq j\leq 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ę 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 $\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.

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 Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle g:[a,b]\toR} , norma ta jest zdefiniowana jako

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

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

$\displaystyle F_{M}^{r}([a,b])\,=\,\{\,f\in C^{(r+1)}([a,b]):\,\|f^{(r+1)}\|_{C([a,b])}\leq 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_{M}^{r}([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

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

$\displaystyle \|v-w_{v}\|_{C([a,b])}\,=\,{\frac {M}{(r+1)!}}\cdot \max _{a\leq x\leq b}|(x-x_{0})\cdots (x-x_{r})|,$ co kończy dowód.

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

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 $\displaystyle e(F_{M}^{r}([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\leq x\leq 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_{M}^{r}([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\leq j\leq 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\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 $\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\geq 1$). Ponadto wielomian $\displaystyle T_{k}$ można dla $\displaystyle |x|\leq 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\leq t\leq \pi$, oraz wzór na sumę cosinusów otrzymujemy dla $\displaystyle k\geq 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 $\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\leq x\leq 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 {j}{k}}\pi {\Big )},\qquad 0\leq j\leq 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\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 $\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\geq 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.

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.

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.