MN11: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przykry (dyskusja | edycje)
mNie podano opisu zmian
 
 
(Nie pokazano 16 wersji utworzonych przez 3 użytkowników)
Linia 1: Linia 1:


==Nadokreślone układy równań liniowych==
<!--
Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki
wprowadzi� przykry@mimuw.edu.pl
-->
=Funkcje sklejane (splajny)=
 
{{powrot |Metody numeryczne | do strony głównej
przedmiotu <strong>Metody numeryczne</strong>}}
 
Interpolacja  wielomianami interpolacyjnymi, chociaż korzysta z funkcji gładkich
i łatwo reprezentowalnych w komputerze, ma jednak również pewne wady.
Zauważmy, że błąd interpolacji może być bardzo duży ([[MN09#Zjawisko Rungego i dobór węzłów interpolacji|zjawisko Rungego]]), a poza tym interpolacja jest nielokalna: nawet mała zmiana warości funkcji w pojedynczym węźle może powodować dużą zmianę zachowania całego wielomianu interpolacyjnego. Czasem więc lepiej jest zastosować innego rodzaju
interpolację, np. posługując się funkcjami sklejanymi, które tylko ''lokalnie są wielomianami'', sklejonymi w taki sposób, by ''globalnie'' zachować pewien stopień gładkości, tzn. ''różniczkowalność zadaną liczbę razy''. 
 
Tego typu podejście okazało się bardzo owocne m.in. w grafice komputerowej (np. dla wizualizacji scenerii w grach komputerowych), a także np. posłużyło jako narzędzie konstrukcji skalowalnych czcionek komputerowych w Postscripcie (precyzyjniej, korzysta się tam z tzw. <strong>krzywych Beziera</strong> --- pewnych krzywych sklejanych zadanych na płaszczyźnie). Z krzywych Beziera powszechnie korzysta się również w systemach CAD (''Computer Aided Design'').


Zajmiemy się zadaniem wygładzania liniowego,  
Zamiast terminu ''funkcje sklejane'' używa się też często terminów <strong>splajny</strong> (''spline''), albo <strong>funkcje gięte</strong>. Nazwy te biorą się stąd, że zadanie interpolacji naturalnym splajnem kubicznym można interpretować jako model matematyczny aparatu służącego do wytwarzania mebli giętych.
nazywanym też liniowym zadaniem najmniejszych kwadratów.
Jest ono uogólnieniem zadania rozwiązywania kwadratowych układów
równań liniowych do przy\-pa\-dku, gdy układ jest nadokreślony.  


Jest to praktycznie bardzo często pojawiające się zadanie (pewien jego wariant
==Funkcje sklejane==
rozwiązują np. nasze przenośne odbiorniki GPS), a autorem pierwszego rozwiązania
był nie kto inny jak sam wielki Gauss.


===Układ normalny===
W ogólności przez funkcję sklejaną rozumie się każdą
funkcję przedziałami wielomianową. Nas będą jednak
interesować szczególne funkcje tego typu i dlatego termin
<strong>funkcje sklejane</strong> zarezerwujemy dla funkcji
przedziałami wielomianowych i posiadających dodatkowe
własności, które teraz określimy.


Niech <math>\displaystyle A</math> będzie daną macierzą o <math>\displaystyle m</math> wierszach i <math>\displaystyle n</math> kolumnach,
Niech dany będzie przedział skończony <math>[a,b]</math> i węzły
<math>\displaystyle A\inR^{m\times n}</math>, taką, że


<center><math>\displaystyle m\,\ge\,n\,=\, \mbox{rank} (A),  
<center><math>a \,=\, x_0 \,<\, x_1 \,<\, \cdots \,<\, x_n \,=\, b,  
</math></center>
</math></center>


albo równoważnie, taką że jej wektory kolumny są liniowo
przy czym <math>n\ge 1</math>.  
niezależne. Niech także dany będzie wektor <math>\displaystyle  b\inR^m</math>.
Jasne jest, że wtedy układ równań <math>\displaystyle A x= b</math> nie zawsze
ma rozwiązanie - mówimy, że układ jest ''nadokreślony''.  


''Zadanie wygładzania liniowego'' polega na znalezieniu wektora
{{definicja|||
<math>\displaystyle  x^*\inR^n</math>, który minimalizuje ''wektor residualny''
<math>\displaystyle  r= b-A x</math> w normie drugiej, tzn.


<center><math>\displaystyle \| b\,-\,A x^*\|_2\,=\,\min_{ x\inR^n}
Funkcję <math>s:R\to R</math> nazywamy <strong>funkcją sklejaną rzędu <math>r</math></strong> (<math>r\ge 1</math>) odpowiadającą węzłom <math>x_j</math>, <math>0\le j\le n</math>, jeśli spełnione są następujące
    \| b\,-\,A x\|_2.
dwa warunki:
</math></center>
; (i)
:  <math>s</math> jest wielomianem stopnia co najwyżej <math>2r-1</math> na każdym
z przedziałów <math>[x_{j-1},x_j]</math>, tzn.
<math>s_{|[x_{j-1},x_j]}\in\Pi_{2r-1}</math>, <math>1\le j\le n</math>,  
 
; (ii)
:  <math>s</math> jest <math>(2r-2)</math>-krotnie różniczkowalna w sposób
ciągły na całej prostej, tzn. <math>s\in C^{(2r-2)}(R)</math>.
Jeśli ponadto
; (iii)
:  <math>s</math> jest wielomianem stopnia co najwyżej <math>r-1</math> poza
<math>(a,b)</math>, tzn. <math>s_{|(-\infty,a]},s_{|[b,+\infty)}\in\Pi_{r-1}</math>,
to <math>s</math> jest naturalną funkcją sklejaną.
}}
 
Klasę naturalnych funkcji sklejanych rzędu <math>r</math> opartych na
węzłach <math>x_j</math> będziemy oznaczać przez
<math>{\cal S}_r(x_0,\ldots,x_n)</math>, albo po prostu <math>{\cal S}_r</math>,
jeśli węzły są ustalone.
 
Na przykład funkcją sklejaną rzędu pierwszego (<math>r=1</math>)
jest funkcja ciągła i liniowa na poszczególnych
przedziałach <math>[x_{j-1},x_j]</math>. Jest ona naturalna, gdy poza
<math>(a,b)</math> jest funkcją stała. Tego typu funkcje nazywamy
<strong>liniowymi funkcjami sklejanymi</strong>.


<div style="margin-top:1em; padding-top,padding-bottom:1em;">
Najważniejszymi z praktycznego punktu widzenia są jednak
<span  style="font-variant:small-caps;">Przykład</span>  
funkcje sklejane rzędu drugiego odpowiadające <math>r=2</math>. Są
<div class="solution">
to funkcje, które są na <math>R</math> dwa razy różniczkowalne w
sposób ciągły, a na każdym z podprzedziałów są
wielomianami stopnia co najwyżej trzeciego. W tym przypadku
mówimy o <strong>kubicznych funkcjach sklejanych</strong>. Funkcja
sklejana kubiczna <math>s</math> jest naturalna, gdy poza <math>(a,b)</math> jest
wielomianem liniowym, a więc <math>s''(a) = s''(b) = 0</math>.


Przypuśćmy, że dla pewnej funkcji 
==Interpolacja i gładkość==
<math>\displaystyle f:[a,b]\toR</math> obserwujemy jej wartości <math>\displaystyle f_i</math> (dokładne lub
zaburzone) w punktach <math>\displaystyle t_i</math>, <math>\displaystyle 1\le i\le m</math>. Funkcję tą
chcielibyśmy przybliżyć inną funkcją <math>\displaystyle w</math> należącą do
pewnej <math>\displaystyle n</math> wymiarowej przestrzeni liniowej <math>\displaystyle W</math>, np. przestrzeni
wielomianów stopnia mniejszego niż <math>\displaystyle n</math>. Jakość przybliżenia
mierzymy wielkością


<center><math>\displaystyle
Pokażemy najpierw ważny lemat, który okaże się kluczem
  \sum_{i=1}^m (f_i-w(t_i))^2.
do dowodu dalszych twierdzeń.  
</math></center>


Wybierając pewną bazę <math>\displaystyle (w_j)_{j=1}^n</math> w <math>\displaystyle W</math> i rozwijając <math>\displaystyle w</math>  
Niech <math>W^r(a,b)</math> będzie klasą funkcji <math>f:[a,b]\to R</math> takich,
w tej bazie, <math>\displaystyle w(t)=\sum_{j=1}^n c_jw_j(t)</math>, sprowadzamy problem
że <math>f</math> jest <math>(r-1)</math> razy różniczkowalna na <math>[a,b]</math> w sposób
do minimalizacji
ciągły oraz <math>f^{(r)}(x)</math> istnieje prawie wszędzie na <math>[a,b]</math>
i jest całkowalna z kwadratem, tzn.


<center><math>\displaystyle \sum_{i=1}^m\left(f_i-\sum_{j=1}^n c_jw_j(t_i)\right)^2
<center><math>\begin{align} W^r(a,b) &= \{\,f\in C^{(r-1)}([a,b]):\, f^{(r)}(x)
</math></center>
    \mbox{  istnieje p.w. na  }  [a,b] \\
    && \qquad\qquad\qquad\qquad  \mbox{ oraz  }  
      f^{(r)}\in{\cal L}_2(a,b)\,\}.
\end{align}</math></center>


względem <math>\displaystyle c_j</math>, a więc do zadania wygładzania liniowego.
Oczywiście każda funkcja sklejana rzędu <math>r</math> (niekoniecznie
Rzeczywiście, kładąc
naturalna) należy do <math>W^r(a,b)</math>.  
<math>\displaystyle A=(a_{i,j})\inR^{m\times n}</math> z <math>\displaystyle a_{i,j}=w_j(t_i)</math>,  
<math>\displaystyle  b=(f_i)_{i=1}^m</math> i <math>\displaystyle  x=(c_j)_{j=1}^n</math>, wielkość
([[##unorm|Uzupelnic: unorm ]]) jest równa <math>\displaystyle \| b-A x\|_2^2</math>.  
</div></div>


{{lemat|||
{{lemat|||
Zadanie wygładzania liniowego ma jednoznaczne
rozwiązanie <math>\displaystyle  x^*</math>, które spełnia układ równań
Niech <math>f\in W^r(a,b)</math> będzie funkcją
zerującą się w węzłach, tzn.


<center><math>\displaystyle
<center><math>f(x_j)\,=\,0,\qquad 0\le j\le n</math></center>
  A^TA x\,=\,A^T\, b.
 
Wtedy dla dowolnej naturalnej funkcji sklejanej <math>s\in {\cal S}_r</math>
mamy
 
<center><math>\int_a^b f^{(r)}(x)s^{(r)}(x)\,dx\,=\,0.  
</math></center>
</math></center>


}}
}}


Zauważmy, że jeśli macierz <math>\displaystyle A</math> jest kwadratowa, <math>\displaystyle m=n</math>, to
{{dowod|||
rozwiązaniem jest <math>\displaystyle  x^*=A^{-1} b</math> i residuum jest zerem.
Dla <math>r=1</math> funkcja <math>s'</math> jest przedziałami stała.
Zadanie wygładzania liniowego jest więc uogólnieniem
Oznaczając przez <math>a_j</math> jej wartość na <math>[x_{j-1},x_j]</math> dostajemy
rozwiązywania kwadratowych układów równań liniowych.


Równanie ([[##unormal|Uzupelnic: unormal ]]) nazywa się układem ''normalnym''.
<center><math>\begin{align} \int_a^b f'(x)s'(x)\,dx &=
Może ono nam sugerować sposób rozwiązania zadania wygładzania
          \sum_{j=1}^n\int_{t_{j-1}}^{t_j} f'(x)a_j\,dx \\
liniowego. Wystarczy bowiem pomnożyć macierz <math>\displaystyle A^T</math> przez <math>\displaystyle A</math> i
        &= \sum_{j=1}^n a_j(f(x_j)-f(x_{j-1}))\,=\,0,
rozwiązać układ normalny. Zauważmy ponadto, że macierz <math>\displaystyle A^TA</math>
\end{align}</math></center>
jest symetryczna i dodatnio określona, bo <math>\displaystyle (A^TA)^T=A^TA</math> i dla
<math>\displaystyle  x\ne 0</math> mamy
<math>\displaystyle  x^T(A^TA) x=(A x)^T(A x)=\|A x\|_2>0</math>, przy
czym ostatnia nierówność wynika z faktu, że kolumny macierzy <math>\displaystyle A</math>
są liniowo niezależne i dlatego <math>\displaystyle A x\ne 0</math>. Przy mnożeniu
<math>\displaystyle A^T</math> przez <math>\displaystyle A</math> wystarczy więc obliczyć tylko elementy na głównej
przekątnej i pod nią, a do rozwiązania równania z macierzą
<math>\displaystyle A^TA</math> można zastosować algorytm Banachiewicza-Choleskiego
opisany w U. [[##BC|Uzupelnic: BC ]]. Jak łatwo się przekonać, koszt takiego
algorytmu wynosi <math>\displaystyle n^2(k+n/3)</math>, przy czym dominuje koszt mnożenia
obliczenia macierzy <math>\displaystyle A^TA</math>.


Ma on jednak pewne wady. Mnożenie macierzy powoduje w <math>\displaystyle fl_\nu</math>  
ponieważ <math>f</math> zeruje się w <math>t_j</math>.  
powstanie "po drodze" dodatkowych błędów, które mogą nawet
zmienić rząd macierzy. Na przykład, dla macierzy


<center><math>\displaystyle A\,=\,\left(\begin{array} {cccc}
Rozpatrzmy teraz przypadek <math>r\ge 2</math>. Ponieważ
    1  &  1  &  1  &  1  \\
 
  \epsilon  \\
<center><math>(f^{(r-1)}s^{(r)})'\,=\,f^{(r)}s^{(r)}\,+\,f^{(r-1)}s^{(r+1)},
        &\epsilon \\
        &    &\epsilon \\
        &    &      &\epsilon \end{array} \right)
</math></center>
</math></center>


mamy  
to
 
<center><math>\int_a^b f^{(r)}(x)s^{(r)}(x)\,dx\,=\,f^{(r-1)}(x)s^{(r)}(x)
    \Big |_a^b \,-\,\int_a^b f^{(r-1)}(x)s^{(r+1)}(x)\,dx</math></center>
 
Wobec tego, że <math>s</math> jest poza przedziałem <math>(a,b)</math> wielomianem stopnia
co najwyżej <math>r-1</math> oraz <math>s^{(r)}</math> jest ciągła na <math>R</math>, mamy  
<math>s^{(r)}(a)=0=s^{(r)}(b)</math>, a stąd


<center><math>\displaystyle A^TA\,=\,\left(\begin{array} {cccc}
<center><math>f^{(r-1)}(x)s^{(r)}(x)\Big |_a^b\,=\,0.  
    1+\epsilon^2 & 1 & 1 & 1 \\
    1 & 1+\epsilon^2 & 1 & 1 \\
    1 & 1 & 1+\epsilon^2 & 1 \\
    1 & 1 & 1 & 1+\epsilon^2 \end{array} \right).
</math></center>
</math></center>


Jeśli <math>\displaystyle \epsilon^2<\nu</math> to <math>\displaystyle fl_\nu(1+\epsilon^2)=1</math>, co implikuje
Postępując podobnie, tzn. całkując przez części <math>r-1</math> razy,  
<math>\displaystyle  \mbox{rank} (fl_\nu(A^TA))=1</math>, podczs gdy <math>\displaystyle  \mbox{rank} (fl_\nu(A))=4</math>.
otrzymujemy w końcu
 
<center><math>\int_a^b f^{(r)}(x)s^{(r)}(x)\,dx\,=\,(-1)^{r-1}
      \int_a^b f'(x)s^{(2r-1)}(x)\,dx</math></center>


Poniżej przedstawimy inną metodę rozwiązywania zadania
Funkcja <math>s^{(2r-1)}</math> jest przedziałami stała, a więc możemy
wygładzania liniowego, która oparta jest na specjalnych
teraz zastosować ten sam argument jak dla <math>r=1</math>, aby pokazać,
przekształceniach zwanych odbiciami Householdera.  
że ostatnia całka jest równa zeru.}}


===Odbicia Householdera===
Funkcje sklejane chcielibyśmy zastosować do interpolacji funkcji.
Ważne jest więc, aby odpowiednie zadanie interpolacyjne miało
jednoznaczne rozwiązanie.


Dla danego wektora <math>\displaystyle  w\inR^m</math> o normie
{{twierdzenie|O istnieniu i jednoznaczności naturalnego splajnu interpolacyjnego|O istnieniu i jednoznaczności naturalnego splajnu interpolacyjnego|
<math>\displaystyle \| w\|_2=\sqrt{ w^T w}=1</math>,
''odbicie (macierz) Householdera'' zdefiniowane jest jako


<center><math>\displaystyle H\,=\,I\,-\,2 w w^T.
Jeśli <math>n+1\ge r</math>, to dla dowolnej funkcji
<math>f:[a,b]\to R</math> istnieje dokładnie jedna naturalna funkcja sklejana
<math>s_f\in {\cal S}_r</math> interpolująca <math>f</math> w węzłach <math>x_j</math>, tzn. taka,
że
 
<center><math>s_f(x_j)\,=\,f(x_j),\qquad 0\le j\le n</math></center>
 
}}
 
{{dowod|||
Pokażemy najpierw, że jedyną naturalną
funkcją sklejaną interpolującą dane zerowe jest funkcja zerowa.
Rzeczywiście, jeśli <math>s(x_j)=0</math> dla <math>0\le j\le n</math>, to podstawiając
w poprzednim lemacie <math>f = s</math>, otrzymujemy
 
<center><math>\int_a^b \Big(s^{(r)}(x)\Big)^2\,dx\,=\,0</math></center>
 
Stąd <math>s^{(r)}</math> jest funkcją zerową, a więc <math>s</math> jest
wielomianem stopnia co najwyżej <math>r-1</math> zerującym się w co
najmniej <math>n+1</math> punktach <math>x_j</math>. Wobec tego, że <math>n+1>r-1</math>,
otrzymujemy <math>s\equiv 0</math>.
 
Zauważmy teraz, że problem znalezienia naturalnej funkcji sklejanej
<math>s</math> interpolującej <math>f</math> można sprowadzić do rozwiązania układu równań liniowych z macierzą kwadratową. Na każdym przedziale <math>[x_{i-1},x_i]</math>,
<math>1\le i\le n</math>, jest ona postaci
 
<center><math>s(x)\,=\,w_i(x)\,=\,\sum_{j=0}^{2r-1} a_{i,j}x^j,
</math></center>
</math></center>


Zauważmy, że
dla pewnych współczynników <math>a_{i,j}\in R</math>, a na
<math>(-\infty,a]</math> i <math>[b,\infty)</math> mamy odpowiednio


<center><math>\displaystyle H x\,=\, x\,-\,2( w^T x) w,
<center><math>s(x)\,=\,w_0(x)\,=\,\sum_{j=0}^{r-1} a_{0,j}x^j
</math></center>
</math></center>


a ponieważ <math>\displaystyle ( w^T x) w=( x, w)_2 w</math>
i
jest rzutem prostopadłym <math>\displaystyle  x</math> na kierunek wektora <math>\displaystyle  w</math>
(<math>\displaystyle (\cdot,\cdot)_2</math> oznacza iloczyn skalarny), to <math>\displaystyle H x</math> jest
odbiciem lustrzanym wektora <math>\displaystyle  x</math> względem hiperpłaszczyzny
(wymiaru <math>\displaystyle m-1</math>) prostopadłej do <math>\displaystyle  w</math>.


Odbicia Householdera są przekształceniami nieosobliwymi
<center><math>s(x)\,=\,w_{n+1}(x)\,=\,\sum_{j=0}^{r-1} a_{n+1,j}x^j.
spełniającymi
</math></center>


<center><math>\displaystyle H^{-1}\,=\,H\,=\,H^T.
Aby wyznaczyć <math>s</math>, musimy więc znaleźć ogółem <math>2r(n+1)</math>
współczynników <math>a_{i,j}</math>, przy czym są one związane
<math>(2r-1)(n+1)</math> warunkami jednorodnymi wynikającymi z gładkości,
 
<center><math>w_i^{(k)}(x_i)\,=\,w_{i+1}^{(k)}(x_i)
</math></center>
</math></center>


Rzeczywiście, ponieważ <math>\displaystyle  w</math> ma normę jednostkową, mamy
dla <math>0\le i\le n</math> i <math>0\le k\le 2r-2</math>, oraz <math>n+1</math>
niejednorodnymi warunkami interpolacyjnymi,
 
<center><math>w_i(x_i)\,=\,f(x_i)
</math></center>
 
dla <math>0\le i\le n</math>. Otrzymujemy więc układ <math>2r(n+1)</math>
równań liniowych ze względu na <math>2r(n+1)</math> niewiadomych
<math>a_{i,j}</math>.
 
Naturalna funkcja sklejana interpolująca <math>f</math> jest wyznaczona jednoznacznie wtedy i tylko wtedy, gdy układ ten ma jednoznaczne rozwiązanie. To zaś zachodzi, gdy zero jest jedynym rozwiązaniem układu jednorodnego. Rzeczywiście, układ jednorodny odpowiada zerowym warunkom interpolacyjnym, przy których, jak pokazaliśmy wcześniej, zerowa funkcja sklejana (której odpowiada <math>a_{i,j}=0</math>, <math>\forall i,j</math>) jest jedynym rozwiązaniem zadania interpolacyjnego.}}
 
Naturalnych funkcji sklejanych możemy więc używać do
interpolacji funkcji. Pokażemy teraz inną ich własność,
która jest powodem dużego praktycznego zainteresowania
funkcjami sklejanymi.
 
{{twierdzenie|O ekstremalnej własności splajnów naturalnych|O ekstremalnej własności splajnów naturalnych|
 
Niech <math>f\in W^r(a,b)</math> i niech <math>s_f\in {\cal S}_r</math>
będzie naturalną funkcją sklejaną rzędu <math>r</math> interpolującą
<math>f</math> w węzłach <math>x_j</math>, <math>0\le j\le n</math>. Wtedy


<center><math>\displaystyle H^2 \,=\, (I-2 w w^T)^2\,=\,
<center><math>
  I-4 w w^T+4 w( w^T w) w^T \,=\, I,
  \int_a^b \Big( f^{(r)}(x)\Big)^2\,dx\,\ge\,
  \int_a^b \Big(s_f^{(r)}(x)\Big)^2\,dx.
</math></center>
</math></center>
}}
{{dowod|||
Jeśli przedstawimy <math>f</math> w postaci
<math>f=s_f+(f-s_f)</math>, to
<center><math>\begin{align} \int_a^b\Big(f^{(r)}(x)\Big)^2\,dx 
    &= \int_a^b\Big(s_f^{(r)}(x)\Big)^2\,dx\,+\,
        \int_a^b\Big((f-s_f)^{(r)}(x)\Big)^2\,dx \\
    & & \qquad\qquad + 2\,\int_a^b s_f^{(r)}(x)(f-s_f)^{(r)}(x)\,dx.
\end{align}</math></center>
Funkcja <math>f-s_f</math> jest w klasie <math>W^r(a,b)</math> i zeruje się w węzłach
<math>x_j</math>, <math>0\le j\le n</math>. Z lematu wynika więc, że
<math>\int_a^b s_f^{(r)}(x)(f-s_f)^{(r)}(x)dx=0</math>, a stąd wynika teza.}}
Wartość całki <math>\int_a^b(f^{(r)}(x))^2 dx</math> może być w
ogólności uważana za miarę gładkości funkcji.  Dowiedzioną nierówność możemy więc zinterpretować w następujący sposób. Naturalna funkcja sklejana jest w
klasie <math>W^r(a,b)</math> najgładszą funkcją spełniającą dane
warunki interpolacyjne w wybranych węzłach <math>x_j</math>.
Jak już wspomnieliśmy, najczęściej używanymi są kubiczne
funkcje sklejane. Dlatego rozpatrzymy je oddzielnie.
==Kubiczne funkcje sklejane==
Jeśli zdecydowaliśmy się na użycie kubicznych funkcji
sklejanych, powstaje problem wyznaczenia <math>s_f\in{\cal S}_2</math>
interpolującej daną funkcję <math>f</math>, tzn. takiej, że
<math>s_f(x_i)=f(x_i)</math> dla <math>0\le i\le n</math>. W tym celu, na każdym
przedziale <math>[x_i,x_{i+1}]</math> przedstawimy <math>s_f</math> w postaci jej
rozwinięcia w szereg Taylora w punkcie <math>x_i</math>,
<center><math>s_f(x)\,=\,w_i(x)\,=\,a_i+b_i(x-x_i)+
    c_i\frac{(x-x_i)^2}2+d_i\frac{(x-x_i)^3}6</math>,</center>
i podamy algorytm obliczania <math>a_i,b_i,c_i,d_i</math>
dla <math>0\le i\le n-1</math>.
Warunki brzegowe i warunki ciągłości
dla <math>s_f''</math> dają nam <math>w_0''(x_0)=0=w_{n-1}''(x_n)</math> oraz
<math>w_i''(x_{i+1})=w_{i+1}''(x_{i+1})</math>, czyli
<center><math>\begin{align} c_0 &= 0, \\ 
  c_i+d_ih_i &= c_{i+1}, \qquad 0\le i\le n-2, \\
  c_{n-1}+d_{n-1}h_{n-1} &= 0,
\end{align}</math></center>
gdzie <math>h_i=x_{i+1}-x_i</math>. Stąd, przyjmując dodatkowo
<math>c_n=0</math>, otrzymujemy
<center><math>
  d_i\,=\,\frac{c_{i+1}-c_i}{h_i},\qquad 1\le i\le n-1</math></center>
Z warunków ciągłości dla <math>s_f'</math> dostajemy z kolei
<center><math>b_i+c_ih_i+d_i\frac{h_i^2}{2}\,=\,b_{i+1},
    \qquad 0\le i\le n-2</math>,</center>


oraz  
oraz  


<center><math>\displaystyle H^T\,=\,(I-2 w w^T)^T\,=\,I-2( w^T)^T w^T\,=\,I.  
<center><math>
  b_{i+1}\,=\,b_i+h_i\frac{c_{i+1}+c_i}{2},
    \qquad 0\le i\le n-2.  
</math></center>
</math></center>


W szczególności <math>\displaystyle H</math> jest więc przekształceniem ''ortogonalnym'',
Warunki ciągłości <math>s_f</math> dają w końcu
<math>\displaystyle H^{-1}=H^T</math>, czyli nie zmienia długości wektora,


<center><math>\displaystyle \|H x\|_2\,=\,\sqrt{(H x)^T(H x)}\,=\,
<center><math>
    \sqrt{ x^T(H^TH) x}\,=\,\sqrt{ x^T x}\,=\,
  a_i+b_ih_i+c_i\frac{h_i^2}{2}+d_i\frac{h_i^3}{6}\,=\,a_{i+1},  
     \| x\|_2.
     \qquad 0\le i\le n-2.  
</math></center>
</math></center>


Odbicia Householdera zastosujemy do przeprowadzenia danego wektora
Powyższe równania definiują nam na odcinku <math>[a,b]</math> naturalną kubiczną funkcję
<math>\displaystyle  x\ne 0</math> na kierunek innego niezerowego wektora, powiedzmy
sklejaną. Ponieważ poszukiwana funkcja sklejana <math>s_f</math> ma interpolować <math>f</math>, mamy dodatkowych <math>n+1</math> warunków interpolacyjnych <math>w_i(x_i)=f(x_i)</math>, <math>0\le i\le n-1</math>,  
<math>\displaystyle  e</math>, tzn.
oraz <math>w_{n-1}(x_n)=f(x_n)</math>, z których
 
<center><math>
  a_i\,=\,f(x_i), \qquad 0\le i\le n-1</math></center>
 
Teraz możemy warunki ciągłości przepisać jako


<center><math>\displaystyle H x\,=\,(I-2 w w^T) x\,=\,\alpha\, e.
<center><math>f(x_{i+1})\,=\,f(x_i)+b_ih_i+c_ih_i^2+d_i\frac{h_i^3}{6}</math>,</center>
 
przy czym wzór ten zachodzi również dla <math>i=n-1</math>.
Po wyrugowaniu <math>b_i</math> i podstawieniu <math>d_i</math>,
mamy
 
<center><math>
b_i\,=\,f(x_i,x_{i+1})+h_i\frac{c_{i+1}+2c_i}{6}, \qquad 0\le i\le n-1,  
</math></center>
</math></center>


<!--
gdzie <math>f(x_i,x_{i+1})</math> jest odpowiednią różnicą
[[Image:MNhouseholderidea.png|thumb|400px|Odbicie Househodera]]
dzieloną. Możemy teraz powyższe wyrażenie na <math>b_i</math>  
-->
podstawić, aby otrzymać
 
<div class="thumb tright"><div><flash>file=householderidea.swf</flash><div.thumbcaption>Odbicie Househodera</div></div></div>
<center><math>c_i\frac{h_i}{6}+c_{i+1}\frac{h_i+h_{i+1}}{3}+c_{i+1}\frac{h_{i+1}}{6}
  \,=\,f(x_{i+1},x_{i+2})-f(x_i,x_{i+1}).
</math></center>


Załóżmy dla uproszczenia, że <math>\displaystyle \| e\|_2=1</math>.
Wprowadzając oznaczenie
Aby wyznaczyć <math>\displaystyle H</math> zauważmy, że


<center><math>\displaystyle w\,=\,\frac{ x-\alpha e}{2( w^T x)},  
<center><math>
c_i^*\,=\,\frac{c_i}{6},  
</math></center>
</math></center>


a ponieważ <math>\displaystyle \alpha=\pm\| x\|_2</math> i <math>\displaystyle \| w\|_2=1</math> to  
możemy to równanie przepisać jako


<center><math>\displaystyle w\,=\,\frac{ x\mp\| x\|_2 e}
<center><math>\frac{h_i}{h_i+h_{i+1}}c_i^*\,+\,2\,c_{i+1}^*\,+\,
                {\| x\mp\| x\|_2 e\|_2}.
  \frac{h_{i+1}}{h_i+h_{i+1}}c_{i+2}^*\,=\,
  f(x_i,x_{i+1},x_{i+2}),
</math></center>
</math></center>


W szczególności, jeśli <math>\displaystyle  e= e_1</math> jest pierwszym
<math>0\le i\le n-2</math>, albo w postaci macierzowej
wersorem to powyższe wzory dają


<center><math>\displaystyle H\,=\,I\,-\,\frac{ u u^T}{\gamma},  
<center><math>
</math></center>
\left(\begin{array} {cccccc}
    2  & w_1 \\
  u_2 &  2  & w_2 \\
      & u_3 &  2    & w_3 \\
      &    &\ddots &\ddots  &\ddots \\
      &    &      & u_{n-2} &  2  & w_{n-2} \\
      &    &      &        & u_{n-1} &  2
    \end{array} \right)
  \left(\begin{array} {c}
    c_1^*\\ c_2^*\\ c_3^*\\ \vdots\\ c_{n-2}^*\\ c_{n-1}^*
      \end{array} \right) \,=\,
  \left(\begin{array} {c}
    v_1\\ v_2\\ v_3\\ \vdots\\ v_{n-2}\\ v_{n-1}
      \end{array} \right)</math>,</center>


gdzie  
gdzie  


<center><math>\displaystyle u_i\,=\,\left\{\begin{array} {ll}
<center><math>\begin{align} && u_i\,=\,\frac{h_{i-1}}{h_{i-1}+h_i},\qquad
        x_1\mp\| x\|_2  &\quad i=1, \\
    w_i\,=\,\frac{h_i}{h_{i-1}+h_i}, \\
                      x_i &\quad 2\le i\le m,  
  && v_i\,=\,f(x_{i-1},x_i,x_{i+1}).
        \end{array} \right.
\end{align}</math></center>
 
Ostatecznie, aby znaleźć współczynniki <math>a_i,b_i,c_i,d_i</math>
należy najpierw rozwiązać układ równań liniowych,
a potem zastosować wzory definiujące pozostałe współczynniki.
 
Zauważmy, że macierz układu równań liniowych jest trójdiagonalna i ma dominującą przekątną. Układ  można więc rozwiązać kosztem proporcjonalnym do wymiaru <math>n</math> używając [[MN08#Macierze trójdiagonalne|algorytmu przeganiania]]. Koszt znalezienia wszystkich współczynników kubicznej funkcji sklejanej interpolującej <math>f</math> jest więc też proporcjonalny do <math>n</math>.
 
MATLAB i Octave mają wbudowaną funkcję wyznaczającą naturalny kubiczny splajn interpolujący zadane wartości:
 
<div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>s = spline(x,y);
</pre></div>
Aby wyznaczyć wartości takiego splajnu w zadanych punktach <math>X</math>, także musimy użyć specjalnej funkcji,
 
<div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>Y = ppval(s,X);
</pre></div>
Na końcu oszacujemy jeszcze błąd interpolacji naturalnymi
kubicznymi funkcjami sklejanymi na przedziale <math>[a,b]</math>.
Będziemy zakładać, że <math>f</math> jest dwa razy
różniczkowalna w sposób ciągły.
 
{{twierdzenie|O błędzie interpolacji splajnem kubicznym|O błędzie interpolacji splajnem kubicznym|
 
Jeśli <math>f\in F^1_M([a,b])</math> to
 
<center><math>\|f-s_f\|_{C([a,b])}\,\le\,5\,M\,
      \max_{1\le i\le n}(x_i-x_{i-1})^2</math></center>
 
W szczególności, dla podziału równomiernego
<math>x_i=a+i\frac{b-a}{n}</math>, <math>0\le i\le n</math>, mamy
 
<center><math>\|f-s_f\|_{ C([a,b])}\,\le\,
      5\,M\,\Big(\frac{b-a}n\Big)^2.  
</math></center>
</math></center>
}}
{{dowod|||
Wykorzystamy obliczoną wcześniej 
postać interpolującej funkcji sklejanej <math>s_f</math>.
Dla <math>x\in [x_i,x_{i+1}]</math> mamy
<center><math>\begin{align} w_i(x) &= f(x_i)\,+\,\left(\frac{f(x_{i+1})-f(x_i)}{h_i} -h_i(c_{i+1}^*+2c_i^*)\right)(x-x_i) \\ &&\qquad\qquad \,+\, 3c_i^*(x-x_i)^2\,+\,\frac{c_{i+1}^*-c_i^*}{h_i}(x-x_i)^3.
\end{align}</math></center>
Z rozwinięcia <math>f</math> w szereg Taylora w punkcie <math>x_i</math> dostajemy
<math>f(x)=f(x_i)+f'(x_i)(x-x_i)+f''(\xi_1)(x-x_i)^2/2</math> oraz
<math>(f(x_{i+1})-f(x_i))/h_i=f'(x_i)+h_if''(\xi_2)/2</math>. Stąd
<center><math>\begin{align} f(x)-s_f(x) \,=\, f(x)-w_i(x) \\
  &= \frac{f''(\xi_1)}2(x-x_i)^2-\left(\frac{f''(\xi_2)}2
      -(c_{i+1}^*+2c_i^*)\right)h_i(x-x_i) \\
  & & \qquad\qquad\qquad -3c_i^*(x-x_i)^2
      -\frac{c_{i+1}^*-c_i^*}{h_i}(x-x_i)^3,
\end{align}</math></center>


oraz  
oraz  


<center><math>\displaystyle \aligned \gamma &= \frac 12\| u\|_2^2\,=\,
<center><math>\begin{align}
      \frac 1 2\Big((x_1\mp\| x\|_2)^2+\sum_{i=2}^m x_i^2\Big) \\
|f(x)-s_f(x)| &\le &(M+2|c_{i+1}^*|+6|c_i^*|)h_i^2 \\
  &= \frac 1 2 \Big(\sum_{i=1}^m x_i^2\,+\,\| x\|_2^2\,\mp\,
    &= (M+8\max_{1\le i\le n-1}|c_i^*|)h_i^2.  
      2 x_1\|x\|_2\Big) \,=\,\|x\|_2^2\,\mp\,x_1 \|x\|_2.  
\end{align}</math></center>
\endaligned</math></center>


Otrzymaliśmy dwa odbicia Householdera przekształcające dany wektor
Niech teraz <math>\max_{1\le i\le n-1}|c_i^*|=|c_s^*|</math>.  
<math>\displaystyle  x</math> na kierunek pierwszego wersora, w zależności od wybranego
Z postaci układu otrzymujemy
znaku przy <math>\displaystyle \| x\|_2</math>. Ustalimy ten znak na plus gdy <math>\displaystyle x_1\ge 0</math>,
oraz na minus gdy <math>\displaystyle x_1<0</math>, co pozwoli na obliczenie <math>\displaystyle u_1</math> i <math>\displaystyle \gamma</math>
z małym błędem względem w <math>\displaystyle fl_\nu</math>. Wtedy bowiem mamy


<center><math>\displaystyle u_1\,=\,\left\{\begin{array} {ll}
<center><math>\begin{align} |c_s^*| &= 2|c_s^*|-|c_s^*|(u_s+w_s) \,\le\,
       x_1+\|x\|_2 & \quad x_1\ge 0, \\
    |u_sc_{s-1}^*+2c_s^*+w_sc_{s+1}| \\
      x_1-\|x\|_2 & \quad x_1<0, \end{array} \right.
  &= |f(x_{s-1},x_s,x_{s+1})|\,\le\,
       \Big|\frac{f''(\xi_3)}2\Big|\,\le\,\frac 12 M,
\end{align}</math></center>
 
a stąd
 
<center><math>|f(x)-s_f(x)|\,\le\, 5\, M\, h_i^2,
</math></center>
</math></center>


oraz <math>\displaystyle \gamma=\| x\|_2^2+|x_1|\,\| x\|_2</math>, czyli zawsze
co kończy dowód.
dodajemy liczby tych samych znaków. Ponadto pierwsza współrzędna
}}
wektora <math>\displaystyle H x</math> jest równa <math>\displaystyle -\| x\|_2</math> dla <math>\displaystyle x_1\ge 0</math>, oraz
 
<math>\displaystyle +\| x\|_2</math> dla <math>\displaystyle x_1<0</math>.  
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
 
<span  style="font-variant:small-caps;">Przykład</span>
===Rozkład QR===
<div class="solution" style="margin-left,margin-right:3em;">
 
Porównanie interpolacji splajnowej i Lagrange'a.
 
[[Image:MNrunge17czebyspline.png|thumb|550px|center|Interpolacja splajnowa wydaje się lepiej spełniać zadanie odtworzenia kształtu funkcji]]
 
Jak widać, w przeciwieństwie do wielomianu interpolacyjnego, splajn interpolacyjny praktycznie pokrywa się z wykresem funkcji, tutaj: <math>f(x) = \frac{1}{1+x^2}</math>.
 
</div></div>
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Uwaga</span>
<div class="solution" style="margin-left,margin-right:3em;">
 
Niech
 
<center><math>W^r_M(a,b)\,=\,\Big\{\,f\in W^r(a,b):\,
    \int_a^b\left(f^{(r)}(x)\right)^2\,dx\le M\,\Big\}</math></center>
 
Ustalmy węzły <math>a=x_0<\cdots<x_n=b</math>. Dla <math>f\in W^r_M(a,b)</math>, niech <math>s_f</math> będzie naturalną funkcją sklejaną interpolującą <math>f</math> w <math>x_j</math>, <math>0\le j\le n</math>, a <math>a_f</math> dowolną inną aproksymacją korzystającą jedynie z informacji o wartościach <math>f</math> w tych węzłach, tzn.
 
<center><math>a_f\,=\,\phi(f(x_0),\ldots,f(x_n))</math></center>


Odbić Householdera można użyć do rozkładu macierzy
Załóżmy, że błąd aproksymacji mierzymy nie w normie Czebyszewa, ale w normie średniokwadratowej zdefiniowanej jako
<math>\displaystyle A\inR^{m\times n}</math> na iloczyn ortogonalno-trójkątny.


Niech <math>\displaystyle A=( a_1, a_2,\ldots, a_n)</math>, gdzie <math>\displaystyle  a_j</math> są
<center><math>\|g\|_{{\cal L}_2(a,b)}\,=\,\sqrt{\int_a^b (g(x))^2\,dx}.
wektorami-kolumnami macierzy <math>\displaystyle A</math>. Wybierzmy pierwsze odbicie
</math></center>
Householdera <math>\displaystyle H_1=I_m- u_1 u_1^T/\gamma_1</math> tak, aby
przekształcało pierwszy wektor-kolumnę macierzy <math>\displaystyle A</math> na kierunek
<math>\displaystyle  e_1</math>. Efektem pomnożenia macierzy <math>\displaystyle A</math> z lewej strony przez
<math>\displaystyle H_1</math> będzie wtedy macierz


<center><math>\displaystyle A^{(1)}\,=\,( a^{(1)}_1,\ldots, a^{(1)}_n)
Wtedy
  \,=\,(H_1 a_1,\ldots, H_1 a_n),
 
<center><math>\sup_{f\in W^r_M(a,b)}\|f-s_f\|_{{\cal L}_2(a,b)}\,\le\,
  \sup_{f\in W^r_M(a,b)}\|f-a_f\|_{{\cal L}_2(a,b)}.
</math></center>
</math></center>


w której pierwsza kolumna <math>\displaystyle  a^{(1)}_1</math> ma niezerową tylko
Aproksymacja naturalnymi funkcjami sklejanymi jest więc optymalna w klasie <math>W^r_M(a,b)</math>.  
pierwszą współrzędną. W następnym kroku wybieramy drugie
przekształcenie Householdera 
<math>\displaystyle \bar H_2=I_{m-1}- v_2 v_2^T/\gamma_2</math> wymiaru <math>\displaystyle m-1</math> tak,
aby przeprowadzało wektor <math>\displaystyle (a^{(1)}_{i,2})_{i=2}^m</math> na kierunek
pierwszego wersora w <math>\displaystyle R^{m-1}</math>. Rozszerzając <math>\displaystyle  v_2\inR^{m-1}</math>
do wektora <math>\displaystyle  u_2\inR^m</math> przez dodanie zera jako pierwszej
współrzędnej, <math>\displaystyle  u_2=(0, v_2)^T</math>, otrzymujemy
przekształcenie (macierz) Householdera
<math>\displaystyle H_2=I_m- u_2 u_2^T/\gamma_2</math> w <math>\displaystyle R^m</math> postaci


<center><math>\displaystyle H_2\,=\,\left(\begin{array} {cccc}
Można również pokazać, że interpolacja <math>s_f^*</math> naturalnymi funkcjami sklejanymi na węzłach równoodległych <math>x_j=a+(b-a)j/n</math>, <math>0\le j\le n</math>, jest optymalna co do
    1 &  0^T \\
rzędu w klasie <math>W^r_M(a,b)</math>, wśród wszystkich aproksymacji korzystających jedynie z informacji o wartościach funkcji w <math>n+1</math> dowolnych punktach, oraz
      0 & \bar H_2  \end{array} \right).
 
<center><math>\max_{f\in W^r_M(a,b)}\|f-s^*_f\|_{{\cal L}_2(a,b)}
    \,\asymp\,n^{-r}.  
</math></center>
</math></center>


Pomnożenie macierzy <math>\displaystyle A^{(1)}</math> z lewej strony przez <math>\displaystyle H_2</math> spowoduje
</div></div>
teraz wyzerowanie drugiej kolumny macierzy pod elementem
 
<math>\displaystyle a^{(1)}_{2,2}</math>, przy czym pierwszy wiersz i pierwsza kolumna
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
pozostaną niezmienione. Postępując tak dalej <math>\displaystyle n</math> razy
<span  style="font-variant:small-caps;">Uwaga</span>
(albo <math>\displaystyle n-1</math> razy gdy <math>\displaystyle m=n</math>) otrzymujemy
<div class="solution" style="margin-left,margin-right:3em;">
 
Tak jak wielomiany, naturalne funkcje sklejane interpolujące dane funkcje można
reprezentować przez ich współczynniki w różnych bazach. Do tego celu można na przykład użyć bazy kanonicznej <math>K_j</math>, <math>0\le j\le n</math> zdefiniowanej równościami
 
<center><math>K_j(x_i)\,=\,\left\{\begin{array} {ll}
        0 &\quad i\ne j, \\
        1 &\quad i=j, \end{array} \right.</math></center>
 
przy której <math>s_f(x)=\sum_{j=0}^n f(x_j)K_j(x)</math>. Baza kanoniczna jest jednak niewygodna w użyciu, bo funkcje <math>K_j</math> w ogólności nie zerują się na żadnym podprzedziale, a tym samym manipulowanie nimi jest utrudnione.


<center><math>\displaystyle H_nH_{n-1}\cdots H_2H_1A\,=\,R,
Częściej używa się bazy B-sklejanej. Można ją zdefiniować dla splajnów dowolnego rzędu za pomocą wzoru rekurencyjnego (przyjmując, że dla dowolnego <math>i</math>, <math>x_i < x_{i+1}</math>):
 
<center><math>B_i^0(x) = \left\{ \begin{matrix}  1, &  \mbox{ jeśli }  x_i \leq x < x_{i+1},\\
0, &  \mbox{ w przeciwnym przypadku} ;
\end{matrix}
\right.</math></center>
 
<center><math>B_i^r(x) = \frac{x-x_i}{x_{i+r}-x_i} B_i^{r-1}(x) +
\frac{x_{i+r+1}-x}{x_{i+r+1}-x_{i+1}} B_{i+1}^{r-1}(x)</math></center>
 
W przypadku naturalnych splajnów kubicznych, <math>r=2</math>, baza B-sklejana jest jawnie zdefiniowana przez następujące warunki:
 
<center><math>\begin{align} B_j(x_j) &= 1, \qquad \mbox{ dla  }  0\le j\le n, \\
    B_j(x) &= 0,\qquad \mbox{ dla  }  x\le x_{j-2}, j\ge 2,
      \mbox{  oraz  dla  }  x\ge x_{j+2}, j\le n-2. 
\end{align}</math></center>
 
Dla <math>B_0</math> i <math>B_1</math> dodatkowo żądamy, aby
 
<center><math>B_0''(x_0)\,=\,0\,=\,B_1''(x_0),\qquad B_1(x_0)\,=\,0,  
</math></center>
</math></center>


gdzie <math>\displaystyle R\inR^{m\times n}</math> jest uogólnioną macierzą trójkątną
a dla <math>B_{n-1}</math> i <math>B_n</math> podobnie
górną, tzn. <math>\displaystyle r_{i,j}=0</math> dla <math>\displaystyle i>j</math>. Stąd, podstawiając
 
<math>\displaystyle Q=H_1H_2\cdots H_n</math>, dostajemy rozkład macierzy na iloczyn
<center><math>B_{n-1}''(x_n)\,=\,0\,=\,B_n''(x_n),
ortogonalno-trójkątny
      \qquad B_{n-1}(x_{n-1})\,=\,0</math></center>
 
Wtedy <math>B_j</math> nie zeruje się tylko na przedziale
<math>(x_{j-2},x_{j+2})</math>. Wyznaczenie współczynników
rozwinięcia w bazie <math>\{B_i\}_{i=0}^n</math> funkcji sklejanej
interpolującej <math>f</math> wymaga rozwiązania układu liniowego
z macierzą trójdiagonalną <math>\{B_j(x_i)\}_{i,j=0}^n</math>, a
więc koszt obliczenia tych współczynników jest
proporcjonalny do <math>n</math>.
</div></div>


<center><math>\displaystyle
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
   A\,=\,Q\cdot R.  
<span  style="font-variant:small-caps;">Uwaga</span>
<div class="solution" style="margin-left,margin-right:3em;">
 
Oprócz naturalnych funkcji sklejanych często rozpatruje się też <strong>okresowe funkcje sklejane</strong>. Są to funkcje <math>\tilde s:R\to R</math> spełniające  warunki '''(i)''', '''(ii)''' \link{splinedef}{definicji funkcji sklejanej}, oraz warunek:
    
; (iii)'
:  <math>\tilde s^{(i)}</math> jest dla <math>0\le i\le r-1</math>
funkcją okresową o okresie <math>(b-a)</math>, tzn.
<math>\tilde s^{(i)}(x)=\tilde s^{(i)}(x+(b-a))</math>, <math>\forall x</math>.
Klasę okresowych funkcji sklejanych rzędu <math>r</math> oznaczymy przez <math>\widetilde{\cal S}_r</math>. Funkcje te mają podobne własności jak naturalne funkcje sklejane. Dokładniej, niech
 
<center><math>
\tilde W^r(a,b)\,=\,\{\,f\in W^r(a,b):\,  f^{(i)}(a)=f^{(i)}(b),\; 0\le i\le r-1\,\}</math>,</center>
 
tzn. <math>\tilde W^r(a,b)</math> jest klasą funkcji z <math>W^r(a,b)</math>, które można przedłużyć do funkcji, krórych wszystkie pochodne do rzędu <math>r-1</math> włącznie są <math>(b-a)</math>-okresowe na <math>R</math>. Wtedy dla dowolnej funkcji <math>f\in\tilde W^r(a,b)</math>
zerującej się w węzłach <math>x_j</math>, oraz dla dowolnej  <math>\tilde s\in\widetilde{\cal S}_r</math> mamy
 
<center><math>
\int_a^b f^{(r)}(x)\tilde s^{(r)}(x)\,dx\,=\,0.  
</math></center>
</math></center>


Rzeczywiście, macierz <math>\displaystyle Q\inR^{m\times m}</math> jest ortogonalna, bo
Wynika z niego jednoznaczność rozwiązania zadania interpolacyjnego dla okresowych funkcji <math>f</math> (tzn. takich, że <math>f(a)=f(b)</math>), jak również odpowiednia własność minimalizacyjna okresowych funkcji sklejanych. Dokładniej, jeśli <math>f\in\tilde W^r(a,b)</math> oraz <math>\tilde s_f\in\widetilde{\cal S}_r</math> interpoluje <math>f</math> w węzłach <math>x_j</math>, <math>0\le j\le n</math>, to
 
<center><math>\int_a^b \Big(  f^{(r)}(x)\Big)^2\,dx\,\ge\,
  \int_a^b \Big(\tilde s_f^{(r)}(x)\Big)^2\,dx.
</math></center>


<center><math>\displaystyle \aligned Q^{-1} &= (H_1H_2\cdots H_n)^{-1}\,=\,
</div></div>
    H_n^{-1}\cdots H_2^{-1}H_1^{-1} \\
  &= H_n^T\cdots H_2^TH_1^T \,=\,
    (H_1H_2\cdots H_n)^T\,=\,Q^T.
\endaligned</math></center>


Dyspunując rozkładem ([[##orttr|Uzupelnic: orttr ]]) zadanie wygładzania liniowego
==Dygresja o najlepszej aproksymacji==
można rozwiązać następująco. Ponieważ mnożenie przez macierz
ortogonalną nie zmienia normy drugiej wektora, mamy


<center><math>\displaystyle \aligned \| r\|_2 &= \| b-A x\|_2\;=\;\| b-QR x\|_2 \\
Klasyczne zadanie aproksymacyjne w przestrzeniach funkcji definiuje się
    &= \|Q(Q^T b-R x)\|_2 \;=\;\| c-R x\|_2,
w następujący sposób.
\endaligned</math></center>


gdzie <math>\displaystyle  c=Q^T b=H_n\cdots H_2H_1 b</math>.
Niech <math>F</math> będzie pewną przestrzenią liniową funkcji <math>f:[a,b]\to R</math>, w której określona została norma <math>\|\cdot\|</math>. Niech <math>V_n\subset F</math> będzie podprzestrzenią
Rozbijając wektor <math>\displaystyle  c</math> na <math>\displaystyle  c=( c_I, c_{II})^T</math>,
w <math>F</math> wymiaru <math>n</math>. Dla danej <math>f\in F</math>, należy znaleźć funkcję <math>v_f\in F</math> taką, że
gdzie <math>\displaystyle  c_I\inR^n</math> i <math>\displaystyle  c_{II}\inR^{m-n}</math>, oraz macierz
<math>\displaystyle R</math> na


<center><math>\displaystyle R\,=\,\left(\begin{array} {c} R_I \\ 0\end{array} \right),
<center><math>\|f-v_f\|\,=\,\min_{v\in V_n}\|f-v\|.
</math></center>
</math></center>


gdzie <math>\displaystyle R_I\inR^{n\times n}</math> jest macierzą trójkątną górną, a
Okazuje się, że tak postawione zadanie ma rozwiązanie <math>v_f</math>, choć nie zawsze jest ono wyznaczone jednoznacznie.
<math>\displaystyle 0</math> jest macierzą zerową wymiaru <math>\displaystyle (m-n)\times n</math>, otrzymujemy


<center><math>\displaystyle \| r\|_2^2\;=\;\| c_I-R_I x\|_2^2\,+\,
Jako przykład, rozpatrzmy <math>F=W^r(a,b)</math>. Utożsamiając funkcje <math>f_1,f_2\in W^r(a,b)</math> takie, że <math>f_1(x)-f_2(x) \in\Pi_{r-1}</math>, zdefiniujemy w <math>W^r(a,b)</math> normę
    \| c_{II}\|_2^2.  
 
<center><math>\|f\|\,=\,\sqrt{\int_a^b\left(f^{(r)}(x)\right)^2\,dx}.
</math></center>
</math></center>


Rozwiązanie <math>\displaystyle  x^*</math> zadania wygładzania jest więc
Dla ustalonych węzłów <math>a=x_0<\cdots<x_n=b</math>, niech
rozwiązaniem układu liniowego trójkątnego,  


<center><math>\displaystyle x^*\,=\,R_I^{-1} c_I,
<center><math>V_{n+1}\,=\,{\cal S}_r
</math></center>
</math></center>


oraz <math>\displaystyle \| r^*\|_2=\| b-A x^*\|_2=\| c_{II}\|_2</math>.  
będzie podprzestrzenią w <math>W^r(a,b)</math> naturalnych funkcji sklejanych rzędu <math>r</math> opartych węzłach <math>x_j</math>, <math>0\le j\le n</math>. Oczywiście <math>\mbox{dim} {\cal S}_r=n+1</math>, co wynika z jednoznaczności rozwiązania w <math>{\cal S}_r</math> zadania interpolacji. Okazuje się, że wtedy optymalną dla <math>f\in W^r(a,b)</math> jest naturalna funkcja sklejana <math>s_f</math> interpolująca <math>f</math> w węzłach <math>x_j</math>, tzn.  


Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde
<center><math>\|f-s_f\|\,=\,\min_{s\in {\cal S}_r}\|f-s\|.  
z kolejnych przekształceń Householdera <math>\displaystyle H_k</math> wyznaczamy przez
</math></center>
obliczenie <math>\displaystyle \gamma_k</math> oraz współrzędnych wektora <math>\displaystyle  u_k</math>.
Wektor ten ma tylko <math>\displaystyle m-k+1</math> współrzędnych niezerowych, a ponadto
<math>\displaystyle u_{k,i}=a^{(k-1)}_{i,k}</math> dla <math>\displaystyle k+1\le i\le m</math>. Dzięki takiej
reprezentacji <math>\displaystyle H_k</math>, mnożenia <math>\displaystyle H_k x</math> możemy dla dowolnego
<math>\displaystyle  x</math> realizować według wzoru


<center><math>\displaystyle (H_k x)_i\,=\,x_i\,-\,s\,u_{k,i},  
Rzeczywiście, ponieważ norma w przestrzeni <math>W^r(a,b)</math> generowana jest przez iloczyn skalarny
 
<center><math>(f_1,f_2)\,=\,
    \int_a^b f_1^{(r)}(x)f_2^{(r)}(x)\,dx,  
</math></center>
</math></center>


gdzie <math>\displaystyle s= u_k^T x/\gamma_k</math>.
jest to przestrzeń unitarna. Znane twierdzenie mówi, że w przestrzeni unitarnej najbliższą danej <math>f</math> funkcją w dowolnej domkniętej podprzestrzeni <math>V</math> jest rzut
prostopadły <math>f</math> na <math>V</math>, albo równoważnie, taka funkcja <math>v_f\in V_{n+1}</math>, że iloczyn skalarny


Uwzględnizjąc obecność zerowych elementów w <math>\displaystyle  u_k</math>,
<center><math>(f-v_f, v)\,=\,0, \qquad\forall v\in V.  
przejście od macierzy <math>\displaystyle A^{(k-1)}</math> do <math>\displaystyle A^{(k)}</math> kosztuje rzędu
</math></center>
<math>\displaystyle 4(m-k+1)(n-k)</math> operacji arytmetycznych i obliczenie jednego
 
pierwiastka kwadratowego. Cały rozkład <math>\displaystyle A=QR</math> kosztuje więc
W naszym przypadku, ostatnia równość jest równoważna
rzędu (dla dużych <math>\displaystyle m</math> i <math>\displaystyle n</math>)


<center><math>\displaystyle \sum_{k=1}^n 4(m-k+1)(n-k)\,\approx\,\frac 43n^3+2n^2(m-n)
<center><math>\int_a^b(f-v_f)^{(r)}(x)s^{(r)}(x)\,dx\,=\,0,
  \,=\,2n^2(m-n/3)
    \quad\forall s\in{\cal S}_r.
</math></center>
</math></center>


operacji arytmetycznych i <math>\displaystyle n</math> pierwiastków kwadratowych. Zauważmy,
To zaś jest prawdą, gdy <math>v_f</math> interpoluje <math>f</math> w punktach <math>x_j</math>, czyli <math>v_f=s_f</math>.
że w przypadku <math>\displaystyle m=n</math>, a więc dla kwadratowego układu równań,
 
koszt ten wynosi <math>\displaystyle (4/3)n^3</math> i jest dwa razy większy od kosztu
Dodajmy jeszcze, że nie zawsze interpolacja daje najlepszą
eliminacji Gaussa.  
aproksymację w sensie klasycznym.


===Uwarunkowanie===
==Literatura==


===Biblioteki===
W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 6.4</b> w
* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.
Warto także zapoznać się (nieobowiązkowo) z rozdziałami&nbsp;6.5 i&nbsp;6.6 tamże.

Aktualna wersja na dzień 11:16, 12 wrz 2023


Funkcje sklejane (splajny)

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Interpolacja wielomianami interpolacyjnymi, chociaż korzysta z funkcji gładkich i łatwo reprezentowalnych w komputerze, ma jednak również pewne wady. Zauważmy, że błąd interpolacji może być bardzo duży (zjawisko Rungego), a poza tym interpolacja jest nielokalna: nawet mała zmiana warości funkcji w pojedynczym węźle może powodować dużą zmianę zachowania całego wielomianu interpolacyjnego. Czasem więc lepiej jest zastosować innego rodzaju interpolację, np. posługując się funkcjami sklejanymi, które tylko lokalnie są wielomianami, sklejonymi w taki sposób, by globalnie zachować pewien stopień gładkości, tzn. różniczkowalność zadaną liczbę razy.

Tego typu podejście okazało się bardzo owocne m.in. w grafice komputerowej (np. dla wizualizacji scenerii w grach komputerowych), a także np. posłużyło jako narzędzie konstrukcji skalowalnych czcionek komputerowych w Postscripcie (precyzyjniej, korzysta się tam z tzw. krzywych Beziera --- pewnych krzywych sklejanych zadanych na płaszczyźnie). Z krzywych Beziera powszechnie korzysta się również w systemach CAD (Computer Aided Design).

Zamiast terminu funkcje sklejane używa się też często terminów splajny (spline), albo funkcje gięte. Nazwy te biorą się stąd, że zadanie interpolacji naturalnym splajnem kubicznym można interpretować jako model matematyczny aparatu służącego do wytwarzania mebli giętych.

Funkcje sklejane

W ogólności przez funkcję sklejaną rozumie się każdą funkcję przedziałami wielomianową. Nas będą jednak interesować szczególne funkcje tego typu i dlatego termin funkcje sklejane zarezerwujemy dla funkcji przedziałami wielomianowych i posiadających dodatkowe własności, które teraz określimy.

Niech dany będzie przedział skończony [a,b] i węzły

a=x0<x1<<xn=b,

przy czym n1.

Definicja

Funkcję s:RR nazywamy funkcją sklejaną rzędu r (r1) odpowiadającą węzłom xj, 0jn, jeśli spełnione są następujące dwa warunki:

(i)
s jest wielomianem stopnia co najwyżej 2r1 na każdym

z przedziałów [xj1,xj], tzn. s|[xj1,xj]Π2r1, 1jn,

(ii)
s jest (2r2)-krotnie różniczkowalna w sposób

ciągły na całej prostej, tzn. sC(2r2)(R).

Jeśli ponadto

(iii)
s jest wielomianem stopnia co najwyżej r1 poza

(a,b), tzn. s|(,a],s|[b,+)Πr1,

to s jest naturalną funkcją sklejaną.

Klasę naturalnych funkcji sklejanych rzędu r opartych na węzłach xj będziemy oznaczać przez 𝒮r(x0,,xn), albo po prostu 𝒮r, jeśli węzły są ustalone.

Na przykład funkcją sklejaną rzędu pierwszego (r=1) jest funkcja ciągła i liniowa na poszczególnych przedziałach [xj1,xj]. Jest ona naturalna, gdy poza (a,b) jest funkcją stała. Tego typu funkcje nazywamy liniowymi funkcjami sklejanymi.

Najważniejszymi z praktycznego punktu widzenia są jednak funkcje sklejane rzędu drugiego odpowiadające r=2. Są to funkcje, które są na R dwa razy różniczkowalne w sposób ciągły, a na każdym z podprzedziałów są wielomianami stopnia co najwyżej trzeciego. W tym przypadku mówimy o kubicznych funkcjach sklejanych. Funkcja sklejana kubiczna s jest naturalna, gdy poza (a,b) jest wielomianem liniowym, a więc s(a)=s(b)=0.

Interpolacja i gładkość

Pokażemy najpierw ważny lemat, który okaże się kluczem do dowodu dalszych twierdzeń.

Niech Wr(a,b) będzie klasą funkcji f:[a,b]R takich, że f jest (r1) razy różniczkowalna na [a,b] w sposób ciągły oraz f(r)(x) istnieje prawie wszędzie na [a,b] i jest całkowalna z kwadratem, tzn.

Wr(a,b)={fC(r1)([a,b]):f(r)(x) istnieje p.w. na [a,b] oraz f(r)2(a,b)}.

Oczywiście każda funkcja sklejana rzędu r (niekoniecznie naturalna) należy do Wr(a,b).

Lemat

Niech fWr(a,b) będzie funkcją zerującą się w węzłach, tzn.

f(xj)=0,0jn

Wtedy dla dowolnej naturalnej funkcji sklejanej s𝒮r mamy

abf(r)(x)s(r)(x)dx=0.

Dowód

Dla r=1 funkcja s jest przedziałami stała. Oznaczając przez aj jej wartość na [xj1,xj] dostajemy

abf(x)s(x)dx=j=1ntj1tjf(x)ajdx=j=1naj(f(xj)f(xj1))=0,

ponieważ f zeruje się w tj.

Rozpatrzmy teraz przypadek r2. Ponieważ

(f(r1)s(r))=f(r)s(r)+f(r1)s(r+1),

to

abf(r)(x)s(r)(x)dx=f(r1)(x)s(r)(x)|ababf(r1)(x)s(r+1)(x)dx

Wobec tego, że s jest poza przedziałem (a,b) wielomianem stopnia co najwyżej r1 oraz s(r) jest ciągła na R, mamy s(r)(a)=0=s(r)(b), a stąd

f(r1)(x)s(r)(x)|ab=0.

Postępując podobnie, tzn. całkując przez części r1 razy, otrzymujemy w końcu

abf(r)(x)s(r)(x)dx=(1)r1abf(x)s(2r1)(x)dx

Funkcja s(2r1) jest przedziałami stała, a więc możemy teraz zastosować ten sam argument jak dla r=1, aby pokazać,

że ostatnia całka jest równa zeru.

Funkcje sklejane chcielibyśmy zastosować do interpolacji funkcji. Ważne jest więc, aby odpowiednie zadanie interpolacyjne miało jednoznaczne rozwiązanie.

Twierdzenie O istnieniu i jednoznaczności naturalnego splajnu interpolacyjnego

Jeśli n+1r, to dla dowolnej funkcji f:[a,b]R istnieje dokładnie jedna naturalna funkcja sklejana sf𝒮r interpolująca f w węzłach xj, tzn. taka, że

sf(xj)=f(xj),0jn

Dowód

Pokażemy najpierw, że jedyną naturalną funkcją sklejaną interpolującą dane zerowe jest funkcja zerowa. Rzeczywiście, jeśli s(xj)=0 dla 0jn, to podstawiając w poprzednim lemacie f=s, otrzymujemy

ab(s(r)(x))2dx=0

Stąd s(r) jest funkcją zerową, a więc s jest wielomianem stopnia co najwyżej r1 zerującym się w co najmniej n+1 punktach xj. Wobec tego, że n+1>r1, otrzymujemy s0.

Zauważmy teraz, że problem znalezienia naturalnej funkcji sklejanej s interpolującej f można sprowadzić do rozwiązania układu równań liniowych z macierzą kwadratową. Na każdym przedziale [xi1,xi], 1in, jest ona postaci

s(x)=wi(x)=j=02r1ai,jxj,

dla pewnych współczynników ai,jR, a na (,a] i [b,) mamy odpowiednio

s(x)=w0(x)=j=0r1a0,jxj

i

s(x)=wn+1(x)=j=0r1an+1,jxj.

Aby wyznaczyć s, musimy więc znaleźć ogółem 2r(n+1) współczynników ai,j, przy czym są one związane (2r1)(n+1) warunkami jednorodnymi wynikającymi z gładkości,

wi(k)(xi)=wi+1(k)(xi)

dla 0in i 0k2r2, oraz n+1 niejednorodnymi warunkami interpolacyjnymi,

wi(xi)=f(xi)

dla 0in. Otrzymujemy więc układ 2r(n+1) równań liniowych ze względu na 2r(n+1) niewiadomych ai,j.

Naturalna funkcja sklejana interpolująca f jest wyznaczona jednoznacznie wtedy i tylko wtedy, gdy układ ten ma jednoznaczne rozwiązanie. To zaś zachodzi, gdy zero jest jedynym rozwiązaniem układu jednorodnego. Rzeczywiście, układ jednorodny odpowiada zerowym warunkom interpolacyjnym, przy których, jak pokazaliśmy wcześniej, zerowa funkcja sklejana (której odpowiada ai,j=0, i,j) jest jedynym rozwiązaniem zadania interpolacyjnego.

Naturalnych funkcji sklejanych możemy więc używać do interpolacji funkcji. Pokażemy teraz inną ich własność, która jest powodem dużego praktycznego zainteresowania funkcjami sklejanymi.

Twierdzenie O ekstremalnej własności splajnów naturalnych

Niech fWr(a,b) i niech sf𝒮r będzie naturalną funkcją sklejaną rzędu r interpolującą f w węzłach xj, 0jn. Wtedy

ab(f(r)(x))2dxab(sf(r)(x))2dx.

Dowód

Jeśli przedstawimy f w postaci f=sf+(fsf), to

ab(f(r)(x))2dx=ab(sf(r)(x))2dx+ab((fsf)(r)(x))2dx+2absf(r)(x)(fsf)(r)(x)dx.

Funkcja fsf jest w klasie Wr(a,b) i zeruje się w węzłach xj, 0jn. Z lematu wynika więc, że

absf(r)(x)(fsf)(r)(x)dx=0, a stąd wynika teza.

Wartość całki ab(f(r)(x))2dx może być w ogólności uważana za miarę gładkości funkcji. Dowiedzioną nierówność możemy więc zinterpretować w następujący sposób. Naturalna funkcja sklejana jest w klasie Wr(a,b) najgładszą funkcją spełniającą dane warunki interpolacyjne w wybranych węzłach xj.

Jak już wspomnieliśmy, najczęściej używanymi są kubiczne funkcje sklejane. Dlatego rozpatrzymy je oddzielnie.

Kubiczne funkcje sklejane

Jeśli zdecydowaliśmy się na użycie kubicznych funkcji sklejanych, powstaje problem wyznaczenia sf𝒮2 interpolującej daną funkcję f, tzn. takiej, że sf(xi)=f(xi) dla 0in. W tym celu, na każdym przedziale [xi,xi+1] przedstawimy sf w postaci jej rozwinięcia w szereg Taylora w punkcie xi,

sf(x)=wi(x)=ai+bi(xxi)+ci(xxi)22+di(xxi)36,

i podamy algorytm obliczania ai,bi,ci,di dla 0in1.

Warunki brzegowe i warunki ciągłości dla sf dają nam w0(x0)=0=wn1(xn) oraz wi(xi+1)=wi+1(xi+1), czyli

c0=0,ci+dihi=ci+1,0in2,cn1+dn1hn1=0,

gdzie hi=xi+1xi. Stąd, przyjmując dodatkowo cn=0, otrzymujemy

di=ci+1cihi,1in1

Z warunków ciągłości dla sf dostajemy z kolei

bi+cihi+dihi22=bi+1,0in2,

oraz

bi+1=bi+hici+1+ci2,0in2.

Warunki ciągłości sf dają w końcu

ai+bihi+cihi22+dihi36=ai+1,0in2.

Powyższe równania definiują nam na odcinku [a,b] naturalną kubiczną funkcję sklejaną. Ponieważ poszukiwana funkcja sklejana sf ma interpolować f, mamy dodatkowych n+1 warunków interpolacyjnych wi(xi)=f(xi), 0in1, oraz wn1(xn)=f(xn), z których

ai=f(xi),0in1

Teraz możemy warunki ciągłości przepisać jako

f(xi+1)=f(xi)+bihi+cihi2+dihi36,

przy czym wzór ten zachodzi również dla i=n1. Po wyrugowaniu bi i podstawieniu di, mamy

bi=f(xi,xi+1)+hici+1+2ci6,0in1,

gdzie f(xi,xi+1) jest odpowiednią różnicą dzieloną. Możemy teraz powyższe wyrażenie na bi podstawić, aby otrzymać

cihi6+ci+1hi+hi+13+ci+1hi+16=f(xi+1,xi+2)f(xi,xi+1).

Wprowadzając oznaczenie

ci*=ci6,

możemy to równanie przepisać jako

hihi+hi+1ci*+2ci+1*+hi+1hi+hi+1ci+2*=f(xi,xi+1,xi+2),

0in2, albo w postaci macierzowej

(2w1u22w2u32w3un22wn2un12)(c1*c2*c3*cn2*cn1*)=(v1v2v3vn2vn1),

gdzie

ui=hi1hi1+hi,wi=hihi1+hi,vi=f(xi1,xi,xi+1).

Ostatecznie, aby znaleźć współczynniki ai,bi,ci,di należy najpierw rozwiązać układ równań liniowych, a potem zastosować wzory definiujące pozostałe współczynniki.

Zauważmy, że macierz układu równań liniowych jest trójdiagonalna i ma dominującą przekątną. Układ można więc rozwiązać kosztem proporcjonalnym do wymiaru n używając algorytmu przeganiania. Koszt znalezienia wszystkich współczynników kubicznej funkcji sklejanej interpolującej f jest więc też proporcjonalny do n.

MATLAB i Octave mają wbudowaną funkcję wyznaczającą naturalny kubiczny splajn interpolujący zadane wartości:

s = spline(x,y);

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

Y = ppval(s,X);

Na końcu oszacujemy jeszcze błąd interpolacji naturalnymi kubicznymi funkcjami sklejanymi na przedziale [a,b]. Będziemy zakładać, że f jest dwa razy różniczkowalna w sposób ciągły.

Twierdzenie O błędzie interpolacji splajnem kubicznym

Jeśli fFM1([a,b]) to

fsfC([a,b])5Mmax1in(xixi1)2

W szczególności, dla podziału równomiernego xi=a+iban, 0in, mamy

fsfC([a,b])5M(ban)2.

Dowód

Wykorzystamy obliczoną wcześniej postać interpolującej funkcji sklejanej sf. Dla x[xi,xi+1] mamy

wi(x)=f(xi)+(f(xi+1)f(xi)hihi(ci+1*+2ci*))(xxi)+3ci*(xxi)2+ci+1*ci*hi(xxi)3.

Z rozwinięcia f w szereg Taylora w punkcie xi dostajemy f(x)=f(xi)+f(xi)(xxi)+f(ξ1)(xxi)2/2 oraz (f(xi+1)f(xi))/hi=f(xi)+hif(ξ2)/2. Stąd

f(x)sf(x)=f(x)wi(x)=f(ξ1)2(xxi)2(f(ξ2)2(ci+1*+2ci*))hi(xxi)3ci*(xxi)2ci+1*ci*hi(xxi)3,

oraz

|f(x)sf(x)|(M+2|ci+1*|+6|ci*|)hi2=(M+8max1in1|ci*|)hi2.

Niech teraz max1in1|ci*|=|cs*|. Z postaci układu otrzymujemy

|cs*|=2|cs*||cs*|(us+ws)|uscs1*+2cs*+wscs+1|=|f(xs1,xs,xs+1)||f(ξ3)2|12M,

a stąd

|f(x)sf(x)|5Mhi2,

co kończy dowód.

Przykład

Porównanie interpolacji splajnowej i Lagrange'a.

Interpolacja splajnowa wydaje się lepiej spełniać zadanie odtworzenia kształtu funkcji

Jak widać, w przeciwieństwie do wielomianu interpolacyjnego, splajn interpolacyjny praktycznie pokrywa się z wykresem funkcji, tutaj: f(x)=11+x2.

Uwaga

Niech

WMr(a,b)={fWr(a,b):ab(f(r)(x))2dxM}

Ustalmy węzły a=x0<<xn=b. Dla fWMr(a,b), niech sf będzie naturalną funkcją sklejaną interpolującą f w xj, 0jn, a af dowolną inną aproksymacją korzystającą jedynie z informacji o wartościach f w tych węzłach, tzn.

af=ϕ(f(x0),,f(xn))

Załóżmy, że błąd aproksymacji mierzymy nie w normie Czebyszewa, ale w normie średniokwadratowej zdefiniowanej jako

g2(a,b)=ab(g(x))2dx.

Wtedy

supfWMr(a,b)fsf2(a,b)supfWMr(a,b)faf2(a,b).

Aproksymacja naturalnymi funkcjami sklejanymi jest więc optymalna w klasie WMr(a,b).

Można również pokazać, że interpolacja sf* naturalnymi funkcjami sklejanymi na węzłach równoodległych xj=a+(ba)j/n, 0jn, jest optymalna co do rzędu w klasie WMr(a,b), wśród wszystkich aproksymacji korzystających jedynie z informacji o wartościach funkcji w n+1 dowolnych punktach, oraz

maxfWMr(a,b)fsf*2(a,b)nr.

Uwaga

Tak jak wielomiany, naturalne funkcje sklejane interpolujące dane funkcje można reprezentować przez ich współczynniki w różnych bazach. Do tego celu można na przykład użyć bazy kanonicznej Kj, 0jn zdefiniowanej równościami

Kj(xi)={0ij,1i=j,

przy której sf(x)=j=0nf(xj)Kj(x). Baza kanoniczna jest jednak niewygodna w użyciu, bo funkcje Kj w ogólności nie zerują się na żadnym podprzedziale, a tym samym manipulowanie nimi jest utrudnione.

Częściej używa się bazy B-sklejanej. Można ją zdefiniować dla splajnów dowolnego rzędu za pomocą wzoru rekurencyjnego (przyjmując, że dla dowolnego i, xi<xi+1):

Bi0(x)={1, jeśli xix<xi+1,0, w przeciwnym przypadku;
Bir(x)=xxixi+rxiBir1(x)+xi+r+1xxi+r+1xi+1Bi+1r1(x)

W przypadku naturalnych splajnów kubicznych, r=2, baza B-sklejana jest jawnie zdefiniowana przez następujące warunki:

Bj(xj)=1, dla 0jn,Bj(x)=0, dla xxj2,j2, oraz dla xxj+2,jn2.

Dla B0 i B1 dodatkowo żądamy, aby

B0(x0)=0=B1(x0),B1(x0)=0,

a dla Bn1 i Bn podobnie

Bn1(xn)=0=Bn(xn),Bn1(xn1)=0

Wtedy Bj nie zeruje się tylko na przedziale (xj2,xj+2). Wyznaczenie współczynników rozwinięcia w bazie {Bi}i=0n funkcji sklejanej interpolującej f wymaga rozwiązania układu liniowego z macierzą trójdiagonalną {Bj(xi)}i,j=0n, a więc koszt obliczenia tych współczynników jest proporcjonalny do n.

Uwaga

Oprócz naturalnych funkcji sklejanych często rozpatruje się też okresowe funkcje sklejane. Są to funkcje s~:RR spełniające warunki (i), (ii) \link{splinedef}{definicji funkcji sklejanej}, oraz warunek:

(iii)'
s~(i) jest dla 0ir1

funkcją okresową o okresie (ba), tzn. s~(i)(x)=s~(i)(x+(ba)), x.

Klasę okresowych funkcji sklejanych rzędu r oznaczymy przez 𝒮~r. Funkcje te mają podobne własności jak naturalne funkcje sklejane. Dokładniej, niech

W~r(a,b)={fWr(a,b):f(i)(a)=f(i)(b),0ir1},

tzn. W~r(a,b) jest klasą funkcji z Wr(a,b), które można przedłużyć do funkcji, krórych wszystkie pochodne do rzędu r1 włącznie są (ba)-okresowe na R. Wtedy dla dowolnej funkcji fW~r(a,b) zerującej się w węzłach xj, oraz dla dowolnej s~𝒮~r mamy

abf(r)(x)s~(r)(x)dx=0.

Wynika z niego jednoznaczność rozwiązania zadania interpolacyjnego dla okresowych funkcji f (tzn. takich, że f(a)=f(b)), jak również odpowiednia własność minimalizacyjna okresowych funkcji sklejanych. Dokładniej, jeśli fW~r(a,b) oraz s~f𝒮~r interpoluje f w węzłach xj, 0jn, to

ab(f(r)(x))2dxab(s~f(r)(x))2dx.

Dygresja o najlepszej aproksymacji

Klasyczne zadanie aproksymacyjne w przestrzeniach funkcji definiuje się w następujący sposób.

Niech F będzie pewną przestrzenią liniową funkcji f:[a,b]R, w której określona została norma . Niech VnF będzie podprzestrzenią w F wymiaru n. Dla danej fF, należy znaleźć funkcję vfF taką, że

fvf=minvVnfv.

Okazuje się, że tak postawione zadanie ma rozwiązanie vf, choć nie zawsze jest ono wyznaczone jednoznacznie.

Jako przykład, rozpatrzmy F=Wr(a,b). Utożsamiając funkcje f1,f2Wr(a,b) takie, że f1(x)f2(x)Πr1, zdefiniujemy w Wr(a,b) normę

f=ab(f(r)(x))2dx.

Dla ustalonych węzłów a=x0<<xn=b, niech

Vn+1=𝒮r

będzie podprzestrzenią w Wr(a,b) naturalnych funkcji sklejanych rzędu r opartych węzłach xj, 0jn. Oczywiście dim𝒮r=n+1, co wynika z jednoznaczności rozwiązania w 𝒮r zadania interpolacji. Okazuje się, że wtedy optymalną dla fWr(a,b) jest naturalna funkcja sklejana sf interpolująca f w węzłach xj, tzn.

fsf=mins𝒮rfs.

Rzeczywiście, ponieważ norma w przestrzeni Wr(a,b) generowana jest przez iloczyn skalarny

(f1,f2)=abf1(r)(x)f2(r)(x)dx,

jest to przestrzeń unitarna. Znane twierdzenie mówi, że w przestrzeni unitarnej najbliższą danej f funkcją w dowolnej domkniętej podprzestrzeni V jest rzut prostopadły f na V, albo równoważnie, taka funkcja vfVn+1, że iloczyn skalarny

(fvf,v)=0,vV.

W naszym przypadku, ostatnia równość jest równoważna

ab(fvf)(r)(x)s(r)(x)dx=0,s𝒮r.

To zaś jest prawdą, gdy vf interpoluje f w punktach xj, czyli vf=sf.

Dodajmy jeszcze, że nie zawsze interpolacja daje najlepszą aproksymację w sensie klasycznym.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 6.4 w

  • D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Warto także zapoznać się (nieobowiązkowo) z rozdziałami 6.5 i 6.6 tamże.