Cvičení č. 3: Procedurální styl

Sečteme ještě jednu řadu pro \(\pi\). Výslednému programu už dáme jistou strukturu: každou myšlenku zabalíme do samostatné programové jednotky – funkce, a ty sdružíme do modulu. Nazrál čas i na počty, které už lze vizualizovat: Lissajousovy obrazce. Pro kreslení budeme teď i potom používat gnuplot a příležitostně, jako bonbonek, ParaView.

Ramanujanův součet

Pokračujeme ve výpočtech \(\pi\) podle vzorců ze souhrnu v PDF, a to fascinující Ramanujanovou řadou (fascinující je i Ramanujanovo CV):

\[\frac1\pi=\frac{\sqrt8}{99^2}\sum_{n=0}^\infty\frac{(4n)!}{(n!)^4}\frac{1103+26390n}{396^{4n}} =\frac{\sqrt8}{9801}\left(1103+\frac{54986}{2049271488}+\cdots\right)\]

Podobně jako minule u Viètova součinu se hodí vyjádřit n-tý člen řady pomocí (n-1)-ho členu:

\[a_n=a_{n-1}\frac{(4n-3)(4n-2)(4n-1)(4n)}{(396n)^4},\quad a_0=1\]

Budeme tak čelit přetékání v čitateli i jmenovateli (resp. aritmetice velkých integerů v Pythonu) při doslovném počítání podle zadání; nekonfliktní faktor \(1103+26390n\) jsme vyčlenili stranou. I tak se ovšem jistíme přechodem k vyčíslování v reálné aritmetice použitím reálných literálů ve výrazech, např. (4.*n-3) v Pythonu a (4._8*n-3) ve Fortranu.

Funkce

Výpočty \(\pi\) podle jednotlivých předpisů nyní strukturujeme jako programové jednotky nazývané funkce, parametrizované svými argumenty a vracející hodnotu použitelnou jako výraz ve výpisu. Výsledný program bude v Pythonu tvořen dvojicí souborů, a to modulem mPi.py s definicemi funkcí:

# Modul mPi.py: výpočty pí podle http://geo.troja.mff.cuni.cz/~lh/NOFY056/cviceni-2022/_static/piFormulas.pdf
# Procedurální styl

from math import pi,atan,sqrt

# Machin
def piMachin():
  return (4*atan(1/5)-atan(1/239))*4

# Leibniz
def piLeibniz(nmax):
  piLeibniz=0.
  z=1.
  for n in range(1,nmax+1):
    piLeibniz+=z/(2*n-1)
    z=-z
  return piLeibniz*4

# Euler: kumulace od nejmenších členů
def piEuler(nmax):
  piEuler=0.
  for n in range(nmax,0,-1):
    piEuler+=1/(n*n)
  return sqrt(piEuler*6)

# Viete
def piViete(nmax):
  piViete=1.
  a=0.
  for n in range(1,nmax+1):
    a=sqrt(0.5+0.5*a)
    piViete*=a
  return 2/piViete

# Ramanujan
def piRamanujan(nmax):
  a=1.
  piRamanujan=a*1103
  for n in range(1,nmax+1):
    a=a*(4*n-3)*(4*n-2)*(4*n-1)*4*n/(396.*n)**4
    piRamanujan+=a*(1103+26390*n)
  return 1/(piRamanujan*sqrt(8)/9801)

a skriptem, který modul připojí a vytěží:

# Python: volání funkcí z modulu mPi nebo mPiFast
from mPi import pi,piMachin,piLeibniz,piEuler,piViete,piRamanujan
# from mPiFast import pi,piMachin,piLeibniz,piEuler,piViete,piRamanujan
nmax=10_000_000        # vhodné pro mPi
# nmax=1_000_000_000   # vhodné pro mPiFast
format='%25s%20.15f'
print(format%('pi podle math =',pi))                    # printf-style string formatting
# print('{:>25}{:20.15f}'.format('pi podle math =',pi)) # string format() method
# print(f'{"pi podle math =":>25}{pi:20.15f}')          # formatted string literals
print(format%('pi podle Machina =',piMachin()))
print(format%('pi podle Leibnize =',piLeibniz(nmax)))
print(format%('pi podle Eulera =',piEuler(nmax)))
print(format%('pi podle Vieta =',piViete(25)))
print(format%('pi podle Ramanujana =',piRamanujan(2)))

