Wikipedista:Kolarp/Gaussova kvadratura

Porovnání dvoubodové Gaussovy kvadratury s lichoběžníkovou metodou. Modrá křivka je graf funkce ${\textstyle y(x),}$ jejíž určitý integrál na intervalu $\langle -1,1\rangle$ se má vypočítat.
Lichoběžníková metoda aproximuje funkci ${\textstyle y(x)}$ oranžovou čárkovanou úsečkou, která spojuje funkční hodnoty na okrajích intervalu. Z obrázku je patrné, že tato aproximace není dobrá, takže chyba spočítané hodnoty integrálu je velká (lichoběžníková metoda dává hodnotu y(–1) + y(1) = –10, zatímco správná hodnota je 2/3). Aby se lichoběžníkovou metodou dosáhlo přesnějšího výsledku, bylo by nutné rozdělit interval na mnoho podintervalů, funkci ${\textstyle y(x)}$ nahradit v každém podintervalu lineární funkcí a počítat mnoho funkčních hodnot funkce ${\textstyle y(x).}$
Naproti tomu Gaussova metoda vybírá vhodnější body integrace, takže jediná lineární funkce (znázorněná černou čárkovanou úsečkou) aproximuje funkci ${\textstyle y(x)}$ mnohem lépe. Výsledek integrace dvoubodovou Gaussovou metodou je ${\textstyle y\left(-{\sqrt {\frac {1}{3}}}\right)+y\left({\sqrt {\frac {1}{3}}}\right)={\frac {2}{3}}.}$ Tato hodnota je dokonce přesný výsledek, protože zelená oblast má stejnou plochu jako součet ploch červených oblastí. Důvodem je, že integrovaná funkce ${\textstyle y(x)}$ je ve skutečnosti polynom ${\textstyle y(x)=7x^{3}-8x^{2}-3x+3}$ stupně 3.

Numerické metody integrace jsou v numerické matematice metody výpočtu přibližných hodnot určitých integrálů funkcí. Při numerickém výpočtu určitého integrálu nějaké funkce se obvykle používá vážený součet hodnot této funkce v určitých bodech integračního intervalu. Další informace o numerických metodách integrace jsou ve článcích Numerická integrace a Kvadratura (matematika).

Gaussova kvadratura, přesněji n-bodová Gaussova metoda numerické integrace, pojmenovaná po Carlu Friedrichu Gaussovi je metoda numerické integrace, která díky vhodné volbě uzlů x_i a vah w_i pro i = 1, ..., n pro polynomy do stupně 2n − 1 dává přesné výsledky. Moderní formulaci pomocí ortogonálních polynomů vyvinul Carl Gustav Jacob Jacobi v roce 1826. Nejobvyklejším integračním intervalem pro toto pravidlo je $\langle -1,1\rangle$, kdy pravidlo má tvar

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

a poskytuje přesné výsledky pro polynomy do stupně 2n − 1. Toto přesné pravidlo je známé jako Gaussova-Legendreova metoda numerické integrace. Metoda numerické integrace bude přesnou aproximací integrálu uvedeného výše pouze v případě, když funkci f(x) lze na intervalu $\langle -1,1\rangle$ dobře aproximovat polynomem stupně nejvýše 2n − 1.

Gaussova-Legendreova metoda numerické integrace se zpravidla nepoužívá pro integraci funkcí se singularitami v koncových bodech. V takovém případě se snažíme integrand zapsat jako

$f(x)=\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta }g(x),\quad \alpha ,\beta >-1,$

kde g(x) je funkce, kterou lze dobře aproximovat polynomem nízkého stupně. Pak použití alternativních uzlů $x_{i}'$ a vah $w_{i}'$ obvykle dává přesnější metody numerické integrace, které jsou známy jako Gaussova-Jacobiho kvadraturní pravidla:

$\int _{-1}^{1}f(x)\,dx=\int _{-1}^{1}\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta }g(x)\,dx\approx \sum _{i=1}^{n}w_{i}'g\left(x_{i}'\right).$

Jako váhy se často používají hodnoty ${\textstyle {\frac {1}{\sqrt {1-x^{2}}}}}$ (Čebyševova–Gaussova) a ${\sqrt {1-x^{2}}}$. Můžeme také chtít integrovat na intervalech neomezených z jedné (Gaussova-Laguerrova kvadratura) nebo obou stran (Gaussova–Hermitova kvadratura).

Lze ukázat (viz Press, et al. nebo Stoer a Bulirsch), že kvadraturní uzly x_i jsou kořeny polynomu patřícího do třídy ortogonálních polynomů (třída ortogonální vůči váženému vnitřnímu součinu). To je klíčové pozorování pro výpočet uzlů a vah pro Gaussovu kvadraturu.

Gaussova–Legendreova kvadratura

Grafy Legendrových polynomů (do stupně n = 5)

Pro nejjednodušší případ numerické integrace uvedený výše, kdy je funkce f(x) dobře aproximovaná polynomy na intervalu $\langle -1,1\rangle$ související/příslušnými ortogonální polynomy jsou Legendreovy polynomy, označované P_n(x). Je-li n-tý polynom normalizovaný, aby dával P_n(1) = 1, i-tý Gaussův uzel, x_i, je i-tý kořen polynomu P_n a váhy jsou dané vzorcem

