\documentclass[10pt,a4paper]{article}

%% Polonizacja
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}
\usepackage{polski}
\usepackage{indentfirst}

%% Odnośniki
\usepackage{varioref}

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts}

%% Tabele
\input tap.tex
\usepackage{pstricks,pst-node}

%% Fragmenty kodu
\usepackage[!labels]{listings}
\lstloadlanguages{C}
\lstset{
  language=C,
  basicstyle=\small,
  commentstyle=\slshape,
  stringstyle=\ttfamily,
  showstringspaces=false,
  fontadjust=true,
  extendedchars=true,
  inputencoding=latin2,
  mathescape=true,
  }

%% Twierdzenia
\newtheorem{thm}{Twierdzenie}[section]
\newtheorem{df}{Definicja}[section]

\newcommand{\n}[1][n]{$#1$\nobreakdash-\hspace{0pt}}
\newcommand{\abs}[1]{\lvert#1\rvert}


\begin{document}

\newcounter{seria}
\setcounter{seria}{3}
\newcommand*{\preparedby}{Jakub Narębski}
\input{naglowek.tex}

\section{Interpolacja wielomianowa}
\label{sec:interp-wielomian}

\textbf{Zadanie interpolacyjne Lagrange'a} polega na znalezieniu dla
danej funkcji $f$ (dla zadanych liczb $f_i$) wielomianu $L_n$ stopnia
nie wyższego niż $n$, którego wartości w~$n+1$ punktach $x_i$ są takie
same jak wartości funkcji $f$: $L_n (x_i) = f(x_i)$ (są równe zadanym
liczbom $f_i$: $L_n (x_i) = f_i$.  Wielomian $n$--tego stopnia jest
jednoznacznie wyznaczony przez wartości w~$n+1$ punktach.



\subsection{Interpolacja liniowa}
\label{sec:liniowa}

Przez dwa punkty $(x_1, y_1)$ i~$(x_2, y_2)$ możemy przeprowadzić
prostą.  Aby obliczyć wartość $y$ przy zadanym $x$ dla punktu leżącego
na tej prostej najlepiej (numerycznie) posłużyć się wzorem:
\begin{equation}
  \label{eq:liniowa}
  y = A y_1 + B y_2, \text{\ gdzie\ }
  A = \frac{x_2 - x}{x_2 - x_1}, B = \frac{x - x_1}{x_2 - x_1}.
\end{equation}
Łatwo sprawdzić, że $A(x=x_1) = 1$, $B(x=x_1) = 0$, natomiast 
$A(x=x_2) = 0$, $B(x=x_2) = 1$.



\subsection{Wzór Lagrange'a}
\label{sec:wzor-Lagrange}

Naturalnym rozszerzeniem wzoru~\eqref{eq:liniowa}
z~podrozdziału~\ref{sec:liniowa} na interpolację liniową do
interpolacji wielomianem dowolnego stopnia jest klasyczny \textbf{wzór
  Lagrange'a}:
\begin{equation}
  \label{eq:Lagrange-wielom}
  \begin{split}
    P(x) &= 
    \frac{(x-x_2)(x-x_3)\dots(x-x_N)}{(x_1-x_2)(x_1-x_3)\dots(x_1-x_N)} y_1 +
    \frac{(x-x_1)(x-x_3)\dots(x-x_N)}{(x_2-x_1)(x_2-x_3)\dots(x_2-x_N)} y_2 \\
    &\quad\null+ \dots +
    \frac{(x-x_1)(x-x_2)\dots(x-x_{N-1})}{(x_N-x_1)(x_N-x_2)\dots(x_N-x_{N-1})} y_N 
  \end{split}
\end{equation}
Powyższy wzór określa wielomian $N-1$ stopnia, taki że $y_1 = P(x_1)$,
$y_2 = P(y_2)$, \ldots, $y_N = P(x_N)$.  Postać ta służy do
udowodnienia twierdzenia o~jednoznaczności zagadnienia interpolacji
wielomianowej (zagadnienia interpolacyjnego Lagrange'a):

\begin{thm}
  \label{thm:Lagrange}
  Zadanie interpolacyjne Lagrange'a ma jednoznaczne rozwiązanie, które
  można przedstawić np. w~postaci:
  \begin{equation}
    \label{eq:Lagrange-tw}
    L_n(x) = \sum_{i=0}^n f(x_i) l_i(x),
  \end{equation}
  gdzie $l_i(x)$ to funkcje pomocnicze zdefiniowane jako
  \begin{equation}
    \label{eq:l_i}
    l_i(x) = \prod_{\substack{j=0\\j\neq i}}^n\frac{x-x_j}{x_i-x_j}, 
    \quad \text{dla\ }i=0, 1, \dots, n.
  \end{equation}
\end{thm}

Wielomian $L_n$ nazywa się wielomianem interpolacyjnym Lagrange'a
$n+1$--go stopnia funkcji $f$ opartym na węzłach $x_0,x_1,\ldots,x_n$.
(\emph{Uwaga:} oznaczenia węzłów są inne niż
w~równaniu~\eqref{eq:Lagrange-wielom}!).

Wzór Lagrange'a~\eqref{eq:Lagrange-wielom} na wielomian interpolacyjny
nie jest zbyt wygodny w~zastosowaniu.  Nie daje żadnego oszacowania
błędu, ponadto jest trudny do zaprogramowania.  Dlatego do
interpolacji wielomianowej w~punkcie używa się raczej przedstawionego
w~następnym podrozdziale \textbf{algorytmu Neville'a} (czasem zwanego
\textit{metodą Aitkena}).



\subsection{Algorytm Neville'a}
\label{sec:Neville}

Zadanie interpolacyjne można początkowo rozpatrzyć dla mniejszej
ilości punktów węzłowych, a~następnie adoptować otrzymane rozwiązania
dla otrzymania właściwego wielomianu interpolacyjnego.

Zdefiniujmy kolejne wielomiany interpolacyjne i~ich wartości w~zadanym
punkcie $x$.  $p_i$ jest wielomianem zerowego stopnia (czyli stałą)
przechodzącą przez punkt $(x_i, y_i)$, a~$p_i = y_i$ jego wartością
w~punkcie $x$.  $p_{i(i+1)}$ będzie wartością w~punkcie $x$ wielomianu
drugiego stopnia przechodzącego przez punkty $(x_i, y_i)$ i~$(x_{i+1},
y_{i+1})$.  Analogicznie $p_{i(i+1)(i+2)}$ itd.  Dla czterech punktów
\(x_0, x_1, x_2, x_3\) wyglądało by to następująco:
\begin{equation*}
\begin{array}{lllll}
x_0\colon & y_0=p_0 &        &         &          \\
          &         & p_{01} &         &          \\
x_1\colon & y_1=p_1 &        & p_{012} &          \\
          &         & p_{12} &         & p_{0123} \\
x_2\colon & y_2=p_2 &        & p_{123} &          \\
          &         & p_{23} &         &          \\
x_3\colon & y_3=p_3 &        &         &       
\end{array}
\end{equation*}
Ogólnie $p_{k \dots l}(x)$ --- wartość, w~punkcie $x$, wielomianu
\n{}tego stopnia przechodzącego przez $l-k+1$ punktów $(x_i,y_i)$,
\((i=k,\dots,l)\).

Wartości wielomianów tego typu spełniają następujący związek
rekurencyjny:
\begin{equation}
  \label{eq:Neville-recurrency}
  p_{i(i+1)\ldots(i+k)} = 
  \frac{(x-x_{i+k})p_{i(i+1) \ldots (i+k-1)} + 
        (x_i-x)p_{(i+1)(i+2) \ldots (i+k)}}{x_i-x_{i+k}},
\end{equation}
gdzie $k$ odpowiada stopniu wielomianu, $k=1,\dots,n$.

Algorytm ten możemy zapisać w~innej notacji: $t_{i+k,k} \equiv
p_{i(i+1)\ldots(i+k)}$.  Pierwszy indeks w~$t_{i,j}$ oznacza ostatni
z~węzłów na którym oparty jest wielomian, zaś drugi oznacza krok
iteracji (stopień wielomianu).  W~tej notacji związek rekurencyjny
przyjmuje postać:
\begin{equation}
  \label{eq:Neville-t_ij}
  \begin{split}
    t_{kl} 
    &=
    \frac{(x_{k-l}-x)t_{k,l-1} - (x_{k}-x)t_{k-1,l-1}}{x_{k-l}-x_{k}} = \\
    &=
    t_{k,l-1} + (t_{k,l-1} - t_{k-1,l-1}) \frac{x - x_k}{x_k - x_{k-l}} \,.
  \end{split}
\end{equation}

Implementując metodę Neville'a warto zauważyć, że można ją
zaimplementować korzystając wyłącznie z~jednego wektora pomocniczego
o~rozmiarze równym ilości punktów na których rozpinamy interpolację,
$\mathcal{O}(m)$.  Nie należy stosować rekurencji, ponieważ skutkuje
ona wielokrotnym wyliczaniem tych samych wartości pomocniczych
i~wymaga więcej pamięci.  Zakładając, że chcemy aby końcowy wynik
interpolacji znalazł się w~tablicy pod indeksem 0, możemy pożądany
rezultat uzyskać na (przynajmniej) dwa sposoby, opisane
poniżej. Załóżmy, że na początku mamy tablicę
\lstinline!w[0..3]! wypełnioną wartościami $p_i = y_i$.
\deftable\aitkenup
\begintable
  \begintableformat
    \right " &\left 
  \endtableformat
  \B"^ \lstinline!w[0]! |            $y_0 = p_0$  | \Rnode{p01}{$p_{01}$} | \Rnode{p012}{$p_{012}$} | \Rnode{p0123}{$p_{0123}$} \E"
  \B"+ \lstinline!w[1]! | $y_1 = \Rnode{p1}{p_1}$ | \Rnode{p12}{$p_{12}$} | \Rnode{p123}{$p_{123}$} \E"
  \B"+ \lstinline!w[2]! | $y_2 = \Rnode{p2}{p_2}$ | \Rnode{p23}{$p_{23}$} \E"
  \B"_ \lstinline!w[3]! | $y_3 = \Rnode{p3}{p_3}$ \E"
  \psset{nodesep=2pt,offset=0pt,linewidth=0.4pt}
  %%\ncline{-}{p1}{p01}
  %%\ncline{-}{p2}{p12} \ncline{-}{p12}{p012}
  %%\ncline{-}{p3}{p23} \ncline{-}{p23}{p123} \ncline{-}{p123}{p0123} 
\endtable
\begin{center}
\centerline{\begintable \aitkenup \endtable}
\end{center}

Pierwszą z~nich jest wyliczanie wartości kolumna po kolumnie, tzn
w~kolejności $p_{01}, p_{12}, p_{23}; p_{012}, p_{123}; p_{0123}$,
korzystając z~faktu że np. po wyliczeniu $p_{01}$ wartość $p_{0}$ nie
będzie już potrzebna.

Drugą jest wyliczanie wartości skosem (po przekątnej), wyliczając
w~wewnętrznej iteracji wszystkie $p_{\ldots{}i}$ kończące się w~danym
punkcie, tzn.  w~kolejności $p_{01}; p_{12}, p_{012}; p_{23}, p_{123},
p_{0123}$.

Aby przybliżenie interpolacyjne miało tendencję do optymalności i~było
szybko zbieżne należy uporządkować punkty w~kolejności rosnącej
odległości od punktu intepolacji, tzn. %
\(
  \abs{x_i - x} \leq \abs{x_j - x}
  \text{ dla }
  i < j
\).
Ulepszona metoda (patrz następny rozdział) nie ma tego ograniczenia.


\subsection{Ulepszenie metody Neville'a (z Numerical Recipies)}
\label{sec:ulepsz-neville}

Równanie rekurencyjne~\eqref{eq:Neville-recurrency} metody Neville'a
można ulepszyć obliczając małe \emph{różnice} między przodkami
(wynikami cząstkowymi dla \n[k]tego stopnia) i~potomkami (wynikami dla
\n[k+1]go stopnia).  Zdefiniujmy
\begin{subequations}
  \label{eq:def:DOWN-UP}
  \begin{align}
    \label{eq:def:DOWN}
    D_{i,m} &\equiv p_{i\dots{}(i+m)} - p_{i\dots{}(i+m-1)}   \\
    \label{eq:def:UP}
    U_{i,m} &\equiv p_{i\dots{}(i+m)} - p_{(i+1)\dots{}(i+m)} 
  \end{align}
\end{subequations}
Z~wzoru~\eqref{eq:Neville-recurrency} łatwo można wyprowadzić wzory
rekurencyjne na $D_{i,m}$ (ang. \emph{\textbf{d}own}) i~$U_{i,m}$
(ang. \emph{\textbf{u}p}):
\begin{subequations}
   \label{eq:DOWN-UP}
  \begin{align}
    \label{eq:DOWN}
    D_{i,m+1} &= \frac{D_{i+1,m} - U_{i,m}}{x_i - x_{i+m+1}} (x_{i} - x)  \\
    \label{eq:UP}
    U_{i,m+1} &= \frac{D_{i+1,m} - U_{i,m}}{x_i - x_{i+m+1}} (x_{i+m+1} - x)
  \end{align}
\end{subequations}
Na każdym poziomie $m$, wielkości $D_{i,m}$ i~$U_{i,m}$ są poprawkami
potrzebnymi do interpolacji jeden rząd wyżej.  $p_{0\dots{}n}$ będące
końcowym wynikiem otrzymujemy biorąc sumę \emph{dowolnego} $y_i$ oraz
zestawu $D$ i/lub~$U$ które tworzą ścieżkę przez drzewo od wybranego
$y_i$ aż do $p_{0\dots{}n}$.

W~celu minimalizacji błędów najlepiej zacząć od $y_i$ należącego do
pary $(x_i, y_i)$ dla $x_i$ najbliższego punktowi w~którym
interpolujemy.  W~implementacji z~,,Numerical Recipies'' prowadzimy
ścieżkę jak najbardziej środkiem drzewa, wybierając ścieżkę w~dół gdy
jestesmy powyżej środka i~ścieżkę w~górę gdy jesteśmy poniżej środka.
Ostatnia poprawka może służyć jako oszacowanie błędu interpolacji. 

\emph{Uwaga:} notacja jest nieco zmieniona w~porównaniu z~,,Numerical
Recipies in C''.


\section{Znajdowanie punktów, na których chcemy oprzeć interpolację}
\label{sec:punkty}

Często jest tak, że funkcję mamy zadaną za pomocą wielu punktów,
a~interpolację w~zadanym punkcie $x$ chcemy oprzeć tylko o~kilka
najbliższych $x$ punktów.  Przykładem może być interpolacja
wielomianowa: wielomiany wyższego stopnia zazwyczaj silnie oscyluja
i~moga dawać błędne wyniki interpolacji.

Zakładamy, że punkty w~których mamy stablicowaną funkcję są ułożone
rosnąco, t.j. $-\inf < x_0 \leq x_1 \leq \dots \leq x_{n-1} < +\int$.
Mamy zadaną wartość $x$, chcemy znaleźć taki indeks $i$, że $x_i \leq
x \leq x_{i+1}$.  Wartość $i = -1$ oznacza $x_{-1} = -\infty$, wartość
$i = n$ oznacza $x_n = +\infty$.

\begin{itemize}
\item Metoda \textbf{bisekcji}.  W~każdym kroku połowimy przedział
  w~którym musi znajdować się $x$.
  
\item Metoda \textbf{polowania}.  Zaczynamy z~jakiegoś przybliżenia na
  indeks $i^{\ast}$, następnie poruszamy się w~krokach będących
  kolejnymi potęgami dwójki (1, 2, 4, itp.) aż ,,osaczymy'' $x$.
  Następnie wykonujemy bisekcję zaczynając od znalezionego przedziału.
\end{itemize}

Aby przekazać funkcji interpolującej (niech nazywa się ona
\lstinline!polint!) odpowiedni fragment tablicy o~długości $m$ (chcemy
oprzeć interpolację w~punkcie $x$ na $m$ najbliższych mu punktach),
należy znaleźć początek (początkowy indeks) fragmentu oraz jego
długość.  Zakładamy, że znaleźliśmy takie $i$, że $x_i \leq x \leq
x_{i+1}$, tablica jest numerowana od 0 i~ma $n$ elementów.  Początkowy
indeks znajdujemy za pomocą wzoru:
\begin{equation}
  \label{eq:k}
  k = \min(\max(i - \frac{m}{2}, 1), n - m).
\end{equation}
W~ten sposób zabezpieczamy się przed przypadkiem $k < 0$ oraz $k + m >
n$ (wyjściem poza tablicę).  Długość fragmentu to oczywiście 
$\min(m, n-k)$.  Naszą funkcję wywołujemy zatem z~parametrami
\begin{center}
  \lstinline!polint(tab + k, MIN(m, n-k), \dots)!,
\end{center}
gdzie \lstinline!MIN(a, b)! możemy zdefiniować jako
\lstinline!#define MIN(a,b) ((a)>(b)?(a):(b))!

\medskip

\textbf{Zadanie 1.}  Napisać procedurę znajdującą dla zadanej wartości
$x$ taki indeks $i$ w~uporządkowanej rosnąco tablicy \texttt{xx}
(\texttt{xx[0]} $<$ \texttt{xx[1]} $<$ \ldots $<$ \texttt{xx[N-1]},
gdzie \texttt{N} jest rozmiarem tablicy), że $\text{\texttt{xx[i-1]}}
< x \leq \text{\texttt{xx[i]}}$ metodą \textbf{bisekcji}.  Jeśli $x$
jest mniejsze od najmniejszego elementu tablicy (czyli od
\texttt{xx[0]}) należy zwrócić wartość \texttt{0}, jeżeli natomiast
$x$ jest większe od największego elementu tablicy (czyli od
\texttt{xx[N-1]}) należy zwrócić \texttt{N}.

\textbf{Zadanie 2.}  Napisać procedurę, znajdującą dla zadanego
przybliżenia $i^{\ast}$ i~wartości $x$ taki indeks $i$
w~uporządkowanej rosnąco tablicy \texttt{xx}, że
$\text{\texttt{xx[i-1]}} < x \leq \text{\texttt{xx[i]}}$.  Należy
zastosować algorytm \textbf{polowania} zaczynając od indeksu
$i^{\ast}$.  Można skorzystać z~napisanej poprzednio procedury
(funkcji) stosującej algorytm bisekcji.

\section{Zadanie domowe}
\label{sec:zadanie}
\noindent \emph{Rok akademicki 2004/05.}

Napisać program dokonujący interpolacji wielomianem zadanego stopnia za pomocą
metodą Neville'a.  Przy funkcji zadanej na $N >= m + 1$
punktach (gdzie $m$ jest stopniem wielomianu interpolacyjnego),
wyznaczać wartość wielomianu interpolacyjnego w punkcie na podstawie
interpolacji opartej o najbliższe (w sensie odległości albo w sensie
indeksów) $m+1$ punktów.  

Przykładowe funkcje do testowania: 
\begin{align*}
  f(x) &= |x| \,,\\
  f(x) &= \sin(\pi x) \,,\\
  f(x) &= \frac{1}{1 + x^2} \,,\\
  f(x) &= \exp\bigl(-x^2\bigr) \,,\\ 
  f(x) &= \exp\left(\frac{-1}{x^2}\right) \,,\\
  f(x) &= a_0 + a_1 x + a_2 x^2 + \dots + a_n x^n \quad\text{wielomian} \,.
\end{align*}


\noindent \textbf{Wymagania}: 
\begin{enumerate}
\item 
  Napisanie funkcji (procedury, metody) znajdującej wartość interpolacji wielomianowej
  wielomianem określonego stopnia $m$, rozpiętej na $m+1$ punktach, w~zadanym
  punkcie $x$ (70\%).  Za użycie rekurencji bądź macierzy pomocniczej -10\%.
  Za użycie algorytmu poprawek z Numerical Recipies +15\%.

\item
  Napisanie kodu generującego wykres interpolacji wielomianowej, zadanych
  punktów i~porównanie interpolacji z interpolowaną (prawdziwą) funkcją
  (10\%).  Program może zapisywać wyniki do pliku (punkty interpolacji
  oraz wynik interpolacji), a~wykres może być robiony zewnętrznym
  programem np. \texttt{gnuplot}. 

\item 
  Napisanie funkcji znajdującej punkty na których należy rozpiąć
  interpolację, przy zadanym punkcie interpolacji (i~w~przypadku metody
  polowania zgadniętego indeksu) oraz stopniu wielomianu interpolującego.  Użyć
  funkcji pomocniczej do znalezienia indeksu ``środkowego'', tzn
  takiego że \lstinline!xx[i] <= x < xx[i+1]! (20\%).
\end{enumerate}
%%\input{podpis.tex}
\vfill
\noindent
Przygotował: \preparedby

\end{document}
