Splajny kubiczne i liniowe. Interpolacja splajnowa

W tym rozdziale zajmiemy się interpolacją splajnową, czyli interpolowaniem danej funkcji za pomocą splajnów - inaczej funkcji giętych.

Skupimy się na splajnach kubicznych, czyli funkcjach, które są klasy $ C^2 $ na odcinku $ [a,b] $ i dla danego podziału tego odcinka:

$$<br /> a=x_0<x_1\ldots<x_N=b<br /> $$

na pododcinki. Te funkcje obcięte do każdego pododcinka $ [x_i,x_{i+1}] $ są wielomianami kubicznymi.

Zadanie interpolacji splajnami kubicznymi polega na znalezieniu spaljnu kubicznego $ s $ spełniającego:
\begin{eqnarray*}
s(x_0)=y_0\\
s(x_1)=y_1 \\
\vdots \\
s(x_N)=y_N
\end{eqnarray*}
dla zadanych wartości $ y_k $. Okazuje się, że tak postawione zadanie nie jest jednoznaczne; trzeba dodać dwa dodatkowe warunki na $ s $. Zazwyczaj są to odpowiednie warunki brzegowe, tzn. związane z wartościami $ s $, pierwszych lub drugich pochodnych $ s $ w końcach odcinka.

  1. Funkcje octave'a spline() and ppval.

    Zapoznaj się z pomocą do tych funkcji (help spline i help ppval).

    Wykorzystując te funkcje narysuj wykres splajnu kubicznego $ s_1 $ na podziale równomiernym
    odcinka $ [-3,3] $ z węzłami $ \{x_k=k\} $ dla $ k=-3,-2,\ldots,3 $ przyjmującego wartości $ s_1(x_k)=(-1)^k $ w tych węzłach.

    Następnie znajdź współczynniki splajnu kubicznego $ s_2 $ na tym samym podziale odcinka i przyjmującego te same wartości w węzłach co $ s_1 $, ale który dodatkowo przyjmuje wartości pochodnych
    w końcowych węzłach równe zero, tzn. wywołaj
    funkcje spline() podając dwie wartości więcej.

    Następnie narysuj wykresy splajnów $ s_1 $ i $ s_2 $ na tym samym rysunku.

    Czy otrzymaliśmy te same splajny?

    Policz przybliżoną normę maksimum
    różnicy $ s_1-s_2 $ na odcinku $ [-3,3] $.

  2. \textrm{Splajn kubiczny bazowy.}

    Dla danych węzłów równoodległych $ \{k\}_{k=-5,-4,\ldots,5} $ na $ [-5,5] $ narysuj wykres splajnu kubicznego typu not-a-knot (czyli splajnu, którego współczynniki zwróci funkcja spline() przy najprostszym wywołaniu przez podanie wektora węzłów i wektora wartości w tych węzłach, por. help spline)
    takiego, że $ s(0)=1 $ i $ s(k)=0 $ dla węzłów $ k\not=0 $.

    Określ na podstawie wykresu nośnik tego splajnu.

  3. (#)
    \textrm{Splajn kubiczny o minimalnym nośniku.}

    Dla danych węzłów równoodległych $ \{k\}_{k=-5,-4,\ldots,5} $ na $ [-5,5] $ narysuj wykres splajnu kubicznego
    takiego, że $ s(-1)=s(1)=1,s(0)=4 $ i $ s(k)=0 $ dla węzłów $ k\not \in\{-1,0,1\} $ oraz ma pochodne równe zero w węzłach skrajnych, tzn. : $ -5 $ i $ 5 $. Czy poza $ [-2,2] $ ten splajn jest równy zero?

    Policz przybliżone normy maksimum na $ [-5,-2] $ i $ [2,5] $ dla tego splajnu.

  4. (#) Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego z hermitowskimi warunkami brzegowymi.

    Korzystając z funkcji octave'a spline() znajdź współczynniki interpolacyjnego splajnu kubicznego hermitowskiego $ S_N $ na $ N $ węzłach równoodległych dla funkcji $ f(x)=\sin(x) $ na odcinku $ [-\pi,2*\pi] $ dla $ N=2^kN_0 $ dla $ N_0=5 $ i $ k=1,2,3,4,5 $.

    Następnie

    • narysuj wykresy funkcji $ f(x) $ i splajnów $ S_N $ dla różnych $ N $.
    • oblicz dyskretną normę maksimum na siatce równomiernej złożonej z tysiąca punktów na tym odcinku, tzn. $ e_N=\max_k|\sin(x_k)-S_N (x_k)| $ dla $ x_k $ punktów siatki.
    • policz równocześnie
      współczynnik $ \frac{e_{N}}{e_{2N}} $. Czy prawdą jest, że

      $$<br /> \frac{e_{N}}{e_{2N}}\approx 2^p<br /> $$

      dla jakiegoś $ p $ całkowitego np. $ p=4,8 $ lub $ 16 $?

  5. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego bez warunków brzegowych (splajn typu not a knot).

    Powtórz zadanie   [link], ale dla splajnów interpolacyjnych otrzymanych przez spline() bez podawania wartości pochodnych w skrajnych węzłach. Czy współczynniki $ \frac{e_{N}}{e_{2N}} $ są te same? Tzn. czy szybkość zbieżności $ \|\sin(x)-S_N\|_\infty $ jest taka sama?

  6. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego naturalnego (warunek brzegowy - zerowanie się drugich pochodnych w końcach odcinka).

    Powtórz zadanie  [link],
    ale dla splajnów interpolacyjnych naturalnych. Tu trzeba wykorzystać funkcję z octave-forge (czyli rozszerzenia pakietu octave)\\ pp=csape(x,y,'variational') \\- ostatni argument określa to, że splajn będzie naturalny.

    Podajemy link do strony www z pomocą do funkcji csape(): http://octave.sourceforge.net/splines/function/csape.html

  7. \textrm{Przykład Rungego,
    czyli $ f(x)=1/(1+x*x) $ i odcinek $ [-5,5] $, a zbieżność interpolacji splajnami kubicznymi.}

    Przetestuj jak w poprzednich zadaniach, czy splajny interpolacyjne kubiczne z podanymi warunkami na pochodne w końcach odcinka zbiegają w normie supremum do $ f $, tzn. korzystając z funkcji octave'a spline() znajdź współczynniki splajnu interpolacyjnego kubicznego $ S_N $ na $ N $ węzłach równoodległych dla $ f $ na odcinku $ [-5,5] $ dla $ N=2^kN_0 $ dla $ N_0=5 $ i $ k=1,2,3,4,5 $ oraz
    narysuj wykresy $ f $ i tych splajnów dla różnych $ N $.

    Następnie
    oblicz dyskretną normę maksimum na siatce złożonej z tysiąca punktów na tym odcinku, tzn. $ e_N=\max|\sin(x_k)-S_N \sin(x_k)| $ dla $ x_k=-5+k*0.01 $ z $ k=0,\ldots,1000. $

    Policz równocześnie
    współczynnik $ \frac{e_{N}}{e_{2N}} $. Czy
    $ \frac{e_{N}}{e_{2N}}\approx 2^p<br /> $dla jakiegoś $ p $ całkowitego?

  8. Funkcja octave'a mkpp() .
    Zapoznaj się z tą funkcją (help mkpp()). Utwórz przy pomocy mkpp() splajn kubiczny $ s $ na podziale $ -1\leq0\leq1 $ odcinka $ [-1,1] $
    taki, że $ s $ jest wielomianem trzeciego stopnia na całym odcinku $ [-1,1] $ np. $ s(x)=x^2 $ lub $ (x+1)^3 $.

  9. Funkcja octave'a umkpp().

    Zapoznaj się z tą funkcją (help umkpp()).

    Utwórz przy pomocy spline() splajn kubiczny $ s $ na podziale $ -1\leq0\leq1 $ odcinka $ [-1,1] $
    taki, że $ s $ interpoluje wielomian trzeciego stopnia na całym odcinku $ [-1,1] $ np. $ f(x)=(x+1)^3 $.

    Następnie sprawdź współczynniki $ s $ w bazie
    $ \{(x-x_k)^j\}_{j=0,1,2,3} $ za pomocą umkpp() na obu przedziałach, tzn. dla $ x_0=-1 $ na przedziale $ [-1,0] $ i $ x_1=0 $ na $ [0,1] $.

  10. Znajdź za pomocą funkcji octave'a umkpp() współczynniki splajnu z zadania  [link] na wszystkich pododcinkach $ [k,k+1] $ w bazie $ \{(x-k)^j\}_{j=0,\ldots,3} $ dla
    $ k=-2,\ldots,1 $ czyli tam, gdzie ten B-splajn ma nośnik.

    Porównaj z wynikami otrzymanymi teoretycznie.

  11. \textrm{Interpolacja splajnami liniowymi.}

    Dla danego równomiernego podziału odcinka $ [-\pi,2*\pi] $ na $ N $ pododcinków utwórz za pomocą mkpp() strukturę splajnu liniowego $ s_n $ interpolującego funkcję $ \sin(x) $ w węzłach dla $ N=3,6,9,18 $.

    • Narysuj
      wykresy funkcji $ \sin(x) $ oraz tych splajnów liniowych na jednym wykresie
    • Policz przybliżone normy maksimum (na 1000 punktach z tego odcinka) błędu $ e_N=\|\sin -s_n\|_\infty $.
    • Policz
      współczynnik $ \frac{e_N}{e_{2N}} $. Czy widać, że

      $$<br /> \frac{e_N}{e_{2N}}\approx 2^p<br /> $$

      dla jakiegoś $ p $ całkowitego np. $ p=2,4 $ lub $ 8 $?

    W tym zadaniu można wykorzystać funkcję z następnego zadania, tj. zadania  [link].

  12. (#)
    Napisz w octavie funkcję function pp=linspline(x,y), która dla

    • $ x $ wektora $ N+1 $ różnych węzłów
      uszeregowanych ($ a=x_0<x_1\ldots<x_N=b $)
    • $ y $ - wektora $ N+1 $ wartości
      funkcji $ y=f(x) $

    zwróci w strukturze $ pp $ współczynniki
    splajnu liniowego $ s $ dla podziału zadanego węzłami $ x_k $.

    Strukturę należy utworzyć funkcją octave'a mkpp() w taki sposób, aby można było obliczyć wartość
    tego splajnu w punkcie (tablicy punktów) za pomocą funkcji octave'a ppval().

  13. Testowanie rzędu zbieżności interpolacji splajnami liniowymi w zależności od gładkości funkcji.

    Dla funkcji

    $$<br /> f_j(x)=(x)_+^j=\left\{<br /> \begin{array}{lcr}<br /> 0 &\quad &x<0\\<br /> x^j& \quad& x>0<br /> \end{array}<br /> \right. \qquad j=1,2,3<br /> $$

    oraz
    dla podziału odcinka $ [a,b] $ z $ a=-\pi $ i $ b=3 $ na węzły równoodległe

    $$\{x_k=-\pi+k*h\}$$

    dla $ h=(b-a)/N $ i $ N=4,8,16,32,64 $.

    Przetestuj rząd zbieżności splajnu liniowego interpolacyjnego.

    Bardziej szczegółowo:

    • przy pomocy funkcji octave'a z zadania  [link]
      znajdź współczynniki odpowiedniego splajnu liniowego $ s_Nf_j $.
    • policz przybliżoną normę maksimum błędu (na $ 1000 $ równomiernych punktach), tzn. przybliżenie
      $$e_N=\|s_Nf_j -f_j\|_{\infty,[a,b]}.$$
    • policz
      współczynniki $ \frac{e_{N}}{e_{2N}} $. Czy widać, że

      $$<br /> \frac{e_{N}}{e_{2N}}\approx 2^p<br /> $$

      dla jakiegoś $ p $ całkowitego np. $ p=2,4 $ lub $ 8 $?
      Czy widać różnicę dla różnych $ j $?