∫ Całkowanie numeryczne

Całkowanie numeryczne

Nie zawsze możliwe jest wyznaczenie analitycznego wzoru będącego wynikiem całkowania danej funkcji f(x).

Praktycznie zawsze możne jednak wyznaczyć całkę oznaczoną funkcji przy podanych granicach całkowania a i b:

( )

b

a

f x dx

∫

Pamiętajmy, że człon dx oznacza wyłącznie różniczkę (pochodną) zmiennej, względem której zachodzi całkowanie i że człon ten nie wchodzi do obliczeń!

Wiadomo, że całka oznaczona funkcji f(x) w granicach a i b jest równa polu pod wykresem, ograniczonym wartościami a i b:

Oczywiście podanie ogólnego wzoru na pole figury o dowolnym kształcie jest niemożliwe. Problem ten rozwiązuje się poprzez podzielenie go na szereg bardzo wąskich pasków, wyznaczenie ich pól i dodanie do siebie:

Czasem pola te wylicza się zakładając, że paski te są prostokątami (tzw. metoda prostokątów) Jak widzimy jednak, znacznie lepszym przybliżeniem kształtu paska jest „leżący na boku” trapez, o wysokości dx oraz podstawach f(x0) i f(x1). Jeżeli policzymy pole każdego trapezu i dodamy je do siebie, otrzymamy pole powierzchni całego obszaru pod wykresem. Oczywiście, im węższe są trapezy (czyli im jest ich więcej), tym dokładność obliczeń wyższa – jest to związane ze zwykle złym dopasowaniem „ukośnego” boku trapezu do krzywizny funkcji. Przy trapezach o bardzo małej wysokości (czyli małym dx) różnica między rzeczywistym polem figury ograniczonej podstawami trapezu a trapezem maleje, co zwiększa dokładność obliczeń.

W praktyce oznacza to konieczność dzielenia zakresu całkowania na możliwie jak największą liczbę odcinków dx.

Pole pojedynczego (i-tego) trapezu wyznaczymy ze znanego wzoru:

( ) ( )

2

i i

i

f x f x P = dx ⋅ +

natomiast pole całego obszaru od wykresem będzie reprezentowane wzorem:

( ) ( )

1

1

1

2

n

i i

i

f x f x

P dx

− +

=

= ∑ ⋅ +

Dla przykładu policzymy, jaka jest wartość następującej całki:

0

sin x x dx

π

∫ ⋅

Na początku musimy zawsze wpisać do arkusza dolną i górną granicę całkowania, jak również liczbę odcinków, na którą dzielimy obszar całkowania (n) – zwykle liczba ta jest równa przynajmniej 100, jednak tutaj, dla ułatwienia, ustalimy ją równą zaledwie 10 (proponuję przypisać nazwy do zawartości odpowiednich komórek, w celu późniejszego łatwiejszego nimi operowania). Następnie wyznaczam wartość dx (wysokość trapezu), będącą po prostu ilorazem zakresu całkowania (górnej granicy całkowania minus dolnej granicy całkowania) oraz liczby n.

W efekcie otrzymujemy:

Oznacza to, że granice całkowania od 0 do π podzieliliśmy na 10 odcinków, z których każdy ma długość 0,314159265. W dalszej kolejności tworzymy kolumnę kolejnych wartości x, począwszy od dolnej granicy całkowania, a skończywszy na górnej granicy całkowania…

…przy czym każda następna wartość x jest większa od poprzedniej o dx.

Powielanie tej formuły kontynuujemy aż do osiągnięcia wartości określonej przez komórkę górna:

Warto zauważyć, że skoro pierwszy element kolumny x ma adres B10, a zdefiniowaliśmy liczbę n równą 10, to ostatnia komórka musi mieć adres B(10+n), czyli B20. Można to traktować jako wzajemną kontrolę poprawności wykonania tej części zadania – komórka ostatnia musi mieć wartość równą wartości górna i jednocześnie musi mieć właściwy adres.

