Pavel Stránský
Contact: Contact
Česky English
Last updated: 15.7.2019

Použití PC ve fyzice

Akademický rok 2018-2019

Cíl cvičení

Účelem cvičení je pohrát si s jednoduchými fyzikálními problémy a vyzkoušet si různé způsoby jejich řešení (i triviální řešení je řešení) a zobrazení výsledků. Většina příkladů bude demonstrována v populárním programovacím jazyce Python, ale dojde i na základy psaní textů v LaTeXu, software Mathematica nebo na kancelářský Excel. Další informace k předmětu jsou na stránce SISu.

Ve cvičení nebude probíhat systematická výuka žádného programovacího jazyka ani numerických algoritmů. K tomu jsou určeny specializované přednášky na oboru fyzika, například

NEVF523 - Numerické metody počítačové fyziky
NPRF050 - Programování v Pythonu
NEVF107 - C++ pro fyziky
NPRF020 - Úvod do programování v prostředí MATLAB, Octave a Scilab
NPRF006 - Pokročilé metody programování
NTMF048 - Použití systémů počítačové algebry ve fyzice
NJSF081 - Software a zpracování dat ve fyzice částic

a nepřeberně přednášek na oboru informatika.

Zápočet

Tři možné cesty, jak získat zápočet (seřazeno od té nejupřednostňovanější):

1. Naprogramujete řešení nějakého jednoduchého fyzikálního problému (o složitosti srovnatelné s příklady procvičenými během hodin), který si zvolíte podle svých vlastních zájmů a zálib. Hotový program mi představíte a vysvětlíte buď na posledním cvičení semestru, nebo kdykoliv po domluvě. Program by neměl být totožný s těmi, za které jste (byli, budete) klasifikováni v jiných předmětech. Doporučuji vám, abyste se mnou zadání konzultovali, než se pustíte do práce.
2. Pokud vás nic, ale opravdu nic nenapadne, zadám na posledním cvičení semestru jednoduchý problém vycházející z procvičených příkladů, který během dvouhodinovky vyřešíte a předvedete.
3. Pokud se nebudete moci dostavit v zápočtovém týdnu na cvičení (nebo se dostavíte, ale zápočtovou úlohu se vám nepodaří vyřešit), kontaktujete mě a já vám zadám úlohu na doma. Až ji naprogramujete, řešení mi předvedete (nejpozději do poloviny září 2019).

Ačkoliv převážná část předmětu bude demonstrována příklady v programovacím jazyce Python 3.7, rozhodně není nutné, abyste zápočtový program psali v tomto jazyce. Zvolte si sami, jaký jazyk či vývojové prostředí jsou vám sympatické.

Literatura

Na webu existuje nepřeberné množství tutoriálů a návodů k Pythonu a jeho knihovnám, které vám pomůžou s konkrétními problémy. Zde uvádím spíš monografie pro pokročilejší, v nichž naleznete seriózní příklady využití tohoto programovacího jazyka ve vědě.

[1]P.R. Turner, T. Arildsen, K. Kavanagh, Applied Scientific Computing With Python (Springer 2018)
[2]R. Johansson, Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib (Springer 2019)
[3]S. Nagar, Introduction to Python for Engineers and Scientists (Springer 2018)
[4]S. Lynch, Dynamical Systems with Applications using Python (Springer 2018)
[5]B.J. Korites, Python Graphics: A Reference for Creating 2D and 3D images (Springer 2018)

Poznámky a zdrojové kódy

(Pozor: řešení uvedená v přiložených kódech zdaleka nejsou jediná možná a mohou obsahovat chyby. Vždy dejte na vlastní hlavu. Kódy jsou většinou mírně modifikovány v porovnání s tím, co jsme naprogramovali na cvičeních.)

1. cvičení

Stíhací křivka

Křivka, po které se pohybuje pes, když se v každém okamžiku pohybuje směrem k myslivci. Předpokládáme, že velikost rychlosti psa i myslivce jsou konstantní.

