FYKOS

Seriál o počítačové fyzice (21. ročník)

Úvod

Tento ročník se budeme zabývat pro vás nadmíru užitečným tématem – využitím počítače pro řešení nespočetného množství úloh, a to nejen fyzikálních problémů. Letošní seriál nás také naučí jinému pohledu na matematický aparát, který využijeme při výpočtech na počítači. To nám mimo jiné umožní snáze vypočítat derivace, integrály, řešit diferenciální rovnice. Zkrátka téměř cokoliv.

Sami jistě cítíte, že si budeme moci hrát s jakýmkoliv fyzikálním systémem, omezeni budeme jen naší silou a rychlostí procesorů. A i kdyby si kvůli vašim budoucím zběsilým výpočtům vaše mašinky přály raději odejít do křemíkového nebe, můžeme jejich práci rozdělit mezi více procesorů a spustit je na superpočítači, který běží v našem univerzitním výpočetním středisku.

Prvním přírodním jevem, který si vezmeme na mušku, je gravitační působení a jeho účinky. Důvody jsou k tomu dva. Klasický popis je pro každého jednoduše pochopitelný. Avšak již úloha, kdy na sebe gravitačně působí tři volná tělesa, například tři hvězdy, které kolem sebe obíhají, není analyticky řešitelná. Co to znamená? Že například dráhy těchto těles, dobu oběhu apod. můžeme vypočítat pouze použitím nějaké přibližné numerické metody (neděste se prosím vás tohoto sprostého slova, naopak řešit něco numericky, tedy číselně, je principiálně jednodušší postup než obecný analytický, tedy odvozování vzorečků) a zpravidla jen na počítači. Vývoj celého systému je třeba počítačově nasimulovat.

Takto modelovat se dá cokoliv. Bloudění opilého námořníka, požár v pralese, růst buněk, kvantová turbulence v supratekutém héliu (pokud byste to zvládli), časový vývoj naší Sluneční soustavy či crash test automobilu. Většinu času budeme věnovat simulacím čistě fyzikálním. Výběr témat můžete hodně ovlivnit sami, pokud se s námi během roku podělíte o své choutky.

A v čem to budeme dělat? Málokdo tuší, kolik se toho dá simulovat už v tom „blbém“ Excelu. Bez problému nám v něm může rotovat celá hvězdokupa. Navíc díky FYKOSímu textu Úvod do programování, si během hodinky úplně každý, kdo umí na počítači alespoň nainstalovat nějaký program, nainstaluje překladač Pascalu a naučí se nejnutnějším základům vyššího programování.

Těšíme se na vaše řešení. Uvítáme rovněž jakékoliv ohlasy. Napište nám prosím, co vám jde, co vám nejde a co byste se rádi dozvěděli! Uvědomujeme si, že je třeba si osvojit spoustu dovedností s počítačem, které se vám však budou ve fyzice i životě nesmírně hodit. A jako se vším je nejtěžší začít. Proto, abychom vás neodradili, a naopak vám umožnili se spoustu skvělého naučit, tak další díly se pokusíme ušít na míru vašim ohlasům. Jeho převážná většina může být pro ty pokročilejší opakováním, nicméně ti zajisté uznají, že dobře si osvojit esenciální základy programování je důležité. Věříme, že pak už to půjde všem místo po banánové slupce jako po másle.

Chtěl bych zmínit pár dobrých zdrojů. V knihovničce Fyzikální olympiády je k nalezení text Šedivý, P.: Modelování pohybů numerickými metodami. Dále se numerickými metodami zabývá seriál osmého ročníku FYKOSu. Určitě proto neváhejte navštívit náš Archiv ročenek. Dále mohu doporučit Houmpejdž Tomáše Ledvinky, počítačového fyzika z Ústavu teoretické fyziky, jenž na Matematicko-fyzikální fakultě přednáší Programování pro fyziky. Naším tématem se zabývá i devátá kapitola prvního dílu Feynmanových přednášek z fyziky.