W dalszej kolejności obok kolumny x musimy utworzyć kolumnę f(x) (o tej samej długości), której każda komórka zawierała będzie wartość całkowanej funkcji f(x) (w naszym przykładzie x·sin(x)) dla danej (sąsiedniej) wartości x:

W tym momencie pozostaje już tylko policzyć pola kolejnych trapezów, zgodnie z wcześniej podanym wzorem:

UWAGA! Zauważmy, że trapezy powielamy do przedostatniej komórki powstającej tabeli – w komórce C20 znajduje się wartość górnej (prawej) krawędzi ostatniego trapezu, który zaczyna się w komórce C19!

Pozostaje nam już tylko policzyć sumę wszystkich trapezów:

Czyli w rezultacie:

Wyznaczona w ten sposób wartość przykładowej całki jest równa 3,115711487. Ma ona jednak dwie wady:

• Niską dokładność obliczeń – można ją zwiększyć zwiększając liczbę podziałów n

• Nie znamy dokładności otrzymanego wyniku – nie wiemy, ile miejsc znaczących możemy wypisać Sposobem, który częściowo rozwiązuje obydwa problemy, jest obliczenie tej samej całki drugą metodą. Jest to jednocześnie sposób na autokontrolę otrzymanego wyniku – wartości otrzymane obiema metodami powinny być wzajemnie do siebie podobne.

Tą drugą, godną polecenia metodą obliczenia numerycznego całek oznaczonych, jest metoda 1/3 Simpsona.

W metodzie tej pole pod wykresem obliczane jest jako suma pól nie trapezów, lecz figur składających się z pary opisywanych wcześniej, sąsiadujących ze sobą pasków, przy czym bok powstałej w ten sposób figury, przyległy do krzywej funkcyjnej, przybliżony jest parabolą:

W efekcie pole całkowite obszaru pod krzywą funkcyjną opisane będzie wzorem:

( ) ( ) ( )

2

1 2

1,3,5,

1 4

3

n

i i i

i

P dx f x f x f x

−

+ +

=

= ∑ + ⋅ +

…

Proszę zwrócić uwagę na indeksy i zasięg sumowania – „pomijane” są paski parzyste gdyż, jak wspomniano, każdy „pasek Simpsona” składa się z dwóch pasków. Oznacza to, że metoda Simpsona wymaga podzielenia obszaru całkowania na parzystą liczbę pasków, co zwykle nie stanowi problemu.

Pokażę teraz, jak policzyć metodą Simpsona wartość całki, którą wcześniej liczyliśmy metodą trapezów. Dla ułatwienia użyję tabelki stworzonej na potrzeby metody trapezów (trzeba uważać, żeby nie pomylić kolumn!) Dla ułatwienia i zachowania przejrzystości w pola tworzonej kolumny wpisuję jedynie „człon Simpsona” (czyli wyrażenie pod sumą), musze jednak pamiętać o tym, aby potem trzymaną sumę pomnożyć przez 1

3dx.

Następnie zaznaczamy bok komórek składający się z komórki zawierającej wynik obliczeń „członu Simpsona”

i komórki kolejnej, i ten właśnie blok powielamy, aż do pozycji E18 (jest to początek ostatniego „paska Simpsona”):

w efekcie otrzymując:

Aby zakończyć obliczanie wartości całki, należy policzyć sumę wszystkich „Simpsonów” i pomnożyć ją przez 1

3dx :

W rezultacie otrzymujemy:

Całka, której wartość właśnie policzyliśmy numerycznie, posiada rozwiązanie analityczne i jest ono równe π, czyli 3,14159265… Widzimy, że w omawianym przykładzie metoda 1/3 Simpsona dała wynik znacznie bliższy rozwiązaniu analitycznemu. W przypadku ogólnym nieznanej funkcji f(x) nie możemy być jednak pewni, która metoda jest dokładniejsza. Przez porównanie wartości, otrzymanych obiema metodami mogę jednak ustalić zakres zgodności – w naszym przykładzie obydwa wyniki są zgodne do pierwszego miejsca dziesiętnego, czyli mogę jedynie powiedzieć, że rozwiązaniem zadania jest liczba 3,1, pozostałe cyfry są w zasadzie niepewne.

Zwiększyć dokładność mogłoby zwiększenie liczny n, czyli ilości pasków, na który dzielimy obszar podcałkowy.

Obliczenia numeryczne przy bardzo dużej liczbie komórek

Poniżej przedstawię sposób pracy przy liczbie pasków równej 65000 (liczbe tę wybrałem, gdyż maksymalna liczba wierszy w arkuszu wynosi 65536). Policzona będzie ta sama całka, co we wcześniejszych przykładach.

1. Ustalam parametry całkowania numerycznego:

2. Zawartość komórki B11 (formułę) kopiuję do schowka

Skoro pierwsza komórka kolumny x ma adres B10, to ostatnia będzie miała adres B(10+n), czyli B65010. Musimy wkleić „masowo” formułę z komórki B11 do wszystkich komórek od B12 do B65010.

3. W tym celu przechodzimy do komórki B65010, wpisując ten adres do pola nazwy i zatwierdzając poprzez naciśnięcie przycisku Enter.

4. Wciskając kombinację klawiszy Ctrl Shift Home zaznaczamy cały obszar od komórki B65010 do A1.

Trzymając nadal wciśnięty klawisz Shift przy pomocy klawiszów kursora („strzałek”) ustawiamy zaznaczenie tak, aby zaczynało się w komórce B12 i kontynuowało się w dół kolumny B:

5. Puszczamy klawisz Shift i naciskamy klawisz Enter, co powoduje wypełnienie wszystkich zaznaczonych komórek (od B10 do B65010) znajdującą się w schowku formułą:

UWAGA! Kombinacja klawiszy Ctrl Home i Ctrl End umożliwia szybkie poruszanie się po arkuszu w granicach określonych przez wstawione dane.

6. Tak, jak opisano w metodzie trapezów, wyznaczamy wartość y, poczym zaznaczamy ją do schowka.

7. Poprzez naciśnięcie kombinacji klawiszy Ctrl Shift End zaznaczamy wszystkie komórki począwszy od C10 aż do C65010 (czyli do końca danych)

8. Naciskamy klawisz Enter, wklejając zaznaczoną formułę do wszystkich zaznaczonych komórek:

9. Dokładnie w ten sam sposób wypełniamy kolejną kolumnę (trapezy), pamiętając jednak, że zaznaczenie musi się skończyć na komórce przedostatniej, czyli 65009.

10. Sumujemy wszystkie komórki z kolumny trapezy

otrzymując:

co stanowi bardzo dobre przybliżenie wartości wyznaczanej całki (przypominam, że wartość w tym przypadku powinna wynosić π).

Proponuję tę samą procedurę powtórzyć dla metody 1/3 Simpsona. Wynik powinien wyglądać następująco (obliczanie sumy może trochę w tym przypadku potrwać):

Jak można sprawdzić, tym razem otrzymaliśmy zgodność z liczbą π sięgającą 11 miejsc dziesiętnych!

Niestety, dla innych całek takiej kontroli najczęściej nie mamy – gdybyśmy nie znali wyniku analitycznego opisywanej przykładowej funkcji, we wniosku końcowym należałoby napisać, że na podstawie porównania wartości wyznaczonej metodą trapezów i metodą 1/3 Simpsona wartość całki

0

sin 3,141592653 x x dx

π ⋅ =

∫

Z przykładem obliczenia wartości całki oznaczonej można zapoznać się na stronie:

http://vistula.pk.edu.pl/~sciezor/calkowanie.xls