Zatímco analytické řešení je celkem zdlouhavé a lze ho provést pouze v případě, že se myslivec pohybuje po přímce (řešení najdete například v monografii Viktor Trkal, Mechanika hmotných bodů a tuhého tělesa, ČSAV 1956, str. 279), numerické řešení je jednoduché (a to i v případě silně vrávorajícího myslivce).

Objekt Turtle v Pythonu

Objekt Turtle je kurzor, který můžete jednoduchými příkazy posouvat po obrazovce a vykreslovat jím křivky.

Referenční příručka k objektu Turtle: https://docs.python.org/3/library/turtle.html.

Zdrojový kód naprogramované stíhací křivky: StihaciKrivka.py.

Vyzkoušejte si: 1. Různé velikosti časového kroku. Čím je krok kratší, tím získáme přesnější křivku, ale výpočet trvá déle.
2. Myslivec chodí po složitějších křivkách (kruh, sinusovka).
3. Rychlost psa je menší než rychlost myslivce.

2. cvičení

Objekty a obsluha událostí v Pythonu

Stíhací křivka naprogramovaná objektově, přičemž 1) pohyb myslivce lze ovládat kurzorovými šipkami a 2) program lze rozšířit i pro více psů zároveň.

Zdrojový kód: StihaciKrivkaO.py.

Vyzkoušejte si: Více psů na přímce a na kružnici (jako příklad poslouží funkce Primka(nPes, vMyslivec, vPes, dt) a Kruznice(nPes, vMyslivec, vPes, dt), kde nPes je počet psů). Sledujte, jak se vyvíjí "obalová křivka". Pozor, pro velké množství psů a krátký časový krok dt může být výpočet velmi pomalý. Mějte trpělivost.

3. cvičení

Zrychlení psa

Směr rychlosti psa se mění, takže pes se pohybuje se zrychlením. Velikost rychlosti psa se nemění, takže zrychlení je "dostředivé", kolmé v každém okamžiku na rychlost. Zrychlení psa budeme počítat přímo ve funkci Pes.Krok.

Knihovna matplotlib a zobrazení grafů v Pythonu

Zrychlení zobrazíme v grafu pomocí knihovny matplotlib. Zrychlení se bude zobrazovat v reálném čase. Kladné (záporné) znaménko zrychlení značí, že pes zatáčí vlevo (vpravo).

Zdrojový kód: StihaciKrivkaZrychleni.py. Oproti předchozím kódům počítáme s úhlovou proměnnou v radiánech, čehož se docílí voláním metody radians() třídy Turtle. Je nutné upravit všechny části kódu, ve kterých jsme doposud uvažovali úhly ve stupních (například při zatáčení myslivce).

Pozor! Knihovna matplotlib není s instalací Pythonu nainstalována automaticky. K jejímu doinstalování spusťte na svých počítačích v adresáři s Pythonem verze, kterou používáte, příkaz python -m pip install matplotlib.

Omezení velikosti zrychlení psa

Pes má konečnou sílu nebo se pohybuje po povrchu, který klouže. Jak se změní stíhací křivka, pokud zrychlení psa nemůže překročit nějakou maximální velikost?

Zdrojový kód: StihaciKrivkaZrychleniLimit.py.

Vyzkoušejte si: Pokud je aMax velmi malé (představte si myslivce a psa na zamrzlém rybníce), pes myslivce jdoucího rovně nedoběhne, ani když se sám pohybuje výrazně rychleji (velká rychlost psa je pak spíš na překážku). Nalezněte takovou kombinaci parametrů.
Zobrazení okamžitého zrychlení

Zrychlení můžeme zobrazit také pomocí úseček přímo do okna se želvičkami. Zdrojový kód: StihaciKrivkaZrychleniLimitN.py.

4. cvičení

Funkce sinus jako řešení diferenciální rovnice

Jedna z možností, jak vykreslit funkci sinus, je řešením pohybové rovnice pro harmonický oscilátor s odpovídajícími počátečními podmínkami. Nejtriviálnější explicitní Eulerova metoda integrace diferenciální rovnice (funkce SinusEuler) dává přesný výsledek pouze při použití velmi krátkého kroku. Je to metoda prvního řádu, a tak k desetinásobnému zlepšení přesnosti (jedna platná cifra navíc) potřebujeme desetkrát menší krok. Vylepšení Eulerovy metody spočívají v lepším odhadu směru, ve kterém děláme krok (funkce SinusEuler2, což je metoda 2. řádu). Jednokrokové explicitní metody se obecně nazývají Runge-Kuttovy metody (RK) a nejčastěji se používají RK 4. řádu.