$w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)\left[P'_{n}(x_{i})\right]^{2}}}.$

Integrační body a váhy pro Gaussovu–Legendreovu kvadraturu nižších řádů na intervalu $\langle -1,1\rangle$ jsou uvedeny v následující tabulce; v další části jsou vzorce pro integraci na jiných intervalech. Nějaké nejméně významný metody numerické integrace jsou tabelovaný níže (přes

Počet bodů, n	Body, xi		Váhy, wi
1	0		2
2	$\pm {\frac {1}{\sqrt {3}}}$	±0.57735...	1
3	0		${\frac {8}{9}}$	0.888889...
	$\pm {\sqrt {\frac {3}{5}}}$	±0.774597...	${\frac {5}{9}}$	0.555556...
4	$\pm {\sqrt {{\frac {3}{7}}-{\frac {2}{7}}{\sqrt {\frac {6}{5}}}}}$	±0.339981...	${\frac {18+{\sqrt {30}}}{36}}$	0.652145...
	$\pm {\sqrt {{\frac {3}{7}}+{\frac {2}{7}}{\sqrt {\frac {6}{5}}}}}$	±0.861136...	${\frac {18-{\sqrt {30}}}{36}}$	0.347855...
5	0		${\frac {128}{225}}$	0.568889...
	$\pm {\frac {1}{3}}{\sqrt {5-2{\sqrt {\frac {10}{7}}}}}$	±0.538469...	${\frac {322+13{\sqrt {70}}}{900}}$	0.478629...
	$\pm {\frac {1}{3}}{\sqrt {5+2{\sqrt {\frac {10}{7}}}}}$	±0.90618...	${\frac {322-13{\sqrt {70}}}{900}}$	0.236927...

Změna intervalu

Pro výpočet integrálu na obecném intervalu $\langle a,b\rangle$ je třeba před použitím Gaussovy metody numerické integrace úlohu převést na integrál na intervalu $\langle -1,1\rangle$. Tuto změnu intervalu lze provést následujícím způsobem:

$\int _{a}^{b}f(x)\,dx=\int _{-1}^{1}f\left({\frac {b-a}{2}}\xi +{\frac {a+b}{2}}\right)\,{\frac {dx}{d\xi }}d\xi$

kde ${\frac {dx}{d\xi }}={\frac {b-a}{2}}$

Následné použití $n$-bodového Gaussova kvadraturního $(\xi ,w)$ pravidla pak dává následující aproximaci:

$\int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}\xi _{i}+{\frac {a+b}{2}}\right).$

Příklad dvoubodového Gaussova kvadraturního pravidla

Použijeme dvoubodovou Gaussovu metodu numerické integrace pro výpočet přibližné vzdálenosti v metrech uražené raketou v čase $t=8\mathrm {s}$ až $t=30\mathrm {s} ,$ je-li dáno $x=\int _{8}^{30}{\left(2000\ln \left[{\frac {140000}{140000-2100t}}\right]-9.8t\right){dt}}$

Změníme meze integračního intervalu, abychom mohli použít váhy a x-ové souřadnice uvedené v Tabulce 1. Také najdeme absolutní hodnotu relativní skutečné chyby. Přesný výsledek je 11061,34m

Řešení:

Změna mezí integračního intervalu z $\left\langle 8,30\right\rangle$ na $\left\langle -1,1\right\rangle$ dává

${\begin{aligned}\int _{8}^{30}{f(t)dt}&={\frac {30-8}{2}}\int _{-1}^{1}{f\left({\frac {30-8}{2}}x+{\frac {30+8}{2}}\right){dx}}\\&=11\int _{-1}^{1}{f\left(11x+19\right){dx}}\end{aligned}}$

Nyní použijeme váhy a hodnoty argumentu z Tabulky 1 pro dvoubodové pravidlo:

• $c_{1}=1.000000000$
• $x_{1}=-0.577350269$
• $c_{2}=1.000000000$
• $x_{2}=0.577350269$

Nyní můžeme použít Gaussův vzorec pro numerickou integraci

${\begin{aligned}11\int _{-1}^{1}{f\left(11x+19\right){dx}}&\approx 11\left[c_{1}f\left(11x_{1}+19\right)+c_{2}f\left(11x_{2}+19\right)\right]\\&=11\left[f\left(11(-0.5773503)+19\right)+f\left(11(0.5773503)+19\right)\right]\\&=11\left[f(12.64915)+f(25.35085)\right]\\&=11\left[(296.8317)+(708.4811)\right]\\&=11058.44\end{aligned}}$

protože

${\begin{aligned}f(12.64915)&=2000\ln \left[{\frac {140000}{140000-2100(12.64915)}}\right]-9.8(12.64915)\\&=296.8317\end{aligned}}$ ${\begin{aligned}f(25.35085)&=2000\ln \left[{\frac {140000}{140000-2100(25.35085)}}\right]-9.8(25.35085)\\&=708.4811\end{aligned}}$

Pokud přesný výsledek je 11061,34m, absolutní hodnota relativní skutečné chyby, $\left|\varepsilon _{t}\right|$ je $\left|\varepsilon _{t}\right|=\left|{\frac {11061.34-11058.44}{11061.34}}\right|\times 100\%=0.0262\%$

Zobecnění

Problém integrace lze vyjádřit poněkud obecnějším způsobem tak, že do integrandu vneseme kladnou váhovou funkci ω a povolíme jiné integrační intervaly než $\langle -1,1\rangle$. Obecně tedy počítáme

