1. Strona główna
  2. Matematyka Obliczeniowa - zadania
  3. Interpolacja Lagrange'a, bazy wielomianów

Interpolacja Lagrange'a, bazy wielomianów

W tym rozdziale zajmiemy się interpolacją wielomianową.
Zadanie interpolacji wielomianowej polega na znalezieniu wielomianu stopnia nie większego od $ n $, spełniającego $ n+1 $ warunków interpolacyjnych.

W octave'ie istnieje funkcja znajdująca współczynniki wielomianu interpolacyjnego Lagrange'a. Jest to funkcja polyfit(), czyli funkcja obliczająca współczynniki wielomianu interpolacyjnego w bazie potęgowej dla zadanych wartości i węzłów.
Z kolei funkcja polyval() jest funkcją obliczającą wartość wielomianu zadanego poprzez współczynniki w bazie potęgowej w jednym lub równocześnie wielu punktach - czyli wektorowo. Co ważne, obie funkcje są ze sobą zgodne. Trzeba uważać na kolejność współczynników; octave numeruje współczynniki w bazie potęgowej

$$(x^j)_{j=0}^n$$

przestrzeni wielomianów stopnia nie większego od $ n $ w odwrotnej kolejności tzn.

$$(x^n,x^{n-1},\ldots,1)$$