Další třídou metod jsou symplektické metody, které jsou navržené pro řešení pohybových rovnic konzervativních systémů. Jejich výhodou je to, že s velkou přesností zachovávají energii (zatímco při použití jiných metod energie systému roste). Nejznámější symplektická metoda je Verletův algoritmus, který vychází z Eulerovy metody 2. řádu (funkce SinusVerlet). Lze však jednoduchým zásahem docílit, aby i Eulerova metoda 1. řádu byla symplektickým algoritmem (funkce SinusEuler1a a SinusEuler1b).

Zdrojový kód: Sinus.py.

Využití knihovny scipy

Knihovna scipy Pythonu obsahuje dvě možnosti řešení diferenciálních rovnic: integrate.odeint (LSODA řešitel knihovny ODEPACK, který podle typu rovnice volí Adamsovu prediktor-korektor metodu nebo BDF metodu) a integrate.ode (více možných řešitelů včetně RK metod; řešitele lze explicitně vybrat). Ve zdrojovém kódu jsou uvedeny příklady obou dvou postupů (funkce SinusOdeint a SinusOde).

Vyzkoušejte si: 1. Porovnejte si přesnost všech metod. Doporučuji volit parametry dx=0.01, xMax=100. Pro metody vyšších řádů zvolte krok větší.
2. Přesvědčte se, že u algoritmů 1. řádu chyba opravdu závisí přímo úměrně na velikosti kroku, zatímco u algoritmů vyšších řádů klesá se zmenšujícím se krokem rychleji. Zkuste odhadnout řád metod knihovny scipy.

5. cvičení

Vazbová energie deuteronu

Vazbovou energii deuteronu získáme řešením bezčasové Schrödingerovy rovnice pro sféricky symetrickou vlnovou funkci. Schrödingerova rovnice se v tomto případě redukuje na obyčejnou lineární diferenciální rovnici 2. řádu, která má dvě lineárně nezávislá řešení chovající se asymptoticky jako rostoucí a klesající exponenciála. Nás zajímá konvergující řešení, cílem je tedy najít takovou vazebnou energii B (chápejte jsi v tuto chvíli jako proměnný parametr), která vynuluje divergující exponenciálu. Tu budeme hledat minimalizováním hodnoty integrálu kvadrátu vlnové funkce s proměnnou B. Nalezení správné hodnoty nám určí nejen vazebnou energii, ale i vlnovou funkci, jejíž kvadrát (po nanormování) udává hustotu pravděpodobnosti, se kterou nalezneme proton ve vzdálenosti r od neutronu.

Podrobnější pojednání o úloze (doc. Jiří Dolejší): deuteron.pdf.

Zdrojový kód: Deuteron.py.

Vyzkoušejte si: 1. Pohrajte si s přesností Eulerovy metody a s hodnototu rMax. Pokud tuto hodnotu zvolíte příliš velkou, může se stát, že vám hodnoty vlnové funkce pro velká r přetečou.
2. Pokud v konzoli IDLE (nebo přímo v kódu programu) provedete přiřazení DeuteronVF = DeuteronVFOdeint, můžete zopakovat celý výpočet s integrací integrate.odeint. Porovnejte výsledky.
3. Vazbovou konstantu můžete hledat také tak, že namísto integrálu kvadrátu vlnové funkce minimalizujete hodnotu kvadrátu vlnové funkce ve vzdálenosti rMax. Modifikujte program a porovnejte obě řešení z hlediska výsledku, stability výpočtu a doby výpočtu.

6. cvičení

Náhodná procházka

Náhodnou procházku s fixním krokem ve 2D prostoru naprogramujeme pomocí objektu Turtle. Uvědomte si, že 2D náhodnou procházku získáte tak, že náhodně volíte úhel, pod kterým se vydá Brownova částice. To vám zaručí izotropní rozdělení směrů.

