\documentclass[10pt,a4paper]{article}

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

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts}

%% Włączanie kodu
\usepackage[!labels]{listings}
\lstset{
  language=C,
  tabsize=2,
  %indent=\parindent,
  xleftmargin=\parindent,
  basicstyle=\small,
  commentstyle=\slshape,
  identifierstyle=\itshape,
  stringstyle=\ttfamily,
  fontadjust=true,
  %stringspaces=false,
  showstringspaces=false,
  mathescape=true,
  extendedchars=true,
  }

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

\begin{document}

\newcounter{seria}
\setcounter{seria}{4}
% 1 - wstęp, dokładność numeryczna
% 2 - wielomiany
% 3 - interpolacja wielomianowa
% 4 - współczynniki wielomianu interpolującego
% 5 - splines
\newcommand*{\preparedby}{Jakub Narębski}
\input{naglowek.tex}

\section{Znajdowanie współczynników wielomianu interpolującego}
\label{sec:glowna}

\subsection{Metody bezpośrednie}
\label{sec:bezposrednie}


\subsubsection{Korzystając ze wzoru Lagrange'a}
\label{sec:Lagrange}


Wielomian interpolacyjny dany jest wprost za pomocą wzoru 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}
Korzystanie z~tej postaci do obliczeń nie jest jednak zbyt wygodne.

Współczynniki w~postaci naturalnej wielomianu interpolacyjnego można
wyznaczyć rozwiązując układ równań liniowych
\begin{equation}
  \label{eq:naturalna-wielom}
  a_0+a_1 x_i+ \dots +a_n x_i^n = f(x_i), \quad i=0,1,\dots, n.
\end{equation}
Ten układ równań może być źle uwarunkowany (tzw. macierz
Vandermonde'a), i~trzeba (należy) do jego rozwiązywania stosować
specjalne metody.


\subsubsection{Metodą ilorazów różnicowych w postaci Newtona}
\label{sec:Newton}


Łatwiej jest znaleźć współczynniki w~postaci Newtona wielomianu
interpolacyjnego.  Korzystamy z~zależności rekurencyjnej miedzy
wielomianami składowymi we wzorze Legendre'a i~obliczamy tzw. ilorazy
różnicowe (za pomocą odpowiedniego wzoru rekurencyjnego).  Gdy mamy
współczynniki wielomianu interpolacyjnego w~postaci Newtona, można
przejść do postaci naturalnej.

Niech $W_n(x)$ będzie wielomianem interpolacyjnym Lagrange'a
przechodzącym przez punkty $(x_0, y_0),\dots,(x_n,y_n)$:
\begin{equation}
  \label{eq:Lagrange}
  W_n(x) = \sum_{j=0}^{n} y_j \prod_{\substack{i=0\\i\neq j}}^n
  \frac{(x - x_i)}{(x_j - x_i)}
\end{equation}
Zachodzi następująca tożsamość
\begin{equation}
  \label{eq:Abel}
  \begin{split}
    W_n(x) &= W_0(x) 
    + \bigl[W_1(x) - W_0(x)\bigr]
    + \bigl[W_2(x) - W_1(x)\bigr]
    + \dots +\\
    &+ \bigl[W_n(x) - W_{n-1}(x)\bigr]\,.
  \end{split}
\end{equation}
Różnica $W_k(x) - W_{k-1}(x)$ jest wielomianem, który można
przedstawić w postaci
\begin{equation}
  \label{eq:diff}
  W_k(x) - W_{k-1}(x) = A_k\,(x-x_0)(x-x_1)\dots(x-x_{k-1})\,,
\end{equation}
gdzie $A$ jest pewną stałą.  Aby ją obliczyć, podstawiamy $x = x_k$
i~otrzymujemy
\begin{equation}
  \label{eq:A}
  A_k = \sum_{i=0}^{k} \frac{y_i}{\omega_k(x_i)}\,,
\end{equation}
gdzie $\omega_k(x)$ oznacza wielomian $k$--tego stopnia w~bazie
Newtona,
\begin{equation}
  \label{eq:Newton-base}
  \omega_k(x) = (x-x_0)(x-x_1)\dots(x-x_k)\,.
\end{equation}
Można udowodnić (metodą indukcji), że
\begin{equation}
  \begin{split}
  f(x_i;x_{i+1};\dots;x_{i+n}) &=
  \frac{y_i}{(x_i - x_{i+1})(x_i - x_{i+2})\dots(x_i-x_{i+n})} +\\
  &+ \frac{y_{i+1}}{(x_{i+1} - x_i)(x_{i+1} - x_{i+2})\dots(x_{i+1}-x_{i+n})} +
  \dots + \\
  &+ \frac{y_{i+n}}{(x_{i+n} - x_i)(x_{i+n} - x_{i+1})\dots(x_{i+n}-x_{i+n-1})} \\
  &= \sum_{k=i}^{i+n} 
  \frac{f(x_k)}{\displaystyle\prod^{i+n}_{\substack{l=i\\l\neq k}} (x_k - x_l)}
  \end{split}
\end{equation}
czego szczególnym przypadkiem jest
\begin{equation}
  f(x_0;x_1;\dots;x_n) 
  = \sum_{j=0}^{n} 
    \frac{f(x_j)}{\displaystyle\prod^{n}_{\substack{k=0\\l\neq k}} (x_k - x_l)}
  = \sum_{j=0}^{n} \frac{y_j}{\omega_n(x_j)}\,,
\end{equation}
gdzie $f(x_0;x_1;\dots;x_n)$ oznacza \emph{iloraz różnicowy $n$--tego
  rzędu} zdefiniowany wzorem rekurencyjnym
\begin{subequations}
  \begin{align}
    f(x_i) &= f(x_i) = y_i \\
    f(x_i;x_{i+1}) &= 
    \frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i} =
    \frac{y_{i+1} - y_i}{x_{i+1} - x_i}\,, \\
    f(x_i;x_{i+1};\dots;x_{i+n}) &=
    \frac{f(x_{i+1};x_{i+2};\dots;x_{i+n}) -
      f(x_i;x_{i+1};\dots;x_{i+n-1})}
    {x_{i+n} - x_i}\,.
  \end{align}
