Simpsonovo pravidlo

Simpsonovu metodu lze odvodit aproximací integrandu f (x) (modře) kvadratickou interpolační funkcí P(x) (červeně).

Animace ukazující jak Simpsonova metoda aproximuje funkci parabolou, a jak se zmenšuje chyba při zmenšení velikosti kroku

Animace ukazující, jak se aproximace Simpsonovou metodou zlepšuje při jemnějším dělení intervalu.

Simpsonova pravidla jsou vzorce používané pro numerický výpočet určitých integrálů, které jsou pojmenovány po Thomasu Simpsonovi (1710–1761).

Nejjednodušší Simpsonovo pravidlo nazývané také Simpsonovo 1/3 pravidlo je

$\int _{a}^{b}f(x)\,dx\approx {\tfrac {b-a}{6}}\left[f(a)+4f\left({\tfrac {a+b}{2}}\right)+f(b)\right].$

V němčině a některých dalších jazycích je pojmenované po Janu Keplerovi, který jej odvodil již v roce 1615 pro výpočet objemu vinných sudů (německy Keplersche Fassregel). Přibližná rovnost v pravidle přechází na přesnou rovnost, pokud f je polynom nejvýše třetího stupně.

Pokud se interval $\langle a,b\rangle$ rozdělí na n stejných dílů a na každý se aplikuje Simpsonovo pravidlo, dostaneme složené Simpsonovo pravidlo, u něhož se pro výpočet integrálu počítá vážený průměr hodnot funkce se střídajícími se vahami 4/3 a 2/3.

Simpsonovo druhé pravidlo, také nazývané Simpsonovo 3/8 pravidlo vyžaduje o jedno vyhodnocení funkce v integračním rozsahu více, ale poskytuje přesný výsledek pro polynomy do třetího stupně.

1/3 a 3/8 Simpsonova pravidla jsou dvěma speciálními případy uzavřených Newtonových–Cotesových vzorců.

Při konstrukci lodí a při odhadech stability lodí se používá také Simpsonovo třetí pravidlo, které pro numerickou matematiku nemá žádný speciální význam, viz Simpsonova pravidla (stabilita lodí).

Simpsonovo 1/3 pravidlo

Odvození Simpsonova pravidla

Kvadratická interpolace

Simpsonovo pravidlo lze odvodit tak, že nahradíme integrovanou funkci $f(x)$ kvadratickým polynomem $P(x)$ (tj. parabolou), který v koncových bodech integračního intervalu $a$ a $b$ a v jeho středu $m=(a+b)/2$ nabývá stejných hodnot jako $f(x)$. Pro odvození vzorce tohoto polynomu můžeme použít Lagrangeovu polynomiální interpolaci:

$P(x)=f(a){\tfrac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\tfrac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\tfrac {(x-a)(x-m)}{(b-a)(b-m)}}.$

Pomocí substituční metody můžeme ukázat, že

$\int _{a}^{b}P(x)\,dx={\tfrac {b-a}{6}}\left[f(a)+4f\left({\tfrac {a+b}{2}}\right)+f(b)\right].$

Pokud velikost kroků, ve kterých počítáme hodnotu funkce $f(x),$ vyjádříme jako $h=(b-a)/2,$ dostaneme vzorec, ze kterého je patrné, proč je v názvu této metody numerické integrace zlomek 1/3:

$\int _{a}^{b}P(x)\,dx={\tfrac {h}{3}}\left[f(a)+4f\left({\tfrac {a+b}{2}}\right)+f(b)\right]$

(zobecnění je popsáno níže).

Průměrování obdélníkové a lichoběžníkové metody

Simpsonovo pravidlo lze odvodit také ze dvou jednodušších aproximací: obdélníkové metody (též Riemannův středový součet):

$M=(b-a)f\left({\tfrac {a+b}{2}}\right)$

a lichoběžníkové metody:

$T={\tfrac {1}{2}}(b-a)(f(a)+f(b)).$

