Całkowanie numeryczne
Grzegorz Kmiecik
Maciej Kopański
Włodzimierz Kiżewski
WFTiMS IS sem IV 2007
Omówienie matematyczne numerycznych metod całkowania.
Całkowanie jest jedną z ważniejszych zagadnień działu matematycznego nazywającego się metodami numerycznymi. Jest dużo metod opracowanych, jedne są dokładniejsze a inne są mniej dokładne. Ja omówię zaledwie kilka spośród wielu metod. W omówieniu tym będę dotykać zagadnienia niezwykle powierzchownie.
Metoda trapezów.
Metoda trapezów jest metodą wg mnie prostą i jednocześnie dosyć dobrze przybliżającą wartość całkowanej funkcji. Metoda ta zakłada przybliżenia za pomocą trapezu. Jest to metoda bardzo podstawowa. Na czym polega? A no na przybliżeniu wartości za pomocą trapezów:
Na powyższym rysunku widać, że funkcja ta została przybliżona za pomocą trapezów. Wzór zaś opiera się na obliczaniu kolejnych pól trapezów.
Wzór prosty: ∫(od a do b) F(x)dx = h/2(f(a)+f(b)).Jedna daje nam ogromny błąd (bo całą funkcję przybliża za pomocą JEDNEGO trapezu). Wzór nieco ciekawszy i lepszy wygląda tak:
Gdzie n jest wielkością przedziału.
Z wyprowadzenia prostego wzoru na całkowanie numeryczne można wyprowadzić metody takie jak metoda Simpsona, Simpsona 3/8 czy Bool'a. Wyprowadzenia na prosty wzór nie będę przedstawiał.
Metoda Monte Carlo
Metoda ta jest prosta obliczeniowo, ale ma dużą niedokładność porównywalną z metodą prostokątów. A więc jest ona często nie wskazana. Metoda ta opiera się na uśrednianiu wartości funkcji w przedziale pomnożoną przez szerokość przedziału. Średnia wyznaczana jest w sposób pseudolosowy jako suma n wartości funkcji w przypadkowo wybranych punktach przedziału całkowania. Wzór wygląda tak:
Gdzie n=b-a zaś a,b są granicami przedziału.
Metoda Gaussa-Legendre'a
To jest tylko jedna metoda z wielu. Mamy też inne pod metody Gaussa takie jak np: metoda Gaussa-Hermite'a czy Gaussa-Jacobiego. No i rozróżniamy też metody proste i złożone, znaczy opierające się na kwadraturze prostej i złożonej. Ja poruszę tylko metodę prostą.
Metoda ta jest stosunkowo szybko zbieżna jednakże stosunkowo trudno się ją liczy. By tego uniknąć stosuje się tylko niskiego rzędu tę metodę a dokładniej się dzieli przedział (zwiększamy n). Zakładam wagę równą 1. Wzór ogólny wygląda następująco:
Ak = - 2 / ( (N+2) * PN+2(xk) * PN+1(xk) )
xk są zerami wielomianu PN+1
Pn(x)= 1 / (2n * n!) * (dn / dxn (x2-1)n
Przedział całkowania: <a,b>
S(f)=(b-a)/2 * [∑(k=0 do N) Ak * f(tk)]
tk=(a+b)/2 + ((b-a)/2) * xk
Proszę zauważyć że S(f) jest właściwym wzorem na całkowanie. Jak widać nie jest on prosty. Można też łatwo zauważyć że nazwa tej metody wzięła się od tego, że do obliczenia stosuje się wielomian Legrende'a.
Podsumowanie
Z trzech powyższych wzorów najefektywniejsza najszybsza i najdokładniejsza jest metoda Gaussa-Legendre'a zaś najłatwiejsza do obliczeń przy zachowaniu dużej dokładności jest metoda trapezów. Najgorszą z tych wszystkich metod jest metoda Monte Carlo - jest najłatwiejsza do obliczeń, ale w zamian mamy bardzo duży błąd.
Omówienie algorytmów programów.
Całkowanie metodą trapezów można podzielić na kilka kroków przy posiadanych kilku danych początkowych.
Posiadamy początek przedziału i koniec przedziału. Znamy również n (czyli na ile części podzielimy przedział).
KROK 1: suma początkowa trapezów równa jest 0.
KROK 2: do zmiennej dx wrzucamy (b-a)/n.
KROK 3: dla licznika i=1,...,n wrzucaliśmy do sumy całkowej zawartość starej sumy dodanej do f(xp + i·dx) czyli s=s+f(xp + i·dx)
KROK 4: po wykonaniu pętli tak utworzoną sumę dodajemy do (f(xp) + f(xk))/2 i całość dodatkowo mnożymy przez dx co daje nam ostatecznie wynik. Czyli wygląda to tak: s=(s+((f(xp) + f(xk))/2))*dx.
Schemat blokowy dla całkowania metodą trapezów:
Całkowanie metodą Simpsona.
Program w celu obliczenia wartości całki wykonuje następujące kroki:
KROK1: Pobór danych wejściowych tzn: przedziałów a i b oraz liczby punktów podziałowych czyli n .
KROK2: Zerujemy wartości całki s=0 oraz wartości funkcji w punktach środkowych q=0 .
KROK3: Obliczamy dx=(b-a)/n czyli obliczamy odległość między sąsiednimi punktami podziałowymi.
KROK4: W pętli i=1 do i=n wykonujemy następujące operacje w poniższej kolejności:
obliczamy pozycje pkt podziałowego: x=a+(i*dx)
obliczamy q =q+f(x-(dx/2)) gdzie f(x) jest funkcją całkowaną
Sprawdzamy czy i<n. Jeśli tak to: s=s+f(x).
KROK5: Obliczamy wartość całki oznaczonej na przedziale <a,b> funkcji f(x) wg wzoru:
s=( dx/6 ) * ( f(a) + f(b) + 2s + 4q).
Całkowanie metodą Monte Carlo
Kroki programu które wykonujemy żeby wyliczyć całkę.
KROK 1: Pobór danych wejściowych czyli pobór przedziału oraz liczby punktów losowych.
Inicjujemy losowanie.
KROK 2: Za wartość całki podstawiamy 0 oraz obliczamy długość przedziału całkowania dx.
Czyli dx=b-a oraz s=0.
KROK 3: W pętli od i=1 do n wykonujemy dwie operacje - tzn, wpierw losujemy x-ksa zmiennoprzecinkowego pojedynczej precyzji a następnie ustalamy nową wartość sumy f(los). Czyli losujemy los oraz sumujemy poszczególne wartość f(los) do zmiennej s.
KROK 4: Po wykonaniu pętli wykonujemy ostatnią operację tzn wyliczamy właściwą wartość całki. Do nowej wartości s wrzucamy stare s (będące sumą f(los)) pomnożone dodatkowo przez dx*n.
Schemacik blokowy:
START
Wczytujemy a,b,n.
s=0;
dx=(b-a)/n;
i=1
i<=n
s=s+f(a+dx*i);
i=i+1;
NIE
s=(s+(f(a)+f(b))/2)*dx
Wypisz s
Koniec
START
Wczytujemy a,b,n
s=0
dx=b-a
i<=n
Losujemy los z <a,b>
s=s+f(los)
i++
NIE
s=s/n*dx
Wypisz s
KONIEC
START
Wczytujemy a,b,n.
s=0;
q=0;
dx=(b-a)/n;
i=1
x=a+(i*dx);
q=q+f(x-(dx/2));
i<n
i<=n
s=s+f(x)
NIE
i++
NIE
s=(dx/6)*(f(a)+f(b)=2s+4q);
Koniec