$\int _{a}^{b}\omega (x)\,f(x)\,dx$

pro určité a, b, a ω. Pro a = −1, b = 1, a ω(x) = 1 jde o stejný problém, jaký je popsán výše. Jiné volby vedou k jiným integračním pravidlům. Některé z nich jsou tabelovány níže. Čísla rovnic ve sloupci „A & S“ jsou podle Abramowitz a Stegun

Interval	ω(x)	Ortogonální polynomy	A & S	Pro více informací, viz ...
$\langle -1,1\rangle$	1	Legendreovy polynomy	25.4.29	Gaussova–Legendreova kvadratura
(−1, 1)	$\left(1-x\right)^{\alpha }\left(1+x\right)^{\beta },\quad \alpha ,\beta >-1$	Jacobiho polynomy	25.4.33 (β = 0)	Gaussova–Jacobiho kvadratura
(−1, 1)	${\frac {1}{\sqrt {1-x^{2}}}}$	Čebyševovy polynomy (prvního druhu)	25.4.38	Čebyševova–Gaussova kvadratura
$\langle -1,1\rangle$	${\sqrt {1-x^{2}}}$	Čebyševovy polynomy (druhého druhu)	25.4.40	Čebyševova–Gaussova kvadratura
$\langle 0,\infty )$	$e^{-x}\,$	Laguerrovy polynomy	25.4.45	Gaussova–Laguerre kvadratura
$\langle 0,\infty )$	$x^{\alpha }e^{-x},\quad \alpha >-1$	Zobecněné Laguerrovy polynomy		Gaussova–Laguerre kvadratura
$(-\infty ,\infty )$	$e^{-x^{2}}$	Hermitovy polynomy	25.4.46	Gaussova–Hermitova kvadratura

Základní věta

Nechť p_n je netriviální polynom stupně n, takový, že

$\int _{a}^{b}\omega (x)\,x^{k}p_{n}(x)\,dx=0,\quad {\text{pro všechna }}k=0,1,\ldots ,n-1.$

Všimněte si, že toto bude splněno pro všechny ortogonální polynomy výše, protože každý polynom p_n je zkonstruován tak, aby byl ortogonální k ostatním polynomům p_j pro j≠n, a x^k je v rozsahu této množiny.

Pokud vybereme n uzlů x_i tak, aby to byly kořeny polynomu p_n, pak existuje n vah w_i, pro které je integrál vypočítaný Gaussovou kvadraturou přesný pro všechny polynomy h(x) stupně nejvýše 2n − 1. Všechny tyto uzly x_i budou navíc ležet v otevřeném intervalu (a, b).

Pro důkaz první části tohoto tvrzení budeme předpokládat, že h(x) je jakýkoli polynom stupně nejvýše 2n − 1. Pokud jej vydělíme ortogonálním polynomem p_n dostaneme

$h(x)=p_{n}(x)\,q(x)+r(x).$

kde q(x) je podíl stupně nejvýše n − 1 (protože součet jeho stupně a dělitele p_n musí být roven stupni dělitele), a r(x) je zbytek také stupně nejvýše n − 1 (protože stupeň zbytku je vždy menší než stupeň dělitele). Protože p_n je podle předpokladu ortogonální ke všem jednočlenům stupně menšího než n, musí být ortogonální také s podílem q(x). Proto

$\int _{a}^{b}\omega (x)\,h(x)\,dx=\int _{a}^{b}\omega (x)\,{\big (}\,p_{n}(x)q(x)+r(x)\,{\big )}\,dx=\int _{a}^{b}\omega (x)\,r(x)\,dx.$

Protože zbytek r(x) má stupeň nejvýše n − 1, můžeme jej interpolovat přesně pomocí n interpolačních bodů pomocí Lagrangeovy interpolace l_i(x), kde

$l_{i}(x)=\prod _{j\neq i}{\frac {x-x_{j}}{x_{i}-x_{j}}}.$

Máme

$r(x)=\sum _{i=1}^{n}l_{i}(x)\,r(x_{i}).$

Jeho integrál pak bude roven

$\int _{a}^{b}\omega (x)\,r(x)\,dx=\int _{a}^{b}\omega (x)\,\sum _{i=1}^{n}l_{i}(x)\,r(x_{i})\,dx=\sum _{i=1}^{n}\,r(x_{i})\,\int _{a}^{b}\omega (x)\,l_{i}(x)\,dx=\sum _{i=1}^{n}\,r(x_{i})\,w_{i},$

kde w_i, váhy přiřazené uzlu x_i, jsou dány tak, aby byly rovny váženému integrálu l_i(x) (jiné vzorce pro váhy jsou uvedeny níže). Protože všechny x_i jsou kořeny polynomu p_n, podle vzorce pro dělení uvedeného výše

$h(x_{i})=p_{n}(x_{i})\,q(x_{i})+r(x_{i})=r(x_{i}),$

pro všechna i. Tedy nakonec máme

$\int _{a}^{b}\omega (x)\,h(x)\,dx=\int _{a}^{b}\omega (x)\,r(x)\,dx=\sum _{i=1}^{n}w_{i}\,r(x_{i})=\sum _{i=1}^{n}w_{i}\,h(x_{i}).$

To dokazuje, že pro jakýkoli polynom h(x) stupně nejvýše 2n − 1 je integrál dán přesně sumou podle Gaussovy kvadratury.

