Numerická integrace
Numerická integrace je často řešenou úlohou v aplikované matematice na výpočet číselné hodnoty určitého integrálu, kterou lze chápat jako výpočet obsahu plochy omezené křivkou definovanou zadanou funkcí. Pro výpočet numerických hodnot určitých integrálů existuje široká rodina algoritmů a jako rozšíření se tento termín také někdy používá pro numerické řešení diferenciálních rovnic. Tento článek se zaměřuje na výpočet určitých integrálů. Termín numerická kvadratura (často zkracovaný na kvadratura) je více nebo méně synonymem pro numerickou integraci, zvláště jednorozměrných integrálů. Někteří autoři používají pro numerický výpočet vícerozměrných integrálů termín kubatura;[1] jiní pod pojem kvadratura zahrnují i vícerozměrnou integraci.
Základním problémem při numerické integraci je výpočet přibližné hodnoty určitého integrálu
s požadovanou přesností. Pokud f(x) je hladká funkce integrovaná přes malý počet rozměrů a definiční obor integrace je omezený, existuje mnoho metod pro aproximaci integrálů s požadovanou přesností.
Historie
Termín „numerická integrace“ se poprvé objevuje v roce 1915 v publikaci Davida Gibba „A Course in Interpolation and Numeric Integration for the Mathematical Laborator“.[2]
Kvadratura je historický matematický termín, který znamená výpočet plochy (obsahu). Problémy kvadratury byly jedním z hlavních zdrojů matematické analýzy. Pythagorejská doktrína chápala jako výpočet plochy jako proces geometrické konstrukce čtverce, který má stejnou plochu jako zadaný geometrický útvar. Proto byl tento proces pojmenován kvadratura. Například kvadratura kruhu, Hippokratových měsíčků, kvadratura paraboly. Tyto konstrukce musejí být prováděny pouze pomocí kružítka a pravítka.
Starověcí Babyloňané používali lichoběžníkovou metodu pro integraci pohybu Jupitera po ekliptice.[3]
Pro kvadraturu obdélníka se stranami a a b je nezbytné zkonstruovat čtverec se stranou (geometrický průměr a a b). Pro toto účel je možné použít následující fakt: jestliže nakreslíme kružnici s průměrem rovným součtu a a b, pak výška BH (z bodu jejich spojení na průsečík s kružnicí) se rovná jejich geometrickému průměru. Podobná geometrická konstrukce řeší problém kvadratury rovnoběžníku a trojúhelníku.
Problémy kvadratury obrazců ohraničených křivkami jsou mnohem obtížnější. I když kvadratura Hippokratových měsíčků byla známa, a kvadratura povrchu koule a Archimédova parabolické výseče se stala největším objevem starověké analýzy, až v 19. století bylo dokázáno, že kvadratura kruhu pomocí kružítka a pravítka není možná.
Archimédovi se připisují objevy, že:
- Plocha povrchu koule se rovná čtyřnásobku obsahu hlavní kružnice této koule.
- Plocha parabolické výseče vyseknuté přímkou je 4/3 plochy trojúhelníka vepsaného do této výseče.
Pro důkaz výsledků Archimedes používal tak zvanou exhaustivní metodu, kterou vyvinul Eudoxos z Knidu.
Ve středověké Evropě znamenalo slovo kvadratura výpočet plochy jakoukoli metodou. Častěji byl používán Cavalieriho princip; byla méně formální, ale jednodušší a výkonnější. S její pomocí Galileo Galilei a Gilles de Roberval nalezli plochu oblouku cykloidy, Grégoire de Saint-Vincent zkoumal plochu pod hyperbolou (Opus Geometricum, 1647) a Alphonse Antonio de Sarasa, de Saint-Vincentův žák a komentátor, zmínil vztahy této plochy s logaritmem.
John Wallis tuto metodu algebraizoval: ve své řadě knih Arithmetica Infinitorum (1656) popsal to, co nyní nazýváme určitý integrál a počítal jeho hodnoty. Isaac Barrow a James Gregory učinili další pokrok: kvadratury pro některé algebraické křivky a spirály. Christiaan Huygens úspěšně prováděl kvadraturu některých rotačních těles.
Kvadratura hyperboly podle Saint-Vincenta a de Sarasa poskytla novou, velmi důležitou funkci přirozený logaritmus.
Objev integrálního počtu přinesl univerzální metody pro výpočet plochy. Následkem toho se termín kvadratura stal zastaralým a místo toho obvykle používá obrat „výpočet jednorozměrného určitého integrálu“.
Důvody numerické integrace
Existuje několik důvodů pro provádění numerické integrace:
- Hodnota integrandu (jak se nazývá integrované funkce) f(x) může být známa pouze v určitých bodech. Některé vestavěné systémy a jiné počítačové aplikace potřebují provádět numerickou integraci funkcí, o nichž je známo pouze několik funkčních hodnot, které byly získány například vzorkováním.
- Vzorec pro integrand může být známý, ale může být obtížné nebo nemožné najít primitivní funkci, která by byla v elementárním tvaru. Příkladem takového integrandu je f(x) = exp(−x2), jejíž primitivní funkci nelze zapsat v elementárním tvaru (funkce se používá v Gaussově chybové funkci).
- Může být možné nelézt primitivní funkci symbolicky, ale výpočet numerické aproximace může být snazší než práce s primitivní funkcí; primitivní funkce může být např. vyjádřena nekonečnou řadou nebo součinem nebo jsou pro její vyhodnocení potřebné speciální funkce, které nejsou dostupné.
Metody pro jednorozměrné integrály
Metody numerické integrace můžeme obecně popsat jako kombinaci vyhodnocování integrandu pro získání aproximace hodnoty integrálu. Integrand se vyčísluje v konečné množině bodů nazývaných integrační body a pro aproximaci integrálu se používá vážený součet těchto hodnot. Integrační body a váhy závisejí na konkrétní použité metodě a požadované přesnosti.
Důležitou částí analýzy jakékoli metody numerické integrace je na studium chování aproximační chyby jako funkce počtu vyhodnocování integrandu. Za nejlepší se obvykle považují metody, které mají malou chybu pro malý počet vyhodnocování. Zmenšení počtu vyhodnocování integrandu snižuje počet použitých aritmetických operací, a proto omezuje celkovou zaokrouhlovací chybu. Každé vyhodnocování také trvá nějakou dobu a integrand může být složitý.
Pokud je integrand rozumný (tj. funkce po částech spojitá a s omezenou variací), lze provést numerickou integraci „hrubou silou“ – vyhodnocováním integrandu po velmi malých krocích.
Kvadraturní pravidla používající interpolační funkce
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é.
Nejjednodušší metodou tohoto typu je použít jako interpolační funkci konstantní funkci (polynom stupně nula), která prochází bodem . 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):
Interpolační funkce může být přímka (afinní zobrazení, tj. polynom stupně 1) procházející body a . Tento přístup je základem tak zvané lichoběžníkové metody nebo lichoběžníkového pravidla:
Výše uvedené metody integrace však v praxi používáme tak, že rozdělíme interval na 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
kde podintervaly mají tvar s a Zde používáme podintervaly stejné délky , ale můžeme také používat intervaly proměnné délky .
Interpolace pomocí polynomů se vyčísluje ve stejně vzdálených (ekvidistantních) bodech v intervalu 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 ve středu intervalu
Zobecněný vzorec pro integraci s použitím hodnoty funkce a jejích derivací ve středu intervalu je
nebo
kde označuje -tou derivaci. Například substitucí a
do zobecněného vzorce pro integraci s použitím hodnoty funkce ve středu intervalu získáme rovnice pro arkus tangens:
kde je imaginární jednotka a
Protože pro každé liché bude čitatel integrandu nulový (), zobecněný vzorec pro integraci s použitím hodnoty funkce ve středu intervalu lze přepsat na tvar
Následující příklad kódu pro program Mathematica generuje graf ukazující rozdíl mezi funkcí arkus tangens a její aproximací pro a :
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 definovaná na intervalu má integrál
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 .
Adaptivní algoritmy
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).
Extrapolační metody
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.
Konzervativní (apriorní) odhad chyby
Nechť mají na omezenou první derivaci, tj. Věta o střední hodnotě pro kde dává
pro nějaké v závislosti na .
Pokud budeme integrovat podle od do na obou stranách a vezmeme absolutní hodnoty, dostaneme
Integrál na pravé straně můžeme dále aproximovat použitím absolutní hodnoty v integrandu a nahrazením členu v horní mezí
|
|
(1) |
kde pro aproximaci bylo použito supremum.
Tudíž, jestliže aproximujeme integrál kvadraturním pravidlem , 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
- .
(Všimněte si, že to je přesně ta chyba, kterou jsme spočítali v příkladu pro .) 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ů.
Integrace na neomezených intervalech
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
a pro intervaly omezené z jedné strany můžeme použít
jako možné transformace.
Množné integrály
Všechna doposud uvedená kvadraturní pravidla byla navržena pro výpočet jednorozměrných integrálů. Jedním z přístupů při výpočtu množných (vícerozměrných) integrálů je převést výpočet množného integrálu na postupný výpočet jednorozměrných integrálů podle Fubiniovy věty (tensorové součinové pravidlo). Tento přístup vyžaduje vyhodnocování funkce ve velkém počtu bodů, který roste exponenciálně se zvyšováním počtu dimenzí. Jsou známé tři metody, jak zabránit tak zvanému prokletí dimenzionality.
Velké množství přídavných technik pro vytváření vícerozměrných integračních pravidel pro množství vážených funkcí je popsáno v monografii od A. H. Strouda.[5]
Monte Carlo
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
Ří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
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]
Spojitost s řešením diferenciálních rovnic
Problém výpočtu integrálu
lze aplikací první části základní věty integrálního počtu převést na počáteční úlohu řešení obyčejné diferenciální rovnice. Zderivováním obou stran výše uvedeného vzorce podle x dostaneme, že pro funkci F platí
Proto lze pro numerickou integraci použít numerické metody vyvinuté pro řešení obyčejných diferenciálních rovnic, jako například Rungeova–Kuttova metoda. Standardní Rungeova–Kuttova metoda čtvrtého řádu aplikovaná na diferenciální rovnici dává Simpsonovo pravidlo uvedené výše.
Díky tomu, že diferenciální rovnice F ' (x) = ƒ(x), která se řeší při integrování, má speciální tvar – pravá strana obsahuje pouze nezávislou proměnnou (zde x) a žádné závislé proměnné (zde F) – teorie i algoritmy se výrazně zjednodušují. Proto je problém výpočtu integrálů vhodné zkoumat samostatně.
Odkazy
Reference
V tomto článku byl použit překlad textu z článku Numerical integration na anglické Wikipedii.
- Cubature [online]. MathWorld [cit. 2020-02-20]. Dostupné online.
- Earliest Known Uses of Some of the Words of Mathematics (Q) [online]. [cit. 2018-03-31]. Dostupné online.
- OSSENDRIJVER, Mathieu. Ancient Babylonian astronomers calculated Jupiter's position from the area under a time-velocity graph. Science. 2016-01-29, roč. 351, čís. 6272, s. 482–484. DOI 10.1126/science.aad8085. PMID 26823423. Bibcode 2016Sci...351..482O.
- LEADER, Jeffery J. Numerical Analysis and Scientific Computation. [s.l.]: Addison Wesley, 2004. Dostupné online. ISBN 978-0-201-73499-7.
- STROUD, A. H. Approximate Calculation of Multiple Integrals. Cliffs, NJ: Prentice-Hall Inc., 1971. Dostupné online.
- BRIOL, François-Xavier; OATES, Chris J.; GIROLAMI, Mark; OSBORNE, Michael A. Frank-Wolfe Bayesian Quadrature: Probabilistic Integration with Theoretical Guarantees. [s.l.]: arXiv.org, 2015-06-08. Dostupné online.
- 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,
Související články
- 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
Externí odkazy
- 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.