Chyby těchto aproximací jsou

${\tfrac {1}{24}}(b-a)^{3}f''(a)+O((b-a)^{4})\quad {\text{a}}\quad -{\tfrac {1}{12}}(b-a)^{3}f''(a)+O((b-a)^{4}),$

kde $O((b-a)^{4})$ označuje člen asymptoticky úměrný $(b-a)^{4}$. Dva členy $O((b-a)^{4})$ si nejsou rovné; pro více detailů viz Landauova notace. Z výše uvedeného vzorce pro chyby obdélníkové a lichoběžníkové metody vyplývá, že pokud použijeme vážený průměr výsledků obou metod a bude chtít, aby měl první chybový člen nulovou hodnotu, musíme použít váhu 2/3 pro výsledek obdélníkové metody a 1/3 pro výsledek lichoběžníkové metody:

${\tfrac {2M+T}{3}}.$

což je právě Simpsonovo pravidlo.

Myšlenku použití váženého průměru můžeme aplikovat na kombinaci výsledků jiných aproximačních metod výpočtu určitého integrálu (například lichoběžníkové metody s dvojnásobným počtem bodů), a odstranit jiný chybový člen. To je princip Rombergovy metody.

Neurčité koeficienty

Třetí odvození Simpsonova pravidla vychází z myšlenky vyjádřit hodnotu integrálu pomocí kombinace hodnot integrované funkce v krajních bodech integračního intervalu a v jeho středu:

${\tfrac {1}{b-a}}\int _{a}^{b}f(x)\,dx\approx \alpha f(a)+\beta f\left({\tfrac {a+b}{2}}\right)+\gamma f(b).$

Koeficienty α, β a γ lze zafixovat tím, že budeme vyžadovat, aby tato aproximace byla přesná pro všechny kvadratické polynomy. Výsledkem je Simpsonovo pravidlo.

Chyba

Chyba při aproximaci integrálu Simpsonovým pravidlem pro $n=2$ je

$-{\tfrac {1}{90}}\left({\tfrac {b-a}{2}}\right)^{5}f^{(4)}(\xi ),$

kde $\xi$ (řecké písmeno ksí) je nějaké číslo z intervalu $\langle a,\ b\rangle .$

Chyba je asymptoticky úměrná $(b-a)^{5}$. Podle výše uvedených odvození by však chyba měla být úměrná $(b-a)^{4}$. Simpsonovo pravidlo však má ve vzorci pro chybu o řád vyšší mocninu, protože body, v nichž se integrand vyčísluje, jsou v intervalu $\langle a,\ b\rangle$ rozložené symetricky.

Když je chybový člen úměrný čtvrté derivaci funkce $f$ v bodě $\xi$, Simpsonovo pravidlo dává přesné výsledky, pokud je funkce $f$ polynomem nejvýše třetího stupně, protože čtvrtá derivace takových polynomů je ve všech bodech nulová.

Pokud existuje druhá derivace $f''$ a v intervalu $(a,\ b)$ je konvexní, platí nerovnost:

$(b-a)f\left({\tfrac {a+b}{2}}\right)+{\tfrac {1}{3}}\left({\tfrac {b-a}{2}}\right)^{3}f''\left({\tfrac {a+b}{2}}\right)\leq \int _{a}^{b}f(x)\,dx\leq {\tfrac {b-a}{6}}\left[f(a)+4f\left({\tfrac {a+b}{2}}\right)+f(b)\right].$

Složené Simpsonovo pravidlo

Pokud interval integrace $\langle a,b\rangle$ je v nějakém smyslu „malý“, pak Simpsonovo pravidlo s $n=2$ podintervaly bude dávat dostatečně přesnou aproximaci hodnoty integrálu. Slovem „malý“ máme na mysli to, že integrovaná funkce je na tomto intervalu $\langle a,b\rangle$ relativně hladká. Pro takové funkce bude hladký kvadratický interpolant, který se používá v Simpsonově pravidle, dávat dobré výsledky.

