1 Wstęp do metod numerycznychWykład 9 Całkowanie numeryczne dr inż. Wojciech Bieniecki Instytut Matematyki i Informatyki

2 Całkowanie funkcji Całkowanie jest działaniem odwrotnym do różniczkowania. Żeby sprawnie liczyć całki, należy wcześniej dobrze opanować liczenie pochodnych. Całkę oznaczamy symbolem który przypomina rozciągniętą literę S. Symbol ten pochodzi od łacińskiego słowa Summa (suma). Poniżej uproszczona definicja całki. Całką funkcji f(x) nazywamy taką funkcję F(x), że: Funkcję F(x) spełniającą powyższy warunek nazywa się funkcją pierwotną. Samą operację całkowania zapiszemy w następujący sposób: 2

3 Przykład całkowania Oblicz całkę funkcji f(x) = 2x + 7.Wykonujemy następujący rachunek: Sprawdzamy rozwiązanie Należy jednak zauważyć, że znaleziona przez nas funkcja F(x) = x2 + 7x nie jest jedynym dobrym rozwiązaniem. Do powyższej funkcji moglibyśmy dodać jeszcze dowolną liczbę i wówczas otrzymalibyśmy inną dobrą funkcję pierwotną dla funkcji f(x) = 2x + 7. Na przykład: Zatem powinniśmy napisać: gdzie C - to dowolna stała liczba rzeczywista (constans). 3

4 Definicje Całką nieoznaczoną funkcji f na przedziale (a, b) nazywamy zbiór wszystkich funkcji pierwotnych funkcji f i oznaczamy symbolem Symbol ten czytamy następująco: "całka f od x (po) dx" Funkcję f nazywamy funkcją podcałkową Zmienną x nazywamy zmienną całkowania Wyrażenie dx nazywamy różniczką zmiennej całkowania . Jeżeli F jest dowolną, ustaloną funkcją pierwotną funkcji f, to przyjmuje się następującą konwencję zapisu: stałą C nazywamy stałą całkowania. Funkcję, dla której na danym przedziale istnieje całka nieoznaczona, nazywamy całkowalną w sensie Newtona ( lub krócej całkowalną ) na tym przedziale.

5 Wzory całkowe Do całkowania prostych funkcji wykorzystujemy wzory całkowe, które są również przydatne przy liczeniu całek bardziej skomplikowanych funkcji. 5

6 Całka oznaczona W interpretacji geometrycznej, obliczenie całki oznaczonej w przedziale [a,b] to obliczenie pola pod wykresem funkcji nad osią OX oraz pod osią OX ze znakiem ujemnym, dokładnie pod fragmentem wykresu dla x zawartych w przedziale ⟨a,b⟩.

7 Całka oznaczona Zgodnie z definicją, liczenie całki oznaczonej funkcji ƒ(x) polega na podzieleniu przedziału ⟨a,b⟩ na fragmenty ⟨a,x1⟩, ⟨x1,x2⟩,…,⟨xn,b⟩ a następnie utworzeniu prostokątów o polach P1=ƒ(a)⋅(x1−a), P2=ƒ(x1)⋅(x2−x1), …, Pn+1=ƒ(xn)⋅(b−xn). Suma pól prostokątów jest przybliżeniem pola pod wykresem funkcji. Im więcej prostokątów o coraz krótszych bokach , tym przybliżona wartość pola jest bliższa rzeczywistej jego wartości. Jeżeli przyjmiemy Δxi=xi−xi-1 to dokładną wartość pola wyraża wzór: Ilość wyrazów dąży do nieskończoności: i→∞ Grubość prostokątów dąży do zera : Δxi→0

