MN09

Interpolacja wielomianowa

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.

• Na podstawie próbki sygnału dźwiękowego (to znaczy: ciągu wartości

amplitud sygnału zmierzonych w kolejnych odstępach czasu) odtwarzanie jego przebiegu.

• Przybliżanie wykresu 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.

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\le j\le n.$

Oznaczmy przez ${\mathrm{\Pi }}_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.

Plik:Lagrange.jpg

Lagrange
Zobacz biografię

Twierdzenie Istnienie i jednoznaczność 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 {\mathrm{\Pi }}_n$ interpolujący $f$ w węzłach $x_j$, $0\le j\le n$.

Dowód

Wybierzmy w ${\mathrm{\Pi }}_n$ dowolną bazę wielomianów ${\phi }_j$, $0\le j\le n$,

${\mathrm{\Pi }}_n=\text{span}\{{\phi }_0,{\phi }_1,\dots ,{\phi }_n\}.$

Wtedy każdy wielomian z ${\mathrm{\Pi }}_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{\phi }_j(\cdot )$ interpolował $f$ jest spełnienie układu $n+1$ równań liniowych

$\sum _{j=0}^n{c}_j{\phi }_j(x_i)=f(x_i),0\le i\le n,$

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

$\left(\begin{array}{cccc}{\phi }_0(x_0)& {\phi }_1(x_0)& \cdots & {\phi }_n(x_0)\\ {\phi }_0(x_1)& {\phi }_1(x_1)& \cdots & {\phi }_n(x_1)\\ ⋮\\ {\phi }_0(x_n)& {\phi }_1(x_n)& \cdots & {\phi }_n(x_n)\end{array}\right)\left(\begin{array}{c}{c}_0\\ c_1\\ ⋮\\ c_n\end{array}\right)=\left(\begin{array}{c}f(x_0)\\ f(x_1)\\ ⋮\\ 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, $f(x_i)=0$, $\mathrm{\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 $({\phi }_j{)}_{j=0}^n$ będzie bazą w przestrzeni ${\mathrm{\Pi }}_n$ wielomianów stopnia co najwyżej $n$. Zadanie interpolacji wielomianowej polega na obliczeniu dla danej funkcji $f$ współ\-czyn\-ni\-ków $c_j$ takich, że wielomian

$w_f(\cdot )=\sum _{j=0}^n{c}_j{\phi }_j(\cdot )$

interpoluje $f$ w punktach $x_j$, $0\le j\le n$.

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żemy pokazać, że jeśli zamiast $f$ rozpatrzyć jej zaburzenie $f+\mathrm{\Delta }f$, gdzie $|\mathrm{\Delta }f|\le ϵ$, to

$|w_f(x)-w_{f+\mathrm{\Delta }f}(x)|\le \text{cond(x,f)}|w_f(x)|ϵ,$

gdzie

$\text{cond(x,f)}=\frac{\sum _{j=0}^n|l_j(x)f(x_j)|}{|p_n(x)|}\ge 1.$

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 $({\phi }_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}_{\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 naturalny sposób powstaje więc problem wyboru "wygodnej" bazy w ${\mathrm{\Pi }}_n$. Rozpatrzymy trzy bazy: Lagrange'a, potęgową i Newtona.

Baza Lagrange'a (kanoniczna)

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

$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 $l_j$ jest stopnia dokładnie $n$ oraz

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

Teraz widać, że wielomiany te stanowią bazę w ${\mathrm{\Pi }}_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)$, $\mathrm{\forall }j$. Wielomian interpolacyjny dla funkcji $f$ można więc zapisać jako

$w_f(\cdot )=\sum _{j=0}^{n}f(x_j)l_j(\cdot ).$

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

Przypuśćmy, że chcielibyśmy obliczyć wartość wielomianu interpolacyjnego $w_f$ w punkcie $x$ różnym od $x_j$, $0\le j\le n$. Podstawiając

$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 $p_n(x)=(x-x_0)\cdots (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 można łatwo zobaczyć 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\left(\begin{array}{c}n\\ j\end{array}\right).$

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

Można pokazać, że wartość $\widetilde{w_f(x)}$ wielomianu iterpolacyjnego obliczona w arytmetyce $f{l}_{\nu }$ według pierwszego algorytmu barycentrycznego spełnia

$\widetilde{w_f(x)}=p_n(x)\sum _{j=0}^n\frac{w_j}{x-x_j}f(x_j)(1+{ϵ}_j),$

gdzie $|{ϵ}_j|\le 5(n+1)$, a więc jest to algorytm numerycznie poprawny. Zachowanie drugiej postaci wzoru barycentrycznego w arytmetyce $f{l}_{\nu }$ jest nieco bardziej skomplikowane w typowych zadaniach.

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, ${\phi }_j(x)=x^j$, $\mathrm{\forall }j$. Jeśli bowiem

$w_f(x)=a_0+a_{1}x+\cdots +a_n{x}^n,$

to również

$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 $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 znając współczynniki wielomianu obliczający jego dokładną wartość 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 {\mathrm{\Pi }}_j$ jest wielomianem interpolacyjnym dla funkcji $f$ opartym na węzłach $x_0,x_1,\dots ,x_j$, $0\le j\le n$, to $w_{f,0}=b_0$ oraz

$w_{f,j}=w_{f,j-1}+b_j{p}_j,1\le j\le 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\ge 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\le j\le n.$

Dowód

Dla $0\le i\le j\le n$, oznaczmy przez $w_{i,j}$ wielomian z ${\mathrm{\Pi }}_{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},\mathrm{\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\le s\le j$. Rzeczywiście, jeśli $i+1\le s\le 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 ${\mathrm{\Pi }}_{j-i}$ interpolującym $f$ w węzłach $x_s$, $i\le s\le 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\ge 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\le j\le 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żemy łatwo obliczyć na podstawie wartości $f(x_j)$, $0\le j\le n$, budując następującą tabelkę:

$\begin{array}{ll}{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)\\ ⋮& ⋮& ⋮& ⋮& \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,\dots ,x_n).\end{array}$

<flash>file=Interpolacja.swf</flash><div.thumbcaption>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\le i<j\le 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</flash><div.thumbcaption>Wyznaczenie tego samego wielomianu $w$, interpolującego zestaw punktów $(0,2)(1,5)(-1,7)$ algorytmem różnic dzielonych --- wykonanym in situ.

Okazuje się, że przy realizacji w $f{l}_{\nu }$ algorytmu różnic dzielonych istotną rolę odgrywa porządek węzłów. Można pokazać, że algorytm liczenia $f(t_0,\dots ,t_n)$ jest numerycznie poprawny ze względu na dane interpolacyjne $f^{(i)}(t_j)$, o ile węzły są uporządkowane nierosnąco lub niemalejąco.

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\le j\le k$, przy czym $\sum _{j=0}^k{n}_j=n+1$. Należy skonstruować wielomian $w_f\in {\mathrm{\Pi }}_n$ taki, że

$w_f^{(i)}(x_j)=f^{(i)}(x_j)\text{ dla }0\le i\le n_j-1,0\le j\le k.$

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

Lemat

Zadanie interpolacji Hermite'a ma jednoznaczne rozwiązanie.

Plik:Hermite.jpg

Hermite
Zobacz biografię

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

$({\overline{x}}_0,{\overline{x}}_1,\dots ,{\overline{x}}_n)=(\underset{n_0}{\underbrace{x_0,\dots ,x_0}},\underset{n_1}{\underbrace{x_1,\dots ,x_1}},\dots ,\underset{n_k}{\underbrace{x_k,\dots ,x_k}})$

i zdefiniujemy uogólnioną bazę Newtona w ${\mathrm{\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

$f({\overline{x}}_i,{\overline{x}}_{i+1},\dots ,{\overline{x}}_j)=\frac{f^{(j-i)}({\overline{x}}_i)}{(j-i)!}$

dla ${\overline{x}}_i={\overline{x}}_{i+1}=\cdots ={\overline{x}}_j$, oraz

$f({\overline{x}}_i,{\overline{x}}_{i+1},\dots ,{\overline{x}}_j)=\frac{f({\overline{x}}_{i+1},\dots ,{\overline{x}}_j)-f({\overline{x}}_i,\dots ,x_{j-1})}{{\overline{x}}_j-{\overline{x}}_i}$

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

Twierdzenie

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({\overline{x}}_0,{\overline{x}}_1,\dots ,{\overline{x}}_j),0\le j\le n.$

Dowód

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

Zauważmy najpierw, że dla ${\overline{x}}_i\ne {\overline{x}}_j$ zachodzi znany nam już wzór,

$w_{i,j}(x)=\frac{(x-{\overline{x}}_i)w_{i+1,j}(x)-(x-{\overline{x}}_j)w_{i,j-1}(x)}{{\overline{x}}_j-{\overline{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 ${\overline{x}}_i,\dots {\overline{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-{\overline{x}}_i)w_{i+1,j}^{(k)}(x_s)-(x_s-{\overline{x}}_j)w_{i,j-1}^{(k)}(x_s)}{{\overline{x}}_j-{\overline{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\ge 1$ wystarczy pokazać, że $b_n=f({\overline{x}}_0,{\overline{x}}_1,\dots ,{\overline{x}}_n)$. W tym celu rozpatrzymy dwa przypadki.

Jeśli ${\overline{x}}_0={\overline{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}{\underbrace{x_0,\dots ,x_0}})$. Jeśli zaś ${\overline{x}}_0\ne {\overline{x}}_j$, to równość $b_n=f({\overline{x}}_0,{\overline{x}}_1,\dots ,{\overline{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 {\mathrm{\Pi }}_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 (przybliżaniu) 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.

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 $\overline{x}\in D$ błąd interpolacji w $\overline{x}$ wyraża się wzorem

$f(\overline{x})-w_f(\overline{x})=(\overline{x}-x_0)(\overline{x}-x_1)\cdots (\overline{x}-x_n)f(x_0,x_1,\dots ,x_n,\overline{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(\overline{x})-w_f(\overline{x})=(\overline{x}-x_0)(\overline{x}-x_1)\cdots (\overline{x}-x_n)\frac{f^{(n+1)}(\xi )}{(n+1)!},$

gdzie $\xi =\xi (\overline{x})$ jest pewnym punktem należącym do najmniejszego przedziału zawierającego punkty $x_0,x_1,\dots ,x_n,\overline{x}$.

Dowód

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

${\overline{w}}_f(x)=w_f(x)+(x-x_0)(x-x_1)\cdots (x-x_n)f(x_0,x_1,\dots ,x_n,\overline{x}),$

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

$0={\psi }^{(n+1)}(\xi )=f^{(n+1)}(\xi )-(n+1)!f(x_0,\dots ,x_n,\overline{x}).$

Stąd

$f(x_0,x_1,\dots ,x_n,\overline{x})=\frac{f^{(n+1)}(\xi )}{(n+1)!},$ co kończy dowód.

Zwykle interesuje nas nie tyle błąd w ustalonym punkcie $\overline{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 $-\mathrm{\infty }<a<b<+\mathrm{\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

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

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

$F_M^r([a,b])=\{f\in C^{(r+1)}([a,b]):\Vert f^{(r+1)}{\Vert }_{C([a,b])}\le M\},$

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

Twierdzenie

Załóżmy, że każdą funkcję $f\in F_M^r([a,b])$ aproksymujemy jej wielomianem interpolacyjnym $w_f\in {\mathrm{\Pi }}_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{x}^{r+1}//(r+1)!$ mamy $v\in F_M^r([a,b])$ oraz

$\Vert v-w_v{\Vert }_{C([a,b])}=\frac{M}{(r+1)!}\cdot \underset{a\le x\le b}{\max }|(x-x_0)\cdots (x-x_r)|,$ co kończy dowód.

Przykład: 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

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 $\underset{a\le x\le b}{\max }|(x-x_0)\cdots (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,\cdots ,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 \left(\frac{2j+1}{2r+2}\pi \right)+\frac{a+b}{2},0\le j\le r.$

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

$e(F_M^r([a,b]);x_0^{*},\dots ,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 $\{T_k{\}}_{k\ge 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}

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\ge 1$). Ponadto wielomian $T_k$ można dla $|x|\le 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\le t\le \pi$, oraz wzór na sumę cosinusów otrzymujemy dla $k\ge 1$

$\cos ((k+1)t)=2\cdot \cos t\cos (kt)-\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

$\Vert T_k{\Vert }_{C([-1,1])}=\underset{-1\le x\le 1}{\max }|T_k(x)|=1$

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

$y_j=\cos \left(\frac{j}{k}\pi \right),0\le j\le 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 \left(\frac{2j+1}{2r}\pi \right),0\le j\le k-1.$

Konsekwencją wymienionych własności jest następująca własność ekstremalna wielomianów Czebyszewa.

Przez ${\bar{\mathrm{\Pi }}}_k$ oznaczymy klasę wielomianów stopnia $k$ o współczynniku wiodącym równym $1$, tzn.

${\bar{\mathrm{\Pi }}}_k=\{w\in {\mathrm{\Pi }}_k:w(x)=x^k+\cdots \}.$

Twierdzenie o minimaksie

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

$\underset{w\in {\bar{\mathrm{\Pi }}}_k}{\min }\Vert w{\Vert }_{C([-1,1])}=\Vert w^{*}{\Vert }_{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.

Miejsca zerowe wielomianu Czebyszewa będziemy nazywać węzłami Czebyszewa.

Dowód Twierdzenia o optymalnym doborze węzłów

Dowód wynika teraz bezpośrednio z twierdzenia o minimaksie. Zauważmy bowiem, że wielomian $(x-x_0)(x-x_1)\cdots (x-x_r)$ jest w klasie ${\bar{\mathrm{\Pi }}}_{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)\cdots (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 ${\bar{\mathrm{\Pi }}}_{r+1}$ minimalną normę Czebyszewa na przedziale $[a,b]$ ma wielomian

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

Stąd

$\Vert w_{a,b}^{*}{\Vert }_{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

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

$||f-w_f|{|}_{C[-1,1]}\le (2+\frac{2}{\pi }\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\le 5$, to wielomian oparty na węzłach Czebyszewa jest co najwyżej 3.02 razy, a gdy $n\le 20$, maksymalnie 4 razy gorszy od optymalnego. Można więc powiedzieć, że jest prawie optymalny.