Často však funkce, kterou chceme integrovat, není na uvažovaném intervalu hladká. Obvykle to znamená, že buď silně osciluje nebo nemá v určitých bodech derivace. V těchto případech může Simpsonovo pravidlo dávat velmi špatné výsledky. Obvyklým způsobem, jak se vypořádat s tímto problémem je rozdělení intervalu $\langle a,b\rangle$ na několik menších podintervalů. Simpsonovo pravidlo se pak aplikuje na každý podinterval a výsledky se sečtou, aby se získala aproximace integrálu na celém intervalu. Tento druh přístupu se nazývá složené Simpsonovo pravidlo.

Předpokládejme, že, že interval $\langle a,b\rangle$ je rozdělují na $n$ podintervalů, kde $n$ je sudé číslo. Složené Simpsonovo pravidlo je pak dáno vzorcem

${\begin{aligned}\int _{a}^{b}f(x)\,dx&\approx {\frac {h}{3}}\sum _{j=1}^{n/2}{\bigg [}f(x_{2j-2})+4f(x_{2j-1})+f(x_{2j}){\bigg ]}\\{}&={\frac {h}{3}}{\bigg [}f(x_{0})+2\sum _{j=1}^{n/2-1}f(x_{2j})+4\sum _{j=1}^{n/2}f(x_{2j-1})+f(x_{n}){\bigg ]},\end{aligned}}$

kde $x_{j}=a+jh$ pro $j=0,1,...,n-1,n$ s $h=(b-a)/n$; konkrétně, $x_{0}=a$ a $x_{n}=b$. Toto složené pravidlo je pro $n=2$ totožné s normálním Simpsonovým pravidlem popsaným v předchozí části.

Chyba složeného Simpsonova pravidla je

$-{\frac {h^{4}}{180}}(b-a)f^{(4)}(\xi ),$

kde $\xi$ je nějaké číslo mezi $a$ a $b$ a $h=(b-a)/n$ je „délka kroku“. Absolutní hodnota chyby je omezena výrazem

${\tfrac {h^{4}}{180}}(b-a)\max _{\xi \in \langle a,b\rangle }|f^{(4)}(\xi )|.$

Toto formulace rozděluje interval $\langle a,b\rangle$ na podintervaly stejné délky. V praxi je často výhodné používat podintervaly různých délek, a podrobněji se zabývat místy, kde je integrand méně rozumný. To vede k adaptivní Simpsonově metodě.

Simpsonovo 3/8 pravidlo

Simpsonovo 3/8 pravidlo, také nazývaný Simpsonovo druhé pravidlo, je druhá metoda pro numerickou integraci, kterou navrhl Thomas Simpson. Používá kubickou interpolaci místo kvadratické. Simpsonovo 3/8 pravidlo popisuje vzorec:

$\int _{a}^{b}f(x)\,dx\approx {\tfrac {3h}{8}}\left[f(a)+3f\left({\tfrac {2a+b}{3}}\right)+3f\left({\tfrac {a+2b}{3}}\right)+f(b)\right]={\tfrac {(b-a)}{8}}\left[f(a)+3f\left({\tfrac {2a+b}{3}}\right)+3f\left({\tfrac {a+2b}{3}}\right)+f(b)\right],$

kde b − a = 3h. Chyba této metody je:

$-{\tfrac {(b-a)^{5}}{6480}}f^{(4)}(\xi ),$

kde $\xi$ je nějaké číslo mezi $a$ a $b$. 3/8 pravidlo má tedy dvojnásobnou přesnost než standardní metoda, ale používá o jednu funkční hodnotu více. Existuje také složené 3/8 pravidlo, obdobně jako výše.

Dalším zobecněním tohoto konceptu na interpolaci pomocí polynomů libovolného stupně jsou Newtonovy–Cotesovy vzorce.

Složené Simpsonovo 3/8 pravidlo