Pro důkaz druhé části tvrzení použijeme polynom ve tvaru součinu kořenových činitelů p_n. Komplexně sdružené kořeny tvoří kvadratické členy, které jsou na celé reálné ose buď kladné nebo záporné. Členy oidpovídající kořenům mimo interval $\langle a,b\rangle ,$ na tomto intervalu nemění znaménko. Členy odpovídající kořenům x_i z intervalu $\langle a,b\rangle ,$ které jsou liché násobnosti, vícenásobně p_n jedním více faktor vytvořit nový polynom

$p_{n}(x)\,\prod _{i}(x-x_{i}).$

Tento polynom nemůže na intervalu $\langle a,b\rangle$ měnit znaménko, protože všechny jeho kořeny mají sudou násobnost. Tak integrál

$\int _{a}^{b}p_{n}(x)\,\left(\prod _{i}(x-x_{i})\right)\,\omega (x)\,dx\neq 0,$

protože váhová funkce ω(x) je vždy nezáporná. Ale p_n je ortogonální ke všem polynomům stupně nejvýše n-1, takže stupeň součinu

$\prod _{i}(x-x_{i})$

musí být alespoň n. Proto p_n má n různých reálných kořenů v intervalu $\langle a,b\rangle .$

Obecný vzorec pro váhy

Váhy je možné vyjádřit jako

$w_{i}={\frac {a_{n}}{a_{n-1}}}{\frac {\int _{a}^{b}\omega (x)p_{n-1}\left(x\right)^{2}dx}{p'_{n}(x_{i})p_{n-1}(x_{i})}}$

(1)

kde $a_{k}$ je koeficient u $x^{k}$ v $p_{k}(x)$. Pro důkaz toto, si všimneme, že při použití Lagrangeovy interpolace můžeme vyjádřit r(x) pomocí $r(x_{i})$ jako

$r(x)=\sum _{i=1}^{n}r(x_{i})\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}$

protože r(x) má stupeň menší než n a je proto určen hodnotami, kterých nabývá v n různých bodech. Znásobením obou strany ω(x) a jejich zintegrováním $\langle a,b\rangle ,$ dostaneme

$\int _{a}^{b}\omega (x)r(x)dx=\sum _{i=1}^{n}r(x_{i})\int _{a}^{b}\omega (x)\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}dx$

Váhy w_i jsou tedy dány vzorcem

$w_{i}=\int _{a}^{b}\omega (x)\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {x-x_{j}}{x_{i}-x_{j}}}dx$

Tento integrální výraz pro $w_{i}$ lze vyjádřit pomocí ortogonálních polynomů $p_{n}(x)$ a $p_{n-1}(x)$ takto.

Můžeme psát

$\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\left(x-x_{j}\right)={\frac {\prod _{1\leq j\leq n}\left(x-x_{j}\right)}{x-x_{i}}}={\frac {p_{n}(x)}{a_{n}\left(x-x_{i}\right)}}$

kde $a_{n}$ je koeficient u $x^{n}$ v $p_{n}(x)$. V limitě x jdoucí k $x_{i}$ dostaneme pomocí L'Hôpitalova pravidla