Text seriálu

  • Kapitola 1: Gravitace, strana 3
  • Kapitola 2: Parciální diferenciální rovnice, strana 2
  • Kapitola 3: Náhodné jevy, strana 16
  • Kapitola 4: Kvantová mechanika, strana 27
  • Kapitola 5: Numerické integrování, strana 29
  • Kapitola 6: To nejlepší nakonec, strana 23

Soubory ke stažení

  • Stáhněte si free pascal.
  • Učte se programovat s naším studijním textem Úvod do programování. Přiložené soubory:
    • rovnice_regula.pas – řešení rovnic metodou regula falsi,
    • minimum_zlatyrez.pas – hledání minima funkce metodou zlatého řezu,
    • konjgrad.pas – hledání minima funkce dvou proměnných metodou konjugovaného gradientu,
    • zihaniTSP.pas – řešení problému obchodního cestujícího metodou simulovaného žíhání a soubor mesta.txt s polohami měst.
  • Vytvořili jsme pro vás v Excelu dva programy laplace.xls a poisson.xls vyšetřující elektrický potenciál jednoduchého dvojrozměrného problému se zadanými okrajovými podmínkami (první z nich v nepřítomnosti nábojů, řeší tedy Laplaceovu rovnici; druhý pak pro zadané rozložení náboje řeší Poissonovu rovnici).
  • Krátký text Nabliáda o trojúhelníčcích.
  • Několik poznatků ze statistické fyziky (Maxwell-Boltzmannovo rozdělení, Boltzmannův vzorec, entropie).
  • Užitečným nástrojem pro numerickou integraci je takzvaná Rombergova metoda. Algoritmus realizující Rombergovu integraci je pro vaši inspiraci k dispozici v demonstračním programu romberg.pas.
  • Přikládáme dva animované GIFy (vizualizace1.gif a vizualizace2.gif) k metodě raytracingu – Vzhled trojrozměrné scény, jak ji vidí relativistický pozorovatel v závislosti na rychlosti v/c, resp. veličině arctg(v/c).
  • Prohlédněte si animované GIFY (animace 1a, animace 1b, animace 2a, animace 2b, animace 3a, animace 3b, animace 4a, animace 4b) s  různými animacemi černé díry v blízkosti Země.
  • Prohlédněte si animované GIFY (animace 1, animace 2, animace 3) s kmitajícím dvojkyvadlem.

Přílohy k řešení úloh

  • K řešení úlohy I.S: pascalovské programy gravitace.pasTriHvezdy2.pas.
  • K řešení úlohy II.S: programy uran.xls, drat.pas a krychle.pas.
  • K zadání úlohy III.S: využijte jednotku OdporovaSit.pas pro výpočet odporu elektrického obvodu. Její použití demonstrujeme v programu vypOdporu.pas, který počítá odpor krychle.
  • K řešení úlohy III.S: pí-obvod – zihani_odporu.pas, epidemie v Praze – epidemie.zip.
  • K řešení úlohy IV.S: Program schrodinger.pas pro řešení časové Schrödingerovy rovnice popisující pohyb vlnového balíku v kvadratickém potenciálu (harmonický oscilátor).
  • K řešení úlohy V.S: Programy elpoz.pas (model termodynamické rovnováhy fotonového, elektronového a pozitronového plynu), hvezda.pas (model hvězdy tvořené degenerovaným elektronovým plynem – bílý trpaslík).
  • K řešení úlohy VI.S: Programy raytracing_console.pas, chaos_dvojkyv.pas (simulace pohybu dvojkyvadla a výpočet Lyapunovova koeficientu).
Facebook icon

Úloha VI . S … na přání

Pokuste se o řešení libovolného problému z šesté kapitoly seriálu.
Řešení: PS PDF
Facebook icon

