P-Programowanie

Całkowanie numeryczne C++

Ostatania modyfikacja: 28 września 2017, kategoria: C++

Całkowanie numeryczne jest to zagadnienie z metod elementów skończonych (MES). Korzystając z całkowania numerycznego możemy obliczyć wartość dowolnej całki jednowymiarowej oznaczonej. Wynik jest zawsze obarczony pewnym błędem, ponieważ wartość całki w C++ liczona jest w przybliżeniu.

Czym jest całka?

Całka oznaczona z dowolnej funkcji f(x) w przedziałach [a,b], jest to pole figury ograniczone wykresem funkcji f(x) oraz osią X w przedziałach [a,b].

Dla przykładu rozważmy funkcję f(x)=-x^2-x+10. Narysuję jej wykres w przedziale [-4,3]. Całkę z tej funkcji \int\limits_{-2}^{1}f(x)dx przedstawia oznaczone na niebiesko pole:

całka C++

Całkowanie numeryczne

W artykule opiszę kilka podstawowych metod całkowania numerycznego i skupię się na przedstawieniu ich implementacji w języku C++. Różne metody całkowania numerycznego dają różną dokładność. Są metody dokładniejsze i mniej dokładne, jednak ich prawdziwa skuteczność zależy od rodzaju funkcji jaką badamy.

Każda metoda całkowania numerycznego opiera się na tej samej koncepcji. Funkcję, którą całkujemy dzielimy na określoną liczbę przedziałów. Ilość przedziałów jest dowolna i zależy od niej dokładność wyniku, który otrzymamy. Wynikiem całki jest suma pól poszczególnych przedziałów.

Metoda prostokątów C++

Metoda prostokątów jest najmniej dokładną metodą całkowania. Polega ona na podzieleniu obszaru całowanego na x prostokątów (przedziałów). Pola prostokątów dodaje się do siebie i otrzymuje przybliżony wynik całki.

Posłużmy się ponownie funkcją f(x)=-x^2-x+10 oraz całką z tej funkcji \int\limits_{-2}^{1}f(x)dx.

metoda prostokątów C++

Obszar całkowania [Xp, Xk] podzieliłem na n równych części, w tym wypadku n=3.

Litera h oznacza krok całkowania czyli wysokość pojedynczego prostokąta. Jej wartość wyznaczamy ze wzoru:

h=\frac{x_{k}-x_{p}}{n} h=\frac{1-(-2)}{3}=1

Pole pojedynczego prostokąta określa wzór:

P_{i}=f(x_{i}) \cdot h, dla i=1..3

Aby poznać długość podstawy f(xi) musimy użyć iteratora. Występuje prosta zależność:

P_{i}=f(x_{p}+i \cdot h) \cdot h

Sprawdźmy obliczenia podstawiając odpowiednie liczby do powyższego wzoru:

P_{1}=f(-2+1 \cdot 1) \cdot 1=f(-1) \cdot 1 = 10 \cdot 1 = 10 P_{2}=f(-2+2 \cdot 1) \cdot 1=f(0) \cdot 1 = 10 \cdot 1 = 10 P_{3}=f(-2+3 \cdot 1)\cdot 1=f(1) \cdot 1 = 8 \cdot 1 = 8

Aby otrzymać ostateczny wynik całki należy zsumować pola prostokątów:

\int\limits_{-2}^{1}f(x)dx=P_{1}+P_{2}+P_{3}=10+10+8 \approx 28

Kod programu w C++ liczący dowolną całkę metodą prostokątów:

 

Metoda trapezów C++

Metoda trapezów daje dość dokładne wyniki, szczególnie przy dużej dokładności N. Polega ona na podzieleniu obszaru całowanego na x trapezów (przedziałów), których pola na końcu sumuje się.

W dalszym ciągu wykorzystajmy funkcję f(x)=-x^2-x+10. Obliczmy z niej całkę: \int\limits_{-2}^{1}f(x)dx.

metoda trapezów C++

Obszar został podzielony na 3 równe części (n=3). W tym wypadku 3 przedziały dają w miarę dokładny wynik (trapezy są bardzo dobrze dopasowane), jednak aby zyskać większą dokładność warto utworzyć przynajmniej 1000 podziałów.

Litera h jest tym razem wysokością trapezu. Ponieważ mamy tylko 3 przedziały, wysokość łatwo odczytać z rysunku. Gdybyśmy chcieli większą dokładność, moglibyśmy wygenerować np. 100 trapezów, wtedy wysokość liczylibyśmy ze wzoru:

h=\frac{x_{k}-x_{p}}{n} h=\frac{1-(-2)}{3}=1

Wykorzystajmy wzór na pole trapezu znany z gimnazjum:

P=\frac{a+b}{2} \cdot h

Przekształćmy wzór dostosowując go do naszego rysunku. Pole i-tego trapezu określa wzór:

P_{i}=\frac{f(x_{i-1})+f(x_{i})}{2} \cdot h

Zauważ, że zamiast podstaw a i b we wzorze pojawiły się wartości funkcji f(xi). Długości podstaw trapezów kończą się w miejscu gdzie przecinają się z wykresem dlatego można je obliczyć wprost z funkcji. Przykładowo aby poznać długość górnej podstawy pierwszego trapezu wystarczy obliczyć wartość funkcji f(x0), a więc do funkcji wstawiamy za x liczbę -2. Otrzymamy wynik 8 co można w tym wypadku odczytać z rysunku.

