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)

 Osobny artykuł: Metoda 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

 Osobny artykuł: Kwadratury 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

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

${\displaystyle \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 ${\displaystyle k=0,\,1,\,\dots ,2n-1}$

${\displaystyle {\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 ${\displaystyle 2n}$ równań

${\displaystyle {\overline {w}}_{1}+{\overline {w}}_{2}+\ldots +{\overline {w}}_{n}=2,}$ ${\displaystyle {\overline {w}}_{1}\xi _{1}+{\overline {w}}_{2}\xi _{2}+\ldots +{\overline {w}}_{n}\xi _{n}=0,}$ ${\displaystyle {\overline {w}}_{1}\xi _{1}^{2}+{\overline {w}}_{2}\xi _{2}^{2}+\ldots +{\overline {w}}_{n}\xi _{n}^{2}={\tfrac {2}{3}},}$ ${\displaystyle {\overline {w}}_{1}\xi _{1}^{3}+{\overline {w}}_{2}\xi _{2}^{3}+\ldots +{\overline {w}}_{n}\xi _{n}^{3}=0,}$ ${\displaystyle .............................}$ ${\displaystyle {\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}}}$ ${\displaystyle .............................}$ ${\displaystyle {\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}},}$ ${\displaystyle {\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 ${\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.}$

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

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

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.}$

${\displaystyle \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 ${\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)

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

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