Kromě funkcí jsme v ukázce nově použili rozšířená přiřazení += a *= a také podtržítka v numerických literálech. A když už máme možnost zapsat čitelně celočíselnou miliardu, zkusme sáhnout po balíčku Numba (instalace: pip install numba), s jehož pomocí lze naše funkce přeložit, zvýšit tak jejich rychlost téměř k fortranské rychlosti, a tedy získat možnost počítat s odezvou v sekundách pro nmax=1_000_000_000. Místo modulu mPi.py připojíme následující modul mPiFast.py:

# Modul mPiFast.py: výpočty pí podle http://geo.troja.mff.cuni.cz/~lh/NOFY056/cviceni-2022/_static/piFormulas.pdf
# Procedurální styl a překlad pomocí balíčku Numba

from math import pi,atan,sqrt
from numba import njit

# Machin ...

# Leibniz
@njit
def piLeibniz(nmax):
  piLeibniz=0.
  z=1.
  for n in range(1,nmax+1):
    piLeibniz+=z/(2*n-1)
    z=-z
  return piLeibniz*4

# Euler: kumulace od nejmenších členů
@njit
def piEuler(nmax):
  piEuler=0.
  for n in range(nmax,0,-1):
    piEuler+=1/(n*n)
  return sqrt(piEuler*6)

# Viete ...

# Ramanujan ...

Vidíme, že v našem případě stačilo importovat z Numby direktivu @njit pro just-in-time compilation v nopython mode a tuto direktivu umístit před definici každé funkce, kterou chceme takto přeložit. Smysl to má zejména pro piLeibniz a piEuler, jejichž pomalá konvergence volá po velkém nmax.

Fortranskou verzi s modulem a funkcemi v něm a programem modul volajícím máme v jediném souboru:

! Fortran: výpočty pí podle http://geo.troja.mff.cuni.cz/~lh/NOFY056/cviceni-2022/_static/piFormulas.pdf
! Procedurální styl, funkce v modulu

module mPi
implicit none
real(8),parameter :: pi=atan(1._8)*4  ! konstanta (atribut parameter)

contains

! Machin
real(8) function piMachin()
piMachin=(4*atan(1._8/5)-atan(1._8/239))*4
end function

! Leibniz
real(8) function piLeibniz(nmax)
integer nmax,n
real(8) z
piLeibniz=0.
z=1.
do n=1,nmax
  piLeibniz=piLeibniz+z/(2*n-1)
  z=-z
end do
piLeibniz=piLeibniz*4
end function

! Euler: kumulace od nejmenších členů
real(8) function piEuler(nmax)
integer nmax,n
real(8) rn
piEuler=0.
do n=nmax,1,-1
  rn=n; piEuler=piEuler+1/(rn*rn)
end do
piEuler=sqrt(piEuler*6)
end function

! Viete
function piViete(nmax) result (result)
real(8) result              ! lokální jméno návratové hodnoty (klauzule result)
integer,intent(in) :: nmax  ! read-only argument (atribut intent(in))
integer n
real(8) a
result=1.
a=0.
do n=1,nmax
  a=sqrt(0.5+0.5*a)
  result=result*a
end do
result=2/result
end function

! Ramanujan
function piRamanujan(nmax) result (result)
real(8) result
integer,intent(in) :: nmax
integer n
real(8) a
a=1.
result=a*1103
do n=1,nmax
  a=a*(4*n-3)*(4*n-2)*(4*n-1)*4*n/(396._8*n)**4
  result=result+a*(1103+26390*n)
enddo
result=1/(result*sqrt(8._8)/9801)
end function

end module

program p3
use mPi       ! připojení modulu
implicit none
integer :: nmax=1000000000
character(100) :: format='(a25,f20.15)'

print format,'pulnocni pi =',pi
print format,'pi podle Machina =',piMachin()
print format,'pi podle Leibnize =',piLeibniz(nmax)
print format,'pi podle Eulera =',piEuler(nmax)
print format,'pi podle Vieta =',piViete(25)
print format,'pi podle Ramanujana =',piRamanujan(2)

end program

DÚ: Lissajousovy obrazce

Počítání jako výše, tedy výpočet jediného čísla – navíc předem přibližně známého – je ve fyzice celkem vzácnost. Typičtější je chrlení množství dat, které je pak třeba dále analyzovat, např. prostou vizualizací. To nyní uděláme.

Budeme generovat vývoj dvou souřadnic popsaný goniometrickými funkcemi a kreslit dráhu, po které se bude pohybovat bod \((x,y)\) v čase \(t\):