Úloha V . S … posloupnosti, horká dutina a bílý trpaslík

  1. Odvoďte Taylorův rozvoj exponenciály a pro x = 1 graficky znázorněte posloupnost částečných součtů řady ∑k = 11 ⁄ k! spolu s posloupností { ( 1 − 1 ⁄ n ) n}n=1,2,....

    Stejným způsobem porovnejte posloupnost { ( 1 − 1 ⁄ n ) n}n=1,2,... a posloupnost částečných součtů řady ∑k = 1xk ⁄ k!, čili posloupnost {∑k = 1nxk ⁄ k!}n=1,2,..., tentokráte však pro x = −1.
  2. Druhý úkol bude určit závislost koncentrace elektronů a pozitronů na teplotě při celkovém náboji Q = 0 v prázdné uzavřené horké dutině. (Bude-li se vám chtít, i při jiných vámi zvolených hodnotách Q.) Dále určete závislost poměru vnitřní energie Ue elektronů a pozitronů ku celkové vnitřní energii systému U (tj. součtu energie elektromagnetického záření a částic) na teplotě a určit hodnoty teploty odpovídající některým význačným hodnotám tohoto poměru (např. 3 ⁄ 4, 1 ⁄ 2, 1 ⁄ 4, …; může tento poměr nabývat všech těchto hodnot?).

    Pokuste se své výsledky pěkně graficky zpracovat ve formě grafů (můžete zkusit i trojrozměrné).

    Při vašem snažení vám může hodně pomoci, pokud si zavedete vhodné bezrozměrné jednotky (např. βE0 místo β apod.).
  3. Řešte soustavu diferenciálních rovnic pro M ( r ) a ρ ( r ) v modelu bílého trpaslíka pro několik vhodně zvolených hodnot ρ ( 0 ) a pro každou z nich sledujte hodnotu, ke které se blíží M ( r )  při r→∞. Ta je zřejmě rovna hmotnosti celé hvězdy. Pokuste se prozkoumat závislost této celkové hmotnosti na ρ ( 0 ) a odhadnout její horní mez. Srovnejte váš výsledek s horní mezí hmotnosti bílého trpaslíka, kterou najdete v literatuře nebo na internetu. Uvažujte, že je hvězda tvořena héliem.
Řešení: PS PDF
Facebook icon

Úloha IV . S … kvantový harmonický oscilátor

Modelujte časový vývoj vlnové funkce částice, kterou umístíme do potenciálu Vx ) = ½ kx² a která je v čase τ = 0 popsána vlnovou funkcí

ψR ( X, 0 ) = exp ( −(( X − X0 )²) ⁄ 4 ) ,
ψ
I ( X, 0 ) = 0 .

Jedná se tedy o vlnový balík se středem mimo počátek. Prozradíme vám, že jde o tzv. koherentní stav harmonického oscilátoru a vlnový balík by měl harmonicky kmitat kolem počátku s úhlovou frekvencí √( k ⁄ ) stejně jako klasická částice.

Pokud se vám toto podaří namodelovat, můžete vyzkoušet, jak se budou chovat vlnové balíky o jiné šířce (tedy se jmenovatelem v exponenciále odlišným od čtyř), případně jak bude situace vypadat při jiném průběhu potenciálu.

Řešení: PS PDF
Facebook icon

Úloha III . S … bloudění námořníka, pí-obvod a epidemie v Praze

Integrál

Integrujte metodou Monte Carlo funkci e − x² na intervalu [ − 100,100 ] . Zkuste také numericky určit hodnotu tohoto integrálu od −∞ do +∞.

Návod: Funkce je symetrická vůči počátku, čili ji stačí integrovat na intervalu [ 0, +∞ ) . Proveďte substituci x = 1 ⁄ t − 1, čímž změníte meze integrálu od 0 do 1.

Bloudění námořníka

Opilý námořník vstoupil na molo dlouhé 50 kroků a široké 20 kroků. Jde směrem k pevnině. Při každém kroku dopředu však zavrávorá zároveň o krok nalevo nebo napravo. Zjistěte, s jakou pravděpodobností námořník dojde až na břeh a s jakou pravděpodobností spadne do moře a utopí se.

Námořník měl štěstí a neutopil se. Druhou noc se však opět vydává opilý z lodi na pevninu. Tentokrát však vane stálý vítr o rychlosti 3 m·s−1, který způsobí to, že na jednu stranu udělá krok s pravděpodobností 0,8 a na druhou stranu s pravděpodobností 0,2. Zjistěte, s jakou pravděpodobností námořník dojde až na břeh a s jakou pravděpodobností spadne do moře a utopí se.

