Laboratorium 2

Laboratorium 2: Interpolacja

Wstęp

Rozpatrzmy takie zadanie: Mamy zestaw par punktów \(\{x_i, y_i\}, i=0,1,2,\ldots,N\). Mogą pochodzić one z pomiarów jakiegoś zjawiska (z eksperymentu) albo z próbkowania jakiejś nieznanej funkcji \(f(x)\) \((y_i=f(x_i)\)). Nie znamy zależności pomiędzy \(x\)\(y\) i poszukujemy takiej funkcji \(w(x)\), że: \(w(x_i)=y_i\).

Co więcej bardzo często żądamy dodatkowo, żeby funkcja pozwalała na łatwe manipulacje (dodawanie, mnożenie funkcji, różniczkowanie), a wyliczanie jej wartości nie było zbyt kosztowne.

Czasami zadanie stawiane jest niej ambitnie — mamy zestaw \(\{x_i, y_i\}, i=0,1,2,\ldots,N\), a interesują nas wartości funkcji „między” tymi punktami.

W przypadku, gdy \(A\le x_1 \le x_2 \le \cdots \le x_n \le B\), a interesują nas wartści funkcji dla argumentów \(x<A\) lub \(x>B\) (czyli spoza zakresu pomiarów) zadanie będziemy nazywać ekstrapolacją.

Wielomiany

Bardzo często jako funkcję \(w(x)\) wybierane są wielomiany. Łatwo je różniczkować, stosunkowo łatwo obliczać wartość a i zadanie znalezienia współczynników wielomianu interpolacyjnego jest stosunkowo proste. \[w_n(x)=\sum_{i=0}^{N} a_i x^i\]

Do wyliczania wartości takiego wielomianu najlepiej stosować schemat Hornera.

Interpolacja Lagrange’a

Interpolacja Lagrange’a polega na znalezieniu dla danego zestawu danych \(\{x_i, y_i\}, i=0,1,2,\ldots,N\) wielomianu \(W_n\) stopnia nie wyższego niż \(n\), którego wartości w \(n+1\) punktach \(x_i\) są takie same jak wartości interpolowanej funkcji, tzn.: \[W_n(x_i)=y_i, \quad \text{dla } i=0,1,2,\ldots,N\] przy czym zakładamy, że gdy \(i\ne j\) to \(x_i \ne x_j\).

Można pokazać, że tak postawione zadanie interpolacyjne ma jednoznaczne rozwiązanie, które można przedstawić w postaci: \[W_n(x)= \sum_{i=0}^N y_i l_i(x)\] gdzie funkcja pomocnicza \(l_i(x_j)=\delta_{ij}= \begin{cases} 1& \text{dla } i=j\\ 0& \text{dla } i\ne j \end{cases}, \) w szczególności: \[l_i(x)=\prod_{\substack{j=0\\j\ne i}}^N \frac{x-x_j}{x_i-x_j}\] Można pokazać że jest to jednoznaczne rozwiązanie problemu interpolacji Lagrange’a.

Niestety, postać Lagrange’a wielomianu interpolacyjnego jest bardzo niewygodna do prowadzenia jakichkolwiek obliczeń. Zazwyczaj wylicza się współczynniki wielomianu klasycznego i używa ich w obliczeniach.

Aby wyliczyć współczynniki wielomianu należy rozwiązać układ \(n+1\) równań liniowych: \[a_0+a_1x_+a_2 x_i^2+\cdots +a_N x_i^N = y_i, \quad i=0,1,\ldots,N\] Z względu na specyficzną postać współczynników (\(1, x, x^2,\ldots,x^N)\) — są one zależne, co może być problemem w przypadku równoodległych węzłów) najlepiej stosować specjalne metody rozwiązywania tego układu równań.

Jedna z nich, oparta na wyliczaniu ilorazów różnicowych opisana zostanie pokrótce niżej. Wspominam o niej z „konikarskiego obowiązku” — to właśnie do wyliczania wartości wielomianów interpolacyjnych Babbage budował swój pierwszy „komputer:” Maszynę różnicową.

Wielomian interpolacyjny \(W_n(x)\) zapisywany jest w alternatywnej postaci: \[W_N(x) = \sum_{i=0}^N b_k p_k(x)\] Wielomiany \(p_k(x)\) opisane są wzorem: \[\begin{split} p_0(x) & = 1\\ p_k(x) & = (x-x_0)(x-x_1)\ldots(x-x_{k-1}), \quad k=1,2,\ldots, N \end{split}\] Współczynniki \(b_k\) dane są wzorem: \[b_k=\sum_{i=0}^k\frac{y_i}{\displaystyle\prod_{\substack{j=0\\j \ne i}}^k(x_i-x_j)}\]

Współczynniki \(b_k\) noszą nazwę ilorazów różnicowych funkcji \(f\) (pamiętamy, że \(y_i=f(x_i)\)) opartej na węzłach \(x_0, x_1,\ldots,x_k\). Ilorazy różnicowe oznaczać będziemy tak: \(f[x_l,x_{l+1},\ldots,x_{l+k}]\) \[f[x_l,x_{l+1},\ldots,x_{l+k}] = \sum_{i=l}^{l+k}\frac{y_i}{\displaystyle\prod_{\substack{j=l\\j \ne i}}^{l+k}(x_i-x_j)}\]