Pokud interval $\langle a,b\rangle$ rozdělíme na $n$ podintervalů délky $h=(b-a)/n$ a budeme počítat hodnoty funkce $f$ v bodech $x_{i}=a+ih$ dostáváme vzorec

${\begin{aligned}\int _{a}^{b}f(x)\,dx&\approx {\tfrac {3h}{8}}\left[f(x_{0})+3f(x_{1})+3f(x_{2})+2f(x_{3})+3f(x_{4})+3f(x_{5})+2f(x_{6})+\cdots +3f(x_{n-2})+3f(x_{n-1})+f(x_{n})\right].\\&={\frac {3h}{8}}\left[f(x_{0})+3\sum _{i\neq 3k}^{n-1}f(x_{i})+2\sum _{j=1}^{n/3-1}f(x_{3j})+f(x_{n})\right]\qquad {\text{For: }}k\in \mathbb {N} _{0}\end{aligned}}$

přičemž zbytek je:

$-{\frac {h^{4}}{80}}(b-a)f^{(4)}(\xi ),$

Pro použití tohoto pravidla musí být $n$ dělitelné třemi.

Alternativní rozšíření Simpsonova pravidla

Složené Simpsonovo pravidlo lze také aplikovat na překrývající se segmenty, což dává:

$\int _{a}^{b}f(x)\,dx\approx {\tfrac {h}{48}}{\bigg [}17f(x_{0})+59f(x_{1})+43f(x_{2})+49f(x_{3})+48\sum _{i=4}^{n-4}f(x_{i})+49f(x_{n-3})+43f(x_{n-2})+59f(x_{n-1})+17f(x_{n}){\bigg ]}.$

Tento vzorec lze získat zkombinováním původního složeného Simpsonova pravidla s pravidlem, kde se použije Simpsonovo 3/8 pravidlo v krajních podintervalech a standardní tříbodové pravidlo ve zbývajících podintervalech. Výsledek získáme použitím průměru obou vzorců.

Simpsonova pravidla v případě úzké špičky

Při zjišťování plochy ohraničené funkcí, která má úzké špičky jsou Simpsonova pravidla mnohem méně efektivní než lichoběžníková metoda. Jmenovitě složené Simpsonovo 1/3 pravidlo vyžaduje 1,8krát více bodů, aby se dosáhlo stejné přesnosti jako při použití lichoběžníkové metody. Složené Simpsonovo 3/8 pravidlo je v tomto případě ještě méně přesné. Integraci Simpsonovým 1/3 pravidlem lze považovat za součet 2/3 výsledku získaného lichoběžníkovou metodou s krokem h a 1/3 aplikace obdélníkového pravidla s krokem 2h. Zprůměrování složených součtů Simpsonova 1/3 pravidla se správně posunutými rámci dává následující vzorce:

$\int _{a}^{b}f(x)\,dx\approx {\tfrac {h}{24}}{\bigg [}-f(x_{-1})+12f(x_{0})+25f(x_{1})+24\sum _{i=2}^{n-2}f(x_{i})+25f(x_{n-1})+12f(x_{n})-f(x_{n+1}){\bigg ]}$

kde se používají dva body mimo integrační oblast a

$\int _{a}^{b}f(x)\,dx\approx {\tfrac {h}{24}}{\bigg [}9f(x_{0})+28f(x_{1})+23f(x_{2})+24\sum _{i=3}^{n-3}f(x_{i})+23f(x_{n-2})+28f(x_{n-1})+9f(x_{n}){\bigg ]}$

kde se používají pouze body uvnitř integrační oblasti. Aplikace druhého pravidla na oblast se 3 body dává Simpsonovo 1/3 pravidlo, na oblast se 4 body Simpsonovo 3/8 pravidlo.