Třetí noc se námořník opět vydává opilý na pevninu. Tentokrát však vane proměnlivý vítr. Vane podle normálního rozdělení se střední hodnotou 0 m·s−1 a disperzí 2 m·s−1. Zjistěte, s jakou pravděpodobností tentokrát námořník dojde až na břeh a s jakou pravděpodobností spadne do moře a utopí se. Můžete uvažovat, že námořník jde pomalu a setrvačnost větru lze zanedbat. Komu by to vadilo, nechť vymyslí, jak by vítr v po sobě jdoucích krocích koreloval.

Pí-obvod

Máme k dispozici 50 rezistorů o odporech 50 Ω a chceme z nich sestavit obvod, jehož celkový odpor v ohmech bude co nejblíže číslu π. Pokuste se metodou simulovaného žíhání najít obvod, který by tomuto požadavku vyhovoval co nejlépe.

Pro určování celkového odporu obvodu si můžete přizpůsobit program, který najdete na našich webových stránkách.
Pokud se na tento úkol necítíte, můžete zkusit zahrnout do problému obchodního cestujícího zakřivení zemského povrchu a pokusit se jej vyřešit pro nějakou konkrétní množinu měst na Zemi (například všechna hlavní města v Evropě, USA atd.).

Epidemie v Praze

Zkoumejte vývoj epidemie v Praze, uvažujte 1 milión obyvatel. Intenzita nákazy β je 0,4 ⁄ 1 000 000 za den, uzdravení γ je ( čtyřidny )1. Na počátku je nakaženo 100 lidí. Porovnejte průběh epidemie při očkování předem dvaceti procent lidí s průběhem epidemie při očkování až během epidemie s rychlostí půl procenta denně. A také s průběhem bez očkování. Konec epidemie vyhlásíme, bude-li méně jak 20 lidí nemocných.

Je spousta údajů, které můžete z počítačové simulace získat. Krom středovaného průběhu epidemie uveďte pro zajímavost též graf, kde ukážete prvních pět náhodných simulací. Dále můžete sledovat fluktuace. Můžete též výsledky porovnat s deterministickým modelem, když neuvažujete náhodnost nakažení. Těžištěm hodnocení bude, kolik různých zajímavých dat dokážete hezky zpracovat.

Řešení: PS PDF
Facebook icon

Úloha II . S … porcování divokých rovin

Skladování uranu

Palčivá otázka jaderné energetiky je skladování vyhořelého radioaktivního paliva. Většinou se skladuje ve válcových článcích ponořených ve vodní lázni, která drží jejich povrch na konstantní teplotě asi 20 °C. Na vás je nyní zjistit, jaké bude rozložení teploty v článcích tvaru kvádru se čtvercovou podstavou o hraně délky 20 cm. Článek bude poměrně vysoký a proto nás zajímá rozložení tepla v příčném řezu. Uran bude zaujímat koncentrický kvádr se čtvercovou podstavou o hraně 5 cm. Ze zkušenosti s válcovými kapslemi víme, že bude mít konstantní teplotu okolo 200 °C.

Zahřívající se drát

Máme velmi dlouhý drát kruhového průřezu o poloměru r z materiálu o tepelné vodivosti λ a měrné elektrické vodivosti σ. Přiložíme na něj konstantní elektrické napětí. Nechť je intenzita elektrického pole (tj. napěťový spád) uvnitř drátu konstantní, rovnoběžná s jeho osou a její velikost budiž E. Pak drátem bude procházet proud o plošné hustotě j = σE a bude se vytvářet Jouleovské teplo s objemovým výkonem p = σE².