Można pokazać, że zachodzi następująca zależność rekurencyjna: \[f[x_l,x_{l+1},\ldots,x_{l+k}]=\frac{f[x_{l+1},x_{l+2},\ldots, x_{l+k}] – f[x_{l},x_{l+1},\ldots, x_{l+k-1}]}{x_{l+k}-x_l}\]

Wyliczenie współczynników \(b_k\) sprowadza się do rekurencyjnego wyliczania ilorazów różnicowych funkcji \(f\). Budujemy tablicę zawierającą wartości \(x_i\), wartości \(f(x_i)\) oraz odpowiednie ilorazy różnicowe: \[\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] \\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots\\ x_n & f(x_n) & f[x_{n-1},x_n] & \ldots & \ldots & f[x_0,\ldots,x_n] \end{array}\] Szukane współczynniki \(b_k\) równe są elementom przekątniowym powyższej tablicy.

Powyższe obliczenia można wykonać bardzo łatwo korzystając z następującego algorytmu:

Ponieważ najwygodniejszą (i najefektywniejszą) metodą liczenia wartości wielomianów jest metoda Hornera — trzeba ze korzystając ze współczynników \(b_k\) wyliczyć współczynniki \(a_i\) wielomianu w postaci naturalnej. Jest to również dosyć proste [1]. Realizuje je następujący algorytm:

Naszkicowany powyżej algorytm jest praktycznie najefektywniejszym algorytmem interpolacji. Jedynie w przypadku interpolacji trygonometrycznej i użycia Szybkiej Transformaty Fouriera, dla dużych N będzie ona tańsza od algorytmu ilorazów różnicowych. Zainteresowanie Babbage’a tym algorytmem nie powinno zatem budzić zdziwienia.

Przydaność algorytmu interpolacji pokazuje przykładowy notatnik Mathematici, gdzie na podstawie siedemnastu wartości funkcji sinus, które każdy zna (powinien znać) budowany jest wielomian interpolacyjny. Dla dużych \(n\) (w tym wypadku 16) obliczenia chwilę trwają, ale wyniki są całkiem dokładne.

Funkcje sklejane

Interpolacja za pomocą wielomianów stopnia \(n\) gdy \(n\) jest duże ma szereg efektów ubocznych — gwarantuje co prawda, że \(w_n(x_i) = f(x_i)\), ale nie gwarantuje żadnego „przyzwoitego zachowania pomiędzy węzłami. Im więcej punktów — tym wyższy stopień wielomianu i tym „gorzej” się on zachowuje.

Funkcję rzeczywistą \(S\) nazywamy funkcją sklejaną stopnia \(m\) z węzłąmi \(a = x_0 < x_1 < \ldots < x_N = b\) jeśli

  1. w każdym przedziale \((x_{i-1},x_i)\) dla \(i = 0, 1, \ldots, N + 1\) \((x_{-1} = -\infty, x_{N+1} = \infty)\) \(S\) jest wielomianem stopnia nie wyższego niż \(m\),
  2. \(S\) i jej pochodne rzędu \(1, 2, \ldots, m-1\) są ciągłe na całej osi rzeczywistej \(S \in C^{m-1}\).

Gdy \(m=1\), funkcja sklejana jest po prostu łamaną. Otrzymana funkcja interpolacyjna będzie ciągła, ale pochodna będzie nieciągła, co znaczy, że funkcja nie jest gładka. Można (podwyższjąc \(m\)) doprowadzić do sytuacje, że funkcja jest ciągłą, ma pochodną, która również jest ciągła. Otrzymana funkcja będzie głądka. (Natomiast może okazać się, że druga pochodna już ciągłą nie jest — co ma swoje konsekwencje.) Zazwyczają ciągłość pierwszej, a czasami drgiej pochodnej jest wystarczajaca.

Funkcje sklejane to specjalna kategoria funkcji interpolacyjnych. Są one skonstruowane z „kawałków” (funkcja może być inna w każdym przedziale interpolacji), ale mają szereg zalet. nadają się zwłąszcza gdy trzeba na podstawie siatki (na przykłąd MES) utworzyć gładką powierzchnię.

Szczególnym przypadkiem są krzywe Beziera (Bezier curve).

Mathematica

Interpolacja wielominowa

Do wyliczania wielomianów interpolacyjnych służy funkcja Interpolation. Można jeż użyć albo tak:
\(\pmb{\text{data}=\{1,2,-3,5,8,3\}}\)
\(\pmb{f=\text{Interpolation}[\text{data},\text{InterpolationOrder}\to 4]}\)
W tym przypadku zakłada się, że nasz zestaw danych to pary: \(\{(1,1)\), \((2,2)\), \((3,-3)\), \((4,5)\), \((5,8)\), \((6,3)\}\).

