\documentclass[10pt,a4paper]{article}

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

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts}

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

\begin{document}

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

\section{Interpolacja funkcjami sklejanymi}
\label{sec:splines}


\subsection{Sformułowanie zadania}
\label{sec:zadanie-interpolacji}

Przy interpolacji Lagrange'a lub wymiernej zakłada się, że dana
funkcja $f$ przybliżana jest \textbf{jednym} wielomianem (funkcją
wymierną) na całym odcinku $[a,b]$.  Przy tym założeniu jedyną
możliwością uzyskania lepszego przybliżenia jest zwiększenie stopnia
wielomianu.  Nie zawsze musi to jednak działać (przykładem może być
funkcja $1/(1+x^2)$ interpolowana za pomocą wielomianu).  W~tej
sytuacji naturalne jest inne sformułowanie zadania.

Niech ($\Delta$) będzie układem $N+1$ punktów $x_0, x_1, \ldots, x_N$
dzielących przedział $[a,b]$ na $N$ części, tzn.
\begin{equation}
(\Delta)\colon \quad a = x_0 < x_1 < \ldots < x_N = b
\end{equation}

W~każdym z~podprzedziałów $[x_i,x_{i+1}]$ przybliżamy daną funkcję $f$
wielomianem ustalonego (raczej niskiego) stopnia.  Zakładamy ponadto,
że funkcja interpolująca była ciągła wraz z~odpowiednimi pochodnymi na
\textbf{całym odcinku} $\mathbf{[a,b]}$.  Funkcjami o~tych
własnościach są \textbf{funkcje sklejane} (ang. spline).  Omówimy tu
własności interpolacyjnych funkcji sklejanych i~przedstawimy algorytm
ich wyznaczania dla przypadku funkcji stopnia 3.


\subsection{Ogólne fakty dotyczące funkcji sklejane} 
\label{sec:funkcje-sklejane}

\begin{df}
  \label{def:spline}
  Funkcję rzeczywistą $S$ nazywamy \textbf{funkcją sklejaną} stopnia
  $m$ z~węzłami $x_0, x_1,\ldots, x_N$, gdzie $x_0 < x_1 < \ldots <
  x_N$, jeżeli spełnia poniższe własności:
  \begin{itemize}
  \item w~każdym przedziale $(x_{i-1},x_i)$ dla $i=-1,0,1,\ldots ,N+1$\\
    $(x_{-1}=-\inf, x_{N+1}=+\inf)$ $S$ jest wielomianem stopnia nie
    wyższego niż $m$,
    
  \item funkcja $S$ i~jej pochodne rzędu $1,2,\ldots,m-1$ są ciągłe na
    całej osi rzeczywistej, tzn. $S\in \mathcal{C}^{m-1}$.
  \end{itemize}  
\end{df}

W~najprostszym przypadku, dla $m=1$, funkcja sklejana jest po prostu
łamaną przechodzącą przez punkty $x_i$.  Z~definicji tej wynika
również, że wielomiany są szczególnym przypadkiem funkcji sklejanych.

\begin{df}
  \label{def:natural-spline}
  Funkcję sklejana stopnia $2 m-1$ z~węzłami $x_0,x_1,\dots,x_N$
  nazywamy \textbf{naturalną} funkcją sklejaną, jeśli w~przedziałach
  $(-\inf,x_0)$ i~$(x_N,+\inf)$ dana jest wielomianami $m-1$ stopnia
  (a~nie $2m-1$ stopnia).
\end{df}

Dalej będziemy zajmowali się funkcjami sklejanymi interpolującymi dane
punkty $(x_i,y_i)$, tzn. takimi, że
\begin{equation}
  \label{eq:interp-cond}
  S(x_i) = y_i \quad \text{dla $i = 0, 1, \ldots, N$}.
\end{equation}
W~całej klasie funkcji sklejanych zadanie to nie jest wyznaczone
jednoznacznie (co zobaczymy także w~rozdziale~\ref{sec:algorytm}).
Wystarczy jednak ograniczyć się do naturalnych funkcji sklejanych, dla
których jest prawdziwe poniższe twierdzenie.

\begin{thm}
  \label{thm:jednozn-interp-spline}
  Jeżeli węzły $x_i$ są różne dla $i=0,1,\ldots,N$, oraz $1 \leq m
  \leq N+1$, to dla dowolnych wartości $y_i$ istnieje dokładnie jedna
  naturalna funkcja sklejana $S$ interpolująca punkty $(x_i,y_i)$,
  tzn. taka, że $S(x_i) = y_i$ dla $i=0,1,\ldots,N$.
