Cvičení č. 6: Numerické integrování

Počínaje dneškem se budeme věnovat numerickým metodám pro řešení standardních matematických problémů, obvyklých stavebních kamenů pro větší výpočty. Začneme vyčíslením určitých integrálů reálných funkcí jedné a dvou reálných proměnných,

\[\int_a^b f(x)\,dx\ , \quad \int_{x_{min}}^{x_{max}}\int_{y_{min}}^{y_{max}} f(x,y)\,dy\,dx\,.\]

Podíváme se do nabídky pythonovských modulů SymPy a SciPy a pak si to zkusíme vlastníma rukama.

SymPy a SciPy

Tak samozřejmě, pokud existuje rozumně složité analytické řešení problému, je to výhra. Pro analytická řešení musíme dávat pozor v matematice anebo se podívat do Computer Algebra Systems, jako je pythonovský modul SymPy. V něm se můžeme doptávat na určité i neurčité integrály, což si zkusíme s polynomickým, goniometrickým a racionálním integrandem:

\[\int_0^1 x^2\,dx=\frac13\,,\quad \int_0^\pi \sin x\,dx=2\,,\quad \int_0^1\frac{16x-16}{x^4-2x^3+4x-4}\,dx=\pi\,.\]

Definujeme x jako symbolickou proměnnou, její pomocí tři symbolické integrandy a pak už se obrátíme na funkci integrate:

# Python: integrování v SymPy
import sympy
def testSympy():
  x=sympy.Symbol('x')
  fpol=x*x
  fgon=sympy.sin(x)
  frat=(16*x-16)/(x**4-2*x**3+4*x-4)
  print(sympy.integrate(fpol,x),'...',sympy.integrate(fpol,(x,0,1)))
  print(sympy.integrate(fgon,x),'...',sympy.integrate(fgon,(x,0,sympy.pi)))
  print(sympy.integrate(frat,x),'...',sympy.integrate(frat,(x,0,1)))
  return None
testSympy()

Pro Sympy zřejmě snadný úkol:

x**3/3 ... 1/3
-cos(x) ... 2
2*log(x**2 - 2) - 2*log(x**2 - 2*x + 2) + 4*atan(x - 1) ... pi

Kde nepomůže analytické SymPy, nadějí je numerické SciPy. Výpočtu plochy pod křivkou se říká kvadratura a to přešlo do numerické matematiky jako synonymum numerického integrování. Zaznamenáme to i ve SciPy, které ve své části integrate pro výpočet určitých integrálů nazývá hlavní integrátor quad:

# Python: integrování v SciPy
import numpy
import scipy
def fpol(x):
  return x*x
def fgon(x):
  return numpy.sin(x)
def frat(x):
  # return (16*x-16)/(x**4-2*x**3+4*x-4)
  return (16*x-16)/(((x-2)*x*x+4)*x-4)
def quad_scipy(f,a,b):
  return scipy.integrate.quad(f,a,b)
print(quad_scipy(fpol,0,1))
print(quad_scipy(fgon,0,numpy.pi))
print(quad_scipy(frat,0,1))

# (0.33333333333333337, 3.700743415417189e-15)
# (2.0, 2.220446049250313e-14)
# (3.141592653589793, 3.449232861469955e-12)

Ve funkci frat jsme si dali záležet a k vyčíslení polynomu ve jmenovateli jsme použili vytýkání před závorku neboli metodu zvanou Hornerovo schéma, jejíž časová složitost je lineární vzhledem k stupni polynomu a tedy úspornější, než kdybychom vyčíslovali x*x, x*x*x a x*x*x*x navzájem nezávisle nebo se uchýlili k volání funkce power.

Lichoběžníkové pravidlo

Přichází náš čas. Chceme numericky vyčíslit \(\int_a^b f(x)\,dx\). Interval \(\langle a,b\rangle\) rozdělíme na \(N\) podintervalů téže délky, \(h=(b-a)/N\), jejich hraniční body označíme \(x_0=a,\ x_1,\ \dots,\ x_N=b\) a hodnoty integrandu v nich \(f_n=f(x_n),\ n=0,\dots,N\). Mezi body \((x_n,f_n)\) vedeme po částech lineární polynom (lomenou čáru). Hledaný integrál pak lze aproximovat plochou pod touto lomenou čarou, v případě podintervalu \(\langle x_n,x_{n+1}\rangle\) plochou lichoběžníka \(\frac h2(f_n+f_{n+1})\), na celém intervalu pak součtem všech těchto dílčích ploch (tzv. složené lichoběžníkové pravidlo neboli trapezoidal rule):

\[\int_a^b f(x)\,dx \approx L_N \equiv h\left(\frac{f_0}2+f_1+f_2+\dots+f_{N-1}+\frac{f_N}2\right).\]

