Gaussova kvadratura

Porovnání dvoubodové Gaussovy kvadratury s lichoběžníkovou metodou. Modrá křivka je graf funkce $y(x),$ jejíž určitý integrál na intervalu $⟨-1,1⟩$ se má vypočítat.
Lichoběžníková metoda aproximuje funkci $y(x)$ oranžovou čárkovanou úsečkou „Trap“, 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 Šablona:Math, 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 integrační interval na mnoho podintervalů, funkci $y(x)$ nahradit v každém podintervalu lineární funkcí a počítat mnoho funkčních hodnot funkce $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 „Gauss 2“) aproximuje funkci $y(x)$ mnohem lépe. Výsledek integrace dvoubodovou Gaussovou metodou je $y(-\sqrt{\frac{1}{3}})+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 $y(x)$ je ve skutečnosti polynom $y(x)=7x^3-8x^2-3x+3$ stupně 3, který je dvoubodovou Gaussovou kvadraturou aproximován přesně.

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 – viz články 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ů Šablona:Mvar a vah Šablona:Mvar pro Šablona:Math dává pro polynomy do stupně Šablona:Math přesné výsledky. Moderní formulaci, v níž se používají ortogonální polynomy vyvinul Carl Gustav Jacob Jacobi v roce 1826. Nejobvyklejším integračním intervalem pro tuto metodu je $⟨-1,1⟩$, na kterém se pro výpočet určitého integrálu používá vzorec

$\underset{-1}{\overset{1}{\int }}f(x)dx\approx \sum _{i=1}^n{w}_{i}f(x_i),$

který poskytuje přesné výsledky pro polynomy do stupně Šablona:Math. Tato numerická metoda se nazývá Gaussova-Legendreova metoda numerické integrace. Metoda numerické integrace bude přesnou aproximací integrálu uvedeného výše pouze v případě, když funkci Šablona:Math lze na intervalu $⟨-1,1⟩$ dobře aproximovat polynomem stupně nejvýše Šablona:Math.

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)={(1-x)}^{\alpha }{(1+x)}^{\beta }g(x),\alpha ,\beta >-1,$

kde Šablona:Math je funkce, kterou lze dobře aproximovat polynomem nízkého stupně. Pak použití alternativních uzlů ${x_i}^{\prime }$ a vah ${w_i}^{\prime }$ obvykle dává přesnější metody numerické integrace, které jsou známy jako Gaussova-Jacobiho kvadraturní pravidla:

$\underset{-1}{\overset{1}{\int }}f(x)dx=\underset{-1}{\overset{1}{\int }}{(1-x)}^{\alpha }{(1+x)}^{\beta }g(x)dx\approx \sum _{i=1}^n{w_i}^{\prime }g\left({x_i}^{\prime }\right).$

Jako váhy se často používají hodnoty $\frac{1}{\sqrt{1-x^2}}$ (Čebyševova–Gaussova kvadratura) 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 Šablona:Mvar 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ě Šablona:Math

Pro nejjednodušší případ numerické integrace uvedený výše je funkce Šablona:Math aproximována polynomy na intervalu $⟨-1,1⟩,$ a přidružené ortogonální polynomy jsou Legendreovy polynomy označované Šablona:Math. Zde je Šablona:Mvar-tý polynom normalizovaný tak, aby Šablona:Math, Šablona:Mvar-tý Gaussův uzel Šablona:Mvar je Šablona:Mvar-tým kořenem polynomu Šablona:Mvar a váhy jsou dané vzorcemŠablona:Sfn