Manuál ke knihovně pro generování náhodných čísel: https://docs.python.org/3/library/random.html.

Zdrojový kód: NahodnaProchazka.py.

Vyzkoušejte si: Střední vzdálenost Brownovy částice od počátku by měla růst jako odmocnina s časem (lze také počítat rozptyl poloh Brownovy částice v závislosti na čase; ten poroste lineárně s časem). Ověřte.

7. cvičení

Minimalizace pomocí náhodné procházky
Náhodná procházka se hodí k jednoduchému algoritmu pro hledání minima funkce: nacházíme-li se na nějaké "vrstevnici" funkce, zkusíme udělat náhodný krok. Pokud krok vede "dolů", tj. směrem k minimu, provedeme ho. Zdrojový kód: Minimalizace.py. Nejprve si vykreslíme zadanou 2D funkci pomocí funkce contourf z knihovny matplotlib. Do ní pak zakreslujeme náhodnou procházku. Minimalizaci ukončíme, pokud stojime krát neprovedeme žádný krok (předpokládáme, že jsme již nalezli minimum).
Metropolisův algoritmus
Pokud má funkce více minim, může se stát, že nás náhodná procházka zavede do lokálního minima, ze kterého se již nedokáže dostat ven. K překonání bariéry oddělující lokální minimum od globálního využijeme Metropolisův algoritmus, který je inspirován termodynamickým Boltzmannovým rozdělením energie: s nenulovou pravděpodobností můžeme při náhodné procházce udělat krok i "do kopce", avšak čím je kopec strmější, tím je pravděpodobnost menší. Minimalizační procedura je rovněž v kódu Minimalizace.py.
Vyzkoušejte si: 1. Pro funkci V1 zkoušejte různé počáteční polohy a různé délky kroku.
2. Pro funkci V2 zkuste nejprve obyčejnou náhodnou procházku a poté náhodnou procházku s Metropolisovým algoritmem. Pokuste se najít optimální hodnotu teploty tak, abychom náhodná procházka skončila v globálním minimumu (přesněji fluktuovala v jeho bezprostředním okolí), avšak už z něj neutekla pryč. V praxi je možné algoritmus modifikovat tak, aby postupně zmenšoval teplotu a procházka nakonec "zamrzla" v globálním minimu.
3. Pohrajte si s Metropolisovým algoritmem v případě funkce V3 (plato na vajíčka). Naimplementujte cyklické okrajové podmínky a nechte algoritmus běžet dlouho. Pro vhodnou velikost teploty vám zvýrazní místa minim.
Náhodná procházka v d-rozměrném prostoru
Izotropní rozdělení směrů v d-rozměrném prostoru je o něco složitější než ve 2D. Možnosti jsou tři:
1. Ve 3D je možné generování pomocí dvou úhlů sférických souřadnic (je nutné uvážit netriviální objemový element pro azimutální úhel). Pro vícerozměrné prostory je tento přístup prakticky nerealizovatelný (vede na problém inverzních funkcí k funkcím daným řadou goniometrických funkcí).
2. Generujeme bod v d-rozměrné hyperkrychli o délce hrany 2 a vybíráme jen ty body, které padnou dovnitř vepsané jednotkové hyperkoule, které následně promítneme na povrch hyperkoule. Pro rostoucí d exponenciálně klesá pravděpodobnost zásahu hyperkoule (pro d=10 se s náhodným bodem z hyperkrychle trefíme do hyperkoule jen jednou ze zhruba 400 pokusů; to znamená, že musíme nagenerovat průměrně 4000 náhodných čísel, abychom mohli udělat jeden krok náhodné procházky).
3. Generujeme d čísel z normálního Gaussovského rozdělení a takto získaný vektor promítneme na povrch d rozměrné hyperkoule. Rozdělení směrů je již izotropní, jak je ukázáno ve článcích
M.E. Muller, A note on a method for generating points uniformly on n-dimensional spheres, Comunications of the Association for Computing Machinery, 2, 19 (1959)
G. Marsaglia, Choosing a Point from the Surface of a Sphere, The Annals of Mathematical Statistics 43, 645 (1972)
Ke generování Gaussovsky rozdělených čísel se dostaneme později; jedna metoda je použít součet dostatečného počtu rovnoměrně rozdělených náhodných čísel.