Počet lichoběžníků \(N\) z nebe nespadne, vhodnější by bylo poručit si pro výsledek požadovanou přesnost \(\varepsilon\). Pak bychom začali s \(N=1\) a následně \(N\) zdvojnásobovali tak dlouho, než by nastalo \(|L_N-L_{N/2}|<\varepsilon\). Jistě je to jen aproximace požadované přesnosti, ale běžná. Současně bychom limitovali růst \(N\) vhodnou horní mezí tak, aby byl program schopen skončit i tehdy, nebude-li požadované přesnosti dosaženo.

Nejprve tedy verze s předepsaným \(N\) neboli nmax:

# Python: integrování lichoběžníkovým pravidlem
import numpy
import time
from numba import njit

@njit
def fpol(x):
  return x*x
@njit
def fgon(x):
  return numpy.sin(x)
@njit
def frat(x):
  return (16*x-16)/(((x-2)*x*x+4)*x-4)

# Lichoběžníkové pravidlo pomocí cyklu
@njit
def quad_for(f,a,b,nmax):
  result=(f(a)+f(b))*0.5
  h=(b-a)/nmax
  x=a
  for n in range(1,nmax):
    x+=h
    result+=f(x)
  return result*h

nmax=1_000_000
tic=time.time()
print(quad_for(fpol,0,1,nmax))
print(quad_for(fgon,0,numpy.pi,nmax))
print(quad_for(frat,0,1,nmax))
toc=time.time(); print('Wallclock time:',f'{toc-tic:0.2f} s')

# 0.33333333333399606
# 2.000000000006493
# 3.1415926535793384
# Wallclock time: 0.53 s

Výsledek potvrzuje to, co o chybě lichoběžníkového pravidla říká Eulerův-Maclaurinův vzorec, totiž že je úměrná \(h^2\). Pro nmax=1_000_000 je chyba přibližně 1e-12 a – wow! – my jsme skutečně dostali 12 platných míst. Ale zdání klame, v cíli ještě nejsme. Zkusme nmax=100_000_000:

# 0.33333333417951294
# 1.9999999964431938
# 3.1415926485222867
# Wallclock time: 1.20 s

My čekali 14 platných míst, a zatím přišel mráz. To, co vidíme, je projevem šíření zaokrouhlovací chyby v dlouhé sumaci. Vinu mají oba dva řádky v cyklu for. První z nich lze napravit snadno, hned třemi prakticky ekvivalentními způsoby:

# Python: robustní varianty průchodu aritmetickou posloupností
for n in range(1,nmax):
  x=a+n*h
  result+=f(x)

# arange(start,stop,step) pro hodnoty start,start+step,... z intervalu <start,stop)
for x in numpy.arange(a+h,b,h):
  result+=f(x)

# linspace(start,stop,num) pro num ekvidistatních hodnot z <start,stop>
for x in numpy.linspace(a+h,b-h,nmax-1):
  result+=f(x)

# 0.3333333333332225
# 2.000000000000248
# 3.141592653590256
# Wallclock time: 1.91 s

Počty jsou takto mírně pomalejší, ale eliminace sumování při vytváření aktuálního x zjevně pomohla přesnosti. Se sumováním f(x) je to těžší. Nápravu přinese nahrazení sekvenční sumace párovanou sumací, kterou si zde sice nenaprogramujeme, ale funkce math.fsum ji nabízí hotovou. Funkce sčítá celé pole, můžeme tedy zrušit cyklus a raději rovnou nahradíme funkci quad_for novým exemplářem:

# Python: lichoběžníkové pravidlo pomocí funkcí sum a math.fsum
import math
#@njit
def quad_fsum(f,a,b,nmax):
  x=numpy.linspace(a,b,nmax+1)
# result=(f(a)+f(b))*0.5+sum(f(x[1:nmax]))
  result=(f(a)+f(b))*0.5+math.fsum(f(x[1:nmax]))
  return result*(b-a)/nmax

# 0.33333333333333337
# 2.0
# 3.1415926535897927
# Wallclock time: 16.51 s

Pokud bychom užili řádek se standardní funkcí sum, shledali bychom výsledky shodné s předchozím případem. Funkce sum tedy sčítá sekvenčně; jejím pozitivem je, že ji snese Numba a doba výpočtu je tak u 2 sekund. Funkci math.fsum ovšem Numba dosud nerozumí, direktivu @njit jsme proto zakomentovali. Pro nmax=100_000_000 jsme pak v 8bytovém reálném typu získali plnou přesnost, doba běhu se však protáhla skoro o řád. To nás dostává k poslední variantě našich pokusů o lichoběžníkové pravidlo – předepíšeme přesnost eps a budeme zvětšovat počet bodů sítě zvětšovat jen tak dlouho, než dosáhneme eps:

# Python: lichoběžníkové pravidlo s předepsanou přesností
#@njit
def quad_eps(f,a,b,nmax):
  eps=1e-14
  n=1
  Lnew=quad_fsum(f,a,b,n)
  while True:
    Lold=Lnew
    n*=2
    if n>nmax: break
    Lnew=quad_fsum(f,a,b,n)
    if abs(Lnew-Lold)<eps*abs(Lnew): break
  return Lnew

