Programování pro fyziky 2021/22 – Numerické recepty č. 3

Náš průvodní text je především v poznámkách k přednášce. Pomáhat si budeme i odkazy do knih Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR), které jsou (v Second Edition) na webu volně k dispozici. Vidět je ovšem lze jen z Adobe Readeru po instalaci pluginu FileOpen. Interní PDF prohlížeče Firefoxu a Chromu nám nepomohou.

Hádanky

  1. Z jednoho N.-C. kvadraturního vzorce získáme všechny ostatní. Který z nich to je?

  2. Jaké jsou výhody a nevýhody Gaussových kvadraturních vzorců proti N.-C. vzorcům?

  3. Totéž pro integrování metodou Monte Carlo?

Numerické integrování

(NR kap. 4) Na polynomické aproximaci integrandu jsou založeny kvadraturní vzorce neboli vážené součty hodnot integrandu v interpolačních uzlech. Newtonovy-Cotesovy kvadraturní vzorce jsou formulovány na ekvidistantních sítích, kde interpolační polynomy vyšších stupňů příliš oscilují; používají se proto složené N.-C. vzorce, odpovídající polynomické interpolaci po částech (pravidlo lichoběžníkové, Simpsonovo ad.). Přesnější Gaussovy kvadraturní vzorce jsou definovány na sítích tvořených kořeny ortogonálních polynomů. Nejen hazardní hráče přitahuje méně přesná, ale účinná metoda Monte Carlo. Octave poskytuje integrátor quad a několik dalších podobného jména.

Obecný kvadraturní vzorec

Odvozuje se analytickým integrováním Lagrangeova polynomu \(L_n(x)\) interpolujícího integrand na \(n+1\) uzlech,

\[\int_{x_0}^{x_n} f(x)\,dx \approx I_n \equiv \int_{x_0}^{x_n} L_n(x)\,dx = \int_{x_0}^{x_n} \sum_{i=0}^n f_i l_n^{(i)}(x)\,dx = \sum_{i=0}^n w_i f_i,\]

kde

\[w_i = \int_{x_0}^{x_n} l_n^{(i)}(x)\,dx,\ f_i=f(x_i).\]

Newtonovy-Cotesovy kvadraturní vzorce

(NR kap. 4.1) Na intervalu \(\langle a,b\rangle\) se zavede ekvidistantní síť \(x_i=x_0+ih\), kde \(i=0,\dots,n\). Pokud krajní body intervalu jsou interpolačními uzly, \(a=x_0\) a \(b=x_n\), je krok v síti \(h=(b-a)/n\) a výsledný kvadraturní vzorec se nazývá uzavřený, jinak jde o vzorec otevřený. Pokud se interpoluje jedním polynomem \(n\)-tého stupně na \(n+1\) uzlech, říká se kvadraturnímu vzorci jednoduchý, pokud se interpoluje po částech více polynomy, jde o vzorec složený. Analytickým integrováním Lagrangeova interpolačního polynomu se získají následující uzavřené Newtonovy-Cotesovy vzorce rostoucí přesnosti (ve smyslu interpolace integrandu polynomy vyšších stupňů):

jednoduché lichoběžníkové pravidlo (\(n=1\)): \(\displaystyle\int_{x_0}^{x_1} f(x)\,dx \approx \frac12 h (f_0+f_1)\),

neboť \(\displaystyle \int_{x_0}^{x_1}\left(f_0\,l_1^{(0)}+f_1\,l_1^{(1)}\right)dx = \int_{x_0}^{x_1}\left(f_0\frac{x-x_1}{x_0-x_1}+f_1\frac{x-x_0}{x_1-x_0}\right)dx =\frac{f_0}h\left[x_1x-\frac{x^2}2\right]_{x_0}^{x_1}+\frac{f_1}h\left[\frac{x^2}2-x_0x\right]_{x_0}^{x_1} =\frac 12 h (f_0+f_1)\)

jednoduché Simpsonovo pravidlo (\(n=2\)): \(\displaystyle\int_{x_0}^{x_2} f(x)\,dx \approx \frac13 h (f_0+4f_1+f_2)\)