czyli np.
wektor współczynników $ (3,2,1) $ odpowiada wielomianowi $ 3x^2+2x+1 $.
Oczywiście stopnień wielomianu, a dokładniej: dla jakiego $ n $ rozpatrujemy bazę potęgową dla tego wektora, jest długością wektora współczynników pomniejszoną o jeden.

  1. Zapoznaj się z pomocą octave'a do funkcji octave'a polyval().
    Oblicz wartość wielomianu $ x^{50}-1 $ dla $ x=-1,0,1 $.
  2. Korzystając z funkcji polyval() narysuj wykres wielomianu
    $ x^3+x-2 $ bez wykorzystania pętli - czyli wektorowo.

  3. Test funkcji polyfit().

    Zapoznaj się z pomocą octave'a do funkcji octave'a polyfit().

    Wykorzystując funkcję polyfit() znajdź wielomian interpolacyjny $ L_n F $ dla funkcji $ F(x)=\sin(x) $ dla węzłów $ -1,0,1 $.
    Policz wartości różnicy $ F-L_n F $ w węzłach oraz korzystając z funkcji plot() narysuj wykresy $ F $ i $ L_nF $ na jednym rysunku.

    Wykonaj powtórnie to zadanie ale dla węzłów $ -1,0,1,10 $.

    Oblicz wartości wielomianu bez użycia pętli z wykorzystaniem funkcji polyval().

  4. (#)
    Interpolacja Lagrange'a - zbieżność ciągu wielomianów interpolacyjnych dla węzłów równoodległych i węzłów Czebyszewa:

    • Wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne $ L_N f $ dla funkcji $ f=\sin() $ dla $ N $ węzłów równoodległych na $ [0,2*\pi] $ dla $ N=4,8,16,32,64 $.
    • Oblicz dyskretną normę maksimum różnicy $ f-L_N f $ na siatce tysiąca równoodległych punktów na tym odcinku, tzn. $ e_N=\max|f(x_k)-L_N f(x_k)| $, gdzie $ x_k $ to punkty siatki. Policz stosunek $ e_N/e_{2N} $ dla $ N<64 $.

      Czy błędy maleją do zera? Jak zachowuje się $ e_N/e_{2N} $?

      Siatkę tysiąca równoodległych punktów na odcinku $ [a,b] $ najprościej utworzyć z wykorzystaniem funkcji octave'a: \\ linspace(a,b,1000).

    • Narysuj na ekranie wykresy $ \sin(x) $ i tych wielomianów dla różnych $ N $ - używając funkcji polyval() i plot().
  5. Powtórz zadanie  [link] dla tej samej funkcji i tego samego odcinka dla węzłów Czebyszewa.
    Węzły Czebyszewa to pierwiastki wielomianu:

    $$<br /> T_{n+1}(t)=\cos((n+1)\mathrm{arccos}(t))<br /> $$

    na $ [-1,1] $ odpowiednio przesunięte i przeskalowane.

    Przetestuj, czy błędy $ e_N $ dla węzłów Czebyszewa są mniejsze niż dla węzłów równoodległych.

  6. Napisz funkcję znajdującą współczynniki w bazie potęgowej wielomianu interpolacyjnego zadanego stopnia dla węzłów równoodległych oraz węzłów Czebyszewa dla danej funkcji, odcinka
    $ [a,b] $, tzn. napisz funkcję octave'a (w m-pliku):

     function [LN,eN]=LagrInterp(FCN,a,b,N,type=0)

    która dla parametrów:

    • zadanego wskaźnika funkcyjnego FCN (ang. function handle) do funkcji jednego argumentu: function y=f(x),
    • $ a,b $ - końców odcinka $ [a,b] $,
    • $ N $ - stopnia wielomianu interpolacyjnego
    • type - typu węzłów $ 0 $ - równoodległych, $ 1 $ - Czebyszewa

    zwróci wektor $ LN $ - współczynniki $ L_N f $ wielomianu interpolującego funkcję w tych węzłach oraz
    $ eN $ przybliżenie dyskretnej normy maksimum różnicy $ L_Nf - f $ na tym odcinku.

    Przybliżenie normy maksimum liczymy na dyskretnej siatce zawierającej tysiąc punktów.

  7. Powtórz zadanie  [link] dla obu typów węzłów, tzn. powtórz znajdowanie wielomianów interpolacyjnych na węzłach równoodległych i węzłach Czebyszewa dla funkcji
    $$f(x)=\log(1+x)$$

    na odcinkach

    • $ [0,1] $,
    • $ [0,10] $.

    Czy dla tej funkcji i obu odcinków błędy w normach maksimum maleją wraz ze wzrostem $ N $? Porównaj wyniki otrzymane w tym zadaniu w obliczeniach z oszacowaniami teoretycznymi błędu interpolacji Lagrange'a.

  8. Interpolacja Lagrange'a - przykład Rungego. Powtórz zadanie  [link] dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równoodległych i węzłach Czebyszewa, ale dla funkcji:

    $$f(x)=1/(1+x*x)$$

    na $ [-5,5] $.

    Czy obliczone wyniki wskazują na to, że obliczone ciągi wielomianów interpolacyjnych zbiegają do $ f $ jednostajnie dla obu typów węzłów?

  9. Napisz funkcję octave'a obliczającą wartość wielomianu zadanego w bazie
    potęgowej tzn.

    $$w(x)=\sum_{k=0}^n a_kx^k$$

    algorytmem Hornera.
    Parametrami funkcji będą

    • wektor współczynników $ a $
    • macierz wartości $ x $.
    • $ N $ - stopnień wielomianu (to może być parametr opcjonalny, domyślnie przyjmujący wartość równą długości wektora $ a $ minus jeden)

    Funkcja ma zwrócić wartości wielomianu dla wartości w $ x $.

    Przetestuj funkcję dla wielomianów $ 1+x^2 $ oraz $ 1-2x+x^2 $ dla węzłów równoodległych na $ [-1,2] $, tzn. policz wartości wielomianów dla kilku wartości oraz narysuj wykresy tych wielomianów z wykorzystaniem tej funkcji.

  10. Napisz funkcję octave'a znajdującą dla danego wielomianu stopnia $ n $:
    $$w(x)=\sum_{k=0}^n a_kx^k$$

    oraz danej liczby $ q $ współczynniki $ (b_k)_{k=0}^n $ wielomianu

    $$p(x)=\sum_{k=0}^{n-1}b_kx^k$$

    oraz wartość $ r $ takie, że

    $$w(x)=(x-q)*p(x)+r,$$

    obliczone z wykorzystaniem algorytmu Hornera.

  11. Napisz funkcję octave'a znajdującą dla danego wielomianu stopnia $ n $:
    $$w(x)=\sum_{k=0}^n a_kx^k$$

    oraz liczby $ q $ wartość $ w(q) $ i pochodnej $ w'(q) $, obliczone z wykorzystaniem algorytmu Hornera.

  12. Algorytm Hornera w bazie Newtona. Różnice dzielone.

    Zaprogramuj w octavie funkcję ze zmodyfikowanym algorytmem Hornera zwracającą wartość wielomianu zadanego w bazie Newtona dla danych węzłów.
    Parametrami będą

    • $ x $ punkt, w którym obliczamy wielomian (ewentualnie tablica punktów, ale wtedy funkcja też musi zwrócić wektor z wartościami wielomianu w tych punktach),
    • $ N $ - stopnień wielomianu,
    • wektor długości $ N+1 $ ze współczynnikami wielomianu w bazie Newtona.

    Przetestuj na kilku prostych przykładach: dla węzłów $ -1,0,1 $ i wielomian $ w(x)=x^2 $, który w bazie Newtona związanej z tymi węzłami ma następującą postać: $ x^2=(x+1)x - (x+1) +1 $.

  13. Napisz funkcję octave'a, która dla danego wielomianu $ w(x) $, którego współczynniki w bazie
    potęgowej znamy, oblicza współczynniki tego wielomianu w bazie Newtona dla zadanych węzłów podanych w wektorze $ y $. Tzn. parametrami funkcji będą:

    • wektor $ a=(a_k)_k $ taki, że

      $$w(x)=\sum_{k=0}^n a_kx^k$$
    • $ y=(y_k)_k $ wektor współczynników bazy Newtona.

    Funkcja powinna zwrócić wektor współczynników $ b_k $ takich, że

    $$<br /> w(x)=\sum_{k=0}^n b_k w_k,<br /> $$

    gdzie

    $$ w_k(x)=\Pi_{j=0}^{k-1}(x-y_j).<br /> $$
  14. Napisz funkcję, która dla danego wielomianu $ w(x) $, którego współczynniki w bazie
    Newtona $ (w_k)_{k=0}^n $ znamy, oblicza współczynniki tego wielomianu w bazie potęgowej $ (1,x,x^2,\ldots,x^n) $. Tzn. parametrami funkcji jest wektor współczynników $ b=(b_k)_k $ takich, że

    $$w(x)=\sum_{k=0}^n b_k w_k$$

    i wektor $ y=(y_k)_k $ węzłów bazy Newtona $ (w_k)_{k=0}^n $
    dla

    $$ w_k(x)=\Pi_{j=0}^{k-1}(x-y_j).<br /> $$

    Funkcja ma zwrócić wektor współczynników $ a=(a_k)_k $ takich, że

    $$<br /> w(x)=\sum_{k=0}^n a_k x^k.<br /> $$
  15. Sprawdź eksperymentalnie ile wynosi dla różnych wartości $ N=4,8,16,32,64,\ldots $ przybliżenie:
    $$<br /> A_{r,N}=\sum_{k=0}^N\|l_k\|_{\infty,[-1,1]}<br /> $$

    dla $ \{l_k\}_{k=0}^N $ bazy Lagrange'a dla węzłów równoodległych na
    $ [-1,1] $, tzn. dla $ x_k=-1+k*h $ dla $ h=2/N $.

    Normę maksimum funkcji $ \|f\|_{\infty,[a,b]} $ liczymy w sposób przybliżony obliczając dyskretną normę maksimum na $ \max\{1000,100*N\} $ równoodległych punktach z odcinka $ [a,b] $.

  16. Sprawdź eksperymentalnie ile wynosi dla różnych $ N $ np.
    $$N=4,8,16,32,64,\ldots$$

    przybliżenie:

    $$<br /> A_{c,N}=\sum_{k=0}^N\|l_k\|_{\infty,[-1,1]}<br /> $$

    dla $ \{l_k\}_{k=0}^N $ bazy Lagrange'a dla węzłów Czebyszewa na
    $ [-1,1] $, tzn. dla $ x_k $ zer wielomianu $ T_{N+1}(x)=\cos((N+1)\mathrm{arccos}(x)) $.

    Normę maksimum funkcji $ \|f\|_{\infty,[a,b]} $ liczymy w sposób przybliżony obliczając dyskretną normę maksimum na $ \max\{1000,100*N\} $ równoodległych punktach z odcinka $ [a,b] $.

  17. Powtórz dwa poprzednie zadania dla węzłów równoodległych i węzłów Czebyszewa na odcinku $ [0,10] $ i dla odpowiedniej normy maksimum na tym odcinku.
  18. Policz iloraz
    $$<br /> \frac{\|L_N f\|_{\infty,[-5,5]}}{\sum_{k=0}^N\|l_k\|_{\infty,[-5,5]}}<br /> $$

    dla $ f=1/(1+x^2) $ i $ L_N f $ wielomianu interpolującego $ f $ w węzłach równoodległych na $ [-5,5] $ dla $ N=10,20,40,80 $. Tutaj $ \{l_k\}_{k=0}^N $ baza Lagrange'a dla tych węzłów.

    Normę maksimum funkcji $ \|f\|_{\infty,[a,b]} $ liczymy w sposób przybliżony obliczając dyskretną normę maksimum na $ \max\{1000,100*N\} $ równoodległych punktach z odcinka $ [a,b] $.

  19. Powtórz poprzednie zadanie dla tych samych funkcji i odcinka dla węzłów Czebyszewa zamiast węzłów równoodległych.