8 Całka oznaczona - definicjaDla danej funkcji ƒ(x) określonej w przedziale ⟨a,b⟩, który dzielimy w dowolny sposób na części wstawiając punkty podziału: a największą z różnic Δxi=xi+1−xi (i=0,1,2…n−1) oznaczając przez δ Całką oznaczoną funkcji ƒ(x) w przedziale ⟨a,b⟩ nazywamy granicę: gdzie ξi jest dowolną liczbą w pomiędzy xi+1 a xi (xi≤ξi≤xi+1). Jest to definicja Riemanna Całka oznaczona przyporządkowuje funkcji pewną liczbę. Dla funkcji ciągłej i określonej w przedziale ⟨a,b⟩ zachodzi wzór: gdzie F(x) jest funkcją pierwotną funkcji ƒ(x). Aby skorzystać z tego wzoru należy obliczyć całkę nieoznaczoną. Jest to definicja Newtona-Leibniza

9 Przykład – całka oznaczonaPrzykład 1: Oblicz pole pod parabolą y=x2 w przedziale ⟨0;1⟩. Obliczamy całkę nieoznaczoną Całka x2 = 1/3 x3 Przykład 2: Oblicz pole pod jednym garbem sinusoidy y=sinx – przedział ⟨0;Π⟩ UWAGA: Z definicji, funkcja musi być określona w przedziale całkowania, dlatego niedopuszczalne jest liczenie całki w poniższy sposób W przedziale całkowania (−1; 1) znajduje się x=0, dla którego funkcja ta jest nieokreślona. Wynik jest absurdem, zwłaszcza, że funkcja nigdy nie przyjmuje wartości ujemnych.

10 Zadanie całkowania numerycznegoNa podstawie znajomości wartości yi funkcji f(x) w punktach xi (narzuconych lub wybieranych) wyznaczyć wartość I całki oznaczonej z funkcji w przedziale : Zastosowania całek • Splot • Wartość średnia, skuteczna • Energia sygnału • Pola powierzchni, objętości • Prawdopodobieństwo na podstawie funkcji gęstości prawdopodobieństwa Wartość przybliżona: odpowiadające im operacje sumacyjne na sygnałach dyskretnych

11 Metody całkowania numerycznegoMetoda prostokątów Metoda trapezów Metody interpolacji wielomianowej Newtona-Cotesa Gaussa-Legendre’a Metoda Simpsona Metoda Monte-Carlo

12 Kwadratura ϕ(x) = f0 N0(x) + f1 N1(x) + · · · + fn Nn(x)Gdy przedział całkowania jest skończony, wówczas numeryczne całkowanie polega na zastąpieniu funkcji podcałkowej f (x) odpowiednim wielomianem interpolacyjnym lub aproksymacyjnym ϕ(x) zbudowanym na zbiorze n + 1 węzłów o współrzędnych xi , i = 0, 1, 2, , n . Wymaga to wówczas całkowania jedynie prostych funkcji bazowych z wykorzystaniem wzoru na I(f ). Metody całkowania numerycznego wykorzystują interpolację lub aproksymację funkcji za pomocą wielomianów algebraicznych. Podstawiając w miejsce funkcji podcałkowej f (x) wielomian algebraiczny ϕ(x) = f0 N0(x) + f1 N1(x) + · · · + fn Nn(x) otrzymamy tzw. wzór kwadraturowy.

13 Kwadratura Wzorem kwadraturowym albo krócej kwadraturą nazywamy:W którym są tzw. współczynnikami wagowymi (wagami). Wartość wi określa wielkość udziału rzędnej fi ≡ f (xi) w wartości całej sumy S(f) . Dokładność kwadratury S(f ) jest tym większa, im mniejsza jest różnica I(f) − S(f) nazywana błędem kwadratury.

14 Kwadratury dla węzłów równoodległychNajprostszymi kwadraturami interpolacyjnymi są kwadratury zbudowane na węzłach równoodległych, o danych współrzędnych xi = x0 + i · h, i = 0, 1, 2, , n . Niewiadome współczynniki wi są obliczane z układu n + 1 liniowych równań algebraicznych, które otrzymamy na podstawie kwadratury zastosowanej dla wielomianów Wk(x) = xk , k = 0, 1, 2, , n, dla których I(Wk ) = S(Wk). skąd

15 Kwadratury interpolacyjne – układ równańRozwiązaniem tego układu algebraicznych równań liniowych są wartości wag wi .