$\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\left(x_{i}-x_{j}\right)={\frac {p'_{n}(x_{i})}{a_{n}}}$

Můžeme tedy psát integrál výraz pro váhy jako

$w_{i}={\frac {1}{p'_{n}(x_{i})}}\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx$

(2)

V integrandu zapisování

${\frac {1}{x-x_{i}}}={\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}+\left({\frac {x}{x_{i}}}\right)^{k}{\frac {1}{x-x_{i}}}$

dává

$\int _{a}^{b}\omega (x){\frac {x^{k}p_{n}(x)}{x-x_{i}}}dx=x_{i}^{k}\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx$

poskytnutý $k\leq n$, protože

${\frac {1-\left({\frac {x}{x_{i}}}\right)^{k}}{x-x_{i}}}$

je polynom stupně k − 1 který pak je ortogonální na $p_{n}(x)$. Takže, pokud q(x) je polynom nejvýše n-tého stupně, máme

$\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {1}{q(x_{i})}}\int _{a}^{b}\omega (x){\frac {q(x)p_{n}(x)}{x-x_{i}}}dx$

Můžeme vyhodnotit integrál na pravé straně pro $q(x)=p_{n-1}(x)$ takto. Protože ${\frac {p_{n}(x)}{x-x_{i}}}$ je polynom stupně n − 1, máme

${\frac {p_{n}(x)}{x-x_{i}}}=a_{n}x^{n-1}+s(x)$

kde s(x) je polynom stupně $n-2$. Protože s(x) je ortogonální na $p_{n-1}(x)$ máme

$\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {a_{n}}{p_{n-1}(x_{i})}}\int _{a}^{b}\omega (x)p_{n-1}(x)x^{n-1}dx$

Pak můžeme psát

$x^{n-1}=\left(x^{n-1}-{\frac {p_{n-1}(x)}{a_{n-1}}}\right)+{\frac {p_{n-1}(x)}{a_{n-1}}}$

Člen v závorkách je polynom stupně $n-2$, který je proto ortogonální s $p_{n-1}(x)$. Integrál tedy můžeme napsat jako

$\int _{a}^{b}\omega (x){\frac {p_{n}(x)}{x-x_{i}}}dx={\frac {a_{n}}{a_{n-1}p_{n-1}(x_{i})}}\int _{a}^{b}\omega (x)p_{n-1}(x)^{2}dx$

Podle rovnice (2), váhy získáme dělením toto výrazem $p'_{n}(x_{i})$ a to dává výraz v rovnici (1).

$w_{i}$ mohou být také vyjádřeny členy ortogonálních polynomů $p_{n}(x)$ a nyní $p_{n+1}(x)$. V trojčlenném rekurentním vztahu $p_{n+1}(x_{i})=(a)p_{n}(x_{i})+(b)p_{n-1}(x_{i})$ člen s $p_{n}(x_{i})$ bude mít nulovou hodnotu, tak $p_{n-1}(x_{i})$ v Eq. (1) lze nahradit ${\textstyle {\frac {1}{b}}p_{n+1}\left(x_{i}\right)}$.

Důkaz, že váhy jsou kladné

Uvažujme následující polynom stupně $2n-2$

$f(x)=\prod _{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}{\frac {\left(x-x_{j}\right)^{2}}{\left(x_{i}-x_{j}\right)^{2}}}$

kde, jak je uvedeno výše, x_j jsou kořeny polynomu $p_{n}(x)$. Jasně $f(x_{j})=\delta _{ij}$. Protože stupeň $f(x)$ je méně než $2n-1$, Gaussův vzorec pro numerickou integraci obsahujícím váhy a uzly získaný z $p_{n}(x)$ se použije. Protože $f(x_{j})=0$ pro j není rovný i, máme

$\int _{a}^{b}\omega (x)f(x)dx=\sum _{j=1}^{N}w_{j}f(x_{j})=\sum _{j=1}^{N}\delta _{ij}w_{j}=w_{i}>0.$

Protože jak $\omega (x)$ tak $f(x)$ jsou nezáporné funkce, z toho plyne, že $w_{i}>0$.

Numerické integrace Gaussovou metodou

Existuje mnoho algoritmů pro výpočet uzlů x_i a vah w_i Gaussovy metody numerické integrace. Nejoblíbenější jsou Golubův-Welschův algoritmus vyžadující O(n²) operací, Newtonova metoda pro řešení $p_{n}(x)=0$ pomocí trojčlenného rekurentního vztahu pro vyhodnocování vyžaduje O(n²) operací, a asymptotické vzorce pro velká n vyžadují O(n) operací.

Rekurentní vztah

Pro ortogonální polynomy $p_{r}$ platí, že $(p_{r},p_{s})=0$ pro $r\neq s,$ kde $(\,\,)$ je skalární součin, stupeň $(p_{r})=r$ a koeficient u nejvyššího členu je jedna (tj. jde o monické ortogonální polynomy) vyhovuje rekurentní vztah

$p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)\cdots -a_{r,0}p_{0}(x)$

a skalární součin definovaný

$(f(x),g(x))=\int _{a}^{b}\omega (x)f(x)g(x)dx$

pro $r=0,1,\ldots ,n-1$ kde n je maximální stupeň, který může být nekonečno, a kde ${\textstyle a_{r,s}={\frac {\left(xp_{r},p_{s}\right)}{\left(p_{s},p_{s}\right)}}}$. Nejdříve polynomy definované vztahem rekurentní vztah začínající $p_{0}(x)=1$ má úvodní koeficient jedna a správný stupeň. Je-li dán počátek $p_{0}$, ortogonalitu $p_{r}$ lze dokázat indukcí. Pro $r=s=0$ máme

$(p_{1},p_{0})=(x-a_{0,0})(p_{0},p_{0})=(xp_{0},p_{0})-a_{0,0}(p_{0},p_{0})=(xp_{0},p_{0})-(xp_{0},p_{0})=0.$

Pokud jsou nyní $p_{0},p_{1},\ldots ,p_{r}$ ortogonální, pak také $p_{r+1}$, protože v

$(p_{r+1},p_{s})=(xp_{r},p_{s})-a_{r,r}(p_{r},p_{s})-a_{r,r-1}(p_{r-1},p_{s})\cdots -a_{r,0}(p_{0},p_{s})$

všechny skalární součiny jsou nulové kromě prvního a toho, kde $p_{s}$ je stejný ortogonální polynom. Proto

$(p_{r+1},p_{s})=(xp_{r},p_{s})-a_{r,s}(p_{s},p_{s})=(xp_{r},p_{s})-(xp_{r},p_{s})=0.$

Pokud však pro skalární součin platí $(xf,g)=(f,xg)$ (což je případ Gaussovy kvadratury), rekurentní vztah se zjednoduší na trojčlenný rekurentní vztah: Pro $s<r-1,xp_{s}$ je polynom nejvýše stupně r − 1. Na druhou stranu, $p_{r}$ je ortogonální ke každému polynomu stupně nejvýše r − 1. Proto máme $(xp_{r},p_{s})=(p_{r},xp_{s})=0$ a $a_{r,s}=0$ pro s < r − 1. Rekurentní vztah se pak zjednoduší na

$p_{r+1}(x)=(x-a_{r,r})p_{r}(x)-a_{r,r-1}p_{r-1}(x)$

nebo

$p_{r+1}(x)=(x-a_{r})p_{r}(x)-b_{r}p_{r-1}(x)$

(s konvencí $p_{-1}(x)\equiv 0$) kde

$a_{r}:={\frac {(xp_{r},p_{r})}{(p_{r},p_{r})}},\qquad b_{r}:={\frac {(xp_{r},p_{r-1})}{(p_{r-1},p_{r-1})}}={\frac {(p_{r},p_{r})}{(p_{r-1},p_{r-1})}}$

(poslední protože $(xp_{r},p_{r-1})=(p_{r},xp_{r-1})=(p_{r},p_{r})$, protože $xp_{r-1}$ se liší od $p_{r}$ stupněm menším než r).

Golubův-Welschův algoritmus

Trojčlenný rekurentní vztah lze zapsat v maticovém tvaru $J{\tilde {P}}=x{\tilde {P}}-p_{n}(x)\times \mathbf {e} _{n}$ kde ${\tilde {P}}={\begin{pmatrix}p_{0}(x)&p_{1}(x)&\ldots &p_{n-1}(x)\end{pmatrix}}^{\mathsf {T}}$, $\mathbf {e} _{n}$ je $n$-tý vektor standardní báze, tj. $\mathbf {e} _{n}={\begin{pmatrix}0&\ldots &0&1\end{pmatrix}}^{\mathsf {T}}$, a J je tak zvaná Jacobiho matice:

$\mathbf {J} ={\begin{pmatrix}a_{0}&1&0&\ldots &\ldots &\ldots \\b_{1}&a_{1}&1&0&\ldots &\ldots \\0&b_{2}&a_{2}&1&0&\ldots \\0&\ldots &\ldots &\ldots &\ldots &0\\\ldots &\ldots &0&b_{n-2}&a_{n-2}&1\\\ldots &\ldots &\ldots &0&b_{n-1}&a_{n-1}\end{pmatrix}}$

Nulové body $x_{j}$ polynomů do stupně n, které se používají jako uzly pro Gaussovu kvadraturu, lze nalézt výpočtem vlastních hodnot této třídiagonální matice. Tento postup je znám jako Golubův–Welschův algoritmus.

Pro výpočet vah a uzlů je vhodnější uvažovat symetrickou třídiagonální matici ${\mathcal {J}}$ s prvky

${\begin{aligned}{\mathcal {J}}_{i,i}=J_{i,i}&=a_{i-1}&&i=1,\ldots ,n\\{\mathcal {J}}_{i-1,i}={\mathcal {J}}_{i,i-1}={\sqrt {J_{i,i-1}J_{i-1,i}}}&={\sqrt {b_{i-1}}}&&i=2,\ldots ,n.\end{aligned}}$

J a ${\mathcal {J}}$ jsou podobné matice a proto mají stejná vlastní čísla (uzly). Váhy lze vypočítat z odpovídajících vlastních vektorů: Pokud $\phi ^{(j)}$ je normalizovaný vlastní vektor (tj. vlastní vektor s Eukleidovskou normou rovné jedné) příslušný k vlastní hodnotě x_j, odpovídající váhu lze vypočítat z první složky tohoto vlastního vektoru, jmenovitě:

$w_{j}=\mu _{0}\left(\phi _{1}^{(j)}\right)^{2}$

kde $\mu _{0}$ je integrál váhové funkce

$\mu _{0}=\int _{a}^{b}\omega (x)dx.$

Další detaily lze nalézt v monografii Numerical Methods for Special Function.

Odhady chyby

Pro integrand, který má 2n spojitých derivací, platí následující odhad chyby Gausovy metody numerické integrace:

$\int _{a}^{b}\omega (x)\,f(x)\,dx-\sum _{i=1}^{n}w_{i}\,f(x_{i})={\frac {f^{(2n)}(\xi )}{(2n)!}}\,(p_{n},p_{n})$

pro nějaké ξ v (a, b), kde p_n je monický (tj. úvodní koeficient je 1) ortogonální polynom stupně n a kde

$(f,g)=\int _{a}^{b}\omega (x)f(x)g(x)\,dx.$

V důležitém speciálním případě, kdy ω(x) = 1, je odhad chyby:

${\frac {\left(b-a\right)^{2n+1}\left(n!\right)^{4}}{(2n+1)\left[\left(2n\right)!\right]^{3}}}f^{(2n)}(\xi ),\qquad a<\xi <b.$

Tento odhad chyby může být v praxi problematický, protože odhad derivace řádu 2n může být obtížný, a skutečná chyba může být navíc mnohem menší než jakou dává tento odhad. Alternativním přístupem může být použití dvou Gaussových metod numerické integrace s různými parametry, a odhadnutí chyby jako rozdílu obou výsledků. K tomuto účelu může být užitečná Gaussova–Kronrodova metoda numerické integrace.

Gaussova–Kronrodova pravidla

Podrobnější informace naleznete na stránce Gaussův–Kronrodův kvadraturní vzorec.

Pokud bychom zjemnili rozdělení intervalu $\langle a,b\rangle$, Gaussovy body vyhodnocování na nových podintervalech se nebudou shodovat s předchozími body (až bodu nula pro liché n), a integrand by se tedy musel vyčíslovat v každém bodě. Gaussova–Kronrodova pravidla jsou rozšířením Gaussovy metody numerické integrace generované přidáním n + 1 bodů k n-bodovému pravidlu takovým způsobem, že výsledné pravidlo má řád 2n + 1 a umožňuje vypočítat odhady vyšších řádů s použitím funkčních hodnot z odhadů nižších řádů. Rozdíl mezi hodnotami získanými Gaussovou metodou numerické integrace a jejím Kronrodovým rozšířením se často používá pro odhad chyby aproximace.

Gaussova–Lobattova pravidla

Metoda známá také jako Lobattova kvadratura pojmenovaná po nizozemském matematikovi Rehuelu Lobattovi se podobá Gaussově kvadratuře s následujícími rozdíly:

1. Integračními body jsou také koncové body integračního intervalu.
2. Je přesná pro polynomy do stupně 2n – 3, kde n je počet integračních bodů.

Lobattova kvadratura funkce f(x) na interval $\langle -1,1\rangle$:

$\int _{-1}^{1}{f(x)\,dx}={\frac {2}{n(n-1)}}[f(1)+f(-1)]+\sum _{i=2}^{n-1}{w_{i}f(x_{i})}+R_{n}.$

X-ové souřadnice: x_i je $(i-1)$-tý kořen $P'_{n-1}(x)$, zde $P_{m}(x)$označuje standardní Legendreův polynom m-tého stupeň a apostrof označuje derivaci.

Váhy:

$w_{i}={\frac {2}{n(n-1)\left[P_{n-1}\left(x_{i}\right)\right]^{2}}},\qquad x_{i}\neq \pm 1.$

Zbytek:

$R_{n}={\frac {-n\left(n-1\right)^{3}2^{2n-1}\left[\left(n-2\right)!\right]^{4}}{(2n-1)\left[\left(2n-2\right)!\right]^{3}}}f^{(2n-2)}(\xi ),\qquad -1<\xi <1.$

Některé váhy jsou:

Počet bodů, n	Body, xi	Váhy, wi
$3$	$0$	${\frac {4}{3}}$
	$\pm 1$	${\frac {1}{3}}$
$4$	$\pm {\sqrt {\frac {1}{5}}}$	${\frac {5}{6}}$
	$\pm 1$	${\frac {1}{6}}$
$5$	$0$	${\frac {32}{45}}$
	$\pm {\sqrt {\frac {3}{7}}}$	${\frac {49}{90}}$
	$\pm 1$	${\frac {1}{10}}$
$6$	$\pm {\sqrt {{\frac {1}{3}}-{\frac {2{\sqrt {7}}}{21}}}}$	${\frac {14+{\sqrt {7}}}{30}}$
	$\pm {\sqrt {{\frac {1}{3}}+{\frac {2{\sqrt {7}}}{21}}}}$	${\frac {14-{\sqrt {7}}}{30}}$
	$\pm 1$	${\frac {1}{15}}$
$7$	$0$	${\frac {256}{525}}$
	$\pm {\sqrt {{\frac {5}{11}}-{\frac {2}{11}}{\sqrt {\frac {5}{3}}}}}$	${\frac {124+7{\sqrt {15}}}{350}}$
	$\pm {\sqrt {{\frac {5}{11}}+{\frac {2}{11}}{\sqrt {\frac {5}{3}}}}}$	${\frac {124-7{\sqrt {15}}}{350}}$
	$\pm 1$	${\frac {1}{21}}$

Adaptivní varianta tohoto algoritmus se dvěma vnitřními uzl je implementována v GNU Octave a MATLAB jako quadl a integrate.

Odkazy

Reference

V tomto článku byl použit překlad textu z článku Gaussian quadrature na anglické Wikipedii.

• KABIR1, H.; MATIKOLAEI, S. A. H. H., 2017. Implementing an Accurate Generalized Gaussian Quadrature Solution to Find the Elastic Field in a Homogeneous Anisotropic Media. Journal of the Serbian Society for Computational Mechanics. Roč. 11, čís. 1, s. 11-19. Dostupné online. DOI 10.24874/jsscm.2017.11.01.02. 
• ABRAMOWITZ, Milton; STEGUN, 1983 [June 1964]. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Applied Mathematics Series. Svazek 55. Washington D.C.; New York: United States: Department of Commerce, National Bureau of Standards; Dover Publications. ISBN 978-0-486-61272-0. 
• ANDERSON, Donald G., 1965. Gaussian quadrature formulae for$\int _{0}^{1}-\ln(x)f(x)dx$. Math. Comp.. Roč. 19, čís. 91, s. 477–481. DOI 10.1090/s0025-5718-1965-0178569-1. 
• GOLUB, Gene H.; WELSCH, John H., 1969. Calculation of Gauss Quadrature Rules. Mathematics of Computation. Roč. 23, čís. 106, s. 221–230. DOI 10.1090/S0025-5718-69-99647-1. JSTOR 2004418. 
• GAUTSCHI, Walter, 1968. Construction of Gauss–Christoffel Quadrature Formulas. Math. Comp.. Roč. 22, čís. 102, s. 251–270. DOI 10.1090/S0025-5718-1968-0228171-0. 
• GAUTSCHI, Walter, 1970. On the construction of Gaussian quadrature rules from modified moments. Math. Comp.. Roč. 24, s. 245–260. DOI 10.1090/S0025-5718-1970-0285117-6. 
• PIESSENS, R., 1971. Gaussian quadrature formulas for the numerical integration of Bromwich's integral and the inversion of the laplace transform. J. Eng. Math.. Roč. 5, čís. 1, s. 1–9. DOI 10.1007/BF01535429. Bibcode 1971JEnMa...5....1P. 
• DANLOY, Bernard, 1973. Numerical construction of Gaussian quadrature formulas for$\int _{0}^{1}(-\log x)x^{\alpha }f(x)dx$and$\int _{0}^{\infty }E_{m}(x)f(x)dx$. Math. Comp.. Roč. 27, čís. 124, s. 861–869. DOI 10.1090/S0025-5718-1973-0331730-X. 
• KAHANER, David; MOLER, Cleve; NASH, Stephen, 1989. Numerical Methods and Software. [s.l.]: Prentice-Hall. Dostupné online. ISBN 978-0-13-627258-8. 
• SAGAR, Robin P., 1991. A Gaussian quadrature for the calculation of generalized Fermi-Dirac integrals. Comput. Phys. Commun.. Roč. 66, čís. 2–3, s. 271–275. DOI 10.1016/0010-4655(91)90076-W. Bibcode 1991CoPhC..66..271S. 
• YAKIMIW, E., 1996. Accurate computation of weights in classical Gauss-Christoffel quadrature rules. J. Comput. Phys.. Roč. 129, čís. 2, s. 406–430. DOI 10.1006/jcph.1996.0258. Bibcode 1996JCoPh.129..406Y. 
• LAURIE, Dirk P., 1999. Accurate recovery of recursion coefficients from Gaussian quadrature formulas. J. Comput. Appl. Math.. Roč. 112, čís. 1–2, s. 165–180. DOI 10.1016/S0377-0427(99)00228-9. 
• LAURIE, Dirk P., 2001. Computation of Gauss-type quadrature formulas. J. Comput. Appl. Math.. Roč. 127, čís. 1–2, s. 201–217. DOI 10.1016/S0377-0427(00)00506-9. Bibcode 2001JCoAM.127..201L. 
• RIENER, Cordian; SCHWEIGHOFER, Markus, 2018. Optimization approaches to quadrature: New characterizations of Gaussian quadrature on the line and quadrature with few nodes on plane algebraic curves, on the plane and in higher dimensions. Journal of Complexity. Roč. 45, s. 22–54. DOI 10.1016/j.jco.2017.10.002. arXiv 1607.08404. 
• STOER, Josef; BULIRSCH, Roland, 2002. Introduction to Numerical Analysis. 2. vyd. [s.l.]: Springer. ISBN 978-0-387-95452-3. .
• TEMME, Nico M., 2010. NIST Handbook of Mathematical Functions. Příprava vydání in Olver, Frank W. J.; Lozier, Daniel M.; Boisvert, Ronald F.; Clark, Charles W.. [s.l.]: Cambridge University Press. ISBN 978-0-521-19225-5. 
• PRESS, WH; TEUKOLSKY, SA; VETTERLING, WT; FLANNERY, BP, 2007. Numerical Recipes: The Art of Scientific Computing. 2. vyd. New York: Cambridge University Press. ISBN 978-0-521-88068-8. Kapitola Section 4.6. Gaussian Quadratures and Orthogonal Polynomials. 
• GIL, Amparo; SEGURA, Javier; TEMME, Nico M., 2007. Numerical Methods for Special Functions. [s.l.]: SIAM. ISBN 978-0-89871-634-4. Kapitola §5.3: Gauss quadrature. 
• QUARTERONI, Alfio; SACCO, Riccardo; SALERI, Fausto. Numerical Mathematics. New York: Springer-Verlag, 2000. ISBN 0-387-98959-5. S. 422, 425. 
• GAUTSCHI, Walter. A Software Repository for Gaussian Quadratures and Christoffel Functions. [s.l.]: SIAM, 2000. ISBN 978-1-611976-34-2. 

Související články

• Numerická integrace
• Newtonovy–Cotesovy vzorce

Literatura

• SEGETHOVÁ, J. Základy numerické matematiky. Praha: MFF UK, 2002. 

Externí odkazy

• Gauss quadrature formula] Encyclopedia of Mathematics, EMS Press, 2001 [1994]
• ALGLIB obsahuje sadu algoritmů pro numerickou integration (v jazycích C# / C++ / Delphi / Visual Basic / atd.)
• GNU Scientific Library — obsahuje implementaci algoritmů QUADPACK pro programovací jazyk C (viz též GNU Scientific Library)
• From Lobatto Quadrature to the Euler constant e
• Gaussian Quadrature Rule of Integration – Notes, PPT, Matlab, Mathematica, Maple, Mathcad at Holistic Numerical Methods Institute
• Kolarp/Gaussova kvadratura v encyklopedii MathWorld (anglicky)
• Gaussian Quadrature autoři: Chris Maes a Anton Antonov, Wolfram Demonstrations Project.
• Tabulované váhy a uzly zdrojovým kódem pro systém Mathematica, high přesností (16 a 256 desítkových míst) Legendreova-Gaussovy kvadraturní váhy a uzly, pro n=2 through n=64, zdrojovým kódem pro systém Mathematica.
• Mathematica source code distributed under the GNU LGPL pro generování uzlů vah pro libovolné váhové funkce W(x), integrační intervaly a přesnosti.
• Gaussova kvadratura v Boost.Math, pro libovolnou přesnost a řád aproximace
• Gaussova-Kronrodova kvadratura v Boost.Math
• Uzly a váhy pro Gaussovu kvadraturu