1. Strona główna
  2. Matematyka Obliczeniowa - zadania
  3. Kwadratury

Kwadratury

Zadania z tego rozdziału służą przetestowaniu najprostszych kwadratur numerycznych, czyli metod przybliżonego obliczania całek po odcinkach.

W tym rozdziale zapoznamy się także z funkcją octave'a quad() służącą przybliżonemu obliczaniu całek jednowymiarowych postaci

$$<br /> \int_a^b f(t) \: dt \qquad -\infty\leq a<b\leq +\infty,<br /> $$

Możemy więc znajdować przybliżenia całek po całej prostej rzeczywistej czy półprostych.

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

    Policz przy pomocy quad() całkę $ \int_a^b f(t)\:dt $ dla następujących funkcji i odcinków:

    • $ x^2 $ na $ [-1,1] $
    • $ x^9 $ na $ [0,1] $
    • $ \sin(x) $ na $ [0,\pi] $
    • $ \cos(100*x) $ na $ [0,\pi] $
    • $ \cos(100*x)*\cos(1000*x) $ na $ [0,\pi] $

    Czy wyniki są zgodne z teorią?

  2. Złożona kwadratura trapezów.
    1. Zaprogramuj w octavie funkcję

      function c=kwadtrapez(FCN,a,b,n)

      obliczającą złożoną kwadraturę trapezów:

      $$<br /> T_n f=h\left(0.5[f(a)+f(b)]+\sum_{k=1}^{n-1} f(a+k*h)\right)<br /> $$

      dla $ h=(b-a)/n $.

      Parametry funkcji:

      • $ FCN $ - wskaźnik do funkcji octave'a function y=f(x)
        obliczającej wartość funkcji podcałkowej
      • $ a $ - lewy koniec odcinka
      • $ b $ - prawy koniec odcinka
      • $ n $ - ilość obliczeń wartości funkcji podcałkowej w kwadraturze trapezów minus jeden (tak, jak we wzorze powyżej)

      Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą powyższego wzoru, tzn. złożonej kwadratury trapezów.

      Funkcja powinna działać również jeśli ją wywołamy tylko z trzema parametrami. Wtedy $ n $ powinno domyślnie przyjąć wartość sto.

    2. Przetestuj funkcję octave'a z poprzedniego podpunktu dla $ N_k=2^kN_0 $ $ k=0,1,2,\ldots $ z ustalonym $ N_0=5 $ licząc:

      $$<br /> \frac{E_N}{E_{2N}}<br /> $$

      dla błędu

      $$ E_N=|\int_a^b f \:dt - T_N f|$$

      dla następujących funkcji, dla których wartości całek znamy:

      • $ x^2 $ na $ [0,1] $,
      • $ \sin(x) $ na $ [0,\pi] $ i $ N_0=5 $ - funkcja analityczna,
      • $ \sin(100*x) $ na $ [0,\pi] $ i $ N_0=5 $ - funkcja analityczna,
        silnie oscylująca - duże wartości drugiej pochodnej,
      • $ x^{j+0.5} $ na $ [0,1] $ dla $ j=0,1,2 $ czyli funkcji w $ C^j([0,1]) $
        i o nieograniczonej w otoczeniu zera $ j+1 $ pochodnej.

      Porównaj wyniki obliczane kwadraturą trapezów z wynikami funkcji quad().
      Można sprawdzić dla jakiego $ n $ błąd obliczony złożoną kwadraturą trapezów jest na poziomie błędu funkcji octave quad() - i porównać ilość wywołań funkcji $ f $ przez obie procedury.

  3. Czy wielomiany Czebyszewa tworzą układ funkcji ortogonalnych w $ L^2 $ na $ [-1,1] $ z wagą $ \frac{1}{\sqrt{1-x^2}} $.

    Policz za pomocą funkcji octave quad():

    $$\int_{-1}^1\frac{1}{\sqrt{1-x*x}}T_n(x)T_m(x)\:d x$$

    dla $ m=2 $ i $ n=3 $. Czy wynik jest zgodny z teorią?

  4. Kwadratura Gaussa-Czebyszewa dla
    $$I(f)=\int^1_{-1}\frac{1}{\sqrt{1-x^2}}f(x)\:dx$$

    czyli

    $$<br /> GC_{n+1}f=\frac{\pi}{n+1}\sum_{k=0}^n<br /> f(x_k)= \frac{\pi}{n+1}\sum_{k=0}^n<br /> f\left(\cos\left(\frac{0.5*\pi+k\pi}{n+1}\right)\right)<br /> $$

    dla $ x_k $ zer $ T_{n+1} $ - $ n+1 $ wielomianu Czebyszewa.

    Zaimplementuj funkcję function c=gaussczeb(FCN,n)
    obliczającą wartość kwadratury $ GC_n f $. Parametry
    - wskaźnik do funkcji function y=f(x) i ilość punktów kwadratury.

    Przetestuj jej działanie:

    • policz przybliżenia całek $ I(T_k^2) $. Czy $ GC_n (T_k^2) $
      przybliża $ \pi/2 $ dla $ k>0 $ dla $ n=20,40,80 $?

    • policz przybliżenie całki z wielomianów różnych stopni dla $ n=4,10,20 $ -
      czy kwadratura jest dokładna dla wielomianów Czebyszewa stopnia $ 0<k<2n $ (tzn. czy zwraca zero)?

    • dla funkcji matematycznej $ f(x)=\exp(\mathrm{arccos}(x)) $, której całkę z tą wagą możemy obliczyć dokładnie.

      Porównaj wartość kwadratury
      z wynikiem dokładnym
      dla $ n=10,20,40,80 $, tzn. policz

      $$e_n=|GC_n(f)-I(f)|$$

      przyjmując za $ I(f) $ wartość dokładną.

    • analogicznie do poprzedniego podpunktu tego zadania,
      porównaj wynik obliczony za pomocą kwadratury Gaussa - Czebyszewa, tzn. $ GC_n(f) $, dla funkcji $ f(x)=\exp(\mathrm{arccos}(x)) $ z wynikiem obliczonym funkcją quad(). Wypisz na ekran $ e_n=|GC_n(f)-c| $ dla $ c $ wartości uzyskanej za pomocą funkcji octave'a quad().
    • Przetestuj rząd zbieżności tej kwadratury dla funkcji
      $ f(x)=\exp(\mathrm{arccos}(x)) $, jak w przypadku kwadratury trapezów, tzn. policz błędy $ e_n $, analogicznie do poprzednich podpunktów i wypisz $ e_n/e_{2n} $.

  5. Złożona kwadratura prostokątów.
    1. Zaprogramuj w octave funkcję

      function c=kwadprost(FCN,a,b,n),

      która zwraca wartość złożonej kwadratury prostokątów dla funkcji $ f $ na odcinku $ [a,b] $ na $ n $ równoodległych punktach:

      $$<br /> P_n f=h \sum_{k=1}^n f(a+ (k-0.5)*h)), \qquad h=(b-a)/n.<br /> $$

      Parametry funkcji:

      • $ FCN $ - wskaźnik do funkcji octave'a postaci
        function y=f(x)

          y=......

        endfunction

        obliczającej wartość funkcji, której całkę chcemy obliczyć, tj. $ f(x) $, w danym punkcie $ x $.

      • $ a,b $ - końce odcinka,
      • $ n $ - ilość węzłów wykorzystywanych przez kwadraturę.

      Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą powyższego wzoru, tzn. złożonej kwadratury prostokątów.

      Funkcja powinna działać również jeśli ją wywołamy bez ostatniego parametru. Wtedy $ n $ powinno domyślnie przyjąć ustaloną wartość, np. dwieście.

    2. Testy dla funkcji, których całki znamy.

      Przetestuj tę funkcję octave'a kwadprost() dla $ N_k=2^kN_0 $ $ k=0,1,2,\ldots $ z ustalonym $ N_0=5 $ licząc:

      $$<br /> \frac{E_N}{E_{2N}}, \qquad E_N=|\int_a^b f \:dt - P_N f|<br /> $$

      dla całek z następujących funkcji po odpowiednich odcinkach:

      • $ x^2 $ na $ [0,1] $
      • $ \sin(x) $ na $ [0,\pi] $ i $ N_0=5 $ - funkcja analityczna,
      • $ \sin(100*x) $ na $ [0,\pi] $ i $ N_0=5 $ - funkcja analityczna
        silnie oscylująca (duże wartości drugiej pochodnej),
      • $ x^{j+0.5} $ na $ [0,1] $ dla $ j=0,1,2 $ czyli funkcji w $ C^j([0,1]) $
        i o nieograniczonej w otoczeniu zera $ j+1 $ pochodnej.

      Porównaj wyniki obliczane kwadraturą prostokątów z wynikami funkcji quad() .

  6. Porównaj wyniki obliczane kwadraturą prostokątów z wynikami funkcji obliczającej całkę za pomocą złożonej kwadratury trapezów.

    Porównaj błędy dla wartości $ N $ i $ f(x) $ z poprzedniego zadania - czyli całek, których wartość teoretycznie znamy - z wynikami dla obu kwadratur:

    • dla tej samej ilości wywołań funkcji $ f $, tzn. porównaj $ P_N f $ z $ T_{N-1} f $,
    • porównaj obie kwadratury dla tego samego $ h $, tzn.
      $ T_Nf $ z $ P_Nf $, aby ocenić błąd w terminach parametru $ h=(b-a)/N $.
  7. \textrm{Kwadratura Romberga.}
    1. Zaprogramuj w octave funkcję

      function c=Romberg(FCN,a,b,p,N0)

      z kwadraturą Romberga $ R_{p,N0} $ obliczającą przybliżenie
      całki $ \int_a^bf(t) \:dt $.

      Parametry funkcji:

      • $ FCN $ wskaźnik do funkcji octave'a function y=f(x)
        obliczającej wartość funkcji podcałkowej
      • $ a $ lewy koniec odcinka
      • $ b $ prawy koniec odcinka
      • $ p $ ilość poziomów w kwadraturze Romberga
      • $ N0 $ startowa ilość punktów w kwadraturze Romberga

      Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą kwadratury Romberga.

      Funkcja powinna działać również jeśli ją wywołamy tylko z trzema albo czterema parametrami. Parametr $ p $ powinien wtedy domyślnie przyjąć wartość $ 5 $,
      a $ N0 $ wartość sto.

    2. Przetestuj kwadraturę Romberga dla całki $ \int_a^bf\:dt $ dla następujących funkcji, odcinków i wartości parametrów $ N0 $ i $ p $:
      • $ \sin(x) $ dla $ [0,\pi] $
      • $ \sin(x) $ dla $ [0,100*\pi] $
      • $ \sin(10*x) $ dla $ [0,\pi] $
      • $ x^{j+0.5} $ na $ [0,1] $ dla $ j=0,1,2,\ldots,10 $,

      dla $ N0=100 $ i dla $ p=2,3,4,5,6 $.

      Czy dla $ \sqrt{x} $ i rosnącego $ p $ błąd maleje?

      Czy stosunek $ E_N/E_{2N} $ maleje tak samo dla wszystkich funkcji?

      Tutaj użyliśmy oznaczenia na błąd $ E_N=|\int_a^bf\:dt - R_{p,N0}| $.