Numerická integrace

Rozsáhlou třídu kvadraturních pravidel lze odvodit zkonstruováním takové interpolační funkce, kterou lze snadno integrovat. Interpolačními funkcemi bývají obvykle polynomy. Protože však polynomy vysokého stupně mají sklon k divokým kmitům, v praxi se používají pouze polynomy nízkého stupně, obvykle lineární a kvadratické.

Ilustrace obdélníkové metody.

Nejjednodušší metodou tohoto typu je použít jako interpolační funkci konstantní funkci (polynom stupně nula), která prochází bodem ${\textstyle \left({\frac {a+b}{2}},f\left({\frac {a+b}{2}}\right)\right)}$. Tento princip je základem jedné z variant obdélníkové metody nebo obdélníkového pravidla (s použitím hodnoty funkce ve středu intervalu):

${\displaystyle \int _{a}^{b}f(x)\,dx\approx (b-a)f\left({\frac {a+b}{2}}\right).}$

Ilustrace lichoběžníkové metody.

Interpolační funkce může být přímka (afinní zobrazení, tj. polynom stupně 1) procházející body ${\displaystyle \left(a,f(a)\right)}$ a ${\displaystyle \left(b,f(b)\right)}$. Tento přístup je základem tak zvané lichoběžníkové metody nebo lichoběžníkového pravidla:

${\displaystyle \int _{a}^{b}f(x)\,dx\approx (b-a)\left({\frac {f(a)+f(b)}{2}}\right).}$

Ilustrace Simpsonova pravidla.

Výše uvedené metody integrace však v praxi používáme tak, že rozdělíme interval ${\displaystyle \langle a,b\rangle }$ na ${\displaystyle n}$ podintervalů, spočítáme aproximaci integrálu na každém podintervalu, a sečteme všechny dílčí výsledky, čímž obvykle dosáhneme přesnějšího výsledku. V angličtině se tyto úpravy označují jako composite rule, extended rule nebo iterated rule, v češtině se obvykle obdélníkovou nebo lichoběžníkovou metodou myslí až tyto složené metody. Vzorec pro lichoběžníkovou metodu integrace pak má tvar

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{n}}\left({f(a) \over 2}+\sum _{k=1}^{n-1}\left(f\left(a+k{\frac {b-a}{n}}\right)\right)+{f(b) \over 2}\right),}$

kde podintervaly mají tvar ${\displaystyle [a+kh,a+(k+1)h]\subset \langle a,b\rangle ,}$ s ${\textstyle h={\frac {b-a}{n}}}$ a ${\displaystyle k=0,\ldots ,n-1.}$ Zde používáme podintervaly stejné délky ${\displaystyle h}$, ale můžeme také používat intervaly proměnné délky ${\displaystyle \left(h_{k}\right)_{k}}$.

Interpolace pomocí polynomů se vyčísluje ve stejně vzdálených (ekvidistantních) bodech v intervalu ${\displaystyle \langle a,b\rangle }$ dává Newtonovy–Cotesovy vzorce, z nichž obdélníková metoda a lichoběžníkové metoda jsou příklady. Simpsonova metoda, která používá polynomy 2. řádu, patří také mezi Newtonovy–Cotesovy vzorce.

Kvadraturní pravidla s ekvidistantními body mají velmi příjemnou vlastnost nazývanou nesting. Příslušná metoda pro každý znovu rozdělený interval využívá všechny původní body integrace, takže lze tyto hodnoty integrandu znovu použít.

Pokud dovolíme, aby se velikost interpolačního intervalu měnila, získáme jinou skupinu kvadraturních vzorců, jako například Gaussovo kvadraturní pravidlo. Gaussovo kvadraturní pravidlo je obvykle přesnější než Newtonovo–Cotesovo pravidlo, které vyžaduje stejný počet vyhodnocování funkce, jestliže integrand je hladký (tj. jestliže je dostatečně derivovatelná). Ke kvadraturním metodám s proměnným intervalem patří Clenshawovu–Curtisovu kvadraturu (také nazývaný Fejérova kvadratura) metody, který umožňuje nesting.

Pro Gaussovo kvadraturní pravidlo nelze nesting použít, ale pro příbuzný Gaussův–Kronrodův kvadraturní vzorec ano.

Zobecněný vzorec pro integraci s použitím hodnoty funkce a jejích derivací ve středu intervalu je

${\displaystyle \int _{0}^{1}{f(x)dx}=\sum _{m=1}^{M}{\sum _{n=0}^{\infty }{{\frac {{{\left({-1}\right)}^{n}}+1}{{{\left({2M}\right)}^{n+1}}\left({n+1}\right)!}}{{\left.{{f^{\left(n\right)}}\left(x\right)}\right|}_{x={\frac {m-1/2}{M}}}}}}}$

