Całkowanie numeryczne

Całkowanie numeryczne – metoda numeryczna^[1] polegająca na przybliżonym obliczaniu całek oznaczonych. Termin kwadratura numeryczna, często po prostu kwadratura, jest synonimem całkowania numerycznego, w szczególności w odniesieniu do całek jednowymiarowych. Dwu- i wielowymiarowe całkowania nazywane są czasami kubaturami, choć nazwa kwadratura odnosi się również do całkowania w wyższych wymiarach.

Proste metody całkowania numerycznego polegają na przybliżeniu całki za pomocą odpowiedniej sumy ważonej wartości całkowanej funkcji w kilku punktach. Aby uzyskać dokładniejsze przybliżenie dzieli się przedział całkowania na niewielkie fragmenty. Ostateczny wynik jest sumą oszacowań całek w poszczególnych podprzedziałach. Najczęściej przedział dzieli się na równe podprzedziały, ale bardziej wyszukane algorytmy potrafią dostosowywać krok do szybkości zmienności funkcji.

Metoda prostokątów

Najprostsza z metod kwadraturowych polega na zastosowaniu wzoru

${\displaystyle \int _{x_{0}}^{x_{n}}\!\!f(x)\approx h\sum _{i=0}^{n-1}f(x_{i}+\alpha h),\quad h={\tfrac {x_{n}-x_{0}}{n}},}$

w którym ${\displaystyle n}$ jest liczbą podprzedziałow o długości ${\displaystyle h.}$

Metoda ta ma trzy warianty:

• lewych prostokątów, gdy ${\displaystyle \alpha =0,}$
• średnich prostokątów, gdy ${\displaystyle \alpha ={\tfrac {1}{2}}}$   – ten wariant daje najlepsze przybliżenie,
• prawych prostokątów, gdy ${\displaystyle \alpha =1.}$

Istnieje oczywiście wariant ogólny, gdy ${\displaystyle \alpha \in [0,\,1].}$

Metoda trapezów

Metoda trapezów polega na tym, że całkowaną funkcję aproksymujemy liniowo w każdym z podprzedziałów ${\displaystyle x_{i+1},\,x_{i},\;i=0,\,1,\,\dots \,n-1}$ o długości ${\displaystyle h={\tfrac {x_{n}-x_{0}}{n}}.}$ Dzięki temu otrzymujemy po wprowadzeniu oznaczenia ${\displaystyle f_{i}\equiv f(x_{i})}$

${\displaystyle \int _{x_{0}}^{x_{n}}\!\!f(x)dx\approx \sum _{i=0}^{n-1}{\tfrac {h}{2}}[f(x_{i+1})+f(x_{i}))]=h({\tfrac {1}{2}}f_{0}+f_{1}+f_{2}+\,\dots \,+f_{n-1}+{\tfrac {1}{2}}f_{n}).}$

Oszacowanie błędu tej metody wynosi

${\displaystyle R_{n}=\left|\int \limits _{a}^{b}f(x)dx-S_{n}\right|\leqslant {\frac {(b-a)^{3}M'{'}}{12n^{2}}},}$

gdzie:

${\displaystyle M'{'}=\max _{\langle a,b\rangle }|f'{'}|.}$

Metoda parabol (Simpsona)

Ta metoda wymaga podzielenia przedziału całkowania na parzystą liczbę ${\displaystyle 2n}$ podprzedziałów, tzn.

${\displaystyle h={\frac {x_{2n}-x_{0}}{2n}},}$

${\displaystyle f_{i}=f(x_{i}),\quad x_{i}=x_{0}+ih,\quad i=0,\,1,\,\dots \,2n.}$

Stosując kwadratową interpolację Lagrange’a na dwu sąsiadujących podprzedziałach, otrzymujemy

${\displaystyle \int _{x_{i}}^{x_{i+2}}f(x)dx\approx {\tfrac {h}{3}}[f_{i}+4f_{i+1}+f_{i+2}],}$