\[\begin{split}x(t) & = a_x \sin(f_xt+\varphi_x) \\ y(t) & = a_y \sin(f_yt+\varphi_y)\end{split}\]

Taková dráha, zvláště je-li podíl frekvencí v x a y směru racionálním číslem, se nazývá Lissajousův obrazec.

Náš pythonský skript lissajous.py i fortranský program lissajous.f90 si nesou vstupní data ve svém těle (hned nahoře, dobře na očích) a výstupní data vypisují volitelně buď na obrazovku nebo do souboru lissajous.dat:

# Python: příprava dat pro Lissajousovy obrazce
from math import pi,sin
nmax=1000   # počet kroků
ax=1        # x amplituda
ay=1        # y amplituda
fx=3        # x frekvence
fy=2        # y frekvence
phx=0       # x fáze
phy=0       # y fáze
tmin=0      # počáteční čas
tmax=pi*2   # koncový čas
output=1    # 1: na obrazovku, 2: do souboru
filename='lissajous.dat'
if output==2:
  f=open(filename,'w')
for n in range(0,nmax+1):
  t=(tmax-tmin)*n/nmax
  x=ax*sin(fx*t+phx)
  y=ay*sin(fy*t+phy)
  match output:
    case 1: print(f'{x:8.4f}{y:8.4f}')
    case 2: print(f'{x:8.4f}{y:8.4f}',file=f)
if output==2:
  f.close()
  print('Vytvořen soubor '+filename+'.')
! Fortran: příprava dat pro Lissajousovy obrazce
program Lissajous
implicit none
real(8),parameter :: pi=atan(1._8)*4
integer nmax,n,output,f
real(8) ax,ay,fx,fy,phx,phy,tmin,tmax,t,x,y
character(20) filename
nmax=1000   ! počet kroků
ax=1        ! x amplituda
ay=1        ! y amplituda
fx=3        ! x frekvence
fy=2        ! y frekvence
phx=0       ! x fáze
phy=0       ! y fáze
tmin=0      ! počáteční čas
tmax=pi*2   ! koncový čas
output=2    ! 1: na obrazovku, 2: do souboru
filename='lissajous.dat'
if (output==2) open (newunit=f,file=filename)
do n=0,nmax
  t=(tmax-tmin)*n/nmax
  x=ax*sin(fx*t+phx)
  y=ay*sin(fy*t+phy)
  select case (output)
  case (1); print '(2f8.4)',x,y
  case (2); write (f,'(2f8.4)') x,y
  end select
enddo
if (output==2) then
  close (f)
  print '(a)','Vytvoren soubor '//trim(filename)//'.'
endif
end program

Jakmile máme data v souboru, můžeme je vykreslit např. v gnuplotu.

Získat vizualizační program (gnuplot) není trnité. Po instalaci (na mysli máme dále Windows, v Linuxu to není moc jiné) se mezi programy nabídky Start objeví gnuplot s položkou gnuplot 5.4.4 (nebo jiná verze), která spustí wgnuplot.exe s mini-GUI nebo console version pro gnuplot.exe. Kde se gnuplot spustil, na to odpoví příkaz pwd (print working directory). Pokud jsou Vaše data (soubor lissajous.dat) v adresáři např. D:\User\Data, je třeba zadat příkaz cd 'D:\User\Data' (change directory). Je samozřejmě také možné nastavit si pracovní adresář v Properties položky gnuplot 5.4.4 v nabídce Start nastavením adresáře Start in, nebo přidat adresář s wgnuplot.exe do proměnné PATH způsobem, který jsme popsali na prvním cvičení. Pak už by nic nemělo stát v cestě příkazu plot 'lissajous.dat' with lines neboli zkráceně p 'lissajous.dat' w l, kterým se otevře vizualizační okno s výsledkem. Z toho se do příkazového okna lze vrátit mj. stisknutím Ctrl+Space, a můžete pokračovat dalším příkazem nebo skončit zadáním quit neboli q.

_images/gnuplot.png

Sekvenci příkazů vedoucích až k obrázku v PNG formátu shrnuje následující skript:

# gnuplot: kreslení dat ze souboru na obrazovku a v PNG formátu
plot 'lissajous.dat' with lines
set terminal pngcairo      # pokrocilejsi terminal nez starsi png
set output 'lissajous.png'
replot
set output

# zkracovaný zápis
# p 'lissajous.dat' w l
# set term pngc
# set o 'lissajous.png'
# rep
# set o

který se buď spustí v gnuplotu příkazem load 'lissajous.gn' nebo se uvedené příkazy zadají v gnuplotu po jednom ručně.

