\(\newcommand{\footnotename}{footnote}\)
\(\def \LWRfootnote {1}\)
\(\newcommand {\footnote }[2][\LWRfootnote ]{{}^{\mathrm {#1}}}\)
\(\newcommand {\footnotemark }[1][\LWRfootnote ]{{}^{\mathrm {#1}}}\)
\(\newcommand \ensuremath [1]{#1}\)
\(\newcommand {\LWRframebox }[2][]{\fbox {#2}} \newcommand {\framebox }[1][]{\LWRframebox } \)
\(\newcommand {\setlength }[2]{}\)
\(\newcommand {\addtolength }[2]{}\)
\(\newcommand {\setcounter }[2]{}\)
\(\newcommand {\addtocounter }[2]{}\)
\(\newcommand {\cline }[1]{}\)
\(\newcommand {\directlua }[1]{\text {(directlua)}}\)
\(\newcommand {\luatexdirectlua }[1]{\text {(directlua)}}\)
\(\newcommand {\protect }{}\)
\(\def \LWRabsorbnumber #1 {}\)
\(\def \LWRabsorbquotenumber "#1 {}\)
\(\def \mathchar {\ifnextchar "\LWRabsorbquotenumber \LWRabsorbnumber }\)
\(\def \mathcode #1={\mathchar }\)
\(\let \delcode \mathcode \)
\(\let \delimiter \mathchar \)
\(\let \LWRref \ref \)
\(\renewcommand {\ref }{\ifstar \LWRref \LWRref }\)
\(\newcommand {\intertext }[1]{\text {#1}\notag \\}\)
\(\newcommand {\mathllap }[2][]{{#1#2}}\)
\(\newcommand {\mathrlap }[2][]{{#1#2}}\)
\(\newcommand {\mathclap }[2][]{{#1#2}}\)
\(\newcommand {\mathmbox }[1]{#1}\)
\(\newcommand {\clap }[1]{#1}\)
\(\newcommand {\LWRmathmakebox }[2][]{#2}\)
\(\newcommand {\mathmakebox }[1][]{\LWRmathmakebox }\)
\(\newcommand {\cramped }[2][]{{#1#2}}\)
\(\newcommand {\crampedllap }[2][]{{#1#2}}\)
\(\newcommand {\crampedrlap }[2][]{{#1#2}}\)
\(\newcommand {\crampedclap }[2][]{{#1#2}}\)
\(\newenvironment {crampedsubarray}[1]{}{}\)
\(\newcommand {\crampedsubstack }{}\)
\(\newcommand {\smashoperator }[2][]{#2\limits }\)
\(\newcommand {\adjustlimits }{}\)
\(\newcommand {\SwapAboveDisplaySkip }{}\)
\(\require {extpfeil}\)
\(\Newextarrow \xleftrightarrow {10,10}{0x2194}\)
\(\Newextarrow \xLeftarrow {10,10}{0x21d0}\)
\(\Newextarrow \xhookleftarrow {10,10}{0x21a9}\)
\(\Newextarrow \xmapsto {10,10}{0x21a6}\)
\(\Newextarrow \xRightarrow {10,10}{0x21d2}\)
\(\Newextarrow \xLeftrightarrow {10,10}{0x21d4}\)
\(\Newextarrow \xhookrightarrow {10,10}{0x21aa}\)
\(\Newextarrow \xrightharpoondown {10,10}{0x21c1}\)
\(\Newextarrow \xleftharpoondown {10,10}{0x21bd}\)
\(\Newextarrow \xrightleftharpoons {10,10}{0x21cc}\)
\(\Newextarrow \xrightharpoonup {10,10}{0x21c0}\)
\(\Newextarrow \xleftharpoonup {10,10}{0x21bc}\)
\(\Newextarrow \xleftrightharpoons {10,10}{0x21cb}\)
\(\newcommand {\LWRdounderbracket }[3]{\underset {#3}{\underline {#1}}}\)
\(\newcommand {\LWRunderbracket }[2][]{\LWRdounderbracket {#2}}\)
\(\newcommand {\underbracket }[1][]{\LWRunderbracket }\)
\(\newcommand {\LWRdooverbracket }[3]{\overset {#3}{\overline {#1}}}\)
\(\newcommand {\LWRoverbracket }[2][]{\LWRdooverbracket {#2}}\)
\(\newcommand {\overbracket }[1][]{\LWRoverbracket }\)
\(\newcommand {\LATEXunderbrace }[1]{\underbrace {#1}}\)
\(\newcommand {\LATEXoverbrace }[1]{\overbrace {#1}}\)
\(\newenvironment {matrix*}[1][]{\begin {matrix}}{\end {matrix}}\)
\(\newenvironment {pmatrix*}[1][]{\begin {pmatrix}}{\end {pmatrix}}\)
\(\newenvironment {bmatrix*}[1][]{\begin {bmatrix}}{\end {bmatrix}}\)
\(\newenvironment {Bmatrix*}[1][]{\begin {Bmatrix}}{\end {Bmatrix}}\)
\(\newenvironment {vmatrix*}[1][]{\begin {vmatrix}}{\end {vmatrix}}\)
\(\newenvironment {Vmatrix*}[1][]{\begin {Vmatrix}}{\end {Vmatrix}}\)
\(\newenvironment {smallmatrix*}[1][]{\begin {matrix}}{\end {matrix}}\)
\(\newenvironment {psmallmatrix*}[1][]{\begin {pmatrix}}{\end {pmatrix}}\)
\(\newenvironment {bsmallmatrix*}[1][]{\begin {bmatrix}}{\end {bmatrix}}\)
\(\newenvironment {Bsmallmatrix*}[1][]{\begin {Bmatrix}}{\end {Bmatrix}}\)
\(\newenvironment {vsmallmatrix*}[1][]{\begin {vmatrix}}{\end {vmatrix}}\)
\(\newenvironment {Vsmallmatrix*}[1][]{\begin {Vmatrix}}{\end {Vmatrix}}\)
\(\newenvironment {psmallmatrix}[1][]{\begin {pmatrix}}{\end {pmatrix}}\)
\(\newenvironment {bsmallmatrix}[1][]{\begin {bmatrix}}{\end {bmatrix}}\)
\(\newenvironment {Bsmallmatrix}[1][]{\begin {Bmatrix}}{\end {Bmatrix}}\)
\(\newenvironment {vsmallmatrix}[1][]{\begin {vmatrix}}{\end {vmatrix}}\)
\(\newenvironment {Vsmallmatrix}[1][]{\begin {Vmatrix}}{\end {Vmatrix}}\)
\(\newcommand {\LWRmultlined }[1][]{\begin {multline*}}\)
\(\newenvironment {multlined}[1][]{\LWRmultlined }{\end {multline*}}\)
\(\let \LWRorigshoveleft \shoveleft \)
\(\renewcommand {\shoveleft }[1][]{\LWRorigshoveleft }\)
\(\let \LWRorigshoveright \shoveright \)
\(\renewcommand {\shoveright }[1][]{\LWRorigshoveright }\)
\(\newenvironment {dcases}{\begin {cases}}{\end {cases}}\)
\(\newenvironment {dcases*}{\begin {cases}}{\end {cases}}\)
\(\newenvironment {rcases}{\begin {cases}}{\end {cases}}\)
\(\newenvironment {rcases*}{\begin {cases}}{\end {cases}}\)
\(\newenvironment {drcases}{\begin {cases}}{\end {cases}}\)
\(\newenvironment {drcases*}{\begin {cases}}{\end {cases}}\)
\(\newenvironment {cases*}{\begin {cases}}{\end {cases}}\)
\(\newcommand {\MoveEqLeft }[1][]{}\)
\(\def \LWRAboxed #1!|!{\fbox {\(#1\)}&\fbox {\(#2\)}} \newcommand {\Aboxed }[1]{\LWRAboxed #1&&!|!} \)
\( \newcommand {\LWRABLines }[1][\Updownarrow ]{#1 \notag \\}\newcommand {\ArrowBetweenLines }{\ifstar \LWRABLines \LWRABLines } \)
\(\newcommand {\shortintertext }[1]{\text {#1}\notag \\}\)
\(\newcommand {\vdotswithin }[1]{\hspace {.5em}\vdots }\)
\(\newcommand {\LWRshortvdotswithinstar }[1]{\vdots \hspace {.5em} & \\}\)
\(\newcommand {\LWRshortvdotswithinnostar }[1]{& \hspace {.5em}\vdots \\}\)
\(\newcommand {\shortvdotswithin }{\ifstar \LWRshortvdotswithinstar \LWRshortvdotswithinnostar }\)
\(\newcommand {\MTFlushSpaceAbove }{}\)
\(\newcommand {\MTFlushSpaceBelow }{\\}\)
\(\newcommand \lparen {(}\)
\(\newcommand \rparen {)}\)
\(\newcommand {\ordinarycolon }{:}\)
\(\newcommand {\vcentcolon }{\mathrel {\mathop \ordinarycolon }}\)
\(\newcommand \dblcolon {\vcentcolon \vcentcolon }\)
\(\newcommand \coloneqq {\vcentcolon =}\)
\(\newcommand \Coloneqq {\dblcolon =}\)
\(\newcommand \coloneq {\vcentcolon {-}}\)
\(\newcommand \Coloneq {\dblcolon {-}}\)
\(\newcommand \eqqcolon {=\vcentcolon }\)
\(\newcommand \Eqqcolon {=\dblcolon }\)
\(\newcommand \eqcolon {\mathrel {-}\vcentcolon }\)
\(\newcommand \Eqcolon {\mathrel {-}\dblcolon }\)
\(\newcommand \colonapprox {\vcentcolon \approx }\)
\(\newcommand \Colonapprox {\dblcolon \approx }\)
\(\newcommand \colonsim {\vcentcolon \sim }\)
\(\newcommand \Colonsim {\dblcolon \sim }\)
\(\newcommand {\nuparrow }{\mathrel {\cancel {\uparrow }}}\)
\(\newcommand {\ndownarrow }{\mathrel {\cancel {\downarrow }}}\)
\(\newcommand {\bigtimes }{\mathop {\Large \times }\limits }\)
\(\newcommand {\prescript }[3]{{}^{#1}_{#2}#3}\)
\(\newenvironment {lgathered}{\begin {gathered}}{\end {gathered}}\)
\(\newenvironment {rgathered}{\begin {gathered}}{\end {gathered}}\)
\(\newcommand {\splitfrac }[2]{{}^{#1}_{#2}}\)
\(\let \splitdfrac \splitfrac \)
\(\newcommand {\LWRoverlaysymbols }[2]{\mathord {\smash {\mathop {#2\strut }\limits ^{\smash {\lower 3ex{#1}}}}\strut }}\)
\(\newcommand{\alphaup}{\unicode{x03B1}}\)
\(\newcommand{\betaup}{\unicode{x03B2}}\)
\(\newcommand{\gammaup}{\unicode{x03B3}}\)
\(\newcommand{\digammaup}{\unicode{x03DD}}\)
\(\newcommand{\deltaup}{\unicode{x03B4}}\)
\(\newcommand{\epsilonup}{\unicode{x03F5}}\)
\(\newcommand{\varepsilonup}{\unicode{x03B5}}\)
\(\newcommand{\zetaup}{\unicode{x03B6}}\)
\(\newcommand{\etaup}{\unicode{x03B7}}\)
\(\newcommand{\thetaup}{\unicode{x03B8}}\)
\(\newcommand{\varthetaup}{\unicode{x03D1}}\)
\(\newcommand{\iotaup}{\unicode{x03B9}}\)
\(\newcommand{\kappaup}{\unicode{x03BA}}\)
\(\newcommand{\varkappaup}{\unicode{x03F0}}\)
\(\newcommand{\lambdaup}{\unicode{x03BB}}\)
\(\newcommand{\muup}{\unicode{x03BC}}\)
\(\newcommand{\nuup}{\unicode{x03BD}}\)
\(\newcommand{\xiup}{\unicode{x03BE}}\)
\(\newcommand{\omicronup}{\unicode{x03BF}}\)
\(\newcommand{\piup}{\unicode{x03C0}}\)
\(\newcommand{\varpiup}{\unicode{x03D6}}\)
\(\newcommand{\rhoup}{\unicode{x03C1}}\)
\(\newcommand{\varrhoup}{\unicode{x03F1}}\)
\(\newcommand{\sigmaup}{\unicode{x03C3}}\)
\(\newcommand{\varsigmaup}{\unicode{x03C2}}\)
\(\newcommand{\tauup}{\unicode{x03C4}}\)
\(\newcommand{\upsilonup}{\unicode{x03C5}}\)
\(\newcommand{\phiup}{\unicode{x03D5}}\)
\(\newcommand{\varphiup}{\unicode{x03C6}}\)
\(\newcommand{\chiup}{\unicode{x03C7}}\)
\(\newcommand{\psiup}{\unicode{x03C8}}\)
\(\newcommand{\omegaup}{\unicode{x03C9}}\)
\(\newcommand{\Alphaup}{\unicode{x0391}}\)
\(\newcommand{\Betaup}{\unicode{x0392}}\)
\(\newcommand{\Gammaup}{\unicode{x0393}}\)
\(\newcommand{\Digammaup}{\unicode{x03DC}}\)
\(\newcommand{\Deltaup}{\unicode{x0394}}\)
\(\newcommand{\Epsilonup}{\unicode{x0395}}\)
\(\newcommand{\Zetaup}{\unicode{x0396}}\)
\(\newcommand{\Etaup}{\unicode{x0397}}\)
\(\newcommand{\Thetaup}{\unicode{x0398}}\)
\(\newcommand{\Varthetaup}{\unicode{x03F4}}\)
\(\newcommand{\Iotaup}{\unicode{x0399}}\)
\(\newcommand{\Kappaup}{\unicode{x039A}}\)
\(\newcommand{\Lambdaup}{\unicode{x039B}}\)
\(\newcommand{\Muup}{\unicode{x039C}}\)
\(\newcommand{\Nuup}{\unicode{x039D}}\)
\(\newcommand{\Xiup}{\unicode{x039E}}\)
\(\newcommand{\Omicronup}{\unicode{x039F}}\)
\(\newcommand{\Piup}{\unicode{x03A0}}\)
\(\newcommand{\Varpiup}{\unicode{x03D6}}\)
\(\newcommand{\Rhoup}{\unicode{x03A1}}\)
\(\newcommand{\Sigmaup}{\unicode{x03A3}}\)
\(\newcommand{\Tauup}{\unicode{x03A4}}\)
\(\newcommand{\Upsilonup}{\unicode{x03A5}}\)
\(\newcommand{\Phiup}{\unicode{x03A6}}\)
\(\newcommand{\Chiup}{\unicode{x03A7}}\)
\(\newcommand{\Psiup}{\unicode{x03A8}}\)
\(\newcommand{\Omegaup}{\unicode{x03A9}}\)
\(\newcommand{\alphait}{\unicode{x1D6FC}}\)
\(\newcommand{\betait}{\unicode{x1D6FD}}\)
\(\newcommand{\gammait}{\unicode{x1D6FE}}\)
\(\newcommand{\digammait}{\mathit{\unicode{x03DD}}}\)
\(\newcommand{\deltait}{\unicode{x1D6FF}}\)
\(\newcommand{\epsilonit}{\unicode{x1D716}}\)
\(\newcommand{\varepsilonit}{\unicode{x1D700}}\)
\(\newcommand{\zetait}{\unicode{x1D701}}\)
\(\newcommand{\etait}{\unicode{x1D702}}\)
\(\newcommand{\thetait}{\unicode{x1D703}}\)
\(\newcommand{\varthetait}{\unicode{x1D717}}\)
\(\newcommand{\iotait}{\unicode{x1D704}}\)
\(\newcommand{\kappait}{\unicode{x1D705}}\)
\(\newcommand{\varkappait}{\unicode{x1D718}}\)
\(\newcommand{\lambdait}{\unicode{x1D706}}\)
\(\newcommand{\muit}{\unicode{x1D707}}\)
\(\newcommand{\nuit}{\unicode{x1D708}}\)
\(\newcommand{\xiit}{\unicode{x1D709}}\)
\(\newcommand{\omicronit}{\unicode{x1D70A}}\)
\(\newcommand{\piit}{\unicode{x1D70B}}\)
\(\newcommand{\varpiit}{\unicode{x1D71B}}\)
\(\newcommand{\rhoit}{\unicode{x1D70C}}\)
\(\newcommand{\varrhoit}{\unicode{x1D71A}}\)
\(\newcommand{\sigmait}{\unicode{x1D70E}}\)
\(\newcommand{\varsigmait}{\unicode{x1D70D}}\)
\(\newcommand{\tauit}{\unicode{x1D70F}}\)
\(\newcommand{\upsilonit}{\unicode{x1D710}}\)
\(\newcommand{\phiit}{\unicode{x1D719}}\)
\(\newcommand{\varphiit}{\unicode{x1D711}}\)
\(\newcommand{\chiit}{\unicode{x1D712}}\)
\(\newcommand{\psiit}{\unicode{x1D713}}\)
\(\newcommand{\omegait}{\unicode{x1D714}}\)
\(\newcommand{\Alphait}{\unicode{x1D6E2}}\)
\(\newcommand{\Betait}{\unicode{x1D6E3}}\)
\(\newcommand{\Gammait}{\unicode{x1D6E4}}\)
\(\newcommand{\Digammait}{\mathit{\unicode{x03DC}}}\)
\(\newcommand{\Deltait}{\unicode{x1D6E5}}\)
\(\newcommand{\Epsilonit}{\unicode{x1D6E6}}\)
\(\newcommand{\Zetait}{\unicode{x1D6E7}}\)
\(\newcommand{\Etait}{\unicode{x1D6E8}}\)
\(\newcommand{\Thetait}{\unicode{x1D6E9}}\)
\(\newcommand{\Varthetait}{\unicode{x1D6F3}}\)
\(\newcommand{\Iotait}{\unicode{x1D6EA}}\)
\(\newcommand{\Kappait}{\unicode{x1D6EB}}\)
\(\newcommand{\Lambdait}{\unicode{x1D6EC}}\)
\(\newcommand{\Muit}{\unicode{x1D6ED}}\)
\(\newcommand{\Nuit}{\unicode{x1D6EE}}\)
\(\newcommand{\Xiit}{\unicode{x1D6EF}}\)
\(\newcommand{\Omicronit}{\unicode{x1D6F0}}\)
\(\newcommand{\Piit}{\unicode{x1D6F1}}\)
\(\newcommand{\Rhoit}{\unicode{x1D6F2}}\)
\(\newcommand{\Sigmait}{\unicode{x1D6F4}}\)
\(\newcommand{\Tauit}{\unicode{x1D6F5}}\)
\(\newcommand{\Upsilonit}{\unicode{x1D6F6}}\)
\(\newcommand{\Phiit}{\unicode{x1D6F7}}\)
\(\newcommand{\Chiit}{\unicode{x1D6F8}}\)
\(\newcommand{\Psiit}{\unicode{x1D6F9}}\)
\(\newcommand{\Omegait}{\unicode{x1D6FA}}\)
\(\let \digammaup \Digammaup \)
\(\renewcommand {\digammait }{\mathit {\digammaup }}\)
\(\newcommand {\smallin }{\unicode {x220A}}\)
\(\newcommand {\smallowns }{\unicode {x220D}}\)
\(\newcommand {\notsmallin }{\LWRoverlaysymbols {/}{\unicode {x220A}}}\)
\(\newcommand {\notsmallowns }{\LWRoverlaysymbols {/}{\unicode {x220D}}}\)
\(\newcommand {\rightangle }{\unicode {x221F}}\)
\(\newcommand {\intclockwise }{\unicode {x2231}}\)
\(\newcommand {\ointclockwise }{\unicode {x2232}}\)
\(\newcommand {\ointctrclockwise }{\unicode {x2233}}\)
\(\newcommand {\oiint }{\unicode {x222F}}\)
\(\newcommand {\oiiint }{\unicode {x2230}}\)
\(\newcommand {\ddag }{\unicode {x2021}}\)
\(\newcommand {\P }{\unicode {x00B6}}\)
\(\newcommand {\copyright }{\unicode {x00A9}}\)
\(\newcommand {\dag }{\unicode {x2020}}\)
\(\newcommand {\pounds }{\unicode {x00A3}}\)
\(\newcommand {\iddots }{\unicode {x22F0}}\)
\(\newcommand {\utimes }{\overline {\times }}\)
\(\newcommand {\dtimes }{\underline {\times }}\)
\(\newcommand {\udtimes }{\overline {\underline {\times }}}\)
\(\newcommand {\leftwave }{\left \{}\)
\(\newcommand {\rightwave }{\right \}}\)
\(\newcommand {\toprule }[1][]{\hline }\)
\(\let \midrule \toprule \)
\(\let \bottomrule \toprule \)
\(\newcommand {\cmidrule }[2][]{}\)
\(\newcommand {\morecmidrules }{}\)
\(\newcommand {\specialrule }[3]{\hline }\)
\(\newcommand {\addlinespace }[1][]{}\)
\(\newcommand {\LWRsubmultirow }[2][]{#2}\)
\(\newcommand {\LWRmultirow }[2][]{\LWRsubmultirow }\)
\(\newcommand {\multirow }[2][]{\LWRmultirow }\)
\(\newcommand {\mrowcell }{}\)
\(\newcommand {\mcolrowcell }{}\)
\(\newcommand {\STneed }[1]{}\)
\( \newcommand {\multicolumn }[3]{#3}\)
\(\newcommand {\tothe }[1]{^{#1}}\)
\(\newcommand {\raiseto }[2]{{#2}^{#1}}\)
\(\newcommand {\ang }[2][]{(\mathrm {#2})\degree }\)
\(\newcommand {\num }[2][]{\mathrm {#2}}\)
\(\newcommand {\si }[2][]{\mathrm {#2}}\)
\(\newcommand {\LWRSI }[2][]{\mathrm {#1\LWRSInumber \,#2}}\)
\(\newcommand {\SI }[2][]{\def \LWRSInumber {#2}\LWRSI }\)
\(\newcommand {\numlist }[2][]{\mathrm {#2}}\)
\(\newcommand {\numrange }[3][]{\mathrm {#2\,\unicode {x2013}\,#3}}\)
\(\newcommand {\SIlist }[3][]{\mathrm {#2\,#3}}\)
\(\newcommand {\SIrange }[4][]{\mathrm {#2\,#4\,\unicode {x2013}\,#3\,#4}}\)
\(\newcommand {\tablenum }[2][]{\mathrm {#2}}\)
\(\newcommand {\ampere }{\mathrm {A}}\)
\(\newcommand {\candela }{\mathrm {cd}}\)
\(\newcommand {\kelvin }{\mathrm {K}}\)
\(\newcommand {\kilogram }{\mathrm {kg}}\)
\(\newcommand {\metre }{\mathrm {m}}\)
\(\newcommand {\mole }{\mathrm {mol}}\)
\(\newcommand {\second }{\mathrm {s}}\)
\(\newcommand {\becquerel }{\mathrm {Bq}}\)
\(\newcommand {\degreeCelsius }{\unicode {x2103}}\)
\(\newcommand {\coulomb }{\mathrm {C}}\)
\(\newcommand {\farad }{\mathrm {F}}\)
\(\newcommand {\gray }{\mathrm {Gy}}\)
\(\newcommand {\hertz }{\mathrm {Hz}}\)
\(\newcommand {\henry }{\mathrm {H}}\)
\(\newcommand {\joule }{\mathrm {J}}\)
\(\newcommand {\katal }{\mathrm {kat}}\)
\(\newcommand {\lumen }{\mathrm {lm}}\)
\(\newcommand {\lux }{\mathrm {lx}}\)
\(\newcommand {\newton }{\mathrm {N}}\)
\(\newcommand {\ohm }{\mathrm {\Omega }}\)
\(\newcommand {\pascal }{\mathrm {Pa}}\)
\(\newcommand {\radian }{\mathrm {rad}}\)
\(\newcommand {\siemens }{\mathrm {S}}\)
\(\newcommand {\sievert }{\mathrm {Sv}}\)
\(\newcommand {\steradian }{\mathrm {sr}}\)
\(\newcommand {\tesla }{\mathrm {T}}\)
\(\newcommand {\volt }{\mathrm {V}}\)
\(\newcommand {\watt }{\mathrm {W}}\)
\(\newcommand {\weber }{\mathrm {Wb}}\)
\(\newcommand {\day }{\mathrm {d}}\)
\(\newcommand {\degree }{\mathrm {^\circ }}\)
\(\newcommand {\hectare }{\mathrm {ha}}\)
\(\newcommand {\hour }{\mathrm {h}}\)
\(\newcommand {\litre }{\mathrm {l}}\)
\(\newcommand {\liter }{\mathrm {L}}\)
\(\newcommand {\arcminute }{^\prime }\)
\(\newcommand {\minute }{\mathrm {min}}\)
\(\newcommand {\arcsecond }{^{\prime \prime }}\)
\(\newcommand {\tonne }{\mathrm {t}}\)
\(\newcommand {\astronomicalunit }{au}\)
\(\newcommand {\atomicmassunit }{u}\)
\(\newcommand {\bohr }{\mathit {a}_0}\)
\(\newcommand {\clight }{\mathit {c}_0}\)
\(\newcommand {\dalton }{\mathrm {D}_\mathrm {a}}\)
\(\newcommand {\electronmass }{\mathit {m}_{\mathrm {e}}}\)
\(\newcommand {\electronvolt }{\mathrm {eV}}\)
\(\newcommand {\elementarycharge }{\mathit {e}}\)
\(\newcommand {\hartree }{\mathit {E}_{\mathrm {h}}}\)
\(\newcommand {\planckbar }{\mathit {\unicode {x210F}}}\)
\(\newcommand {\angstrom }{\mathrm {\unicode {x212B}}}\)
\(\let \LWRorigbar \bar \)
\(\newcommand {\bar }{\mathrm {bar}}\)
\(\newcommand {\barn }{\mathrm {b}}\)
\(\newcommand {\bel }{\mathrm {B}}\)
\(\newcommand {\decibel }{\mathrm {dB}}\)
\(\newcommand {\knot }{\mathrm {kn}}\)
\(\newcommand {\mmHg }{\mathrm {mmHg}}\)
\(\newcommand {\nauticalmile }{\mathrm {M}}\)
\(\newcommand {\neper }{\mathrm {Np}}\)
\(\newcommand {\yocto }{\mathrm {y}}\)
\(\newcommand {\zepto }{\mathrm {z}}\)
\(\newcommand {\atto }{\mathrm {a}}\)
\(\newcommand {\femto }{\mathrm {f}}\)
\(\newcommand {\pico }{\mathrm {p}}\)
\(\newcommand {\nano }{\mathrm {n}}\)
\(\newcommand {\micro }{\mathrm {\unicode {x00B5}}}\)
\(\newcommand {\milli }{\mathrm {m}}\)
\(\newcommand {\centi }{\mathrm {c}}\)
\(\newcommand {\deci }{\mathrm {d}}\)
\(\newcommand {\deca }{\mathrm {da}}\)
\(\newcommand {\hecto }{\mathrm {h}}\)
\(\newcommand {\kilo }{\mathrm {k}}\)
\(\newcommand {\mega }{\mathrm {M}}\)
\(\newcommand {\giga }{\mathrm {G}}\)
\(\newcommand {\tera }{\mathrm {T}}\)
\(\newcommand {\peta }{\mathrm {P}}\)
\(\newcommand {\exa }{\mathrm {E}}\)
\(\newcommand {\zetta }{\mathrm {Z}}\)
\(\newcommand {\yotta }{\mathrm {Y}}\)
\(\newcommand {\percent }{\mathrm {\%}}\)
\(\newcommand {\meter }{\mathrm {m}}\)
\(\newcommand {\metre }{\mathrm {m}}\)
\(\newcommand {\gram }{\mathrm {g}}\)
\(\newcommand {\kg }{\kilo \gram }\)
\(\newcommand {\of }[1]{_{\mathrm {#1}}}\)
\(\newcommand {\squared }{^2}\)
\(\newcommand {\square }[1]{\mathrm {#1}^2}\)
\(\newcommand {\cubed }{^3}\)
\(\newcommand {\cubic }[1]{\mathrm {#1}^3}\)
\(\newcommand {\per }{/}\)
\(\newcommand {\celsius }{\unicode {x2103}}\)
\(\newcommand {\fg }{\femto \gram }\)
\(\newcommand {\pg }{\pico \gram }\)
\(\newcommand {\ng }{\nano \gram }\)
\(\newcommand {\ug }{\micro \gram }\)
\(\newcommand {\mg }{\milli \gram }\)
\(\newcommand {\g }{\gram }\)
\(\newcommand {\kg }{\kilo \gram }\)
\(\newcommand {\amu }{\mathrm {u}}\)
\(\newcommand {\nm }{\nano \metre }\)
\(\newcommand {\um }{\micro \metre }\)
\(\newcommand {\mm }{\milli \metre }\)
\(\newcommand {\cm }{\centi \metre }\)
\(\newcommand {\dm }{\deci \metre }\)
\(\newcommand {\m }{\metre }\)
\(\newcommand {\km }{\kilo \metre }\)
\(\newcommand {\as }{\atto \second }\)
\(\newcommand {\fs }{\femto \second }\)
\(\newcommand {\ps }{\pico \second }\)
\(\newcommand {\ns }{\nano \second }\)
\(\newcommand {\us }{\micro \second }\)
\(\newcommand {\ms }{\milli \second }\)
\(\newcommand {\s }{\second }\)
\(\newcommand {\fmol }{\femto \mol }\)
\(\newcommand {\pmol }{\pico \mol }\)
\(\newcommand {\nmol }{\nano \mol }\)
\(\newcommand {\umol }{\micro \mol }\)
\(\newcommand {\mmol }{\milli \mol }\)
\(\newcommand {\mol }{\mol }\)
\(\newcommand {\kmol }{\kilo \mol }\)
\(\newcommand {\pA }{\pico \ampere }\)
\(\newcommand {\nA }{\nano \ampere }\)
\(\newcommand {\uA }{\micro \ampere }\)
\(\newcommand {\mA }{\milli \ampere }\)
\(\newcommand {\A }{\ampere }\)
\(\newcommand {\kA }{\kilo \ampere }\)
\(\newcommand {\ul }{\micro \litre }\)
\(\newcommand {\ml }{\milli \litre }\)
\(\newcommand {\l }{\litre }\)
\(\newcommand {\hl }{\hecto \litre }\)
\(\newcommand {\uL }{\micro \liter }\)
\(\newcommand {\mL }{\milli \liter }\)
\(\newcommand {\L }{\liter }\)
\(\newcommand {\hL }{\hecto \liter }\)
\(\newcommand {\mHz }{\milli \hertz }\)
\(\newcommand {\Hz }{\hertz }\)
\(\newcommand {\kHz }{\kilo \hertz }\)
\(\newcommand {\MHz }{\mega \hertz }\)
\(\newcommand {\GHz }{\giga \hertz }\)
\(\newcommand {\THz }{\tera \hertz }\)
\(\newcommand {\mN }{\milli \newton }\)
\(\newcommand {\N }{\newton }\)
\(\newcommand {\kN }{\kilo \newton }\)
\(\newcommand {\MN }{\mega \newton }\)
\(\newcommand {\Pa }{\pascal }\)
\(\newcommand {\kPa }{\kilo \pascal }\)
\(\newcommand {\MPa }{\mega \pascal }\)
\(\newcommand {\GPa }{\giga \pascal }\)
\(\newcommand {\mohm }{\milli \ohm }\)
\(\newcommand {\kohm }{\kilo \ohm }\)
\(\newcommand {\Mohm }{\mega \ohm }\)
\(\newcommand {\pV }{\pico \volt }\)
\(\newcommand {\nV }{\nano \volt }\)
\(\newcommand {\uV }{\micro \volt }\)
\(\newcommand {\mV }{\milli \volt }\)
\(\newcommand {\V }{\volt }\)
\(\newcommand {\kV }{\kilo \volt }\)
\(\newcommand {\W }{\watt }\)
\(\newcommand {\uW }{\micro \watt }\)
\(\newcommand {\mW }{\milli \watt }\)
\(\newcommand {\kW }{\kilo \watt }\)
\(\newcommand {\MW }{\mega \watt }\)
\(\newcommand {\GW }{\giga \watt }\)
\(\newcommand {\J }{\joule }\)
\(\newcommand {\uJ }{\micro \joule }\)
\(\newcommand {\mJ }{\milli \joule }\)
\(\newcommand {\kJ }{\kilo \joule }\)
\(\newcommand {\eV }{\electronvolt }\)
\(\newcommand {\meV }{\milli \electronvolt }\)
\(\newcommand {\keV }{\kilo \electronvolt }\)
\(\newcommand {\MeV }{\mega \electronvolt }\)
\(\newcommand {\GeV }{\giga \electronvolt }\)
\(\newcommand {\TeV }{\tera \electronvolt }\)
\(\newcommand {\kWh }{\kilo \watt \hour }\)
\(\newcommand {\F }{\farad }\)
\(\newcommand {\fF }{\femto \farad }\)
\(\newcommand {\pF }{\pico \farad }\)
\(\newcommand {\K }{\mathrm {K}}\)
\(\newcommand {\dB }{\mathrm {dB}}\)
\(\newcommand {\kibi }{\mathrm {Ki}}\)
\(\newcommand {\mebi }{\mathrm {Mi}}\)
\(\newcommand {\gibi }{\mathrm {Gi}}\)
\(\newcommand {\tebi }{\mathrm {Ti}}\)
\(\newcommand {\pebi }{\mathrm {Pi}}\)
\(\newcommand {\exbi }{\mathrm {Ei}}\)
\(\newcommand {\zebi }{\mathrm {Zi}}\)
\(\newcommand {\yobi }{\mathrm {Yi}}\)
\(\require {mhchem}\)
\(\require {cancel}\)
\(\newcommand {\fint }{âĺŊ}\)
\(\newcommand {\hdots }{\cdots }\)
\(\newcommand {\mathnormal }[1]{#1}\)
\(\newcommand {\vecs }[2]{\vec {#1}_{#2}}\)
\(\renewcommand {\O }{\ensuremath {\mathcal {O}}}\)
\(\newcommand {\vol }{\ensuremath {\operatorname {vol}}}\)
\(\renewcommand {\i }{\ensuremath {\mathrm {i}}}\)
\(\newcommand {\FFT }{\ensuremath {\operatorname {FFT}}}\)
\(\newcommand {\IFFT }{\ensuremath {\operatorname {IFFT}}}\)
\(\newcommand {\diag }{\ensuremath {\operatorname {diag}}}\)
\(\newcommand {\grad }{\ensuremath {\operatorname {grad}}}\)
\(\newcommand {\const }{\ensuremath {\operatorname {const}}}\)
\(\newcommand {\name }[1]{\textsc {#1}}\)
\(\newcommand {\smallpmatrix }[1]{\left (\begin {smallmatrix}#1\end {smallmatrix}\right )}\)
\(\newcommand {\matlab }{{\fontfamily {bch}\scshape \selectfont {}Matlab}}\)
\(\newcommand {\innerproduct }[1]{\left \langle {#1}\right \rangle }\)
\(\newcommand {\norm }[1]{\left \Vert {#1}\right \Vert }\)
\(\renewcommand {\natural }{\mathbb {N}}\)
\(\newcommand {\integer }{\mathbb {Z}}\)
\(\newcommand {\rational }{\mathbb {Q}}\)
\(\newcommand {\real }{\mathbb {R}}\)
\(\newcommand {\complex }{\mathbb {C}}\)
\(\renewcommand {\d }{\mathop {}\!\mathrm {d}}\)
\(\newcommand {\dr }{\d {}r}\)
\(\newcommand {\ds }{\d {}s}\)
\(\newcommand {\dt }{\d {}t}\)
\(\newcommand {\du }{\d {}u}\)
\(\newcommand {\dv }{\d {}v}\)
\(\newcommand {\dw }{\d {}w}\)
\(\newcommand {\dx }{\d {}x}\)
\(\newcommand {\dy }{\d {}y}\)
\(\newcommand {\dz }{\d {}z}\)
\(\newcommand {\dsigma }{\d {}\sigma }\)
\(\newcommand {\dphi }{\d {}\phi }\)
\(\newcommand {\dvarphi }{\d {}\varphi }\)
\(\newcommand {\dtau }{\d {}\tau }\)
\(\newcommand {\dxi }{\d {}\xi }\)
\(\newcommand {\dtheta }{\d {}\theta }\)
\(\newcommand {\tp }{\mathrm {T}}\)
Im Folgenden werden Formeln gesucht, mit denen Integrale der Form \(\int _a^b f(x)\dx \) möglichst gut approximiert werden können. Die vorgestellten Formeln haben alle die Form \(\int _a^b
f(x)\dx \approx \sum _{k=1}^n w_k f(x_k)\). Nun ist die Wahl der Stützstellen \(x_k\) und Gewichte \(w_k\) maßgeblich, d. h. in diesen Parametern unterscheiden sich die verschiedenen
Methoden. Allen gemeinsam sind dabei folgende Forderungen, die man an die Approximationsformeln stellt:
Die Formel sollte exakt für \(f \equiv 1\) sein, d. h. \(b - a= \sum _{k=1}^n w_k\).
Die Gewichte \(w_k > 0\) sollten positiv sein, da man sonst Funktionen definieren könnte, die nur bei einem negativen Gewicht größer Null und sonst überall Null sind. Die Approximation
würde einen negativen Wert ergeben, was offensichtlich sinnlos ist.
Die Formel sollte exakt für Polynome \(f\) von möglichst hohem Grad sein.
Die Gauß-Formel der Ordnung \(n\) approximiert das Integral einer Funktion \(f\) durch das Integral des Interpolationspolynoms an den Nullstellen \(x_1 < \dotsb
< x_n\) des Legendre-Polynoms vom Grad \(n\), d. h.
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _{-1}^1 f(x)\dx \approx \sum _{i=1}^n w_i f(x_i)
\end{align*}
mit \(w_i\) den Integralen der Lagrange-Polynome über das Intervall \([-1, 1]\):
\(\seteqnumber{0}{}{0}\)
\begin{align*}
w_i := \int _{-1}^1 p_i(x)\dx , \quad p_i(x) := \prod _{j \not = i} \frac {x - x_j}{x_i - x_j}.
\end{align*}
Die Formel ist exakt für Polynome vom Grad \(< 2n\) und vor allem für analytische Funktionen sehr genau. Alle Gewichte \(w_i\) sind positiv und die Stützstellen \(x_i\) liegen im
Integrationsintervall \((-1, 1)\).
Gauß-Parameter \(x’\) und \(w’\) für ein beliebiges Integrationsintervall \([a, b]\) erhält man durch lineare Transformation:
\(\seteqnumber{0}{}{0}\)
\begin{align*}
x_k’ = a + \frac {b - a}{2} (x_k + 1), \quad w_k’ = \frac {b - a}{2} w_k.
\end{align*}
Konvergenz der Gauß-Quadratur
Für eine stetige Funktion \(f\) konvergieren die Approximationen der Gauß-Quadratur für ein Integrationsintervall \([a, b]\) mit wachsender Zahl \(n\) der Knoten gegen \(\int _a^b
f(x)\dx \), d. h.
\(\seteqnumber{0}{}{0}\)
\begin{align*}
s_n f := \sum _{i=1}^n w_i^n f(x_i^n) \xrightarrow {n \to \infty } \int _a^b f(x)\dx .
\end{align*}
Fehler der Gauß-Quadratur
Der Fehler der Gauß-Quadratur besitzt die Darstellung
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\sum _{i=1}^n w_i f(x_i) - \int _a^b f(x)\dx = -\gamma _n f^{(2n)}(\xi )(b - a)^{2n + 1}
\end{align*}
mit
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\gamma _n := \frac {(n!)^4}{(2n + 1)((2n)!)^3},\quad \xi \in [a, b].
\end{align*}
Der Kehrwert der Fehlerkonstanten \(\gamma _n\) besitzt für \(n = 1, 2, 3, 4\) die Werte \(24, 4320, 2016000,\) \(1778112000\), d. h. \(\gamma _n\) wird schnell sehr klein.
Gewichtete Gauß-Quadratur
Seien \(x_i\) die Nullstellen des Orthogonalpolynoms vom Grad \(n\) zu einer Gewichtsfunktion \(w\) auf einem Intervall \([a, b]\). Dann ist die auf polynomialer Interpolation basierende Quadraturformel
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _a^b f(x)w(x)\dx \approx \sum _{i=1}^n w_i f(x_i)
\end{align*}
für Polynome vom Grad \(< 2n\) exakt (gewichtete Gauß-Quadratur). Die Gewichte \(w_i\) sind positiv und können als Integrale der
Lagrange-Polynome berechnet werden:
\(\seteqnumber{0}{}{0}\)
\begin{align*}
w_i := \int _a^b p_k(x)w(x)\dx , \quad p_i(x) := \prod _{j \not = i} \frac {x - x_j}{x_i - x_j}
\end{align*}
Der Fehler ist gleich \(\gamma _n f^{(2n)}(\xi )\) für ein \(\xi \in [a, b]\) mit einer von der Gewichtsfunktion abhängigen Konstanten \(\gamma _n\).
Die folgende Tabelle zeigt die Parameter und Gewichtsfunktionen für die klassischen Orthogonalpolynome:
Typ |
\([a, b]\) |
\(w(t)\) |
\(\gamma _n\) |
Legendre |
\([-1, 1]\) |
\(1\) |
\(\frac {2^{2n + 1} (n!)^4}{(2n + 1)((2n)!)^3}\) |
Tschebyscheff |
\([-1, 1]\) |
\(\sqrt {1 - t^2}\) |
\(\frac {\pi }{2^{2n - 1} (2n)!}\), \(n > 0\) |
Jacobi |
\([-1, 1]\) |
\((1 + t)^r (1 - t)^s\) |
\(\frac {2^{2n + r + s + 1} n! \Gamma (n + r + 1) \Gamma (n + s + 1) \Gamma (n + r + s + 1)}{(2n + r + s + 1) \Gamma (2n + r + s + 1)^2
(2n)!}\) |
Laguerre |
\([0, \infty )\) |
\(\exp (-t)\) |
\(\frac {(n!)^2}{(2n)!}\) |
Hermite |
\((-\infty , \infty )\) |
\(\exp (-t^2)\) |
\(\frac {\sqrt {\pi } n!}{2^n (2n)!}\) |
|
|
|
|
Trapezregel
Die Näherung
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _a^b f(x)\dx \approx s_h f := h \left (\frac {f(a)}{2} + f(a + h) + \dotsb + f(b - h) + \frac {f(b)}{2}\right )
\end{align*}
approximiert das Integral durch Summe von Trapezflächen (Trapezregel).
Für eine zweimal stetig differenzierbare Funktion gilt für den Fehler
\(\seteqnumber{0}{}{0}\)
\begin{align*}
s_h f - \int _a^b f(x)\dx = \frac {b - a}{12} f”(\xi ) h^2
\end{align*}
für ein \(\xi \in [a, b]\).
Genauer besitzt der Fehler für glatte Funktionen die asymptotische Entwicklung
\(\seteqnumber{0}{}{0}\)
\begin{align*}
s_h f - \int _a^b f(x)\dx = c_1 (f’(b) - f’(a)) h^2 + c_2 (f”’(b) - f”’(a)) h^4 + \dotsb
\end{align*}
mit von \(f\) und \(h\) unabhängigen Konstanten \(c_j\). Daraus folgt, dass die Trapezregel für \((b - a)\)-periodische Funktionen sehr genau ist. Der Fehler strebt schneller als jede \(h\)-Potenz gegen
Null.
Diese Fehlerformel besagt insbesondere, dass für glatte periodische Integranden schnell eine hohe Genauigkeit erzielt werden kann.
Bernoulli-Polynome
Die normalisierten Bernoulli-Polynome sind definiert durch die Rekursion
\(\seteqnumber{0}{}{0}\)
\begin{align*}
p_i’ := p_{i-1}, \quad p_0(x) := 1 \quad \text {mit} \quad \int _0^1 p_i(x)\dx = 0, \quad i \in \natural .
\end{align*}
Die normalisierten Bernoulli-Polynome sind symmetrisch bzgl. \(t = \frac {1}{2}\), d. h. \(p_{2i-1}(x - \frac {1}{2})\) und \(p_{2i}(x - \frac {1}{2})\) sind ungerade bzw. gerade
Funktionen (\(i \in \natural \)). Für \(i \ge 2\) gilt außerdem
\(\seteqnumber{0}{}{0}\)
\begin{align*}
& p_{2i-1}(0) = p_{2i-1}(1/2) = p_{2i-1}(1) = 0, \\ & p_{2i-1}(t) \not = 0 \quad \text {fÃijr}\quad t \in (0, 1) \setminus \left \{\tfrac
{1}{2}\right \}
\end{align*}
und für \(i \ge 1\) ist
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\gamma _{2i} := p_{2i}(0) = p_{2i}(1)
\end{align*}
entweder ein Minimum oder ein Maximum von \(p_{2i}\) auf \([0, 1]\).
Die Werte \(\gamma _{2i}\) heißen normierte Bernoulli-Zahlen.
Euler-Maclaurin-Entwicklung
Für eine glatte Funktion \(f\) hat der Fehler \(e_h f := s_h f - \int _a^b f(x)\dx \) der Trapezregel die Entwicklung
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\sum _{i=1}^{m-1} \gamma _{2i} (f^{(2i-1)}(b) - f^{(2i-1)}(a)) h^{2i}
\end{align*}
mit dem Restglied
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\gamma _{2m} f^{(2m)}(\xi ) (b - a) h^{2m}
\end{align*}
für ein \(\xi \in [a, b]\) und \(\gamma _{2i}\) den normierten Bernoulli-Zahlen
(Euler-Maclaurin-Entwicklung).
Aus der Entwicklung folgt insbesondere, dass die Trapezregel für unendlich oft differenzierbare \((b - a)\)-periodische Funktionen sehr genau ist. Der Fehler strebt schneller als jede \(h\)-Potenz gegen Null.
Für nicht-periodische Funktionen bildet die Entwicklung die Grundlage für Extrapolationsverfahren, mit denen ebenfalls beliebige Approximationsordnungen erzielt werden können.
Romberg-Algorithmus
Die Genauigkeit der Trapezregel
\(\seteqnumber{0}{}{0}\)
\begin{align*}
s_h^1 := h \left (\frac {f(a)}{2} + f(a + h) + \dotsb + f(b - h) + \frac {f(b)}{2}\right )
\end{align*}
lässt sich durch Extrapolation verbessern (Romberg-Algorithmus).
Die rekursiv definierten Approximationen
\(\seteqnumber{0}{}{0}\)
\begin{align*}
s_h^{j+1} := \frac {4^j s_{h/2}^j - s_h^j}{4^j - 1}
\end{align*}
haben die Fehlerordnung \(\O (h^{2j+2})\) und können in einem Dreiecksschema berechnet werden:
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\begin{array}{ccccc} s_h^1 & \rightarrow & s_h^2 & \rightarrow & s_h^3 \\ & \nearrow & & \nearrow \\ s_{h/2}^1 &
\rightarrow & s_{h/2}^2 \\ & \nearrow \\ s_{h/4}^1 \end {array}
\end{align*}
Es werden solange sukzessive Diagonalen \(s_{2^{-m}}^1, s_{2^{1-m}}^2, \dotsc , s_h^{m+1}\) hinzugefügt, bis mit dem zuletzt generierten Wert die gewünschte Genauigkeit erreicht ist. Bei
den Trapezsummen können bereits berechnete Funktionswerte genutzt werden:
\(\seteqnumber{0}{}{0}\)
\begin{align*}
s_{h/2}^1 = \frac {1}{2} \left (s_h^1 + h \left (f(a + \tfrac {h}{2}) + \dotsb + f\left (b - \tfrac {h}{2}\right )\right )\right ).
\end{align*}
Numerische Integration mit MATLAB
Das Integral \(s = \int _a^b f(x)\dx \) kann in M ATLAB mit dem Befehl s = quad(f, a, b, tol); berechnet werden, wobei die zu integrierende
Funktion f als Funktionshandle, Funktionsname (String) oder Inlinefunktion übergeben wird. tol ist optional (standardmäßig \(10^{-6}\)) und gibt die absolute Genauigkeit vor.
quad basiert dabei auf der Simpson-Regel (abschnittsweise Interpolation mit quadratischen Polynomen und anschließende Integration) mit adaptiver Unterteilung des
Integrationsintervalls, d. h. die Intervalllängen bzw. die Anzahl an Funktionsauswertungen werden an die lokale Komplexität der Funktion angepasst.
Alternativ zu quad gibt es noch die Befehle quadl und dblquad. quadl erzielt eine höhere Approximationsordnung als quad, sollte aber nur bei hohen Genauigkeiten und glatten Integranden verwendet werden. Bei niedrigen Genauigkeiten oder nicht-glatten Integranden empfiehlt sich daher die
Verwendung von quad. dblquad berechnet ein bivariates Integral, d. h. s = dblquad(f, xmin, xmax, ymin, ymax, tol); berechnet das Integral der bivariaten Funktion f
über das Rechteck \([\)xmin\(,\;\)xmax\(] \times [\)ymin\(,\;\)ymax\(]\).
Mehrfachintegrale
Integrationsformeln für Rechteck-Gebiete
\(\seteqnumber{0}{}{0}\)
\begin{align*}
Q := [a_1, b_1] \times \dotsb \times [a_m, b_m]
\end{align*}
erhält man durch Bilden von Tensorprodukten eindimensionaler Quadraturformeln.
Sind die Formeln \(\sum _k w_{k,\nu } f(t_{k,\nu })\) zur Approximation von \(\int _{a_\nu }^{b_\nu } f\) exakt für Polynome vom Grad \(\le n_\nu \), \(\nu = 1, \dotsc , m\),
so ist die Produktformel
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _Q f \approx \sum _{k_1} \dotsb \sum _{k_m} (w_{k_1,1} \dotsb w_{k_m,m}) f(t_{k_1,1}, \dotsc , t_{k_m,m})
\end{align*}
exakt für Polynome vom Koordinatengrad \(\le (n_1, \dotsc , n_m)\).
Für die Trapezregel mit Schrittweite \(h = \frac {b - a}{n}\) sind die Gewichte an den inneren Knoten gleich \(h\) und an den äußeren Knoten \(a\) und \(b\) gleich \(\frac {h}{2}\).
Für ein Rechteck \([a_1, b_1] \times [a_2, b_2]\) und den Schrittweiten \(h_1\) und \(h_2\) erhält man somit drei verschiedene Gewichte: \(h_1 h_2\) für innere Knoten, \(\frac
{h_1 h_2}{2}\) für Knoten auf den Kanten (außer den Ecken) und \(\frac {h_1 h_2}{4}\) für die Ecken.
Allgemein gilt für ein \(m\)-dimensionales Rechteck mit Schrittweiten \(h_\nu \), \(\nu = 1, \dotsc , m\), dass \(w_{k_1, \dotsc , k_m} = h^m 2^{-\alpha _1} \dotsm 2^{-\alpha
_m}\) mit \(\alpha _\nu = 0\) für einen inneren Index bzw. \(\alpha _\nu = 1\) für einen äußeren Index \(k_\nu \), \(\nu = 1, \dotsc , m\).
Der Fehler der multivariaten Trapezregel besitzt eine quadratische Entwicklung, sodass wie für die univariate Formel die Romberg-Extrapolation anwendbar ist.
Für eine bijektive, stetig differenzierbare Transformation \(g\) eines regulären Bereiches \(U \subset \real ^n\) mit \(\det g’(x) \not = 0\) für \(x \in U\) gilt für stetige
Funktionen \(f\)
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _U f \circ g |\det g’| dU = \int _V f dV, \quad V = g(U),
\end{align*}
wobei \(\det g’\) als Funktionaldeterminante der Transformation bezeichnet wird.
Eine Integrationsformel
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _D f \approx \sum _k w_k f(x_k)
\end{align*}
kann man durch Variablensubstitution auf andere Gebiete transformieren.
Ist \(\varphi \colon D \rightarrow \widetilde {D}\) eine bijektive Abbildung, erhält man durch
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\widetilde {w}_k = \left |\det \varphi ’(x_k)\right | w_k, \quad \widetilde {x}_k = \varphi (x_k)
\end{align*}
Gewichte und Punkte für eine Integrationsformel auf \(\widetilde {D}\). Dabei ist \(\varphi ’(x_k)\) die Jacobi-Matrix von \(\varphi \) im Punkt \(x_k\). \(\det \varphi ’(x_k)\) heißt
Funktionaldeterminante der Transformation.
Speziell gilt bei einer affinen Abbildung \(\varphi (x) = Ax + b\) für die Gewichte \(\widetilde {w}_i = |\det A| w_i\).
Die Konvergenz der Integrationsformeln bleibt bei glatten Transformationen \(\varphi \) erhalten. Allerdings werden Polynome auf \(\widetilde {D}\) nicht mehr exakt integriert, da der Integrand die
Funktionaldeterminante der Transformation enthält.
Beispielsweise ist ein durch eine Funktion \(h\) berandeter Integrationsbereich
\(D = \{(x, y) \in \real ^2 \;|\; x \in [a, b],\; y \in [0, h(x)]\}\) als Bild des Rechtecks \(Q = [a, b] \times [0, 1]\) unter der Abbildung \(\varphi \colon Q \rightarrow
D\), \((u, v) \mapsto (u, v h(u))\) mit der Funktionaldeterminante
\(\det \varphi ’(u, v) = \det \begin {pmatrix}1 & 0\\vh’(u) & h(u)\end {pmatrix} = h(u)\) darstellbar. Damit transformieren sich Punkte und Gewichte einer Produktformel
für \(Q\) gemäß \((u_i, v_j) \rightarrow (u_i, v_j h(u_i))\) und \(w_{i,j} \rightarrow w_{i,j} h(u_i)\).
Für einen \(m\)-dimensionalen Simplex \(S\) mit den Ecken \(v_0, \dotsc , v_m\) lässt sich eine Integrationsformel durch Interpolation mit Polynomen vom totalen Grad (d. h. die
Summe der Koordinatengrade) \(\le n\) konstruieren. Das interpolierende Polynom ist eindeutig durch die Werte an den Punkten
\(\seteqnumber{0}{}{0}\)
\begin{align*}
x_k = \frac {1}{n} \sum _{\nu =0}^n k_\nu v_\nu , \quad k_0 + \dotsb + k_n = n
\end{align*}
bestimmt, die ein regelmäßiges Gitter bilden.
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _S f \approx \vol S \sum _{k_0 + \dotsb + k_n = n} w_k f(x_k)
\end{align*}
ist eine Approximation der Ordnung \(n + 1\), d. h. die Integrationsformel ist exakt für alle Polynome vom totalen Grad \(\le n\) (Integrationsformel für Simplizes). Dabei sind \(\vol
S \cdot w_k\) die Integrale über die Lagrange-Polynome zu \(x_k\).
Für allgemeine Gebiete kann die Integrationsformel auf den Simplizes einer Triangulierung angewendet werden. Der Fehler hat dann die Ordnung \(\O (h^{n+1})\), wobei \(h\) den maximalen Durchmesser der
Teilsimplizes bezeichnet.
Monte-Carlo-Verfahren
Lineare Kongruenzmethode
Die lineare Kongruenzmethode definiert durch
\(\seteqnumber{0}{}{0}\)
\begin{align*}
n_\ell & := \alpha n_{\ell -1} \mod \beta \\ x_\ell & := n_\ell / \beta
\end{align*}
mit \(\alpha \in \natural \), \(1 < \alpha < \beta \) und \(\beta \) einer sehr großen Primzahl kann bei geeigneter Wahl der Parameter zur numerischen Simulation von Zufallszahlen \(x_\ell \in
\left [0, 1\right )\) benutzt werden.
Eine minimale Anforderung ist, dass die maximale Periode \(\beta - 1\) erreicht wird. Darüber hinaus soll die Folge \(x_0, x_1, \dotsc \) bei möglichst vielen statistischen Tests gute Ergebnisse
liefern.
Satz von Fermat
Für jede Primzahl \(\beta \) und \(\alpha \not = 0 \mod \beta \) gilt der Satz von Fermat
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\alpha ^{\beta -1} = 1 \mod \beta .
\end{align*}
Maximale Periode bei der linearen Kongruenzmethode
Für eine Primzahl \(\beta \) hat die Folge \(\alpha ^\ell \mod \beta \), \(\ell = 0, 1, \dotsc \) keine kleinere Periode als \(\beta - 1\) genau dann, wenn
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\alpha ^{(\beta -1)/m} \not = 1 \mod \beta
\end{align*}
für alle Primteiler \(m\) von \(\beta - 1\). (Die Periode ist dabei die Anzahl der vorkommenden verschiedenen Zahlen.)
Mithilfe dieses Kriteriums lassen sich geeignete Multiplikatoren \(\alpha \) für die Simulation von Zufallszahlen mit der linearen Kongruenzmethode bestimmen.
Gebräuliche Parameter bei der Generierung von Pseudo-Zufallszahlen mit der linearen Kongruenzmethode sind die Mersenne-Primzahl \(\beta = 2147483647 =
2^{31} - 1\) und der Multiplikator \(\alpha = 16807\). Aufgrund der Größe dieser Zahlen ist das Testen der Periodenlänge nicht ganz einfach.
Aufgrund der Primfaktorzerlegung \(\beta - 1 = 2 \cdot 3^2 \cdot 7 \cdot 11 \cdot 31 \cdot 151 \cdot 331\) erhält man als mögliche Perioden \(m = (\beta - 1)/p_k\):
\(1073741823, 715827882, 306783378, 195225786, 69273666, 14221746,\) \(6487866\).
Zur Berechnung der Potenzen \(\alpha ^m\) geht man zur Dualdarstellung \(m = m_0 + 2m_1 + 4m_2 + \dotsb \), \(m_k \in \{0, 1\}\) über. Man kann dann zunächst rekursiv die
Potenzen \(\alpha _k := \alpha ^{2^k} \mod \beta \) berechnen, indem man \(\alpha _{k+1} \mod \beta = \alpha _k^2 \mod \beta \) rechnet. Damit ist \(\alpha ^m \mod \beta =
(\prod _{m_k=1} \alpha _k) \mod \beta \).
Zum Beispiel ist \(195225786 = (001011101000101110100010111010)_2\) und man erhält
\(16807^{195225786} \mod \beta = 997852928 \not = 1\). Ebenso sind alle sechs anderen zu testenden Potenzen ungleich \(1\) modulo \(\beta \). Damit ist die Periode des Mersenne-Generators maximal.
Spektraltest für die lineare Kongruenzmethode
Für eine Primzahl \(\beta \) lässt sich die durch
\(\seteqnumber{0}{}{0}\)
\begin{align*}
u_\ell := (\alpha ^{\ell m} a \mod \beta ) / \beta , \qquad a := (1, \alpha , \dotsc , \alpha ^{m-1}) \in \real ^m
\end{align*}
definierte Folge von Vektoren \(u_\ell \) durch parallele Hyperebenen im Abstand
\(\seteqnumber{0}{}{0}\)
\begin{align*}
d := (\min \{\norm {n}_2 \;|\; n \in \integer ^m,\; \norm {n}_2 \not = 0,\; a^t n = 0 \mod \beta \})^{-1}
\end{align*}
überdecken.
Der Abstand \(d\) dient zur Beurteilung der Güte der Folge der Pseudo-Zufallsvektoren \(u_\ell \) (Spektraltest). Je kleiner \(d\) ist, um so besser sind im Allgemeinen die statistischen Eigenschaften der
Folge.
Gleichverteilte Folgen
Eine Folge \(x_0, x_1, \dotsc \) in einem Quader
\(\seteqnumber{0}{}{0}\)
\begin{align*}
Q := [a_1, b_1] \times \dotsb \times [a_n, b_n]
\end{align*}
heißt gleichverteilt, falls
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\lim _{\ell \to \infty } \frac {|\{x_k \in Q’ \;|\; k < \ell \}|}{\ell } = \frac {\vol Q’}{\vol Q}
\end{align*}
für alle Teilquader \(Q’ \subseteq Q\).
Allgemeiner heißt die Folge \(m\)-verteilt, falls
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\lim _{\ell \to \infty } \frac {|\{x_{k+1} \in Q’_1, \dotsc , x_{k+m} \in Q’_m \;|\; k < \ell \}|}{\ell } = \frac {\vol Q’_1 \dotsm \vol Q’_m}{\vol Q}
\end{align*}
für alle Teilquader \(Q’_1, \dotsc , Q’_m \subseteq Q\).
Eine Folge, die für alle \(m \in \natural \) \(m\)-verteilt ist, heißt \(\infty \)-verteilt.
Franklin hat gezeigt, dass die Folge
\(\seteqnumber{0}{}{0}\)
\begin{align*}
x_\ell := r^\ell \mod 1,\quad \ell = 0, 1, \dotsc
\end{align*}
für fast alle \(r > 1\) \(\infty \)-verteilt in \([0, 1)\) ist.
Konvergenz der Monte-Carlo-Integration
Für eine gleichverteilte Folge \(x_0, x_1, \dotsc \) in \([0, 1)\) gilt
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\lim _{\ell \to \infty } \frac {1}{\ell } \sum _{k<\ell } f(x_k) = \int _0^1 f
\end{align*}
für jede Riemann-integrierbare Funktion \(f\). Das entsprechende Approximationsverfahren wird aufgrund der quasi zufälligen Wahl der Punkte \(x_k\) als Monte-Carlo-Integration bezeichnet.
Ausgehend von einer gleichverteilten Folge \(x_0, x_1, \dotsc \) in dem Standardintervall \([0, 1)\) können gleichverteilte Zahlen und Vektorfolgen auf allgemeineren Mengen konstruiert werden.
Für \(k = 0, 1, \dotsc \) ist
\(a + (b - a)x_k\) eine Folge in \([a, b)\),
\(\lfloor m + (n - m)x_k \rfloor \) eine Folge in \(\{m, \dotsc , n - 1\}\),
\(U_k := (x_{mk}, \dotsc , x_{mk + m - 1})\) eine Folge in \([0, 1)^m\) und
\(A U_k + b\) mit einer \(m \times m\)-Matrix \(A\) und einem \(m\)-Vektor \(b\) eine Folge in \(A[0,1)^m + b\).
Außerdem ist für eine gleichverteilte Folge in einem Quader \(Q\) die in einer Teilmenge \(D \subseteq Q\) liegende Teilfolge gleichverteilt in \(D\).
Sei \(x_0, x_1, \dotsc \) eine gleichverteilte Folge in \([0, 1)\).
Zur Simulation von Würfeln kann die Folge \(n_k = 1 + 6 \lfloor x_k \rfloor \) verwendet werden.
Eine Folge im abgeschlossenen Einheitskreis \(D\colon u^2 + v^2 \le 1\) erhält man durch die Transformation \((u_k, v_k) = 2(x_{2k}, x_{2k+1}) - (1, 1)\) und Auswahl der Teilfolge, die in
\(D\) liegt.
Gleichverteilte Permutationen von \(\{1, \dotsc , n\}\), etwa zum Simulieren des Mischens von Karten, können durch die Reihenfolge der Indizes beim Sortieren der Komponenten gleichverteilter \(n\)-Vektoren
generiert werden.
Multivariate Monte-Carlo-Integration
Ist \(u_0, u_1, \dotsc \) eine gleichverteilte Folge in einem Quader \(Q\), so lässt sich ein Integral über ein Gebiet \(D \subseteq Q\) durch
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\int _D f = \lim _{\ell \to \infty } \frac {\vol Q}{\ell } \sum _{k < \ell ,\; u_k \in D} f(u_k)
\end{align*}
approximieren (multivariate Monte-Carlo-Integration).
Im Spezialfall \(f = 1\) erhält man ein Verfahren zur Volumenbestimmung:
\(\seteqnumber{0}{}{0}\)
\begin{align*}
\vol D = \lim _{\ell \to \infty } \frac {\vol Q}{\ell } |\{u_k \in D \;|\; k < \ell \}|.
\end{align*}
Der Vorteil multivariater Monte-Carlo-Integrationen ist die weitgehende Unabhängigkeit der Konvergenzrate von der Dimension. Monte-Carlo-Verfahren sind daher besonders in hohen Dimension gut geeignet (z. B.
besser als Trapezregel).