# 0.3333333333333339
# 1.9999999999999942
# 3.1415926535897887
# Wallclock time: 4.12 s

Fortranské řešení, včetně vlastnoruční párované sumace, si můžeme prohlédnout na stránce předmětu Programování prakticky.

Pro vizuální ověření našich počtů vytvoříme v gnuplotu obrázek o třech panelech, v nichž vykreslíme integrandy zadané analyticky a odhadneme plochu pod nimi:

_images/integral3.png
# gnuplot: multiplot and plot with boxes
fpol(x)=x**2; fgon(x)=sin(x); fpi(x)=(16*x-16)/(((x-2)*x*x+4)*x-4)
set terminal pngcairo enhanced crop font 'Sans,10' size 1024,800
set output 'integral3.png'
set samples 1000
set style fill solid
set multiplot
set size .35,.35
set origin 0,0
set xtics 0,.5; set ytics 0,1
set key left top
set title 'plocha zeleně: 1/3'
plot [0:1] fpol(x) with boxes linecolor 'green' title 'x^2'
set origin .33,0
set xtics 0,1; set ytics 0,1
set key right top
set title 'plocha zeleně: 2'
plot [0:pi] fgon(x) with boxes linecolor 'green' title 'sin x'
set origin .66,0
set xtics 0,.5; set ytics 0,4
set key right top
set title 'plocha zeleně: {/Symbol p}'
plot [0:1] fpi(x) with boxes linecolor 'green' title 'fpi(x)'
unset multiplot
set output
unset terminal

Více dimenzí

Kvadraturní vzorce lze zobecnit pro integrování na dvourozměrném (2D) obdélníku i třírozměrném (3D) kvádru – při integraci přes \(y\) se pro každé \(y\) provede 1D integrace přes \(x\), při integraci přes \(z\) se pro každé \(z\) provede 2D integrace přes \(x\) a \(y\). Zvolme si jako ukázku integrál z chobotnice,

\[\int_0^\pi \int_0^\pi \sin x \sin y\,dx\,dy=4\,,\]

na kterou se můžeme podívat opět v gnuplotu:

_images/chobotnice.png
# gnuplot: chobotnice
chobotnice(x,y)=sin(x)*sin(y)
set terminal pngcairo crop size 800,600
set output 'chobotnice.png'
set xtics 0,.5; set ytics 0,.5; set ztics 0.25,.25
set xyplane at 0
set grid
unset colorbox
unset key
set multiplot
set samples 200
set isosamples 200
splot [0:pi][0:pi][0:1] chobotnice(x,y) with pm3d
set isosamples 20
splot [0:pi][0:pi][0:1] chobotnice(x,y) with lines
unset multiplot
set output

Sepíšeme 2D verzi lichoběžníkového pravidla na obdélníku \(\langle x_\min{},x_\max{}\rangle\times\langle y_\min{},y_\max{}\rangle\):

# Python: 2D integrování lichoběžníkovým pravidlem
import numpy
import time
from numba import njit

@njit
def fchobotnice(x,y):
  return numpy.sin(x)*numpy.sin(y)

@njit
def dblquad(f,xmin,xmax,ymin,ymax,nmaxx,nmaxy):
  result=0
  hx=(xmax-xmin)/nmaxx
  hy=(ymax-ymin)/nmaxy
  t=0.5*(f(xmin,ymin)+f(xmax,ymin))
  for x in numpy.linspace(xmin+hx,xmax-hx,nmaxx-1): t+=f(x,ymin)
  result+=0.5*t
  for y in numpy.linspace(ymin+hy,ymax-hy,nmaxy-1):
    t=0.5*(f(xmin,y)+f(xmax,y))
    for x in numpy.linspace(xmin+hx,xmax-hx,nmaxx-1): t+=f(x,y)
    result+=t
  t=0.5*(f(xmin,ymax)+f(xmax,ymax))
  for x in numpy.linspace(xmin+hx,xmax-hx,nmaxx-1): t+=f(x,ymax)
  result+=0.5*t
  return result*hx*hy

nmaxx=1_000
nmaxy=1_000
tic=time.time()
print(dblquad(fchobotnice,0,numpy.pi,0,numpy.pi,nmaxx,nmaxy))
toc=time.time(); print('Wallclock time:',f'{toc-tic:0.2f} s')

# 3.9999934202653558
# Wallclock time: 0.61 s

DÚ: Kruh a koule

Za domácí úkol prosím vytěžte předložené funkce ke spočtení obsahu jednotkového kruhu a 3/4 objemu jednotkové koule v kartézských souřadnicích. Patrně bude vhodné připravit si funkci pulkruh pro horní polovinu jednotkové kružnice na intervalu \(\langle-1,1\rangle\) a funkci polokoule pro horní polovinu povrchu jednotkové koule na čtverci \(\langle-1,1\rangle\times\langle-1,1\rangle\), voláním některé varianty quad*, resp. dblquad spočítat plochu, resp. objem pod nimi a vynásobit dvěma.