${\displaystyle \int _{x_{0}}^{x_{2n}}f(x)dx\approx {\tfrac {h}{3}}[f_{0}+f_{2n}+4(f_{1}+f_{3}+\,\dots \,+f_{2n-1})+2(f_{2}+f_{4}+\,\dots \,f_{2n-2})].}$

Przy całkowaniu funkcji ${\displaystyle f(x)\in C^{(4)}}$ na przedziale ${\displaystyle [a,\,b]}$ błąd metody wynosi

${\displaystyle R=-{\tfrac {(b-a)h^{4}}{180}}f^{IV}(x),\quad x\in [a,\,b].}$

Metoda Gaussa

Istotne podwyższenie dokładności wzorów kwadraturowych można uzyskać, stosując metodę Gaussa^[1]. Jej istota polega na minimalizacji błędu kwadratury przez optymalny wybór położenia węzłów ${\displaystyle \xi _{i}}$ oraz wartości wag ${\displaystyle {\overline {w}}_{i}}$ we wzorze kwadraturowym

w którym ${\displaystyle F(\xi )=f({\tfrac {a+b}{2}}+{\tfrac {b-a}{2}}\xi ),\quad dx={\tfrac {b-a}{2}}d\xi .}$

Dzięki zastosowanemu we wzorze (a) odwzorowaniu ${\displaystyle x={\tfrac {a+b}{2}}+{\tfrac {b-a}{2}}\xi }$ dowolnego przedziału ${\displaystyle (a,\,b)}$ na standardowy przedział ${\displaystyle (-1,\,1),}$ wzór ten jest uniwersalny ponieważ unikalne wartości ${\displaystyle {\overline {w}}_{i},\,\xi _{i}}$ można stablicować raz na zawsze.

Obliczenie wartości ${\displaystyle {\overline {w}}_{i},\,\xi _{i},\;i=1,\,2,\,\dots ,\,n}$ można przeprowadzić żądając, aby całkowanie procedurą Gaussa jednomianów ${\displaystyle x^{k},\;k=0,\,1,\,\dots ,\,2n-1}$ dawało wyniki ścisłe

to znaczy aby było dla ${\displaystyle k=0,\,1,\,\dots ,2n-1}$

${\displaystyle {\overline {w}}_{1}\xi _{1}^{k}+{\overline {w}}_{2}\xi _{2}^{k}+\dots +{\overline {w}}_{n}\xi _{n}^{k}=p_{k},\quad p_{k}={\tfrac {1-(-1)^{k+1}}{k+1}}.}$

Po rozpisaniu otrzymujemy, trudny do rozwiązania, nieliniowy układ ${\displaystyle 2n}$ równań

określający ${\displaystyle 2n}$ wartości ${\displaystyle {\overline {w}}_{i},\,\xi _{i}.}$

W poniższej tabeli zestawiono^[2] obliczone wartości parametrów ${\displaystyle {\overline {w}}_{i},\,\xi _{i}}$ dla wielomianów stopni ${\displaystyle n\leqslant 8.}$

Cytowany w literaturze^[1] sposób rozwiązania układu równań (c) polega na spostrzeżeniu, że dla dowolnie przyjętych wartości liczbowych ${\displaystyle \xi _{1}<\xi _{2}<\,\dots \,<\xi _{n}}$ macierz współczynników ${\displaystyle \xi _{i}^{k},}$ pierwszych ${\displaystyle n}$ równań tego układu, jest macierzą Vandermonde’a. Dzięki temu wiadomo, że istnieje jednoznacznie określone rozwiązanie ${\displaystyle w_{i},\;i=1,\,2,\,\dots \,n.}$ Kwestią otwartą pozostaje jednak wyznaczenie optymalnych wartości ${\displaystyle \xi _{i}.}$

W tym celu warunek (b) modyfikuje się do postaci obowiązującej dla wielomianów stopnia ${\displaystyle N=n+k,\;k=0,\,1,\,\dots ,\,n-1.}$

gdzie ${\displaystyle P_{n}(x)}$ jest wielomianem stopnia ${\displaystyle n.}$