jednoduché tříosminové pravidlo (\(n=3\)): \(\displaystyle\int_{x_0}^{x_3} f(x)\,dx \approx \frac38 h (f_0+3f_1+3f_2+f_3)\)

V případě lichoběžníkového pravidla je integrand aproximován lineárně a plocha pod integrandem je obsahem lichoběžníka, Simpsonovo pravidlo integrand aproximuje kvadraticky, tříosminové pravidlo kubicky.

Složením \(N/n\) jednoduchých vzorců pro \(x_0,\dots,x_N\) (\(N\) sudé pro Simpsonovo pravidlo, dělitelné \(3\) pro \(\frac38\) pravidlo):

složené lichoběžníkové pravidlo: \(L_N\approx h(\frac12 f_0+f_1+f_2+\dots+f_{N-1}+\frac12 f_N)\)

složené Simpsonovo pravidlo: \(S_N\approx \frac13 h\ (f_0+4f_1+2f_2+4f_3+\dots+4f_{N-1}+f_N)\)

složené tříosminové pravidlo: \(T_N\approx \frac38 h\ (f_0+3f_1+3f_2+2f_3+3f_4+\dots+3f_{N-1}+f_N)\)

Jméno má i kvadraturní vzorec založený na (zde pravostranné) aproximaci integrandu konstantním polynomem:

obdélníkové pravidlo: \(\displaystyle\int_{x_0}^{x_1} f(x)\,dx \approx h f_1,\quad O_N \approx h(f_1+f_2+\dots+f_{N-1}+f_N)\)

Složené obdélníkové i lichoběžníkové pravidlo lze snadno použít i na neekvidistantních sítích, \(h_i=x_i-x_{i-1}\): \(O'_N\approx \sum_{i=1}^N h_if_i\,,\) \(L'_N\approx \frac12 \sum_{i=1}^N h_i(f_{i-1}+f_i)\,.\)

Rombergova metoda (NR kap. 4.3) je praktickým postřehem založeným na vzorci pro rozvoj chyby složeného lichoběžníkového pravidla (NR 4.2.1), který praví, že složené Simpsonovo pravidlo na \(N+1\) uzlech lze vyčíslit pomocí lichoběžníkových pravidel na \(N+1\) a \(N/2+1\) uzlech: \(S_N=\frac43 L_N-\frac13 L_{N/2},\) jak lze ověřit dosazením:

\[\frac43L_N-\frac13L_{N/2}=\frac43h\left(\frac{f_0}2+f_1+f_2+\dots+\frac{f_N}2\right) -\frac13 2h\left(\frac{f_0}2+f_2+\dots+\frac{f_N}2\right) =\frac h3\left(f_0+4f_1+2f_2+\dots+f_N\right) =S_N\,.\]

Vyčíslováním \(L_N\) pro \(N=1,\ 2,\ 4,\ 8,\ \dots\ \) tak lze touto a obdobnými kombinacemi výsledky levně zpřesňovat.

Gaussovy kvadraturní vzorce

(NR kap. 4.5) Kvadraturní vzorec \(I_n=\sum w_if_i\) má při cca \(2n\) stupních volnosti (váhy \(w_i\) a uzly \(x_i\)) potenciální řád přesnosti (tj. stupeň polynomu integrovatelného přesně) rovněž cca \(2n\). Triviální volbou cca \(n\) uzlů (ekvidistantní síť N.-C. vzorců) se řád přesnosti sníží o cca \(n\), odpovědnou volbou uzlů se řád přesnosti může uchovat. Jinak řečeno: jednoduché N.-C. vzorce s aproximačním polynomem stupně \(m\) dovolují při \(n+1\) ekvidistantních uzlech integrovat přesně polynomy do stupně \(n\). Gaussovy vzorce s aproximačním polynomem téhož stupně, ale vyčísleným na \(n\) speciálně volených uzlech dovolují integrovat přesně polynomy do stupně \(2n-1\). Vhodnými uzly jsou kořeny ortogonálních polynomů \(p_n(x)\), splňujících požadavek ortogonality s vahou \(w(x)\), \(\int w(x)p_i(x)p_j(x)\,dx=0\) pro \(i \ne j\). Gaussův kvadraturní vzorec je pak \(I_n=\int w(x)f(x)\,dx=\sum w_if_i\), kde váhy \(w_i\) a uzly \(x_i\) pro vyčíslení \(f_i\) jsou dány v učebnicích nebo je lze odvodit. Gaussovy vzorce se mohou zvláště vyplatit, když součást integrandu \(w(x)\) je vahou klasických ortogonálních polynomů a zbytek integrandu \(f(x)\) je polynomem aproximovatelným přesně.