nebo

${\displaystyle \int _{0}^{1}{f(x)dx}=\lim _{M\to \infty }\sum _{m=1}^{M}{\sum _{n=0}^{N}{{\frac {{{\left({-1}\right)}^{n}}+1}{{{\left({2M}\right)}^{n+1}}\left({n+1}\right)!}}{{\left.{{f^{\left(n\right)}}\left(x\right)}\right|}_{x={\frac {m-1/2}{M}}}}}},}$

kde ${\displaystyle f^{(n)}(x)}$ označuje ${\displaystyle n}$-tou derivaci. Například substitucí ${\displaystyle M=1}$ a

${\displaystyle f(x)={\frac {\theta }{1+\theta ^{2}x^{2}}}}$

do zobecněného vzorce pro integraci s použitím hodnoty funkce ve středu intervalu získáme rovnice pro arkus tangens:

${\displaystyle \operatorname {arctg} (\theta )=i\sum _{n=1}^{\infty }{\frac {1}{2n-1}}\left({\frac {1}{\left(1+2i/\theta \right)^{2n-1}}}-{\frac {1}{\left(1-2i/\theta \right)^{2n-1}}}\right)=2\sum _{n=1}^{\infty }{{\frac {1}{2n-1}}{\frac {{{a}_{n}}\left(\theta \right)}{a_{n}^{2}\left(\theta \right)+b_{n}^{2}\left(\theta \right)}}},}$

kde ${\displaystyle i={\sqrt {-1}}}$ je imaginární jednotka a

${\displaystyle {\begin{aligned}a_{1}(\theta )&={\frac {2}{\theta }},\\b_{1}(\theta )&=1,\\a_{n}(\theta )&=\left(1-{\frac {4}{\theta ^{2}}}\right)\,a_{n-1}(\theta )+{\frac {4}{\theta }}\,b_{n-1}(\theta ),\\b_{n}(\theta )&=\left(1-{\frac {4}{\theta ^{2}}}\right)\,b_{n-1}(\theta )-{\frac {4}{\theta }}\,a_{n-1}(\theta ).\end{aligned}}}$

Protože pro každé liché ${\displaystyle n}$ bude čitatel integrandu nulový (${\displaystyle (-1)^{n}+1=0}$), zobecněný vzorec pro integraci s použitím hodnoty funkce ve středu intervalu lze přepsat na tvar

${\displaystyle \int _{0}^{1}{f(x)dx}=2\sum _{m=1}^{M}{\sum _{n=0}^{\infty }{{\frac {1}{{{\left({2M}\right)}^{2n+1}}\left({2n+1}\right)!}}{{\left.{{f^{(2n)}}(x)}\right|}_{x={\frac {m-1/2}{M}}}}}}\,\,.}$

Následující příklad kódu pro program Mathematica generuje graf ukazující rozdíl mezi funkcí arkus tangens a její aproximací pro ${\displaystyle M=5}$ a ${\displaystyle N=10}$:

f[theta_, x_] := theta/(1 + theta^2*x^2);

aTan[theta_, M_, nMax_] :=
    2*Sum[(Function[x, Evaluate[D[f[theta, x], {x, 2*n}]]][(m - 1/2)/
        M])/((2*n + 1)!*(2*M)^(2*n + 1)), {m, 1, M}, {n, 0, nMax}];