Całka we wzorze (d) zeruje się, gdy wielomiany ${\displaystyle P_{n}(x)}$ są ortogonalne z jednomianami ${\displaystyle x^{k}}$ dla ${\displaystyle k<n,}$ przy ${\displaystyle x\in (-1,\,1).}$ Dokładnie taką własność^[3] mają wielomiany Legendre’a. Dla nich mamy wtedy zamiast (d)

Ten warunek jest spełniony tożsamościowo dla dowolnych wartości ${\displaystyle w_{i},}$ gdy ${\displaystyle \xi _{i}}$ są pierwiastkami wielomianu Legendre’a ${\displaystyle n}$-tego stopnia, wtedy bowiem ${\displaystyle P_{n}(\xi _{i})\equiv 0,\;i=1,\,2,\,\dots \,n.}$

Metody losowe

Do przybliżonego obliczania całki oznaczonej można również wykorzystać metody probabilistyczne. Należy pamiętać jednak, że wynik takiego całkowania jest też zmienną losową.

• Monte Carlo

Idea opiera się na policzeniu pola pod wykresem funkcji dla ${\displaystyle f(x)>0}$ i odjęciu pola nad wykresem dla ${\displaystyle f(x)<0}$

• probabilistyczna

${\displaystyle \int \limits _{a}^{b}f(x)dx\approx {\frac {b-a}{n}}\sum _{i=1}^{n}f(x_{i}),}$

${\displaystyle x_{i}}$ jest losowo wybierane z przedziału ${\displaystyle <a,b>,}$

${\displaystyle n}$ określa liczność próbki.

Przykłady

Przykład – metoda prostokątów

Spróbujmy scałkować funkcję ${\displaystyle \cos(x)}$ na przedziale od 0 do 1. Ponieważ da się ją scałkować analitycznie, znamy dokładny wynik i możemy łatwo obliczać błąd przybliżenia różnych metod całkowania. Z dokładnością do 10 miejsc dziesiętnych prawidłowy wynik wynosi:

${\displaystyle \int \limits _{0}^{1}\cos(x)dx=\sin(1)-\sin(0)=0{,}8414709848.}$

Całkowanie numeryczne za pomocą zasady punktu środkowego da nam wynik:

${\displaystyle \int \limits _{0}^{1}\cos(x)dx\approx (1-0)\cos \left({\frac {1}{2}}\right)=0{,}8775825619,}$

co daje błąd 0,0361115771 (błąd względny 4,3%) – niewielki jak na tak prostą metodę, jednak oczywiście niezadowalający do wielu zastosowań.

Żeby uzyskać lepsze przybliżenia możemy podzielić przedział całkowania:

${\displaystyle \int \limits _{0}^{1}\cos(x)dx=\int \limits _{0}^{1/2}\cos(x)dx+\int \limits _{1/2}^{1}\cos(x)dx\approx }$

${\displaystyle {}\quad \approx \left({\frac {1}{2}}-0\right)\cos \left({\frac {1}{4}}\right)+\left(1-{\frac {1}{2}}\right)\cos \left({\frac {3}{4}}\right)=0{,}8503006452}$

z błędem bezwzględnym 0,0088296604 lub względnym 1%.

Dzieląc przedział całkowania na więcej fragmentów, możemy uzyskać lepsze przybliżenie:

Przykład 2

Całkowanie numeryczne przebiegów czasowych. Spróbujmy scałkować spróbkowany przebieg ${\displaystyle \sin(t)}$ na przedziale od 0 do ${\displaystyle 4\cdot \pi }$ [s]. Oznaczmy częstotliwość próbkowania przebiegu przez ${\displaystyle f_{p}}$ [Hz].

Do obliczeń wykorzystamy metodę prostokątów. Średnica podziału ${\displaystyle t_{p}={\frac {1}{f_{p}}}=t_{i+1}-t_{i}}$ wynosi 1. Niech ${\displaystyle X_{i}(t)}$ oznacza próbkę po całkowaniu. Każdy wyraz ${\displaystyle X_{i}}$ można obliczyć jako sumę częściową:

${\displaystyle X_{i}=\sum _{n=0}^{i}x_{i}(t)t_{p}.}$

Zobacz też

• kwadratury Gaussa
• metody Newtona-Cotesa

Przypisy