Rząd interpolacji (InterpolationOrder, standardowo 3) jest bardzo ważnym parametrem. Interpolacja odbywa się w sposób opisany wcześniej tylko wtedy gdy rząd interpolacji równa się \(n-1\) (gdzie \(n\) to liczność zbioru danych). Gdy rząd jest mniejszy niż \(n-1\) — interpolacja dokonywana jest za pomocą jakiejś formy funkcji sklejanych. W szczególności można zażądać, żeby użyta była metoda funkcji sklejanych dodając parametr Method\(\to\)“Spline”.

Wynik interpolacji (w tym wypadku \(f\)) jest funkcją i może być używane tak jak każda funkcja: \(\pmb{p=\text{Plot}[f[x],\{x,1,6\}]}\) dając w efekcie piękny wykres, albo tak:
\(\pmb{y=f[2.5]}\)
\(-2.41406\)

image

Albo gdy mamy zestaw danych \(\pmb{\text{pp}}\), który jest tablicą dwuwymiarową zawierającą współrzędne węzłów w postaci \((x,y)\):
\(\pmb{\text{pp}}\) = \(\left\{\{0,1\},\left\{\frac{1}{10},2\right\},\left\{\frac{1}{5},-3\right\},\left\{\frac{3}{10},5\right\},\left\{\frac{2}{5},8\right\},\left\{\frac{1}{2},3\right\}\right\}\)
\(\pmb{\text{ff}=\text{Interpolation}[\text{pp},\text{InterpolationOrder}\to 1];}\)
\(\pmb{\text{Show}[\text{Plot}[\text{ff}[x],\{x,0,0.5\}],\text{ListPlot}[\text{pp},\text{PlotStyle}\to \text{Red}]]}\)

Przykładowy notatnik Mathematici pokazujący zastosowanie interpolacji.

Interpolacja funkcjami sklejanymi

Do interpolacji funkcjami sklejanymi służyć mogą funkcje Interpolate lub ListInterpolate gdy podamy jako metodę interpolacji “Spline”. (W zasadzie nie widzę między nimi różnicy).

Przykładowy notatnik do porównania interpolacji wielomianowej i splajnów.

Matlab

Zakładam, że wszyscy Państwo znakomicie znacie Matlaba. Więc opisy są minimalne. Sugerują zabawy z Mathematicą.

Interpolacja wielominowa

Funkcja interp1 może być użyta do prostych zadań interpolacji. Korzysta z metod: ’nearest’, ’linear’,’spline’,’pchip’, or ’cubic’ i nie interpoluje w sposób opisany w części teoretycznej.

Użycie:

gdzie x tablica współrzędnych x węzłów interpolacji, v — tablica współrzędnych y węzłów, xq tablica zawierająca wartości punktów, w których chcemy wyliczyć wartości interpolowanej funkcji. Zakłada się, że węzły interpolacji są zapisane w kolejności rosnącej.

Istnieje również możliwość interpolacji na podstawie danych, które nie są w żaden sposób uporządkowane: na przykład mamy zestaw punktów w przestrzeni trójwymiarowej i chcemy narysować poziomice. Używa się wówczas metod opisanych jako interpolating scattered data.

Zadania

  1. Przygotować przykład pokazujący zachowanie zachowanie wielomianu interpolacyjnego pomiędzy węzłami interpolacji.
  2. Użyć interpolacji wielomianowej i za pomocą funkcji sklejanych do przebiegu temperatury z czujnika podczerwonego (plik ma w nazwie temperatura_ir). W wyniku powinniśmy dostać dwa różne „przybliżenia” funkcji ale o tych samych wartościach w węzłach interpolacji. Należy policzyć pochodne obu funkcji interpolacyjnych. Wyniki porównać. Wyciągnąć wnioski. (Zamiast interpolacji wielomianowej — jeżeli trudno lub nie da się jej wyliczyć — można użyć jakiegoś innego sposobu interpolacji.)

    Hint: w przypadku Mathematici, aby uzyskać pochodną funkcji \(f\) wystarczy napisać \(f’\).

Sprawozdanie

Uwaga: Te zajęcia powinny skończyć się krótkim sprawozdaniem opisującym zależności między czasem obserwacji, okresem próbkowania i wynikami otrzymanymi ze stosowania szybkiej transformaty Fouriera.

Instrukcja w postaci jednego pliku…

…jest również dostępna.

  • J. M. i Jankowscy, Przegląd metod i algorytmów numerycznych. Cz. 1, Wydawnictwa Naukowo-Techniczne, 1988.
    [BibTeX] [Online]
    @Book{jankowscy-przeglad-1988,
    Title = {Przegląd metod i algorytmów numerycznych. {Cz}. 1},
    Author = {Jankowscy, Janina i Michał},
    Publisher = {Wydawnictwa Naukowo-Techniczne},
    Year = {1988},
    Copyright = {Biblioteka Politechniki Poznańskiej},
    Keywords = {algorytmy, analiza numeryczna, informatyka},
    Language = {pol},
    Owner = {myszka},
    Timestamp = {2015.09.04},
    Url = {http://www.wbc.poznan.pl/dlibra/docmetadata?id=88897&from=&dirids=1&ver_id=&lp=2&QI=}
    }