16 Metoda prostokątów Skorzystamy z definicji całki oznaczonej Riemanna , w której wartość całki interpretowana jest jako suma pól obszarów pod wykresem krzywej w zadanym przedziale całkowania . Sumę tę przybliżamy przy pomocy sumy pól odpowiednio dobranych prostokątów. Sposób postępowania jest następujący: Przedział całkowania  dzielimy na n równo odległych punktów x1,x2,...,xn. Punkty te wyznaczamy w prosty sposób wg wzoru: odległość między dwoma sąsiednimi punktami, będzie podstawą każdego prostokąta:

17 Metoda prostokątów Dla każdego wyznaczonego w ten sposób punktu obliczamy wartość funkcji f(x) w tym punkcie: Obliczamy sumę iloczynów wyznaczonych wartości funkcji przez odległość dx między dwoma sąsiednimi punktami - da to sumę pól poszczególnych prostokątów ograniczonych wykresem funkcji: S = f1 dx + f2 dx fn dx lub S = dx (f1 + f2  fn) Otrzymana suma jest przybliżoną wartością całki oznaczonej funkcji f(x) w przedziale 

18 Wzór prostokątów - kwadraturaWe wzorze kwadraturowym podstawmy ϕ(x) = f (x0) = const. Otrzymujemy przy czym w0 = b − a .

19 Przykład Obliczymy ręcznie przybliżoną wartość całki oznaczonej z funkcji f(x) = sin(x) w przedziale <0,π>. Przedział podzielimy na n+1 = 5  punktów: Odległość między dwoma sąsiednimi punktami wynosi: Dla każdego z wyznaczonych punktów obliczamy wartość funkcji f(x) = sin(x): Obliczamy sumę pól prostokątów: S = dx (f1 + f2 + f3 + f4) S = 0,7854 (0, , , ,0000) S = 0,7854 × 2,4142 S = 1,8961 Dokładna wartość takiej całki oznaczonej wynosi 2 – popełniliśmy błąd około 5%

20 Algorytm Po odczytaniu informacji o krańcach xp i xk przedziału całkowania ustawiamy sumę s na 0, obliczamy odległość dx pomiędzy sąsiednimi punktami podziałowymi i ustawiamy ich licznik na 1. Rozpoczynamy pętlę iteracyjną, która wykona się n-razy. Wewnątrz pętli obliczamy i-ty punkt podziałowy oraz wartość funkcji w tym punkcie, którą dodajemy do sumy s. Po zakończeniu pętli sumę s musimy jeszcze pomnożyć przez szerokość podstawy prostokątów, czyli odległość dx. Po tej operacji s zawiera wartość przybliżoną całki. Zwracamy ją użytkownikowi i kończymy algorytm.

21 Wybór położenia węzła W zależności od wyboru położenia węzła x0 otrzymujemy wzory: (a) lewych prostokątów, gdy x0 = a

22 Wybór położenia węzła (b) środkowych prostokątów, gdy x0 = (a + b)/2

23 Wybór położenia węzła (c) prawych prostokątów, gdy x0 = b

24 Metoda trapezów Metoda prostokątów nie jest zbyt dokładna, ponieważ pola użytych w niej prostokątów źle odwzorowują pole pod krzywą. Dużo lepszym rozwiązaniem jest zastosowanie zamiast nich trapezów o wysokości dx i podstawach równych odpowiednio wartości funkcji w punktach krańcowych. Sama zasada nie zmienia się. Przedział całkowania  dzielimy na n+1 równo odległych punktów x0, x1, x2, ... , xn. Punkty te wyznaczamy w prosty sposób wg wzoru:

25 Metoda trapezów Obliczamy odległość między dwoma sąsiednimi punktami - będzie to wysokość każdego trapezu: Dla każdego wyznaczonego w ten sposób punktu obliczamy wartość funkcji f(x) w tym punkcie: Pole pod wykresem funkcji przybliżane jest polami n trapezów. Pole i-tego trapezu obliczamy wg wzoru: Przybliżona wartość całki jest sumą pól wszystkich otrzymanych w ten sposób trapezów:

26 Wzór trapezów - kwadraturaJeśli do interpolacji funkcji f (x) zastosujemy interpolację za pomocą wielomianu liniowego zbudowanego na bazie Lagrange’a, to otrzymamy wzór kwadraturowy, nazywany wzorem trapezów. Wagi wi , i = 0, 1 występujące we wzorze trapezów można wyznaczyć rozwiązując układ równań:

27 Przykład Porównamy wynik z metodą prostokątów dla tego samego zadania, czyli Obliczamy sumę pól trapezów: Otrzymaliśmy wynik identyczny jak w przypadku metody prostokątów  i zbyt mały – powinno wyjść 2.

28 Algorytm Po odczytaniu informacji o krańcach xp i xkprzedziału całkowania ustawiamy sumę s na 0, obliczamy odległość dx pomiędzy sąsiednimi punktami podziałowymi i ustawiamy ich licznik na 1. Rozpoczynamy pętlę iteracyjną, która wykona się (n-1)-razy. Wewnątrz pętli obliczamy i-ty punkt podziałowy oraz wartość funkcji w tym punkcie, którą dodajemy do sumy s. Po zakończeniu pętli do sumy s dodajemy średnią wartość funkcji na krańcach przedziału i całość przemnażamy przez dx . Po tej operacji s zawiera wartość przybliżoną całki. Zwracamy ją użytkownikowi i kończymy algorytm.

29 Metoda Simpsona W metodzie Simpsona fragment funkcji pomiędzy węzłami xi i xi+1 będziemy przybliżali parabolą (w metodzie trapezów była to linia prosta przechodząca przez węzły, a w metodzie prostokątów – linia pozioma). Dla paraboli wartość całki można łatwo policzyć analitycznie   Przedział całkowania  dzielimy na n równo odległych punktów x1,x2,...,xn. Punkty te wyznaczamy wg wzoru: Dla każdych dwóch sąsiednich punktów wyznaczamy punkt środkowy ti wg wzoru: odległość między dwoma sąsiednimi punktami:

30 Metoda Simpsona Dla węzłów i punktów środkowych wyznaczamy wartości funkcji f: W każdym podprzedziale  przybliżamy funkcję za pomocą paraboli g(x) o następującej postaci: dla i = 1,2,...,n gi(x) = aix2 + bix + ci,  x Pole pod parabolą w przedziale  będzie równe całce oznaczonej:

31 Metoda Simpsona Przekształcając to wyrażenie możemy pozbyć się współczynników a,b,c, a więc nie trzeba ich będzie wyliczać, jak to było w przypadku różniczkowania Wartość całki dla przedziału [xp,xk] jest przybliżona sumą całek po parabolach

32 Kwadratura dla metody SimpsonaZastosowanie kwadratowej interpolacji Lagrange’a prowadzi do wzoru kwadraturowego: Po przekształceniach Dla 3 węzłów równoodległych, tzn. gdy c = 0.5 · (a + b) współczynniki wagowe wi , i = 0, 1, 2 we wzorze Simpsona oblicza się z układu równań Co daje