Plot[{ArcTan[theta] - aTan[theta, 5, 10]}, {theta, -Pi, Pi},
 PlotRange -> All

Funkce ${\displaystyle g(t)}$ definovaná na intervalu ${\displaystyle (a,b)}$ má integrál

${\displaystyle \int _{a}^{b}{g(t)dt}=\int _{0}^{b-a}{g(\tau +a)d\tau }=(b-a)\int _{0}^{1}{g((b-a)x+a)dx}.}$

Díky tomu můžeme použít zobecněný vzorec pro numerický výpočet integrálu využívající hodnot funkce ve středech intervalů, který je uveden výše, pokud budeme brát, že ${\displaystyle f(x)=(b-a)\,g((b-a)x+a)}$.

Pokud f(x) nemá dost derivací ve všech integračních bodech, nebo pokud jsou tyto derivace příliš velké, pak Gaussovo kvadraturní pravidlo často nedostačuje. V tomto případě bude fungovat lépe algoritmus podobný následujícímu:

def calculate_definite_integral_of_f(f, initial_step_size):
    """
    Tento algoritmus počítá určitý integrál funkce
    od 0 do 1, adaptivně, použitím menších kroků blízko
    problematických bodů.
    """
    x = 0.0
    h = initial_step_size
    accumulator = 0.0
    while x < 1.0:
        if x + h > 1.0:
            h = 1.0 - x  # Na konci jednotkového intervalu nastavit poslední krok na konec v 1.
        if error_too_big_in_quadrature_of_f_over_range(f, [x, x + h]):
            h = make_h_smaller(h)
        else:
            accumulator += quadrature_of_f_over_range(f, [x, x + h])
            x += h
            if error_too_small_in_quadrature_of_over_range(f, [x, x + h]):
                h = make_h_larger(h)  # Neplýtvat časem pro malé kroky.
    return accumulator

Některé detaily algoritmu je nutné pečlivě zvážit. V mnoha případech není zjevný odhad chyby kvadratury funkce f(x) na určitém intervalu. Pak je oblíbeným řešením používat dva různé algoritmy kvadratury a jejich rozdíl použít jako odhad chyby kvadratury. Jiným problémem je rozhodování, co znamená „příliš velká“ nebo „velmi malá“. Lokální kritérium pro „příliš velká“ je, že kvadraturní chyba nesmí být větší než t · h kde reálné číslo t je tolerance, kterou chceme nastavit pro globální chybu. Pak opět pokud h už je velmi malé, může být zbytečné je dále zmenšovat i kdyby kvadraturní chyba byla zjevně velká. Globálním kritériem je, aby součet chyb na všech intervalech byl menší než t. Tento typ analýzy chyb se obvykle nazývá „aposteriorní“, protože chybu počítáme po výpočtu integrálu.

Heuristiky pro adaptivní kvadraturu jsou diskutovaný v Forsythe et al. (Section 5.4).

Přesnost kvadraturního pravidla Newtonova-Cotesova typu je obecně funkcí počtu integračních bodů. Výsledek je obvykle přesnější pro více integračních bodů, nebo, ekvivalentně, když se zmenšuje velikost integračních intervalů. Je přirozené se ptát, jaký bude výsledek, pokud by se velikost kroku blížila k nule. Na to lze odpovědět extrapolací výsledku ze dvou nebo více nenulových velikostí kroku, pomocí metod zrychlení konvergence řady, jako například Richardsonova extrapolace. Extrapolační funkce může být polynom nebo racionální funkce. Extrapolační metody podrobněji popisuje Stoer a Bulirsch (Section 3.4) a jsou implementovány v mnoha funkcích v knihovně QUADPACK.

Nechť ${\displaystyle f}$ mají na ${\displaystyle \langle a,b\rangle }$ omezenou první derivaci, tj. ${\displaystyle f\in C^{1}(\langle a,b\rangle ).}$ Věta o střední hodnotě pro ${\displaystyle f,}$ kde ${\displaystyle x\in \langle a,b),}$ dává