Tato pravidla se velmi podobají Pressovu alternativnímu rozšíření Simpsonova pravidla. Koeficienty ve velké části oblasti, kde se počítá integrál, jsou rovny jedné, rozdíly jsou pouze na jejích okrajích. Tato tři pravidla lze spojit s Eulerovým-MacLaurinovým vzorcem s první derivací, který se nazývá Eulerova-MacLaurinova integrační pravidla prvního řádu. Dvě výše uvedená pravidla se liší pouze způsobem, jak se počítá první derivace na konci oblasti. Člen s první derivací v Eulerově–MacLaurinově integračním pravidle zodpovídá za integrál druhé derivace, který se rovná rozdílu prvních derivací na okrajích integrační oblasti. Eulerova-Maclaurinova pravidla vyššího řádu je možné generovat přidáním rozdílu třetí, páté, atd. derivace s koeficienty, které jsou definovány Eulerovým–MacLaurinovým vzorcem.

Složené Simpsonovo pravidlo pro nepravidelně vzdálené body

V případech, kdy je potřeba, aby byl integrační interval $I=\langle a,b\rangle$ rozdělen na nestejné intervaly – například kvůli nerovnoměrnému vzorkování dat nebo chybějícím nebo poškozeným datům. Předpokládejme, že rozdělíme interval $I$ na sudé počet $N$ podintervalů se šířkami $h_{k}$. Pak složené Simpsonovo pravidlo popisuje vztah

$\int _{a}^{b}f(x)\,\mathrm {d} x=\sum _{i=0}^{N/2-1}\left(\alpha _{i}f_{2i+2}+\beta _{i}f_{2i+1}+\eta _{i}f_{2i}\right),$

kde $f_{k}=f\left(a+\sum _{i=0}^{k-1}h_{i}\right)$ jsou funkce hodnoty v $k$-tém bodě na intervalu $I$, a koeficienty $\alpha _{i},\,\beta _{i},$ a $\eta _{i}$ jsou

$\alpha _{i}={\frac {2h_{2i+1}^{3}-h_{2i}^{3}+3h_{2i}h_{2i+1}^{2}}{6h_{2i+1}(h_{2i+1}+h_{2i})}},$

$\beta _{i}={\frac {h_{2i+1}^{3}+h_{2i}^{3}+3h_{2i+1}h_{2i}(h_{2i+1}+h_{2i})}{6h_{2i+1}h_{2i}}},{\text{and}}$

$\eta _{i}={\frac {2h_{2i}^{3}-h_{2i+1}^{3}+3h_{2i+1}h_{2i}^{2}}{6h_{2i}(h_{2i+1}+h_{2i})}}.$

Pro lichý počet $N$ podintervalů se výše uvedený vzorec používá až na druhý až poslední interval, a poslední interval je ošetřen odděleně přidáním následujícího členu k výsledku:

$\alpha f_{N}+\beta f_{N-1}-\eta f_{N-2},$

kde

$\alpha ={\frac {2h_{N-1}^{2}+3h_{N-1}h_{N-2}}{6(h_{N-2}+h_{N-1})}},$

$\beta ={\frac {h_{N-1}^{2}+3h_{N-1}h_{N-2}}{6h_{N-2}}},\,{\text{and}}$

$\eta ={\frac {h_{N-1}^{3}}{6h_{N-2}(h_{N-2}+h_{N-1})}}.$

Ukázka implementace Simpsonovy metody v Pythonu

from collections.abc import Sequence

def simpson_nonuniform(x: Sequence[float], f: Sequence[float]) -> float:
    """
    Integrace funkce Simpsonovým pravidlem pro nestajně vzdálené body.

    :param x: Seznam bodů, v nichž jsou známy hodnoty integrované funkce
    :param f: Seznam hodnot integrované funkce v bodech zadaných v x

    :return: approximation for the integral

    Viz ``scipy.integrate.simpson`` a ``_basic_simpson``
    pro efektivnější implementaci používající broadcast z modulu numpy.
    """
    N = len(x) - 1
    h = [x[i + 1] - x[i] for i in range(0, N)]
    assert N > 0

    result = 0.0
    for i in range(1, N, 2):
        h0, h1 = h[i - 1], h[i]
        hph, hdh, hmh = h1 + h0, h1 / h0, h1 * h0
        result += (hph / 6) * (
            (2 - hdh) * f[i - 1] + (hph**2 / hmh) * f[i] + (2 - 1 / hdh) * f[i + 1]
        )

    if N % 2 == 1:
        h0, h1 = h[N - 2], h[N - 1]
        result += f[N]     * (2 * h1 ** 2 + 3 * h0 * h1) / (6 * (h0 + h1))
        result += f[N - 1] * (h1 ** 2 + 3 * h1 * h0)     / (6 * h0)
        result -= f[N - 2] * h1 ** 3                     / (6 * h0 * (h0 + h1))
    return result

