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

Niech $D\subset R$ i niech $F$ będzie pewnym zbiorem funkcji $f:D\to R$. 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,

${\mathrm{\Pi }}_n=\{w(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_{1}x+a_0:;a_j\in R,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 $f:D\to R$ 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ółczynnikó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$.

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 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 ${\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(x)=\sum _{j=0}^{n}f(x_j)l_j(x)$

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

Wzory barycentryczne

Przypuśćmy, że chcielibyśmy obliczyć wartość wielomianu interpolacyjnego $w_f$ w punkcie $x$ różnym od $x_j$, $0\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 łatwo widzieć, rozpatrując zadanie interpolacji funkcji $f\equiv 1$. Drugi wzór barycentryczny jest korzystniejszy w implementacji.

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

$w_j=(-1{)}^j\left(\begin{array}{c}n\\ j\end{array}\right)$

Również dla 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 wzoru 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.

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

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

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

Baza Newtona

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

$\begin{array}{rl}{p}_0(x)& =1,\\ p_j(x)& =(x-x_0)(x-x_1)\cdots (x-x_{j-1}),1\le j\le n.\end{array}$

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>v_n = b_n;</math>
for (j=n-1; j >= 0 ; j--)
	<math>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żna ł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|width=550|height=300</flash>

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

Zauważmy przy tym, że "po drodze" obliczamy $f(x_i,x_{i+1},\dots ,x_j)$ dla wszystkich $0\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>b_j</math> = <math>f(x_j)</math>; 
for (j = 1; j <= n; j++)
	for (k = n; k >= j; k--)
		<math>b_k</math> = <math>(b_k-b_{k-1})/(x_k - x_{k-j})</math>;

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

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

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

Okazuje się, że przy realizacji w $f{l}_{\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 $f(t_0,\dots ,t_n)$ jest numerycznie poprawny ze względu na dane interpolacyjne $f(t_j)$, a cały algorytm różnic dzielonych daje w arytmetyce $f{l}_{\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 $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$

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

Przykład

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

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

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

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

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

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

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

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

Y = polyval(c,X);

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

Przykład

Interpolujemy tabelkę

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

wielomianem stopnia co najwyżej 2.

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

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

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

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

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

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

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

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

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

Przypadek węzłów wielokrotnych

Uogólnieniem rozpatrzonego zadania interpolacji jest zadanie interpolacji Hermite'a. Zakładamy, że oprócz (różnych) węzłów $x_j$ dane są również ich krotności $n_j$, $0\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.

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

$\begin{array}{rl}{p}_0(x)& =1,\\ p_j(x)& =(x-{\overline{x}}_0)(x-{\overline{x}}_1)\cdots (x-{\overline{x}}_{j-1}),1\le j\le n.\end{array}$

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 O różnicach dzielonych dla interpolacji Hermite'a

Współczynniki $b_j$ wielomianu interpolacyjnego Hermite'a w bazie Newtona,

$w_f(\cdot )=\sum _{j=0}^n{b}_j{p}_j(\cdot )$,

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

$b_j=f({\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ż

$\begin{array}{rl}{v}^{(k)}(x)& =\frac{k(w_{i+1,j}^{(k-1)}(x)-w_{i,j-1}^{(k-1)}(x))}{{\overline{x}}_j-{\overline{x}}_i}\\ & & +\frac{(x-{\overline{x}}_i)w_{i+1,j}^{(k)}(x)-(x-{\overline{x}}_j)w_{i,j-1}^{(k)}(x)}{{\overline{x}}_j-{\overline{x}}_i},\end{array}$

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 funkcji. Funkcję musimy również aproksymać wtedy, gdy nie jesteśmy w stanie uzyskać pełnej o niej informacji. Na przykład, gdy funkcja reprezentuje pewien proces fizyczny, często zdarza się, że dysponujemy jedynie ciągiem próbek, czyli wartościami tej funkcji w pewnych punktach. Jasne jest, że chcielibyśmy przy tym, aby błąd aproksymacji był możliwie mały.

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

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

Niech $x_0,x_1,\dots ,x_n$ będą (niekoniecznie różnymi) węzłami należącymi do pewnego (być może nieskończonego) przedziału $D\subset R$. Dla danej funkcji $f:D\to R$, 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ę $\psi :D\to R$,

$\begin{array}{r}\psi (x)=f(x)-{\overline{w}}_f(x)\\ & =f(x)-w_f(x)-(x-x_0)(x-x_1)\cdots (x-x_n)f(x_0,\dots ,x_n,\overline{x}).\end{array}$

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 $g:[a,b]\to R$, 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 O najgorszym możliwym błędzie interpolacji w klasie

Załóżmy, że każdą funkcję $f\in F_M^r([a,b])$ aproksymujemy jej wielomianem interpolacyjnym $w_f\in {\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

$\begin{array}{rl}e(F_M^r([a,b]);x_0,x_1,\dots ,x_r)& =\underset{f\in F_M^r([a,b])}{\max }\Vert f-w_f{\Vert }_{C([a,b])}\\ & =\frac{M}{(r+1)!}\cdot \underset{a\le x\le b}{\max }|(x-x_0)\cdots (x-x_r)|.\end{array}$

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

$\begin{array}{rl}\Vert f-w_f{\Vert }_{C([a,b])}& =\underset{a\le x\le b}{\max }|f(x)-w_f(x)|\\ & =\underset{a\le x\le b}{\max }|(x-x_0)\cdots (x-x_r)|\frac{|f^{(r+1)}(\xi (x))|}{(r+1)!}\\ & \le & \frac{M}{(r+1)!}\underset{x\in D}{\max }|(x-x_0)\cdots (x-x_r)|.\end{array}$

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

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

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

Rozważmy zadanie interpolacji funkcji

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

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

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

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

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

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

Zjawisko Rungego: interpolacja w węzłach Czebyszewa

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

Zauważmy, że błąd aproksymacji $e(F_M^r([a,b]);x_0,\dots ,x_r)$ w istotny sposób zależy od wyboru węzłów $x_j$. Naturalne jest więc teraz następujące pytanie: w których punktach $x_j$ przedziału $[a,b]$ należy obliczać wartości funkcji, aby błąd był minimalny? Problem ten sprowadza się oczywiście do minimalizacji wielkości $\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,\dots ,x_r)$ jest minimalny gdy węzły interpolacji są zadane jako węzły Czebyszewa na $(a,b)$, tzn.

$x_j^{*}=\frac{b-a}{2}\cdot \cos \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

$\begin{array}{rl}{T}_0(x)& =1,\\ T_1(x)& =x,\\ T_{k+1}(x)& =2x{T}_k(x)-T_{k-1}(x),\text{ dla }k\ge 1.\end{array}$

Pafnutij Lwowicz Czebyszew
Zobacz biografię

Zauważmy, że $T_k$ jest wielomianem stopnia dokładnie $k$ o współczynniku przy $x^k$ równym $2^{k-1}$ ($k\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$

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

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

Dowód

Dowód wynika teraz bezpośrednio z twierdzenia o minimaksie. Zauważmy bowiem, że wielomian $(x-x_0)(x-x_1)\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, o prawie optymalnej interpolacji w węzłach Czebyszewa

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

$||f-w_f|{|}_{C[-1,1]}\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.

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.