Protože materiál drátu má nenulovou tepelnou vodivost, vytvoří se v něm jisté rovnovážné rozložení teploty, které – jak víme – splňuje Poissonovu rovnici λΔT = −p. Předpokládáme, že okraj drátu udržujeme na dané teplotě T0. Tím máme dánu okrajovou podmínku potřebnou k vyřešení rovnice. Vzhledem k symetrii problému se můžeme omezit na její řešení pouze ve dvou rozměrech – na průřezu vodiče (teplota jistě nebude záviset na posunutí podél osy vodiče). Nyní by již bylo jednoduché problém vyřešit popsanými metodami.

My si však situaci maličko zkomplikujeme a budeme předpokládat (zcela oprávněně), že měrná elektrická vodivost σ závisí na teplotě. Budeme tedy mít rovnici typu ΔT = f(T).

Pokuste se tuto rovnici numericky vyřešit pro nějakou danou závislost vodivosti na teplotě (můžete si ji najít v literatuře, na internetu, nebo si klidně nějakou vymyslet) a najít tak rozložení teploty na průřezu drátu. Můžete se pokusit měnit intenzitu elektrického pole E a nakreslit voltampérovou charakteristiku drátu, vyzkoušet více druhů závislostí σ(T) (třeba pro polovodič, jehož vodivost s rostoucí teplotou na rozdíl od obyčejného kovu roste) atd.

Vaší iniciativě samozřejmě meze neklademe a těšíme na pěkné nápady.

Kapacita krychle

Vypočítejte kapacitu dokonale vodivé krychle o straně délky 2a. Pokud se budete nudit, můžete zkusit kvádr (a třeba závislost kapacity na délkách jednotlivých stran), případně jiné geometrické objekty.

Nápověda. Kapacita je poměr náboje na krychli rozmístěného ku potenciálu povrchu krychle (za předpokladu, že potenciál v nekonečnu je nulový). Problém tedy lze řešit tak, že si zvolíme libovolně potenciál krychle, vyřešíme Laplaceovu rovnici Δφ = 0 vně krychle a vypočítáme celkový náboj na krychli užitím Gaussova zákona (tj. určením intenzity elektrického pole derivováním potenciálu a výpočtem jeho toku vhodně zvolenou plochou obklopující krychli).

Úplným řešením je vymyšlení vhodného fyzikálního modelu, návrh jeho numerického řešení a realizace této úlohy na počítači. Bodově ohodnotíme, pokud úlohu fyzikálně rozvážíte a okomentujete. Nějaký bod by se našel i za návrh algoritmu, který byste rádi počítači předložili.

Řešení: PS PDF
Facebook icon

Úloha I . S … gravitace

Krom programu Planeta.xls prosím stáhněte též DveHvezdy.xlsTriHvezdy.xls. Může se vám hodit přečíst již zmíněný dokument Úvod do programování. Díky němu budete určitě schopni stáhnout a přeložit program Planeta.pas, DveHvezdy.pasTriHvezdy.pas. Cvičením a přípravou na úlohy je tyto programy pochopit a lehce si s nimi pohrát, zkusit si do nich „zašťourat“ či je lehce poupravit.

  1. Úkolem prvním je obohatit alespoň dva programy o něco svého nebo je upravit podle svého. Například nejrůzněji pozměňte počáteční podmínky a hmotnosti. Nebo k systému přidejte další planetu či další hvězdu. Také můžete vyzkoušet pozměnit gravitační zákon a počítat se silou = A ⁄ R² + B ⁄ R³, kde AB jsou pevně zvolené konstanty apod.
  2. Uvažujte dvě stejně těžké hvězdy, které kolem sebe obíhají po kružnici. Po ose této kružnice se k nim začne náhle přibližovat hvězda třetí, která má na začátku stejnou rychlost, jakou se pohybují hvězdy obíhající, a rovněž sdílí i jejich hmotnost. Počítačově nasimulujte, co se bude dít.

Jako řešení obou úloh nám prosím zašlete obrázky s bohatým komentářem. Uveďte alespoň stručné vysvětlení, co to na těch obrázcích je. Dále, jakým způsobem jste při výpočtu postupovali a pomocí kterého výpočetního systému jste jej provedli.

Řešení: PS PDF
©FYKOS – webmaster@fykos.cz