Odkazy

Reference

V tomto článku byl použit překlad textu z článku Simpson's rule na anglické Wikipedii.

• ATKINSON, Kendall E., 1989. An Introduction to Numerical Analysis. 2. vyd. [s.l.]: John Wiley & Sons. ISBN 0-471-50023-2. 
• BURDEN, Richard L.; FAIRES, J. Douglas, 2000. Numerical Analysis. 7. vyd. [s.l.]: Brooks/Cole. Dostupné online. ISBN 0-534-38216-9. 
• MATTHEWS, John H., 2004. Simpson's 3/8 Rule for Numerical Integration [online]. California State University, Fullerton, 2004 [cit. 2008-11-11]. Dostupné v archivu pořízeném z originálu dne 2008-12-04. 
• PRESS, William H.; FLANNERY, Brian P.; VETTERLING, William T.; TEUKOLSKY, Saul A., 1989. Numerical Recipes in Pascal: The Art of Scientific Computing. [s.l.]: Cambridge University Press. Dostupné online. ISBN 0-521-37516-9. 
• SÜLI, Endre; MAYERS, David, 2003. An Introduction to Numerical Analysis. [s.l.]: Cambridge University Press. ISBN 0-521-00794-1. 
• KAW, Autar; KALU, Egwu; NGUYEN, Duc, 2008. Numerical Methods with Applications [online]. 2008. Dostupné online. 
• WEISSTEIN, Eric W., 2010. Newton-Cotes Formulas [online]. MathWorld, 2010 [cit. 2010-08-02]. Dostupné online. 
• CARTWRIGHT, Kenneth V., 2016. Simpson's Rule Integration with MS Excel and Irregularly-spaced Data. Journal of Mathematical Science and Mathematics Education. Roč. 11, čís. 2, s. 34–42. Dostupné online. 
• KALAMBET, Yuri; KOZMIN, Yuri; SAMOKHIN, Andrey. Comparison of integration rules in the case of very narrow chromatographic peaks. Chemometrics and Intelligent Laboratory Systems. 2018, roč. 179, s. 22–30. ISSN 0169-7439. DOI 10.1016/j.chemolab.2018.06.001. 
• KYLÄNPÄÄ, Ilkka, 2019. Computational Physics course. In: Tampere University: [s.n.].

Související články

• Gaussovo kvadraturní pravidlo

Externí odkazy

• Obrázky, zvuky či videa k tématu Simpsonovo pravidlo na Wikimedia Commons
• Simpson formula na Encyclopedia of Mathematics, EMS Press, 2001 [1994]
• "Simpson's Rule" MathWorld. Weisstein, Eric W.
• Aplikace Simpsonova pravidla – Earthwork Excavation (Note: Vzorec popsaný v toto stránka je správný ale existují chyby v výpočet který musí dávají výsledek 569m3 a not 623m3 jako uvedeny)
• Simpson's 1/3. pravidlo integrace – Poznámky/zmiňuje, PPT, Mathcad, Matlab, Mathematica, Maple v Numerický Metody pro KMEN univerzitní
• Podrobný popis počítačové implementace popsal Dorai Sitaram v Teach Yourself Scheme v Fixnum Dnů, Appendix C^{[nedostupný zdroj]}