$w_i=\frac{2}{(1-x_i^2){[P{\text{'}}_n(x_i)]}^2}.$

Integrační body a váhy pro Gaussovu–Legendreovu kvadraturu několika nižších řádů na intervalu $⟨-1,1⟩$ jsou uvedeny v následující tabulce:

Počet bodů, n	Body, Šablona:Mvar		Váhy, Šablona:Mvar
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...

Jiné integrační intervaly

Pro výpočet integrálu na obecném intervalu $⟨a,b⟩$ je třeba před použitím Gaussovy metody numerické integrace úlohu převést na integrál na intervalu $⟨-1,1⟩$ pomocí vzorce

$\underset{a}{\overset{b}{\int }}f(x)dx=\underset{-1}{\overset{1}{\int }}f(\frac{b-a}{2}\xi +\frac{a+b}{2})\frac{dx}{d\xi }d\xi$

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

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

$\underset{a}{\overset{b}{\int }}f(x)dx\approx \frac{b-a}{2}\sum _{i=1}^n{w}_{i}f(\frac{b-a}{2}{\xi }_i+\frac{a+b}{2}).$

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

Použijeme dvoubodovou Gaussovu metodu numerické integrace pro výpočet vzdálenosti, kterou uletí raketa v čase $t=8\mathrm{s}$ až $t=30\mathrm{s},$ s použitím vzorce $x=\underset{8}{\overset{30}{\int }}(2000\ln \left[\frac{140000}{140000-2100t}\right]-9.8t)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,34 metrů

Řešení:

provedeme změnu mezí integračního intervalu z $⟨8,30⟩$ na $⟨-1,1⟩$:

$\begin{array}{cc}\underset{8}{\overset{30}{\int }}f(t)dt& =\frac{30-8}{2}\underset{-1}{\overset{1}{\int }}f(\frac{30-8}{2}x+\frac{30+8}{2})dx\\ & =11\underset{-1}{\overset{1}{\int }}f(11x+19)dx\end{array}$

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$

Načež můžeme použít Gaussův vzorec pro numerickou integraci

$\begin{array}{cc}11\underset{-1}{\overset{1}{\int }}f(11x+19)dx& \approx 11[c_{1}f(11x_1+19)+c_{2}f(11x_2+19)]\\ & =11[f(11(-0.5773503)+19)+f(11(0.5773503)+19)]\\ & =11[f(12.64915)+f(25.35085)]\\ & =11[(296.8317)+(708.4811)]\\ & =11058.44\end{array}$

dostáváme

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

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

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 Šablona:Mvar a povolíme jiné integrační intervaly než $⟨-1,1⟩$. Obecně tedy počítáme

$\underset{a}{\overset{b}{\int }}\omega (x)f(x)dx$

pro určité Šablona:Mvar, Šablona:Mvar a Šablona:Mvar. Pro Šablona:Math, Šablona:Math a Šablona:Math jde o stejný problém, jaký je popsán výše. Jiné volby vedou k jiným integračním pravidlům, z nichž některá jsou tabelována níže. Čísla rovnic ve sloupci „A & S“ jsou podle Abramowitz a Stegun

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

Základní věta

Nechť Šablona:Mvar je netriviální polynom stupně Šablona:Mvar, takový, že

$\underset{a}{\overset{b}{\int }}\omega (x)x^k{p}_n(x)dx=0,\text{pro všechna }k=0,1,\dots ,n-1.$

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

Pokud vybereme Šablona:Mvar uzlů Šablona:Mvar tak, aby to byly kořeny polynomu Šablona:Mvar, pak existuje Šablona:Mvar vah Šablona:Mvar, pro které je integrál vypočítaný Gaussovou kvadraturou přesný pro všechny polynomy Šablona:Math stupně nejvýše Šablona:Math. Všechny tyto uzly Šablona:Mvar budou navíc ležet v otevřeném intervalu Šablona:MathŠablona:Sfn.

Pro důkaz první části tohoto tvrzení budeme předpokládat, že Šablona:Math je jakýkoli polynom stupně nejvýše Šablona:Math. Pokud jej vydělíme ortogonálním polynomem Šablona:Mvar dostaneme

$h(x)=p_n(x)q(x)+r(x).$

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

$\underset{a}{\overset{b}{\int }}\omega (x)h(x)dx=\underset{a}{\overset{b}{\int }}\omega (x)(p_n(x)q(x)+r(x))dx=\underset{a}{\overset{b}{\int }}\omega (x)r(x)dx.$

Protože zbytek Šablona:Math má stupeň nejvýše Šablona:Math, můžeme jej přesně interpolovat pomocí Lagrangeovy interpolace Šablona:Math s Šablona:Mvar interpolačními body, kde

$l_i(x)=\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}.$

tedy

$r(x)=\sum _{i=1}^n{l}_i(x)r(x_i).$

Jeho integrál pak bude roven

$\underset{a}{\overset{b}{\int }}\omega (x)r(x)dx=\underset{a}{\overset{b}{\int }}\omega (x)\sum _{i=1}^n{l}_i(x)r(x_i)dx=\sum _{i=1}^{n}r(x_i)\underset{a}{\overset{b}{\int }}\omega (x)l_i(x)dx=\sum _{i=1}^{n}r(x_i)w_i,$

kde váha Šablona:Math přiřazená uzlu Šablona:Math je definována, aby vyrovnala vážený integrál Šablona:Math (další vzorce pro váhy jsou uvedeny níže). Protože však všechny Šablona:Mvar jsou kořeny polynomu Šablona:Mvar, podle vzorce pro dělení uvedeného výše platí

$h(x_i)=p_n(x_i)q(x_i)+r(x_i)=r(x_i),$

pro všechna Šablona:Mvar. Tím dostáváme

$\underset{a}{\overset{b}{\int }}\omega (x)h(x)dx=\underset{a}{\overset{b}{\int }}\omega (x)r(x)dx=\sum _{i=1}^n{w}_{i}r(x_i)=\sum _{i=1}^n{w}_{i}h(x_i).$

Tím jsme dokázali, že pro jakýkoli polynom Šablona:Math stupně nejvýše Šablona:Math je integrál přesně roven hodnotě spočítané Gaussovou kvadraturou.

Pro důkaz druhé části tvrzení použijeme polynom ve tvaru součinu kořenových činitelů Šablona:Math. Komplexně sdružené kořeny tvoří kvadratické členy, které jsou na celé reálné ose buď kladné nebo záporné. Členy odpovídající kořenům ležícím mimo interval $⟨a,b⟩,$ nemění na tomto intervalu znaménko. A členy odpovídající kořenům Šablona:Mvar z intervalu $⟨a,b⟩,$ které jsou liché násobnosti, násobí Šablona:Math jedním dalším členem, který tvoří nový polynom

$p_n(x)\prod _i(x-x_i).$

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

$\underset{a}{\overset{b}{\int }}{p}_n(x)(\prod _i(x-x_i))\omega (x)dx\ne 0,$

protože váhová funkce Šablona:Math je vždy nezáporná. Protože však Šablona:Math je ortogonální ke všem polynomům stupně nejvýše Šablona:Math, bude stupeň součinu

$\prod _i(x-x_i)$

alespoň Šablona:Mvar. Proto Šablona:Math má Šablona:Mvar různých reálných kořenů v intervalu $⟨a,b⟩.$

Obecný vzorec pro váhy

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

kde $a_k$ je koeficient u $x^k$ v $p_k(x)$. Pro důkaz je třeba si všimnout, že při použití Lagrangeovy interpolace můžeme Šablona:Math vyjádřit pomocí $r(x_i)$ jako

$r(x)=\sum _{i=1}^{n}r(x_i)\prod _{\begin{array}{c}1\le j\le n\\ j\ne i\end{array}}\frac{x-x_j}{x_i-x_j}$

protože Šablona:Math má stupeň menší než Šablona:Mvar a je proto určen hodnotami, kterých nabývá v Šablona:Mvar různých bodech. Když znásobíme obě strany Šablona:Math a zintegrujeme je na intervalu $⟨a,b⟩,$ dostaneme

$\underset{a}{\overset{b}{\int }}\omega (x)r(x)dx=\sum _{i=1}^{n}r(x_i)\underset{a}{\overset{b}{\int }}\omega (x)\prod _{\begin{array}{c}1\le j\le n\\ j\ne i\end{array}}\frac{x-x_j}{x_i-x_j}dx$

Váhy Šablona:Mvar jsou tedy dány vzorcem

$w_i=\underset{a}{\overset{b}{\int }}\omega (x)\prod _{\begin{array}{c}1\le j\le n\\ j\ne i\end{array}}\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{array}{c}1\le j\le n\\ j\ne i\end{array}}(x-x_j)=\frac{\prod _{1\le j\le n}(x-x_j)}{x-x_i}=\frac{p_n(x)}{a_n(x-x_i)}$

kde $a_n$ je koeficient u $x^n$ v $p_n(x)$. Pokud se Šablona:Mvar limitně blíží $x_i,$ dostaneme pomocí L'Hôpitalova pravidla

$\prod _{\begin{array}{c}1\le j\le n\\ j\ne i\end{array}}(x_i-x_j)=\frac{p{\text{'}}_n(x_i)}{a_n}$

Můžeme tedy napsat výraz pro váhy s použitím integrálu ve tvaru

Jmenovatel v integrandu můžeme zapsat

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

což dává

$\underset{a}{\overset{b}{\int }}\omega (x)\frac{x^k{p}_n(x)}{x-x_i}dx=x_i^k\underset{a}{\overset{b}{\int }}\omega (x)\frac{p_n(x)}{x-x_i}dx$

pokud $k\le n$, protože

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

je polynom stupně Šablona:Math, který je tedy ortogonální k $p_n(x)$. Pokud tedy Šablona:Math je polynom nejvýše n-tého stupně, můžeme napsat

$\underset{a}{\overset{b}{\int }}\omega (x)\frac{p_n(x)}{x-x_i}dx=\frac{1}{q(x_i)}\underset{a}{\overset{b}{\int }}\omega (x)\frac{q(x)p_n(x)}{x-x_i}dx$

Integrál na pravé straně můžeme pro $q(x)=p_{n-1}(x)$ vyhodnotit takto: Protože $\frac{p_n(x)}{x-x_i}$ je polynom stupně Šablona:Math, máme

$\frac{p_n(x)}{x-x_i}=a_n{x}^{n-1}+s(x)$

kde Šablona:Math je polynom stupně $n-2$. Protože Šablona:Math je ortogonální s $p_{n-1}(x),$ platí

$\underset{a}{\overset{b}{\int }}\omega (x)\frac{p_n(x)}{x-x_i}dx=\frac{a_n}{p_{n-1}(x_i)}\underset{a}{\overset{b}{\int }}\omega (x)p_{n-1}(x)x^{n-1}dx$

Pak můžeme psát

$x^{n-1}=(x^{n-1}-\frac{p_{n-1}(x)}{a_{n-1}})+\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 vyjádřit jako

$\underset{a}{\overset{b}{\int }}\omega (x)\frac{p_n(x)}{x-x_i}dx=\frac{a_n}{a_{n-1}{p}_{n-1}(x_i)}\underset{a}{\overset{b}{\int }}\omega (x)p_{n-1}(x{)}^{2}dx$

Podle rovnice (Šablona:Odkaz na vzorec) se váhy získají vydělením tohoto výrazu výrazem $p{\text{'}}_n(x_i)$, což dává výraz v rovnici (Šablona:Odkaz na vzorec).

Váhy $w_i$ mohou být také vyjádřeny pomocí ortogonálních polynomů $p_n(x)$ a $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)$ bude mít člen s $p_n(x_i)$ nulovou hodnotu, takže $p_{n-1}(x_i)$ v (Šablona:Odkaz na vzorec) lze nahradit $\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{array}{c}1\le j\le n\\ j\ne i\end{array}}\frac{{(x-x_j)}^2}{{(x_i-x_j)}^2}$

kde, jak je uvedeno výše, Šablona:Mvar jsou kořeny polynomu $p_n(x)$. Je zřejmé, že $f(x_j)={\delta }_{ij}$. Protože stupeň $f(x)$ je menší než $2n-1$, platí Gaussův vzorec pro numerickou integraci obsahující váhy a uzly získané z $p_n(x)$. Protože $f(x_j)=0$, pokud j není rovno i, máme

$\underset{a}{\overset{b}{\int }}\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. Odtud plyne, že $w_i>0$.

Numerická integrace Gaussovou metodou

Existuje mnoho algoritmů pro výpočet uzlů Šablona:Mvar a vah Šablona:Mvar v Gaussově metodě numerické integrace. Nejoblíbenější jsou Golubův-Welschův algoritmus vyžadující Šablona:Math operací, Newtonova metoda pro řešení $p_n(x)=0$ s pomocí trojčlenného rekurentního vztahu pro vyhodnocování, která vyžaduje Šablona:Math operací, a asymptotické vzorce pro velká n, které vyžadují Šablona:Math operací.

Rekurentní vztah

Ortogonální polynomy $p_r$, pro které platí, že $(p_r,p_s)=0$ pro $r\ne 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) vyhovují rekurentnímu 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ý vzorcem

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

pro $r=0,1,\dots ,n-1$, kde Šablona:Mvar je nejvyšší stupeň, který může jít k nekonečnu, a kde $a_{r,s}=\frac{(x{p}_r,p_s)}{(p_s,p_s)}$. Především polynomy, definované rekurentním vztahem začínajícím $p_0(x)=1$, mají úvodní koeficient jedna a jsou správného stupně. Je-li dán počáteční bod $p_0$, lze ortogonalitu $p_r$ dokázat indukcí. Pro $r=s=0$ máme

$(p_1,p_0)=(x-a_{0,0})(p_0,p_0)=(x{p}_0,p_0)-a_{0,0}(p_0,p_0)=(x{p}_0,p_0)-(x{p}_0,p_0)=0.$

Pokud jsou $p_0,p_1,\dots ,p_r$ ortogonální, pak také $p_{r+1}$, protože ve výrazu

$(p_{r+1},p_s)=(x{p}_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)$

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

$(p_{r+1},p_s)=(x{p}_r,p_s)-a_{r,s}(p_s,p_s)=(x{p}_r,p_s)-(x{p}_r,p_s)=0.$

Pokud však pro skalární součin platí $(xf,g)=(f,xg)$ (což je v případě Gaussovy kvadratury splněno), rekurentní vztah lze zjednodušit na vztah obsahující tři členy: Pro $s<r-1$ je $x{p}_s$ polynom stupně nejvýše Šablona:Math. Na druhou stranu, $p_r$ je ortogonální ke každému polynomu stupně nejvýše Šablona:Math. Proto máme $(x{p}_r,p_s)=(p_r,x{p}_s)=0$ a $a_{r,s}=0$ pro Šablona:Math. 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{(x{p}_r,p_r)}{(p_r,p_r)},b_r:=\frac{(x{p}_r,p_{r-1})}{(p_{r-1},p_{r-1})}=\frac{(p_r,p_r)}{(p_{r-1},p_{r-1})}$

(poslední vztah platí, protože $(x{p}_r,p_{r-1})=(p_r,x{p}_{r-1})=(p_r,p_r)$, protože $x{p}_{r-1}$ se liší od $p_r$ ve stupni menším než Šablona:Mvar).

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 {𝐞}_n$, kde $\tilde{P}={\left(\begin{array}{cccc}{p}_0(x)& p_1(x)& \dots & p_{n-1}(x)\end{array}\right)}^{𝖳}$, ${𝐞}_n$ je $n$-tý vektor standardní báze, tj. ${𝐞}_n={\left(\begin{array}{cccc}0& \dots & 0& 1\end{array}\right)}^{𝖳}$, a Šablona:Mvar je tak zvaná Jacobiho matice:

$𝐉=\left(\begin{array}{cccccc}{a}_0& 1& 0& \dots & \dots & \dots \\ b_1& a_1& 1& 0& \dots & \dots \\ 0& b_2& a_2& 1& 0& \dots \\ 0& \dots & \dots & \dots & \dots & 0\\ \dots & \dots & 0& b_{n-2}& a_{n-2}& 1\\ \dots & \dots & \dots & 0& b_{n-1}& a_{n-1}\end{array}\right)$

Nulové body $x_j$ polynomů do stupně Šablona:Mvar, 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 $𝒥$ s prvky

$\begin{array}{cccc}{𝒥}_{i,i}=J_{i,i}& =a_{i-1}& & i=1,\dots ,n\\ {𝒥}_{i-1,i}={𝒥}_{i,i-1}=\sqrt{J_{i,i-1}{J}_{i-1,i}}& =\sqrt{b_{i-1}}& & i=2,\dots ,n.\end{array}$

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

$w_j={\mu }_0{\left({\varphi }_1^{(j)}\right)}^2$

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

${\mu }_0=\underset{a}{\overset{b}{\int }}\omega (x)dx.$

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

Odhady chyby

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

$\underset{a}{\overset{b}{\int }}\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é Šablona:Mvar z intervalu Šablona:Math, kde Šablona:Mvar je monický (tj. má úvodní koeficient roven Šablona:Math) ortogonální polynom stupně Šablona:Mvar a kde

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

V důležitém speciálním případě, kdy Šablona:Math, je odhad chyby:Šablona:Sfn

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

Tento odhad chyby může být v praxi problematický, protože odhad derivace řádu Šablona:Math může být obtížný, a skutečná chyba může být navíc mnohem menší než jakou dává tento odhad.Šablona:Sfn 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

Pokud bychom zjemnili rozdělení intervalu $⟨a,b⟩$, 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 Šablona:Math bodů k Šablona:Mvar-bodovému pravidlu takovým způsobem, že výsledné pravidlo má řád Šablona:Math 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Šablona:Sfn 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ě Šablona:Math, kde Šablona:Mvar je počet integračních bodůŠablona:Sfn.

Lobattova kvadratura funkce Šablona:Math na intervalu $⟨-1,1⟩$:

$\underset{-1}{\overset{1}{\int }}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: Šablona:Mvar je $(i-1)$-tý kořen $P{\text{'}}_{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},x_i\ne \pm 1.$

Zbytek:

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

Integrační body a váhy pro několik nižších řádů jsou:

Počet bodů, n	Body, Šablona:Mvar	Váhy, Šablona:Mvar
$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 algoritmu se dvěma vnitřními uzly je implementována v GNU Octave a MATLAB jako quadl a integrate.

Odkazy

Reference

• Šablona:Citace periodika
• Šablona:Citace monografie
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace monografie
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace periodika
• Šablona:Citace monografie.
• Šablona:Citace monografie
• Šablona:Citace monografie
• Šablona:Citace monografie
• Šablona:Citace monografie
• Šablona:Citace monografie
• Šablona:Citace periodika

Související články

• Numerická integrace
• Newtonovy–Cotesovy vzorce

Literatura

• Šablona:Citace monografie

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
• Šablona:MathWorld
• 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 Šablona:Wayback