Całkowanie numeryczne
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.
Spis treści
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 zielono pole:
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\).
Obszar całkowania \([x_{p}, x_{k}]\) 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(x_{i})\) 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:
#include <iostream>
#include <cstdlib>
using namespace std;
// funkcja do scalkowania
double f1(double x) { return -x*x-x+10; }
int main()
{
float xp, xk, h, calka;
int n;
// przedzialy
xp = -2;
xk = 1;
// im wieksze n tym wieksza dokladnosc (tym wiecej prostokątow)
n = 3;
h = (xk - xp) / (float)n;
cout << "krok: h=" << h << endl;
calka = 0;
for (int i=1; i<=n; i++)
{
calka += f1(xp + i*h)*h;
}
cout << "Wynik calkowania metoda prostokatow: " << calka << endl;
system("PAUSE");
return 0;
}
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\).
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ść liczyli byś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(x_{i})\). 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(x_{0})\), 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 (\(P_{c}\)) 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 \(x_{p}\) i \(x_{k}\) (czyli \(x_{0}\) i \(x_{3}\)):
\(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:
#include <iostream>
#include <cstdlib>
using namespace std;
// funkcja do scalkowania
double f1(double x) { return -x*x-x+10; }
int main()
{
float xp, xk, h, calka;
int n;
// przedzialy
xp = -2;
xk = 1;
// im wieksze n tym wieksza dokladnośc (np. n=1000)
n = 3;
h = (xk - xp) / (float)n;
cout << "krok: h=" << h << endl;
calka = 0;
for (int i=1; i<n; i++)
{
calka += f1(xp + i * h);
}
calka += f1(xp) / 2;
calka += f1(xk) / 2;
calka *= h;
cout << "Wynik calkowania: " << calka << endl;
system("PAUSE");
return 0;
}
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ą.