Wartość całki (Pc) w przybliżeniu równa jest sumie pól wszystkich trapezów. W naszym przypadku posiadamy 3 trapezy stąd:

P_{c}=P_{1}+P_{2}+P_{3}

Rozpisując wzór dla 3 trapezów dojdziemy do postaci:

P_{c}=\frac{f(x_{0})+f(x_{1})}{2} \cdot h + \frac{f(x_{1})+f(x_{2})}{2} \cdot h + \frac{f(x_{2})+f(x_{3})}{2} \cdot h

Wzór można uprościć, wyłączając przed nawias część wspólną poszczególnych sum:

P_{c}=\frac{h}{2} \cdot (f(x_{0}) + f(x_{1})+f(x_{1})+f(x_{2})+f(x_{2})+f(x_{3})) P_{c}=\frac{h}{2} \cdot (f(x_{0}) + 2 \cdot f(x_{1})+2 \cdot f(x_{2})+f(x_{3})) P_{c}=h \cdot (\frac{f(x_{0})}{2} + f(x_{1})+ f(x_{2})+\frac{f(x_{3})}{2})

Wyprowadzony wzór podaje nam w przybliżeniu wartość liczonej przez nas całki, przy podziale na 3 trapezy.

Zwróć uwagę: mając określoną liczbę trapezów n, sumujemy ich wszystkie podstawy pamiętając aby pierwszą i ostatnią podzielić przez 2. Na końcu mnożymy wynik sumy przez wysokość pojedynczego trapezu.

Całość bardzo łatwo zaimplementować w C++ posługując się pętlą for. W pętli sumujemy wszystkie podstawy, pomijając graniczne xp i xk (czyli x0 i x3):

Podst_{i}=f(x_{p} +i \cdot h), dla i=1..2

Po wyjściu z pętli liczymy długości podstaw brzegowych zgodnie z wyprowadzonym wzorem i dzielimy je przez 2:

Podst_{0}=\frac{f(x_{p})}{2} Podst_{3}=\frac{f(x_{k})}{2}

Na samym końcu sumujemy poszczególne podstawy i mnożymy przez krok H (wysokość):

P_{c}=1 \cdot (\frac{8}{2} +10+ 10+\frac{8}{2})= 4+10+10+4 \approx 28

Kod programu w C++ liczący dowolną całkę metodą trapezów:

Podsumowanie

Istnieje wiele metod całkowania numerycznego, które dają lepsze wyniki:

  • metoda prostokątów
  • metoda trapezów
  • metoda Simpsona (metoda parabol)
  • metoda Monte Carlo
  • reguła 3/8
  • cały zbiór kwadratur Gaussa
  • itd.

Każda z użytych metod (nawet metoda prostokątów) może dać bardzo dokładny wynik obliczonej całki (wyobraź sobie, że dzielisz przedział [0, 1] na 10000 małych prostokątów). Im większe przybliżenie, tym więcej operacji musi wykonać program.

Różne metody całkowania numerycznego są używane do rozwiązywania różnych problemów. Bardzo często metody są łączone ze sobą.

Użytkownik Eryk napisał:

15 grudnia 2013


Artykuł przydatny ale po wpisaniu tego wzroy prz biblliotece iostream wyskakuje błąd ;)

Użytkownik Karol napisał:

15 grudnia 2013


Musisz pozmieniać biblioteki zależnie od tego na jakim kompilatorze kompilujesz kod.

Użytkownik Eryk napisał:

21 grudnia 2013


Programuję w Code::Blocks ale myślę, że zmiana biblioteki rozwiąże problem ;)

Użytkownik Mikołaj napisał:

08 stycznia 2015


Jeśli znamy trochę matematyki, można to rozwiązać jeszcze inaczej i dużo dokładniej. Wystarczy policzyć ręcznie funkcję pierwotną (taką, której pochodną będzie funkcja wyjściowa) i odjąć jej wyniki dla końców przedziałów. Na przykład: http://pastie.org/9820433

Użytkownik keraj napisał:

21 stycznia 2015


@Mikołaj

Pewnie, że można. Ale po co wtedy pisać program?

Użytkownik Filip napisał:

28 czerwca 2016


Co w przypadku liczenia całki „online”? np mamy na wyjściu innej funkcji sygnał prędkości i chcemy obliczyć z niego całkę. Nie znamy czasu końcowego (chyba, że jest nim czas aktualny?) całkowania. Sygnał prędkości obliczany jest również online.

Użytkownik Gura napisał:

03 kwietnia 2017


No fajny program, dobry jest działa sobie całkiem fajnie liczy całke dobry.

Użytkownik jan napisał:

24 lutego 2018


metoda prostokątów jest trochę inna f((xi + x(i+1))/2) i czasami jest dokładniejsza od metody trapezów zależności od rozkładu funkcji

Użytkownik bartosz napisał:

09 stycznia 2019


fajny programik ale mi nie dziala

Zachęcam Cię do zostawienia komentarza!

Ilość znaków: 0