Octave: kvadraturní integrátor quad(f,a,b,tol), kde f je funkce skalárního argumentu, též quadgk, quadcc, trapz ad.

f=@(x) sin(x); a=0; b=pi; tol=1e-10; quad(f,a,b,tol)                           % vrátí 2
format long; g=@(x) polyval([16 -16],x)./polyval([1 -2 0 4 -4],x); quad(g,0,1) % vrátí pi

Integrování metodou Monte Carlo

(NR kap. 7.6) Na vstupu metod Monte Carlo se objevují (pseudo)náhodná čísla. Monte Carlo výpočet vícerozměrného určitého integrálu vyčíslováním funkce \(f(x)\) v náhodných bodech \(x_i,\ i=1,\dots,N,\) je popsán vzorcem (NR 7.6.1):

\[\int f\,dV\approx V\langle f\rangle \pm V\sqrt\frac{\langle f^2\rangle-\langle f\rangle^2}N\ ,\]

kde \(V\) je integrační oblast a úhlové závorky značí aritmetický průměr, \(\langle f\rangle=\frac1N\sum_{i=1}^N f_i,\ f_i=f(x_i)\). Chyba je popsána ve statistickém smyslu a je jí vystižena pomalá konvergence metody: počet platných míst výsledku roste s \(\sqrt{N}\). (Chyba složeného lichoběžníkového pravidla \(L_N\) klesá s \(\frac1{N^2}\).)

Pro jednorozměrný integrál nabude Monte Carlo vzorec tvaru (podobného vzorci pro obdélníkové pravidlo \(O_N\)):

\[\int_a^b f(x)\,dx\approx (b-a)\langle f\rangle=\frac{b-a}N\sum_{i=1}^N f_i\,.\]

Oblíbenou variantou je obklopit \(m\)-rozměrný integrand \((m+1)\)-rozměrnou oblastí, tu pokrýt náhodnými body a nasčítat ty z nich, které leží „pod integrandem“. V našem případě tak na intervalu \(\langle a,b\rangle\) odhadneme \(y_{\min}\le f(x)\) a \(y_{\max}\ge f(x)\), zavedeme funkci \(F(x,y)=1\) pro \(y\le f(x),\) \(F(x,y)=0\) pro \(y>f(x),\) a metodou Monte Carlo vyčíslíme dvourozměrný integrál na obdélníku,

\[\int_a^b\int_{y_{\min}}^{y_{\max}} F(x,y)\,dydx\approx S\langle F\rangle=SK/N,\]

kde \(S=(b-a)(y_{\max}-y_{\min})\) a \(K=\sum_{i=1, y_i\le f_i}^N 1\,.\) Je totiž

\[\int_a^b\int_{y_{\min}}^{y_{\max}} F(x,y)\,dydx=\int_a^b\int_{y_{\min}}^{f(x)}1\,dydx=\int_a^bf(x)\,dx-(b-a)y_{\min}\,.\]
format long; f=@(x) sin(x); a=0; b=pi; ymin=0; ymax=1;
N=1000000; x=a+(b-a)*rand(1,N); I1D=(b-a)*mean(f(x))         % vrátí (např.) 1.999
y=ymin+(ymax-ymin)*rand(1,N); K=sum(y<=f(x)); I2D=(b-a)*ymin+(b-a)*(ymax-ymin)*K/N
N=1000; plot(x(1:N),f(x(1:N)),'o',x(1:N),y(1:N),'.'); axis([a b ymin ymax])