8. cvičení

Histogram

Histogram udává četnosti číselných hodnot nějaké množiny na předem zadaných intervalech. Pokud je v zadané množině výběr z nějakého náhodného rozdělení a pokud histogram správně nanormujeme, aproximujeme jím hustotu pravděpodobnosti tohoto rozdělení.

Rovnoměrné rozdělení

Rovnoměrné rozdělení na intervalu <0,1) (výběr získáme pomocí funkce random()) má střední hodnotu 1/2 a rozptyl 1/12.

Pro rozdělení náhodné veličiny dané součtem více nezávislých náhodných veličin platí, že jeho střední hodnota je dána součtem středních hodnot a rozptyl je dán součtem rozptylů jednotlivých náhodných veličin. Hustota pravděpodobnosti je pak dáno konvolucí dílčích hustot. To znamená, že součtem dvou rovnoměrných rozdělení na <0,1) získáme rozdělení se střední hodnotou 1, rozptylem 1/6 a hustotou pravděpodobnosti ve tvaru rovnoramenného trojúhelníku na intervalu <0,2).

Generování výběru z normálního rozdělení

Z centrální limitní věty plyne, že pravděpodobnostní rozdělení součtu více náhodných veličin se blíží rozdělení Gaussovskému. V praxi je vhodné zvolit součet 12 rovnoměrných rozdělení na <0,1), jejichž střední hodnota je 6 a rozptyl 1. Pokud tedy sečteme 12 hodnot vybraných z rovnoměrného rozdělení a odečteme 6, získáme číslo, které lze s velkou přesností považovat za výběr z normálního rozdělení N(0,1).

Zdrojový kód: Histogram.py.

9. cvičení

Poissonovo rozdělení

(Rozdělení fluktuací počtu statisticky nezávislých událostí v nějakém časovém nebo prostorovém intervalu, například počtu rozpadů radioaktivního prvku v jednotce času nebo počtu zásahů nějaké oblasti při využití metody Monte Carlo.) Poissonovo rozdělení lze studovat například na fluktuacích počtu hodnot v intervalech histogramu rovnoměrného rozdělení.

Integrace Monte Carlo
Paralelizace Monte Carla

Monte Carlo lze triviálně paralelizovat: počítáme tutéž úlohu ve více vláknech a nakonec všechny výpočty spojíme dohromady. Například v případě integrace Monte Carlo spočítáme aritmetický průměr výsledků ze všech vláken, který odpovídá jednomu Monte Carlo výpočtu s počtem hodů, který je rovnen součtu počtu hodů ze všech vláken. Jednotlivá vlákna výpočtu se sebou nemusí nijak komunikovat (v anglické terminologii se pro takovýto typ problémů používá označení "embarassingly parallel"). Paralelizace na více výpočetních jádrech je v Pythonu implementována v knihovně multiprocessing (pozor, existuje ještě knihovna threading, která rovněž umožňuje paralelní výpočty ve více vláknech, avšak tato knihovna spouští všechna vlákna na jednom jediném výpočetním jádře, a tedy nevede k urychlení výpočtu).

Zdrojový kód: MonteCarlo.py.

10. cvičení

Úvod do psaní textů v LaTeXu

Příklad: dokument.tex + obrázek kubik.pdf (pro pdfLaTeX) nebo kubik.eps.

11. cvičení

Symbolické výpočty

Při symbolických výpočtech namísto s číselnými hodnotami pracujeme s celými výrazy. To zahrnuje kromě základní úpravy výrazů rovněž derivování, integrování a řešení algebraických i diferenciálních rovnic. Nejčastěji používané programy jsou Mathematica (návod na získání a instalaci fakultní licence) a Maple. Pro symbolické matipulace v Pythonu existuje knihovna Sympy.

Zdrojový kód v programu Mathematica: ZakladyMathematica.nb.