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; 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 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“.
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.
náhled|vlevo|220px|Starověká metoda konstrukce geometrického průměru 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.
náhled|200px|Plocha parabolické výseče 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é nalé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é.
vpravo|náhled|300px|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 . 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):
vpravo|náhled|300px|Ilustrace lichoběžníkové metody. 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:
vpravo|náhled|300px|Ilustrace Simpsonova 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 , nebo , 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 .
Vzhledem k tomu, že koncový bod jednoho podintervalu lichoběžníkové metody je počátečním bodem dalšího podintervalu a pro každý podinterval se do průměru počítá polovina výšky tohoto bodu, je tím pádem plná výška tohoto bodu v součtu započtena celá a tak je tato metoda v podstatě podobně přesná jako obdélníková metoda. Liší se pouze ve způsobu započtení krajních bodů celého intervalu, což při rozdělení na vyšší počet podintervalů je v podstatě zanedbatelný rozdíl a stačí použít obdélníkovou metodu.
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 . 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