\end{subequations}
Kolejne ilorazy różnicowe możemy wyliczać w~analogiczny sposób jak
wyniki cząstkowe metody Neville'a.  Kolejne współczynniki w~bazie
Newtona będą zatem równe kolejnym ilorazom różnicowym:
\begin{equation}
  \label{eq:Newton-intwiel}
  A_k = f(x_0;x_1;\dots;x_k)\quad\text{dla $k=0,1,\dots,n$}\;.
\end{equation}


\subsubsection{Algorytm Hornera dla wielomianu w bazie Newtona}
\label{sec:Horner}


Niech dany będzie wielomian w~postaci Newtona
\begin{equation}
  \label{eq:Newtona}
  W(x) 
  = b_0 + \sum_{i=1}^n b_i (x-x_0) \dots (x-x_{i-1})\,.
\end{equation}
Do obliczania wartości tego wielomianu możemy zastosować uogólniony
algorytm Hornera:
\begin{equation}
  \label{eq:uog-Horner}
  W(x) = (\ldots(b_n (x-x_{n-1}) + b_{n-1}) (x-x_{n-2} + 
         \dots + a_1) (x-x_0) + a_0.
\end{equation}
Algorytm dany jest wzorem
\begin{subequations}
\label{eq:uog-Horner:alg}
  \begin{align}
    w_n(x) &= b_n,\\
    w_i(x) &= w_{n+1}(x) (x-x_i) + b_i,\quad\text{dla\ }i = n-1,n-2,\dots,0\,\\[0.5em]
    \dots\notag
    \intertext{Wynikiem jest}
    W(x) &= w_0(x)
  \end{align}
\end{subequations}


\subsubsection{Zamiana współczynników w bazie Newtona na naturalne (w bazie jednomianów)}
\label{sec:zamiana}


Każdy z~wielomianów pomocniczych $w_i(x)$ w~metodzie Hornera zapiszmy
w~postaci naturalnej:
\begin{equation}
  \label{eq:Newton-Horner-do-potegowej}
  w_i(x) = \sum_{k=i}^{n} a_k^{(i)} x^{k-i},\quad i=n,n-1,\dots,0
\end{equation}
gdzie dostajemy na współczynniki $a_k^{(i)}$ następujący układ równań
rekurencyjnych
\begin{subequations}
  \label{eq:Newton-do-potegowej}
  \begin{align}
    \label{eq:N2p:ann}
    a_n^{(i)} &= a_n^{(i+1)} = \dots = a_n^{(n)} = b_n,\\
    \label{eq:N2p:aii}
    a_i^{(i)} &= b_i - a_{i+1}^{(i+1)} x_i,\\
    \label{eq:N2p:next}
    a_k^{(i)} &= a_k^{(i+1)} - a_{k+1}^{(i+1)} x_i, 
    \quad \text{dla\ }k = i+1, i+2, \ldots, n-1.
  \end{align}
\end{subequations}
\begin{table}[htbp!]
  \centering
  \begin{tabular}{cccc}
    $i=n$ & $i=n-1$ & \dots & $i=0$ \\
    \hline
    $a_n^{(n)}$ & $a_{n-1}^{(n-1)}$ & & $a_0^{(0)}$ \\
                & $a_{n}^{(n-1)}$   & & $a_1^{(0)}$ \\
    \multicolumn{4}{c}{$\cdots$} \\
                &                   & & $a_n^{0}$ 
  \end{tabular}
  \caption{Przejście z bazy Newtona do bazy jednomianów}
  \label{tab:Newton-to-std}
\end{table}
Szukane współczynniki postaci naturalnej wielomianu $W(x)$ to
oczywiście $a_k^{(0)}$ ($k=0,1,\dots,n$).  

Zapis tego algorytmu ma następującą postać
\begin{lstlisting}[language=Algol,morekeywords={step,until}]{newton2natural}
  a[n] := b[n];
  for i := n-1 step -1 until 0 do
  begin
    xi := x[i];
    a[i] := b[i];
    for k := i step 1 until n-1 do
      a[k] := a[k] - xi $\times$ a[k+1];
  end;
\end{lstlisting}
Koszt tego algorytmu to $n (n+1)/2$ mnożeń i~odejmowań.  Jest on
numerycznie stabilny.


\subsection{Zastosowanie algorytmu Neville'a}
\label{sec:neville-wspolcz}

Inną metodą jest użycie algorytmu Neville'a wyznaczania wartości
wielomianu intepolującego w~punkcie.  Niech wielomian interpolujący
będzie dany wzorem
\begin{equation}
  \label{eq:wielom-interp}
  P(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_m x^m.
\end{equation}
Jeśli dokonamy interpolacji (ekstrapolacji) wielomianowej w~punkcie 
$x = 0$, to otrzymana wartość będzie oczywiście równa zerowemu
współczynnikowi wielomianu, czyli $a_0$.  Jeśli teraz odejmiemy od
wielomianu otrzymane $a_0$ i~podzielimy go przez $x$ zmniejszając jego
stopień, to w~ten sam sposób możemy obliczyć następny współczynnik
itd.  

Na język algorytmu Neville'a tłumaczy się to w~następujący sposób.
Najpierw znajdujemy za pomocą algorytmu Neville'a wartość w~$x = 0$
(oznaczmy ją przez $a_k$).  Od każdej stabelaryzowanej wartości $y_i$
odejmujemy otrzymana wartość w~zerze [przesunięcie].  Następnie każdą
stabelaryzowaną wartość $y_i$ dzielimy przez odpowiadające jej $x_i$,
i~odrzucamy jeden z~punktów (np. o~najmniejszej wartości bezwzględnej
$x_i$) [zmniejszenie stopnia wielomianu].
\begin{equation}
  \label{eq:krok}
  y_i := \frac{y_i - c_k}{x_i}.
\end{equation}
Czynność powtarzamy dopóki nie uzyskamy pełnej postaci wielomianu.

Jeżeli $x = 0$ leży daleko poza zakresem stabelaryzowanych wartości,
$x_i$, bądź jedno z~$x_i$ jest równe lub bardzo bliskie 0, należy
rozważyć przesunięcie zagadnienia (obliczanie $P(x - a)$).
Współczynniki wielomianu po przesunięciu możemy wyliczyć np. używając
algorytmu Hornera do wyliczania znormalizowanych pochodnych
$w^{(j)}(x)/j!$ dzieląc wielomian przez $(x-a)$.


\section{Zadanie domowe}
\label{sec:zadanie}

Napisać program znajdujący współczynniki wielomianu interpolacyjnego,
albo w~bazie jednomianów przy zastosowaniu (napisanego wcześniej)
algorytmu Neville'a, albo w~bazie Newtona obliczając rekurencyjnie
odpowiednie ilorazy różnicowe.  Można założyć, że nie ma potrzeby
przesuwania wielomianu, tzn. 0 leży niezbyt daleko od wartości $x_i$,
a~jednocześnie żadne z~$x_i$ nie jest niebezpiecznie bliskie 0.

Sprawdzić, czy algorytm ten zastosowany do wielomianu (odpowiednio
zadanego w~postaci naturalnej lub w~bazie Newtona) odtwarza jego
współczynniki (także dla ilości punktów większej od ilości
współczynników wielomianu); warto sprawdzić, czy algorytm Neville'a
odtwarza wartości wielomianu (czy pierwszy krok daje poprawne wyniki).

Zastosować algorytm do wybranych funkcji.  (W~przeciwieństwie do
zadania z~algorytmem Neville'a nie potrzeba dokonywać interpolacji na
większej ilości punktów niż stopień wielomianu plus jeden).  Wartość
interpolacji w~punkcie należy znajdować obliczając znaleziony
wielomian interpolacyjny w~zadanym punkcie np. za pomocą odpowiedniego
do postaci wielomianu algorytmu Hornera.  Narysować (np. zapisując
wyniki do plików tekstowych i~używając programu \texttt{gnuplot}:
polecenia dla niego też można zapisywać do pliku w~naszym programie)
wykres funkcji, wykres wielomianu interpolacyjnego oraz zaznaczyć
(narysować za pomocą punktów) punkty na których oparta jest
interpolacja.

\textbf{Ogólnie:} należy przetestować algorytm znajdowania
współczynników wielomianu interpolacyjnego.  Powyższe zagadnienia to
plan optymalny; można zrobić mniej, można zrobić więcej.
%%
\\
\bigskip
\\
\input{podpis.tex}

\end{document}