${\displaystyle (x-a)f'(\xi _{x})=f(x)-f(a),}$

pro nějaké ${\displaystyle \xi _{x}\in (a,x\rangle }$ v závislosti na ${\displaystyle x}$.

Pokud budeme integrovat podle ${\displaystyle x}$ od ${\displaystyle a}$ do ${\displaystyle b}$ na obou stranách a vezmeme absolutní hodnoty, dostaneme

${\displaystyle \left|\int _{a}^{b}f(x)\,dx-(b-a)f(a)\right|=\left|\int _{a}^{b}(x-a)f'(\xi _{x})\,dx\right|.}$

Integrál na pravé straně můžeme dále aproximovat použitím absolutní hodnoty v integrandu a nahrazením členu v ${\displaystyle f'}$ horní mezí

kde pro aproximaci bylo použito supremum.

Tudíž, jestliže aproximujeme integrál ${\textstyle \int _{a}^{b}f(x)\,dx}$ kvadraturním pravidlem ${\displaystyle (b-a)f(a)}$, chyba nebude větší než pravá strana vzorce (1). Tento vztah můžeme upravit pro analýzu chyb Riemannova součtu (*), což nám dá horní mez chybového členu této konkrétní aproximace

${\displaystyle {n^{-1} \over 2}\sup _{0\leq x\leq 1}\left|f'(x)\right|}$.

(Všimněte si, že to je přesně ta chyba, kterou jsme spočítali v příkladu pro ${\displaystyle f(x)=x}$.) Použitím derivace vyšších řádů a úpravou kvadratury můžeme provést podobnou analýzu chyb pomocí Taylorovy řady (pomocí částečných součtů a zbytkového členu) funkce f. Tato analýza chyb dává při známých derivacích funkce f přesnou horní mez chyby.

Tuto metodu integrace lze kombinovat s intervalovou aritmetikou pro vytváření počítačových důkazů a ověřených výpočtů.

Existuje několik metod přibližné integrace na neomezených intervalech. Standardní technika používá speciálně odvozená kvadraturní pravidla, jako například Gaussova-Hermitova kvadratura pro integrály na celé reálné ose a Gaussova-Laguerreova kvadratura pro integrály na kladných reálných číslech.[4] Je také možné používat metody Monte Carlo nebo transformaci proměnné na konečný interval; pro výpočet integrálu na celé reálné ose můžeme například použít

${\displaystyle \int _{-\infty }^{\infty }f(x)\,dx=\int _{-1}^{+1}f\left({\frac {t}{1-t^{2}}}\right){\frac {1+t^{2}}{(1-t^{2})^{2}}}\,dt,}$

a pro intervaly omezené z jedné strany můžeme použít

${\displaystyle {\begin{aligned}\int _{a}^{\infty }f(x)\,dx&=\int _{0}^{1}f\left(a+{\frac {t}{1-t}}\right){\frac {dt}{(1-t)^{2}}},\\\int _{-\infty }^{a}f(x)\,dx&=\int _{0}^{1}f\left(a-{\frac {1-t}{t}}\right){\frac {dt}{t^{2}}},\end{aligned}}}$

jako možné transformace.

Na vícerozměrné integrály lze snadno aplikovat Metody Monte Carlo a Kvazi-Monte Carlo. Mohou pro stejný počet vyhodnocení funkce dávat větší přesnost než opakovaná integrace pomocí jednorozměrných metod.

Velkou třídou užitečných metod Monte Carlo integrování jsou algoritmy Markov chain Monte Carlo (MCMC), ke kterým patří Metropolisův–Hastingsův algoritmus a Gibbsovo vzorkování.

Řídké mřížky vyvinul Sergej Abramovič Smoljak pro kvadraturu funkcí vysoké dimenze. Metoda využívá jednorozměrné kvadraturní pravidlo a provádí rafinované kombinace jednorozměrných výsledků. Zatímco tensorové součinové pravidlo zaručuje, že pokud byly váhy kvadraturních bodů kladné, budou kladné i váhy všech bodů kubatury, Smoljakovo pravidlo toto nezaručuje.

Bayesovská kvadratura je statistický přístup k numerickému problému výpočtu integrálů a spadá do oblasti pravděpodobnostních numerických algoritmů. Může poskytovat úplné zpracovávání nejistoty řešení integrálu vyjádřeného jako aposteriorní rozptyl Gaussovského procesu. Poskytuje velmi rychlou konvergenci, která může být až exponenciální funkcí počtu kvadraturních bodů n.[6]

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

• Philip J. Davis a Philip Rabinowitz, Methods of Numerical Integration.
• George E. Forsythe, Michael A. Malcolm a Cleve B. Moler, Computer Methods for Mathematical Computation. Englewood Cliffs, NJ: Prentice-Hall, 1977. (Viz Chapter 5.)
• PRESS, William H.; TEUKOLSKY, Saul A.; VETTERLING, W.T.; FLANNERY, B.P. Numerical Recipes: The Art of Scientific Computing. 3. vyd. New York: Cambridge University Press, 2007. ISBN 978-0-521-88068-8. Kapitola 4. Integrace Funkce.
• Josef Stoer a Roland Bulirsch, Introduction to Numerical Analysis. New York: Springer-Verlag, 1980. (Viz Chapter 3.)
• Carl Benjamin Boyer, A History of Mathematics, 2. ed. rev. by Uta C. Merzbach, New York: Wiley, 1989 isbn = 0-471-09763-2 (1991 pbk ed. isbn = 0-471-54397-7).
• Howard Eves, An Introduction to the History of Mathematics, Saunders, 1990, isbn = 0-03-029558-0,

• Numerické řešení obyčejných diferenciálních rovnic
• Truncation chyba (numerická integrace)
• Clenshaw–Curtis kvadratura
• Gaussova-Kronrodova kvadratura
• Riemannova suma nebo Riemannův integrál
• Lichoběžníková metoda
• Rombergova metoda
• Tanh-sinh kvadratura

• Integration: Background, Simulations, etc. na webu Holistic Numerical Methods Institute
• Lobatto Quadrature na webu Wolfram Mathworld
• Lobatto quadrature formula na webu Encyclopedia of Mathematics
• Implementace mnoha integračních algoritmů na svobodném webu Tracker Component Library.