Całkowanie numeryczne

Całkowanie numeryczne – metoda numeryczna 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

$\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 $n$ jest liczbą podprzedziałow o długości $h.$

Metoda ta ma trzy warianty:

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

Istnieje oczywiście wariant ogólny, gdy $\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 $x_{i+1},\,x_{i},\;i=0,\,1,\,\dots \,n-1$ o długości $h={\tfrac {x_{n}-x_{0}}{n}}.$ Dzięki temu otrzymujemy po wprowadzeniu oznaczenia $f_{i}\equiv f(x_{i})$

$\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

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

gdzie:

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

Metoda parabol (Simpsona)

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

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

$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

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

$\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 $f(x)\in C^{(4)}$ na przedziale $[a,\,b]$ błąd metody wynosi

$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. Jej istota polega na minimalizacji błędu kwadratury przez optymalny wybór położenia węzłów $\xi _{i}$ oraz wartości wag ${\overline {w}}_{i}$ we wzorze kwadraturowym

$I(f)=\int _{a}^{b}\!\!f(x)dx={\tfrac {b-a}{2}}\int _{-1}^{1}F(\xi )d\xi ={\tfrac {b-a}{2}}\sum _{i=1}^{n}{\overline {w}}_{i}F(\xi _{i}),$

(a)

w którym $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 $x={\tfrac {a+b}{2}}+{\tfrac {b-a}{2}}\xi$ dowolnego przedziału $(a,\,b)$ na standardowy przedział $(-1,\,1),$ wzór ten jest uniwersalny ponieważ unikalne wartości ${\overline {w}}_{i},\,\xi _{i}$ można stablicować raz na zawsze.

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

$\sum _{i=1}^{n}{\overline {w}}_{i}\xi _{i}^{k}\equiv \int _{-1}^{1}\!\!\xi ^{k}d\xi$

(b)

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

${\overline {w}}_{1}\xi _{1}^{k}+{\overline {w}}_{2}\xi _{2}^{k}+\ldots +{\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 $2n$ równań

${\overline {w}}_{1}+{\overline {w}}_{2}+\ldots +{\overline {w}}_{n}=2,$ ${\overline {w}}_{1}\xi _{1}+{\overline {w}}_{2}\xi _{2}+\ldots +{\overline {w}}_{n}\xi _{n}=0,$ ${\overline {w}}_{1}\xi _{1}^{2}+{\overline {w}}_{2}\xi _{2}^{2}+\ldots +{\overline {w}}_{n}\xi _{n}^{2}={\tfrac {2}{3}},$ ${\overline {w}}_{1}\xi _{1}^{3}+{\overline {w}}_{2}\xi _{2}^{3}+\ldots +{\overline {w}}_{n}\xi _{n}^{3}=0,$ $.............................$ ${\overline {w}}_{1}\xi _{1}^{k}+{\overline {w}}_{2}\xi _{2}^{k}+\ldots +{\overline {w}}_{n}\xi _{n}^{k}={\begin{cases}{\tfrac {2}{k+1}}\;\;{\text{gdy}}\;\;k\;\;{\text{jest parzyste}}\\\;\;0\quad {\text{gdy}}\;\;k\;\;{\text{jest nieparzyste}}\end{cases}}$ $.............................$ ${\overline {w}}_{1}\xi _{1}^{2n-2}+{\overline {w}}_{2}\xi _{2}^{2n-2}+\ldots +{\overline {w}}_{n}\xi _{n}^{2n-2}={\tfrac {2}{2n-1}},$ ${\overline {w}}_{1}\xi _{1}^{2n-1}+{\overline {w}}_{2}\xi _{2}^{2n-1}+\ldots +{\overline {w}}_{n}\xi _{n}^{2n-1}=0,$

(c)

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

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

${}\;n\;{}$	$\xi _{1};\,\xi _{n}$	$\xi _{2};\,\xi _{n-1}$	$\xi _{3};\,\xi _{n-2}$	$\xi _{4};\,\xi _{n-3}$
1	0
2	$\mp 0{,}57735027$
3	$\mp 0{,}77459667$	0
4	$\mp 0{,}86113631$	$\mp 0{,}33998104$
5	$\mp 0{,}90617985$	$\mp 0{,}53846931$	0
6	$\mp 0{,}93246951$	$\mp 0{,}66120939$	$\mp 0{,}23861919$
7	$\mp 0{,}94910791$	$\mp 0{,}74153119$	$\mp 0{,}40584515$	0
8	$\mp 0{,}96028986$	$\mp 0{,}79666648$	$\mp 0{,}52553242$	$\mp 0{,}18343464$

${}\;n\;{}$	${\overline {w}}_{1};\,{\overline {w}}_{n}$	${\overline {w}}_{2};\,{\overline {w}}_{n-1}$	${\overline {w}}_{3};\,{\overline {w}}_{n-2}$	${\overline {w}}_{4};\,{\overline {w}}_{n-3}$
1	$2$
2	$1$
3	$0{,}55555556$	$0{,}88888889$
4	$0{,}34785484$	$0{,}65214516$
5	$0{,}23692688$	$0{,}47862868$	$0{,}56888889$
6	$0{,}17132450$	$0{,}36076158$	$0{,}46791394$
7	$0{,}12948496$	$0{,}27970540$	$0{,}38183006$	$0{,}41795918$
8	$0{,}10122854$	$0{,}22238104$	$0{,}31370664$	$0{,}36268378$

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

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

$\sum _{i=1}^{n}w_{i}\xi _{i}^{k}P_{n}(\xi _{i})\equiv \int _{-1}^{1}\!\!\xi ^{k}P_{n}(\xi )d\xi ,\quad k=0,\,1,\,\dots \,n-1,$

(d)

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

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

$\sum _{i=1}^{n}w_{i}\xi _{i}^{k}P_{n}(\xi _{i})\equiv 0,\quad k=0,\,1,\,\dots \,n-1.$

(e)

Ten warunek jest spełniony tożsamościowo dla dowolnych wartości $w_{i},$ gdy $\xi _{i}$ są pierwiastkami wielomianu Legendre’a $n$-tego stopnia, wtedy bowiem $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 $f(x)>0$ i odjęciu pola nad wykresem dla $f(x)<0$

• probabilistyczna

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

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

$n$ określa liczność próbki.

Przykłady

Przykład – metoda prostokątów

Spróbujmy scałkować funkcję $\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:

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

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

$\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:

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

${}\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:

Liczba części	Wynik	Błąd
		bezwzględny	względny
1	0,8775825619	0,0361115771	4,29%
2	0,8503006452	0,0088296604	1,05%
4	0,8436663168	0,0021953320	0,26%
8	0,8420190672	0,0005480824	0,07%
	0,8414709848	0	0%

Przykład 2

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

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

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

Zobacz też

• kwadratury Gaussa
• metody Newtona-Cotesa