_images/lissajous-gnuplot.png

Lissajousův obrazec pro poměr frekvencí v x a y směru 3:2. Kresleno v gnuplotu.

A jistě, kdybychom si zrovna neprocvičovali programování a zápis dat do souboru, mohli bychom Lissajousovy obrazce kreslit v gnuplotu rovnou. Gnuplot rozumí funkcím zadaným analyticky a umí i s proměnnými, takže např. omega=1; plot sin(omega*x) vykreslí goniometrickou funkci v default intervalu <-10,10>. (x je default symbol gnuplotu pro nezávisle proměnnou.) Gnuplot rozumí i parametrickému zápisu funkcí, takže výše předložený Lissajousův obrazec lze získat i na zelené louce tímto skriptem (při parametrickém zápisu je default symbolem pro parametr t):

# gnuplot: kreslení funkce zadané analyticky, zde navíc parametricky
set parametric
set samples 1000
fx=3
fy=2
plot sin(fx*t),sin(fy*t)

Za domácí úkol upravte zdrojový text ze cvičení tak, aby počítal trojrozměrné Lissajousovy obrazce (Wikipedia: Lissajousovy uzly). Pro jejich kreslení můžete opět využít gnuplot: splot 'lissajous.dat'. V interaktivním režimu gnuplotu si 3D obrázkem můžete pomocí myši otáčet. Alternativně se můžete podívat na výsledek ve směru jednotlivých os pomocí průmětů, tedy výběrem různých párů z vypočtených tří sloupců: plot 'lissajous.dat' using 1:2 with lines; p 'lissajous.dat' u 1:3 w l; p 'lissajous.dat' u 2:3 w l. Vyberte nějaký úhledný pohled a ve formátu PNG mi jej také pošlete.

Můžete si také prohlédnout dvoustránkový výtah z dokumentace gnuplotu.

P. S. ParaView

Gnuplot je silný v dvourozměrných (2D) obrázcích, slabší (co do kvality výstupu) ve třech dimenzích (3D). Pro 3D vizualizaci je samozřejmě na výběr více cest. Jednou z nejhodnotnějších (a přitom volně dostupných, a pro Windows/Linux/Mac) je ParaView. Tak proč si 3D Lissajousův obrazec neprohlédnout v ParaView?

Vstupní data má ParaView nejraději ve VTK formátu. Pro nakreslení série 3D bodů nám půjde o variantu zvanou polygonální síť alias polydata. Vlastně bude stačit vložit na počátek datového souboru pro gnuplot, který máte získat jako domácí úkol, hlavičku:

# vtk DataFile Version 3.0
comment (max. 1023 chars)
ASCII
DATASET POLYDATA
POINTS 1001 float

a souboru dát příponu vtk. (Budete-li generovat jiný počet bodů než 1001, upravte údaj v hlavičce.) V ParaView pak můžete postupovat třeba takto:

  • File–Open– lissajous.vtk, v panelu Properties vlevo uprostřed tlačítko Apply

  • Filters–Alphabetical–Glyph (nebo nahoře ve třetím řádku ikonek sedmá ikonka zleva), Apply

  • a dál v panelu PropertiesGlyph Source–Glyph Type: Sphere, Apply

  • (pro lepší rozlišení koulí) Glyph Source–Theta/Phi Resolution, Apply (pokud nutno, přepnout předtím Toggle advanced properties kolečkem pod otazníkem vpravo nahoře v Properties)

  • (pro zvětšení koulí) Scale–Scale Factor, Apply

  • (pro pravidelný rozestup koulí) Masking–Glyph Mode–All Points, Apply

  • (pro obarvení koulí) Coloring–Edit a (pro barvu pozadí) Background v Properties dole

a získáte třeba toto:

_images/lissajous-paraview.png

3D Lissajousův obrazec pro poměr frekvencí v x, y a z směru 3:5:7. Kresleno v ParaView.

S obrázkem si v okně ParaView můžete točit, posouvat atd., můžete ho uložit jako bitmapu (File–Save Screenshot) a můžete si také třeba nastavit View–Animation View, pak u modrého křížku dole Camera–Orbit–klik na modrý křížek–OK, ještě Mode–Sequence–No. Frames: 250, nahoře zelenou šipkou Play a dál už si založit ruce a sledovat. Nebo ještě File–Save Animation, Files of type: avi, Frame Rate: 25 a získat toto.

K ParaView a kreslení ve třech rozměrech se ještě vrátíme.