\end{thm}


\subsection{Algorytm wyznaczania interpolującej funkcji sklejanej}
\label{sec:algorytm}


\subsubsection{Algorytm ,,wprost''}
\label{sec:algorytm-wprost}

Zajmijmy się algorytmem wyznaczania naturalnej interpolacyjnej funkcji
sklejanej $S$ stopnia 3.  Najprostszą metodą jest przedstawienie
szukanej funkcji w~każdym z~podprzedziałów jako wielomian stopnia 3.
w~postaci naturalnej:
\begin{equation}
  \label{eq:spline-xn}
  S(x) = a_i + b_i t + c_i t^2 + d_i t^3
\end{equation}
gdzie
\begin{math}
  t = x-x_i \quad \text{dla\  } x \in [x_i,x_{i+1}), \quad i=0,1,\ldots,N-1.  
\end{math}
Oczywiście $a_i = S(x_i) = y_i$. Do wyznaczenia pozostaje zatem
jeszcze $3N$ współczynników. Ponieważ $S, S', S''$ mają być ciągłe
w~węzłach $x_i(i=1,2,\ldots,N-1)$, otrzymamy stąd $3N-3$ równań
(granica z~lewej równa jest granicy prawej).  Brakujące warunki
uzyskamy z~założenia, że $S$ jest funkcja sklejana naturalną, czyli że
drugie pochodne na zewnątrz przedziału $[a, b]$ się zerują.

Nie będziemy tutaj przedstawiać postaci tych równań.  Pozostawiamy do
sprawdzenia, że układ równań na współczynniki wielomianów składowych
funkcji sklejanej jest układem trójdiagonalnym.  Widać, że ta postać
nie jest zbyt wygodna do wyznaczania wartości funkcji sklejanej
w~punkcie (musimy pamiętać dodatkowo trzy współczynniki).  Przydaje
się jednak ona do całkowania/różniczkowania numerycznego.


\subsubsection{Wyznaczanie funkcji sklejanej przy znanych drugich pochodnych}
\label{sec:drugie-pochodne}

Mamy zadaną funkcję w~postaci stabularyzowanej $y_i = f(x_i)$, $i = 0,
1, \dots, N$.  Przypomnijmy wzór na interpolację liniową na przedziale
$(x_j, x_{j+1})$ w~punkcie:
\begin{equation}
  \label{eq:interp-liniowa}
  y = A y_{j} + B y_{j+1}, 
\end{equation}
gdzie
\begin{align}
  \label{eq:interp-liniowa-ozn}
  A &\equiv \frac{x_{j+1} - x}{x_{j+1} - x_{j}},
  &
  B &\equiv 1 - A = \frac{x - x_{j}}{x_{j+1} - x_{j}}.
\end{align}
Jest to szczególna postać wzoru interpolacyjnego Lagrange'a.

Załóżmy, że oprócz stabularyzowanych wartości funkcji $y_i = f(x_i)$
mamy także stabularyzowane wartości drugich pochodnych $y''_i =
f''(x_i)$.  (Później podamy sposób wyznaczania tych wartości).  Do
równania~\eqref{eq:interp-liniowa} możemy dodać wielomian trzeciego
stopnia, którego drugie pochodne w~punktach $x_{j}$ i~$x_{j+1}$ są
równe odpowiednio $y''_{j}$ i~$y''_{j+1}$, zaś wartości w~tych
punktach są równe zeru.  Odpowiednikiem
wzoru~\eqref{eq:interp-liniowa} jest
\begin{equation}
  \label{eq:interp-spline}
  y = A y_{j} + B y_{j+1} + C y''_{j} + D y''_{j+1}, 
\end{equation}
gdzie
\begin{equation}
  \label{eq:interp-spline-ozn}
  \begin{aligned}
    A &\equiv \frac{x_{j+1} - x}{x_{j+1} - x_{j}},         &
    B &\equiv 1 - A = \frac{x - x_{j}}{x_{j+1} - x_{j}},   \\
    C &\equiv \frac{1}{6}(A^3 - A)(x_{j+1} - x_{j})^2,    &
    D &\equiv \frac{1}{6}(B^3 - B)(x_{j+1} - x_{j})^2.   
  \end{aligned}
\end{equation}
Łatwo zauważyć, że powyższy wzór daje dobre wartości w~punkcie.
Pozostawiamy do sprawdzenia, że daje dobre wartości drugich
pochodnych.



\subsubsection{Obliczanie drugich pochodnych z warunku ciągłości}
\label{sec:oblicz-drug-poch}

