Úč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
a nepřeberně přednášek na oboru informatika.
Tři možné způsoby, jak získat zápočet (seřazeno od toho nejupřednostňovanějšího):
| 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, 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é.
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) |
(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.)
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 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. |
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. |
|---|
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.
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.
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ů. |
|---|
Zrychlení můžeme zobrazit také pomocí úseček přímo do okna se želvičkami. Zdrojový kód: StihaciKrivkaZrychleniLimitN.py.
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.
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. |
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. |
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. |
|---|
Počkejte prosím...
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).
| 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. |
| 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. |
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í 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).
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.
(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í.
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.
Příklad: dokument.tex + obrázek kubik.pdf (pro pdfLaTeX) nebo kubik.eps.
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.