33 Przykład Porównamy wynik z wcześniejszymi metodami dla tego samego zadania, czyli Obliczamy sumę pól pod parabolami: S = (x2-x0)/(6∙2)(f0 + f2 +2 f1 + 4(ft1+ft2) S = 3,1415/12 ( ∙ 0,7071) S = 0,2618 × (2+5,6569) S = 2,004

34 Algorytm Odczytujemy krańce przedziału całkowania . Do obliczenia całki musimy zliczyć dwie sumy – s dla wartości funkcji w punktach podziałowych xi oraz wst dla wartości funkcji w punktach środkowych przedziałów ti. Wyznaczamy odległość pomiędzy dwoma sąsiednimi punktami podziałowymi dx i rozpoczynamy pętlę, w której zmienna i pełni rolę numeru punktu podziałowego i środkowego. Pętla ta wykonuje się n-razy od i=1 do i=n włącznie. W pętli wyznaczamy xi i umieszczamy ją w zmiennej x. Następnie obliczamy wartość funkcji w  ti, który jest odległy o połowę dx od wyznaczonego wcześniej punktu xi. Wartość tę dodajemy do sumy st. Drugą sumę tworzymy w zmiennej s. Jednakże powinna ona zawierać jedynie wartości funkcji dla punktów podziałowych od x1 do xn-1.

35 Wzory Newtona-Cotesa Zastosowanie wielomianów ϕ(x) coraz to wyższych stopni we wzorze kwadraturowym prowadzi do tzw. wzorów Newtona-Cotesa. Dla wielomianów ϕ(x) kolejnych stopni n wartości współczynników wagowych wi otrzymuje się z rozwiązania układu równań. W poniższej tablicy są zestawione wartości współczynników wagowych wi : Wartości wag występujące we wzorze N-C obliczane są według wzoru

36 Metoda Monte-Carlo Metoda Crude Monte CarloMonte Carlo jest jedną z najczęściej używanych metod symulacyjnych. Jedną z najczęstszych aplikacji metody Monte Carlo jest całkowanie numeryczne Metoda Crude Monte Carlo Chcemy policzyć całkę z funkcji f(x) W przedziale [a, b] losujemy N punktów sortujemy ich wartości i przybliżamy Wygląda to tak samo jak całka Riemanna.

37 Algorytm Obliczenia rozpoczynamy od pobrania informacji o przedziale całkowania oraz o ilości punktów losowych, które należy wygenerować w celu obliczenia wartości średniej funkcji w tym przedziale. Dokładność metody rośnie wraz ze wzrostem n. W zmiennej s będziemy obliczać sumy wartości funkcji. Zmienna ta posłuży później do wyliczenia średniej oraz samej całki oznaczonej. Na początku obliczeń ustawiamy ją na 0. W zmiennej dx zapamiętujemy szerokość przedziału całkowania. Wartość ta jest wykorzystywana zwykle przy generacji liczby losowej oraz na końcu przy obliczeniu wartości całki. Rozpoczynamy pętlę iteracyjną kontrolowaną przez zmienną i. Pętla ta wykona się n-razy. Wewnątrz pętli generujemy liczbę pseudolosową xlos leżącą w przedziale .

38 M-C – metoda akceptacji – odrzuceńJest to najczęściej stosowana metoda Monte Carlo Losujemy punkty z prostokąta (może być wielowymiarowy) [a, b] ∗ [0, d]. Liczymy proporcje punktów leżących pod wykresem do ilości wszystkich punktów lub leżących nad wykresem. k – liczba punktów pod wykresem, N – liczba losowanych punktów. Przykład: policzmy całkę Losujemy N punktów z kwadratu [0, 1] i patrzymy ile z nich leży pod krzywą.

39 M-C– metoda akceptacji – odrzuceń

40 M-C– metoda akceptacji – odrzuceń

41 Zalety i wady metod M-C W metodzie akceptacji – odrzuceń często trudno jest wyznaczyć wysokość prostokąta. Musimy znać wartość maksymalną i minimalną funkcji w przedziale całkowania. Pojawiają się również kłopoty w przypadku, gdy funkcja zmienia znak. Wady – ogólnie: eksperymenty dla jednak skończonej liczby prób wyniki zawsze będą przybliżeniem wyniki zależą od jakości generatora liczb pseudolosowych Zalety możliwość rozwiązania trudnych problemów prosta forma zastąpienia rozwiązań analitycznych rosnąca moc obliczeniowa komputerów uwalniają użytkownika od skomplikowanej teorii i wzorów, pozwalając skupić się na istocie pytania, na które statystyka ma odpowiedzieć

42 Literatura 1. http://matematyka.pisz.pl/strona/2110.html2. 3. .http://www.matemaks.pl/wprowadzenie-do-calek.html 4. 5. 6. 7. 42