Wyznaczona przez nas funkcja interpolująca jest ciągła i~ma
(z~definicji) ciągłe drugie pochodne.  Nie znamy jednak zazwyczaj
drugich pochodnych funkcji.  Możemy jednak wyznaczyć drugie pochodne
interpolującej funkcji sklejanej przez narzucenie warunku ciągłości
\emph{pierwszych} pochodnych funkcji interpolującej.

Żądane równania uzyskujemy żądając by pochodna
równania~\eqref{eq:interp-spline} obliczona w~punkcie $x_{j}$ na
przedziale $(x_{j}, x_{j+1})$ była równa pochodnej w~tym samym punkcie
ale obliczona na przedziale $(x_{j-1}, x_{j})$.  Po odpowiednich
przekształceniach uzyskujemy równanie
\begin{equation}
  \label{eq:spline-tridiag}
  \frac{x_{j}   - x_{j-1}}{6} y''_{j-1} +
  \frac{x_{j+1} - x_{j-1}}{3} y''_{j}   +
  \frac{x_{j+1} - x_{j}  }{6} y''_{j+1} =
  \frac{y_{j+1} - y_{j}}{x_{j+1} - x_{j}} -
  \frac{y_{j} - y_{j-1}}{x_{j} - x_{j-1}}.
\end{equation}
Jest to $(N+1) - 2$ równań liniowych na $N+1$ zmiennych \(y''_i, i~= 0,
1, \dots, N\).  Brakujące równania (na $y''_0$ i~$y''_N$) możemy
uzyskać na jeden z~poniższych sposobów:
\begin{itemize}
\item ustalić jedną lub obie $y''_0$ i~$y''_N$ na zero, otrzymując
  tzw. \emph{naturalne funkcje sklejane} (patrz
  definicja~\ref{def:natural-spline}),

\item ustalić jedną lub obie $y''_0$ i~$y''_N$ na zadaną wartość,

\item obliczyć jedną lub obie $y''_0$ i~$y''_N$ na podstawie zadanej
  pierwszej pochodnej w~odpowiednim punkcie

\item  obliczyć $y''_0$ i~$y''_N$ narzucając cykliczne warunki
  brzegowe (czyli m. in. pierwsze pochodne w~$x_0$ i~$x_N$ są sobie
  równe).
\end{itemize}
Wszystkie z~powyższych warunków brzegowych z~wyjątkiem ostatniego dają
układ równań w~postaci \emph{trójdiagonalnej}, która jest szczególnie
łatwa do rozwiązania.  Ostatni warunek brzegowy daje układ równań
w~postaci prawie trójdiagonalnej.



\subsubsection{Algorytm ,,dwustopniowy''}
\label{sec:algorytm-dwustopniowy}

Wygodnie zatem zastosować następujący algorytm.  Najpierw wyliczamy
drugie pochodne funkcji sklejanych dla zadanych punktów, rozwiązując
w~tym celu układ równań~\eqref{eq:spline-tridiag}.  Następnie
stosujemy wzory~\eqref{eq:interp-spline}, \eqref{eq:interp-spline-ozn}
do obliczania wartości interpolującej funkcji sklejanej w~punkcie.

W~porównaniu z~algorytmem ,,wprost'' musimy zapamiętać mniej wartości
(jedną na punkt a~nie trzy na punkt danych).  Algorytm ma też lepsze
własności numeryczne (wyznaczanie współczynników wielomianu w~postaci
naturalnej jest często źle uwarunkowane).  Nie daje on nam jednak
równania (postaci analitycznej) funkcji sklejanej.



\subsection{Zadanie}
\label{sec:zadanie}

Napisać program dokonujący interpolację funkcjami sklejanymi trzeciego
stopnia.  Wymagane jest zastosowanie tylko \emph{jednego} z~opisanych
algorytmów.  Należy przetestować interpolację funkcjami sklejanymi
wielomianów, funkcji będących kawałkami wielomianami (np. $\lvert x
\rvert$), oraz jakiejś ,,trudnej'' funkcji, np. $1/(1+x^2)$.  Grafika
wskazana.

\textbf{Uwaga:} funkcje sklejane \emph{naturalne} nie nadają się do
\textbf{ekstrapolacji}!


\appendix
\section{Rozwiązywanie układu równań trójliniowych}
\label{sec:trojliniowe}

Układ równań liniowych typu:
\begin{equation}
  \label{eq:trojliniowe}
  \begin{bmatrix}
    b_1 & c_1 & 0   & \dots & \dots & 0   \\
    a_2 & b_2 & c_2 & 0   & \dots & \dots \\
    0   & a_3 & b_3 & c_3 & 0 & \ldots   \\
    \hdotsfor[1.5]{6} \\
    \dots & \dots  & 0     & a_{N-1} & b_{N-1} & c_{N-1} \\
    0     & \dots  & \dots & 0       & a_N     & c_N     
  \end{bmatrix}
  \begin{bmatrix}
    x_1 \\
    x_2 \\
    x_3 \\
    \hfill\vdots\hfill \\
    x_{N-1} \\
    x_N
  \end{bmatrix}
  =
  \begin{bmatrix}
    v_1 \\
    v_2 \\
    v_3 \\
    \hfill\vdots\hfill \\
    v_{N-1} \\
    v_N
  \end{bmatrix}  
\end{equation}
nazywamy układem w~postaci \textbf{trójliniowej}.

Taki układ można rozwiązać w~bardzo prosty sposób poprzez kolejne
odejmowanie wierszy w~celu wyeliminowania kolejnych zmiennych.  To
znaczy: od drugiego wiersza należy odjąć pierwszy przemnożony przez
$\frac{b_2}{a_1}$.  Pozbywamy się w~ten sposób zmiennej $x_1$
z~drugiego równania.  Powtarzając te czynność w~ostatnim równaniu
otrzymamy równanie na $x_N$.  Następnie idąc do góry będziemy mogli
obliczyć wszystkie $x_i$ na podstawie znanej wartości $x_{i+1}$.

Wybór elementu głównego nie jest zazwyczaj potrzebny przy
rozwiązywaniu układu równań trojliniowych.  Zagadnienia które prowadzą
do równań trójliniowych dają zazwyczaj układy \emph{diagonalnie
  dominujące}, tzn. takie że $|b_j| > |a_j| + |c_j|$.

Poniżej podane są algorytmy rozwiązywania układu równań trójliniowych
(być może niektóre z nich są równoważne):

\textbf{Zmodyfikowany algorytm Gaussa}
\nopagebreak
\begin{subequations}
  \label{eq:jankowscy}
  \begin{align}
    y_N &= \beta_{N+1}, &
    y_i &= \alpha_{i+1} y_{i+1} + \beta_{i+1}                    \quad i = N-1,\dots,1 \\*
    \intertext{gdzie}
    \alpha_2     &= \frac{b_1}{c_1}, &
    \alpha_{i+1} &= \frac{b_i}{c_i - a_i \alpha_i}               \quad i = 2,\dots,N-1 \\*
    \beta_2      &= \frac{v_1}{c_1}, &
    \beta_{i+1}  &= \frac{v_i - a_i \beta_i}{c_i - a_i \alpha_i} \quad i = 2,\dots,N 
  \end{align}
\end{subequations}

\bigskip

\textbf{Algorytm Gaussa}
\nopagebreak
\begin{subequations}
  \label{eq:fortuna:gauss}
  \begin{align}
    && u_1 &= b_1, & y_1 &= v_1,  \notag \\
    l &= \frac{a_i}{u_{i-1}}, & u_i &= b_i - l c_i, & y_i &= v_i - l y_{i-1}, 
  \end{align}
  dla $i = 2, 3, \dots, N$, oraz
  \begin{align}
    x_N &= \frac{y_N}{u_N}, && \notag \\
    x_i &= \frac{y_i - c_i x_{i+1}}{u_i}, &&
  \end{align}
  dla $i = N-1, N-2, \dots, 1$.
\end{subequations}

\bigskip

\textbf{Algorytm Crouta}
\nopagebreak
\begin{subequations}
  \label{eq:fortuna:gauss}
  \begin{align}
    && l &= b_1, & y_1 &= v_1/l,  \notag \\
    u_i &= \frac{c_i}{l} & l &= b_{i+1} - u_i a_{i+1}, & y_i &= \frac{v_{i+1} -  a_{i+1} y_i}{l}, 
  \end{align}
  dla $i = 1, 2, \dots, N-1$, oraz
  \begin{align}
    x_N &= y_N, && \notag \\
    x_i &= y_i - u_i x_{i+1}, &&
  \end{align}
  dla $i = N-1, N-2, \dots, 1$.
\end{subequations}

\bigskip

Należy wybrać jeden z algorytmów, zaimplementować go (przekazując do
procedury \textbf{wyłącznie} niezerowe elementy tablicy, ewentualnie
wektory pomocnicze) i~go przetestować (sprawdzić czy znalezione
rozwiązanie spełnia podane równanie).
%%
\input{podpis.tex}

\end{document}

