Zvýšení kvality vzdělávání na Slezské univerzitě v Opavě ve vazbě na potřeby Moravskoslezského kraje
CZ.02.2.69/0.0/0.0/18_058/0010238
Gabriel Török et al
Studijní opora pro předmět
Numerické metody v přírodních
vědách
Fyzikální ústav
Slezská univerzita v Opavě
Opava
20. 6. 2020
Numerické metody v přírodních vědách
doc. RNDr. Gabriel Török, Ph. D
Mgr. Debora Lančová, RNDr. Jiří Horák, Ph.D, Mgr. Kateřina Klimovičová,
Mgr. Gabriela Urbancová, Mgr. Martin Urbanec, Ph.D, Mgr. Adam Hofer,
Prof. RNDr. Vladimír Karas, DrSc. a RNDr. Martin Blaschke, Ph.D
Dostupné na: physics.cz
Tato studijní opora vznikla v rámci projektu
Zvýšení kvality vzdělávání na Slezské univerzitě v Opavě
ve vazbě na potřeby Moravskoslezského kraje,
registrační číslo CZ.02.2.69/0.0/0.0/18-058/0010238
Fyzikální ústav
Slezská univerzita v Opavě
Bezručovo nám. 13, Opava
Autor šablony: RNDr. Šárka Vavrečková, Ph.D
Sázeno v systému LATEX
Předmluva
Co najdeme v těchto skriptech
Tato studijní opora je určena studentům předmětu Numerické modelování v přírodních vědách. Cílem tohoto
předmětu je ukázat použití různých numerických metod v praxi, zejména při řešení reálných fyzikálních
problémů. Studijní opora je sestavena z 9 kapitol, kdy každá kapitola se věnuje jiné oblasti z fyziky.
Předmět numerické metody v přírodních vědách má za cíl ukázat použití numerických metod v praxi, sestává
tedy, stejně jako tato opora, z jednotlivých příkladů, a spíše než teoretické znalosti zde najdete použití
různých metod v praxi. K studijní opoře je tedy nutné přistupovat jako ke sbírce řešených i neřešených
příkladů. Ve většině kapitol naleznete také ukázku řešení části problémů v jazyce Python, případně Wolfram
Mathematica.
Vzhledem ke komplexnosti zadání sestává opora z úloh, které připravil široký autorský kolektiv. Některé
z nich jsou dílem studentů doktorského studia Fyzikálního ústavu, jiné vytvořili odborní asistenti zde
působící. Další kapitoly byly napsány vědeckými pracovníky z Astronomického ústavu Akademie věd ČR.
Proto je každá kapitola svým charakterem do jisté míry odlišná. Autoři jednotlivých kapitol budou dále
nápomocní při řešení úloh, případně poskytnou další zadání pro semestrální práce.
Nedílnou součástí studijní opory jsou Jupyter a Mathematica notebooky, které Vám poskytne vyučující
předmětu.
Značení
Ve skriptech se používají následující barevné ikony:
• Rychlý náhled (skript, kapitoly), ve kterém se dozvíme, o čem to bude.
• Klíčová slova kapitoly.
• ºº Cíle studia pro kapitolu nám řeknou, co nového se v dané kapitole naučíme.
•  Nové pojmy, značení apod. jsou označeny modrým symbolem, který vidíme zde vlevo. Tuto ikonu
(stejně jako následující) najdeme na začátku odstavce, ve kterém je nový pojem zaváděn.
• ££ Konkrétní postupy a nástroje, způsoby řešení různých situací, do kterých se může správce počítačového
vybavení dostat, atd. jsou značeny také modrou ikonou.
iii
iv
•  Některé části textu jsou označeny fialovou ikonou, což znamená, že jde o nepovinné úseky, které
nejsou probírány (většinou; studenti si je mohou podle zájmu vyžádat nebo sami prostudovat). Jejich
účelem je dobrovolné rozšíření znalostí studentů o pokročilá témata, na která obvykle při výuce
nezbývá moc času.
• Žlutou ikonou jsou označeny odkazy, na kterých lze získat další informace o tématu. Nejčastěji
u této ikony najdeme webové odkazy na stránky, kde se dané tématice jejich autoři věnují podrobněji.
• Červená je ikona pro upozornění a poznámky.
Pokud je množství textu patřícího k určité ikoně větší, je celý blok ohraničen prostředím s ikonami na
začátku i konci, například pro definování nového pojmu:
 Definice
V takovém prostředí definujeme pojem či vysvětlujeme sice relativně známý, ale komplexní pojem s více
významy či vlastnostmi.

Podobně může vypadat prostředí pro delší postup nebo delší poznámku či více odkazů na další informace.
Mohou být použita také jiná prostředí:
M Příklad
Takto vypadá prostředí s příkladem, obvykle nějakého postupu. Příklady jsou obvykle komentovány, aby
byl jasný postup jejich řešení.
M
C Úkol
Otázky a úkoly, náměty na vyzkoušení, které se doporučuje při procvičování učiva provádět, jsou uzavřeny
v tomto prostředí. Pokud je v prostředí více úkolů, jsou číslovány.
C
Souhrn kapitoly
Na konci každé kapitoly je souhrn, ve kterém si ujasníte, o čem daná kapitola byla a co všechno jste se
naučili.
Obsah
Předmluva iii
1 Zákon lomu a Duhová funkce 1
1.1 Jak vzniká duha - fyzikální základ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Lom světla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 Úplný odraz, mezní úhel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.3 Barvy duhy, disperze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.4 Průchod paprsku vodní kapkou . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Sekundární duha a Alexandrův tmavý pás . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Implementace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3.1 Vykreslení duhové funkce 1. a 2. řádu . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3.2 Numerické derivace, hledání minima a maxima funkce . . . . . . . . . . . . . . . . . 9
1.3.3 Hledání kořene - metoda půlení intervalů . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Modelování dat 13
2.1 Metoda nejmenších čtverců . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.1 Lineární regresní funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.2 Příklad z reálného světa - Určení stáří vesmíru . . . . . . . . . . . . . . . . . . . . . 17
2.2 Příklady v Pythonu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Interpolace a extrapolace 22
3.1 Co je to interpolace a extrapolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Polynomiální interpolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.1 Lineární interpolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.2 Lagrangeův interpolační polynom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Interpolace pomocí splajnu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3.1 Lineární interpolační splajn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3.2 Kubický splajn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 Interpolace v programu Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 Numerická řešení rovnic, hledání kořenů 33
4.1 Hledání kořenů kvadratické funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1.1 Hledání kořenů metodou půlení intervalů . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.1.2 Hledání kořenů Newtonovou – Raphsonovou metodou . . . . . . . . . . . . . . . . . 36
4.2 Pohyb částic v poli dvou kolem sebe obíhajících objektů . . . . . . . . . . . . . . . . . . . . 37
v
vi
4.2.1 Teorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.2.2 Implementace v Pythonu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5 Statistický popis dat, distribuční funkce 42
5.1 Základní charakteristiky popisné statistiky . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.1.1 Minimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.2 Maximum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.3 Počet hodnot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.4 Aritmetický průměr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.5 Medián . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.6 Modus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.7 Kvantily . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.1.8 Rozpětí . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.1.9 Mezikvartilové rozpětí a odchylka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.1.10 Rozptyl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.1.11 Směrodatná odchylka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.1.12 Variační koeficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2 Distribuce pravděpodobnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.1 Náhodná veličina . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.2 Pravděpodobnostní funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.3 Distribuční funkce diskrétní veličiny . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.4 Distribuční funkce spojité veličiny . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.2.5 Hustota pravděpodobnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2.6 Centrální limitní věta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2.7 Normální rozdělení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6 Řešení okrajové úlohy relaxační metodou 50
6.1 Okrajová úloha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.1.1 Diskretizace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.1.2 Iterativní řešení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.2 Implementace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Příklad: Gorilla problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.3 Numerické nebo automatické derivace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
7 Bondiho akrece 59
7.1 Akreční procesy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.2 Sférická akrece . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.2.1 Bezrozměrná verze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.2.2 Kritický bod a podmínka regularity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7.3 Analytické řešení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
7.3.1 Asymptotické chování řešení: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.3.2 Grafické znázornění analytických výsledků . . . . . . . . . . . . . . . . . . . . . . . . 63
7.4 Řešení s použitím relaxační metody, řešení problému vlastních hodnot na neznámé oblasti
řešení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7.4.1 Shrnutí problému . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7.4.1.1 Úprava problému vlastních hodnot . . . . . . . . . . . . . . . . . . . . . . . 66
7.4.1.2 Co udělat s vnitřní okrajovou podmínkou? . . . . . . . . . . . . . . . . . . 67
7.4.1.3 Obecná poznámka k síti souřadnic . . . . . . . . . . . . . . . . . . . . . . . 68
vii
7.5 Implementace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
7.5.1 Zkoumání prostoru parametrů . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
8 Některé aproximační metody řešení diferenciálních rovnic 75
8.1 Řešení v podobě Talorovy řady generované v Mathematice . . . . . . . . . . . . . . . . . . . 75
Řešený příklad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.2 Poissonova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
Řešený příklad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.3 Lindstedtova metoda (Poincarého-Lindstedtova metoda) . . . . . . . . . . . . . . . . . . . . 79
Řešený příklad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
9 Spektrální profily 83
Spektrální profil emisní čáry vznikající v rotujícím prstenci . . . . . . . . . . . . . . . . . . . . . 83
Literatura 87
Kapitola 1
Zákon lomu a Duhová funkce
Debora Lančová, Gabriel Török
Rychlý náhled: Tato kapitola se zabývá odvozením rovnice primární a sekundární duhu, jejich vykreslením
a nalezením maxima a minima. Jsou zde popsány základní metody numerické derivace a hledání
kořenů funkce.
Klíčová slova: Numerická derivace, hledání kořenů funkce, metoda půlení intervalů, bisekce, duha,
zákon lomu, duhová funkce, difrakce
ºº Cíle studia: Po prostudování kapitoly budete schopní v Pythonu naprogramovat jednoduchou numerickou
derivaci funkce jedné proměnné, použít a implementovat metodu půlení intervalů (bisekce) –
společně budete schopni najít extrémy funkcí. Seznámíte se také s tím, jak se vykreslují jednoduché grafy
pomocí knihovny matplotlib.
Duha je optický jev, který vzniká odrazem slunečního světla na kapkách deště. Vzniká jen za specifických
podmínek, kapky musí mít správnou velikost, Slunce musí být nízko nad obzorem, a jiné.
1.1 Jak vzniká duha - fyzikální základ
Abychom mohli nastínit základní princip vzniku duhy, musíme nejprve zavést některá zjednodušení - například,
že dešťové kapky jsou kulové. To víceméně platí pro malé kapky, o velikosti 0.5 − 0.7 mm, jejichž tvar
je výrazně ovlivněn povrchovým napětím. Předpokládejme tedy, že Sluneční paprsky dopadají na kulatou
kapku vody a při dopadu se lámou. Bude přitom používat pouze zákony geometrické optiky.
1.1.1 Lom světla
 Zákon lomu, jinak také Snellův zákon, říká, jakým způsobem se mění chod optického paprsku při
přechodu mezi optickými prostředími. Optické vlastnosti prostředí jsou vyjádřeny indexem lomu n, který
říká, o kolik pomaleji se v něm šíří světlo oproti rychlosti šíření ve vakuu
n =
c
v
, (1.1)
kde v je rychlost šíření světla v daném prostředí. Indexy lomu v různých prostředích jsou shrnuty v tabulce
1.1.
1
Kapitola 1 Zákon lomu a Duhová funkce 2
Obrázek 1.1: Dvojitá duha viditelná z budovy Fyzikálního ústavu v Opavě, včetně výrazného Alexandrova
tmavého pásu
Tabulka 1.1: Indexy lomu v různém prostředí.
Prostředí Vzduch Voda Sklo Diamant
Index lomu 1.00 1.33 1.5 2.42
 Podle Fermatova principu nejkratšího času se světlo šíří po takové dráze, která odpovídá nejkratšímu
času. Přechází-li tedy z jednoho prostředí do druhého, ve kterých se liší rychlost šíření, bude se řídit
tímto principem. Zákon lomu má tvar
n1 sin α = n2 sin β, (1.2)
kde nesmíme zapomenout, že úhly α a β se měří od kolmice (jak je nakresleno na obrázku 1.2).
1.1.2 Úplný odraz, mezní úhel
Mezní úhel dopadu je takový úhel, při němž se světlo láme pod úhlem 90◦ od kolmice. Pře jeho překročení
se tedy všechny dopadající paprsky odrážejí - dochází k úplnému odrazu.
Odvození mezního úhlu, který značíme zpravidla αm je jednoduché. Chceme-li, aby β = 90◦, platí
sin β = 1 a pak
sin αm = n2/n1, (1.3)
C Úkol
Spočítejte mezní úhly pro úplný odraz paprsku, který přechází ze vzduchu do vody a naopak.
C
Kapitola 1 Zákon lomu a Duhová funkce 3
Obrázek 1.2: Snellův zákon – zákon lomu
Tabulka 1.2: Indexy lomu světla ve vodě v závislosti na vlnové délce.
Vlnová délka Barva Index lomu ve vodě Index lomu ve skle
396,85 tmavě fialová 1.343 1,471
430,78 fialová 1.341 1,468
486,14 modrá 1.337 1,464
527,00 zelená 1.335 1,461
589,30 žlutá 1.333 1,459
656,28 oranžová 1.331 1,457
686,72 červená 1.330 1,456
760,82 tmavě červená 1.329 1,454
1.1.3 Barvy duhy, disperze
Disperze znamená rozklad bílého světla na jeho jednotlivé barevné složky. K tomu dochází při lomu nebo
odrazu díky tomu, že pro každou vlnovou délku je hodnota indexu lomu jiná - viz tabulka 1.2. Platí, že
pro rostoucí frekvenci klesá rychlost šíření v daném prostředí, a tím pádem také index lomu. Nejznámějším
příkladem disperze je lom na optickém hranolu, kde dochází k rozkladu světla hned dvakrát - obrázek 1.3.
1.1.4 Průchod paprsku vodní kapkou
Nyní si tedy můžeme představit, jak vypadá situace, kdy na kulovou vodní kapku dopadne paprsek světla
(obrázek 1.4). Mějme vodní kapku o poloměru r, do které v bodě A vstupuje paprsek ve vzdálenosti h od
osy procházející středem S (veličině h také říkáme impaktní parametr). Paprsek dopadá pod úhlem α od
kolmice (přímky spojující body S a A), láme se na povrchu kapky pod úhlem β. Poté se odrazí od protější
stěny kapky v bodě B, dopadá opět na povrch v bodě C, láme se a opouští kapku.
Zajímá nás, pod jakým úhlem se paprsek od kapky odrazí - tedy úhel γ. Když se pečlivě podíváte na
Kapitola 1 Zákon lomu a Duhová funkce 4
Obrázek 1.3: Rozklad světla na optickém hranolu [Wikimedia]
Obrázek 1.4: Průchod paprsku vodní kapkou a vznik duhy
obrázek, uvidíte, že platí
γ = 4β − 2α. (1.4)
C Úkol
Odvoďte rovnici (1.4) z obrázku 1.4 tak, že si z ∆ABS spočítáte velikost úhlu δ
a pak z ∆APS úhel γ.
C
Úhly α a β můžeme zapsat následovně
sin α =
h
r
(1.5)
sin β =
1
n
sin α =
1
n
h
r
, (1.6)
kde vztah pro úhel β vyplývá ze zákonu lomu.
Kapitola 1 Zákon lomu a Duhová funkce 5
Obrázek 1.5: Průchod paprsku vodní kapkou, detailní označení úhlů a bodů.
 Definice
Po dosazení rovnic (1.5) a (1.6) do rovnice (1.4) získáme finální vztah
γ = 4arcsin
h
nr
− 2arcsin
h
r
, (1.7)
který označujeme jako Duhová funkce 1. řádu, která popisuje závislost úhlu, pod kterým se paprsek od
vodní kapky odrazí.

Obrázek 1.6: Princip vzniku duhy, jak je popsal ve své knize Discours de la méthode z roku 1637 René
Descartes
Kapitola 1 Zákon lomu a Duhová funkce 6
1.2 Sekundární duha a Alexandrův tmavý pás
Často můžeme pozorovat duhy hned dvě, přičemž druhá duha bude mít převrácené barvy a bude mnohem
méně výrazná (jako na obrázku 1.1). Sekundární duha vniká po dvojitém odrazu paprsku uvnitř kapky,
přičemž paprsek do ní vstupuje pod osou a vystupuje nad osou (obrázek 1.6).
 Definice
Duhová funkce pro sekundární duhu má tvar
γ2 = π − 6arcsin
h
nr
+ 2arcsin
h
r
. (1.8)

C Úkol
Zkuste si vzorec pro sekundární duhu odvodit sami, stejným postupem jako v případě primární duhy.
C
 Mezi primární a sekundární duhou se nachází tzv. Alexandrův tmavý pás, který je mnohem tmavší
než zbytek oblohy. Proč tomu tak je, uvidíte, až si vykreslíme primární a sekundární duhovou funkci.
1.3 Implementace
Naším úkolem bude implementovat funkce pro primární a sekundární duhu v Pythonu, vykreslit si je pro
různé vlnové délky a nalezení maxima a minima pomocí derivace.
1.3.1 Vykreslení duhové funkce 1. a 2. řádu
Nejprve si pouze vykreslíme obě duhové funkce, abychom viděli jejich průběhy a mohli si tak lépe představit
vznik duhy.
Začněme importováním knihoven Numpy pro numerické výpočty a Matplotlib pro tvorbu grafů:
1 import numpy as np
2 from matplotlib import pyplot as plt
Nyní můžee definovat a vykreslit funkci pro jednu vlnovou délku, tedy jeden index lomu. Definujme
také průměr kapky, typicky d = 0.5 mm.
1 #index lomu světla
2 n = 1.3333
3 #typicky prumer destove kapky je 0.5 - 0.7 mm
4 r = 2.5 * 1e-4
5 #duhova funkce
6 def gama(x,r):
7 foo = 4.0 * np.arcsin(1/n * x/r)
8 bar = 2.0 * np.arcsin(x/r)
9 return foo - bar
Kapitola 1 Zákon lomu a Duhová funkce 7
Chceme vykreslit graf závislosti γ na impaktním parametru h, musíme tedy vytvořit array s hodnotami h
od 0 (horizontální osa kapky) do h = r, využijeme funci linspace z knihovny numpy. Pak do další array
gamy uložíme hodnoty funkce gama(h,r) pro všechny prvky h:
1 # impaktni parametry od nuly do r, celkem 100 hodnot
2 h = np.linspace (0,r, num = 100)
3 gamy = gama(h,r)
Nyní už můžeme primární duhovou funkci vykreslit. Avšak, nesmíme zapomenout na převody jednotek.
h jsme si zaveli v metrech, knihovna numpy počítá gonieometrické funkce v radiánech - převedeme si je do
stupňů:
1 plt.plot(h*1000,gamy*180/np.pi)
2 plt.xlabel('h (mm)')
3 plt.ylabel(r'$\gamma$')
4 plt.show()
Obrázek 1.7: Duhová funkce primární duhy pro index lomu n = 1.33
Výsledek můžete vidět na obrázku 1.7. Nyní už jen podle pohledu můžete vidět, že nejvíce paprsků se
od kapky odrazí pod úhlem přibližně 40◦ - tam právě dochází ke vzniku primární duhy. Jaká je konkrétní
hodnota tohoto maxima najdeme pomocí numerické derivace v další kapitole. Nyní pojďme do grafu přidat
barvy a funkci pro sekundární duhu (všimněte si, že jsme převod z radiánů do stupňů přidali již přímo do
funkce pro γ).
1 def gama_n(x,r,n):
2 foo = 4.0 * np.arcsin(1/n * x/r)
3 bar = 2. * np.arcsin(x/r)
4 return (foo - bar) * 180/np.pi
5
Kapitola 1 Zákon lomu a Duhová funkce 8
6 def gama_2(x,r,n):
7 foo = 6.0 * np.arcsin(x/(n*r))
8 bar = 2.0 * np.arcsin(x/r)
9 return (np.pi - foo + bar) * 180/np.pi
10
11 indexy = [1.343,1.341,1.337,1.335,1.333,1.331,1.330,1.329]
12 barvy = ['indigo', 'mediumorchid', 'blue', 'lime', 'yellow', 'orange', 'red', 'maroon']
13
14 for i in range(0,len(indexy)):
15 plt.plot(h*1000,gama_n(h,r,indexy[i]),color = barvy[i])
16 plt.plot(h*1000,gama_2(h,r,indexy[i]),color = barvy[i])
17
18 plt.xlabel('h (mm)')
19 plt.ylabel(r'$\gamma$')
20
21 plt.show()
Všimněte si, že jsme pro vykreslení grafu použili smyčky for. Také jsme zvolili barvy tak, aby víceméně
odpovídali barvě duhy – knihovna Matplotlib nabízí celou řadu barev, přehledně je naleznete na jejich
webových stránkách1. Výsledný graf je na obrázku 1.8, včetně výřezu nejzajímavější části kolem maxima
primární a minima sekundární funkce. Nyní už pravděpodobně chápete, jak vzniká Alexandrův tmavý pás
– mezi primární a sekundární duhou je oblast, do které se odráží žádné paprsky – podle zákonů geometrické
optiky by tedy měl být úplně černý. Vidíte také, proč má sekundární duha otočené barvy oproti primární
duze.
Obrázek 1.8: Duhová funkce pro primární a sekundární duhu – celé (vlevo) a přiblížení na oblast maxima
a minima (vpravo).
C Úkol
Jak se změní vlastnosti duhové funkce 1. a 2. řádu pro jinou velikost kapky?
C
1
https://matplotlib.org/3.1.0/gallery/color/named_colors.html
Kapitola 1 Zákon lomu a Duhová funkce 9
1.3.2 Numerické derivace, hledání minima a maxima funkce
I když můžeme použít vestavěné funkce a knihovny, pojďme si funkci, která bude počítat derivaci, implementovat
sami.
 Definice
Derivace funkce f(x) podle x je limita
df(x)
dx
= lim
δx→0
f(x + δx) − f(x)
δx
(1.9)

Derivaci můžeme použít k nalezení extrémů funkce. Jelikož derivace funkce je rovna směrnici tečny
funkce v daném bodě, pro extrém funkce bude platit, že derivace v tomto místě je rovná nule. Toho
využijeme při hledání úhlu, pod kterým se od kapky odrazí nejvíce paprsků. Při pouhém pohledu na graf
na obrázku 1.8 vidíme, že tato hodnota se bude pohybovat okolo 40◦, ale pojďme najít přesnou hodnotu.
Použijeme nejjednodušší metodu, zvolíme h 1 a řekneme, že derivace funkce f(x) v bodě x0 je
d(x0) =
f(x0 + h) − f(x0)
h
. (1.10)
1 def derivative(f,x,h=0.01):
2 return (f(x + h) - f(x))/h
Spočítáme derivace funkce (1.7) podle h a podíváme se na její průběh. Budeme hledat maximum
třeba pro žlutou barvu (n = 1.333). Nadefinujeme si ale duhovou funkci znova, pomocí lambda funkce.
Ta je užitečná, pokud máme velmi krátké funkce, může obsahovat pouze 1 výraz. Jejich výhoda spočívá
v anonymnosti – nemají název, hodí se, pokud chcete rychle vytvořit jednoduchou funkci, krátce jí používat
a poté zase zapomenout.
1 n = 1.333
2 r = 0.25
3 hs = np.linspace(0,r,100)
4 gammafce = lambda x: (4.0 * np.arcsin(1/n * x/r) - (2. * np.arcsin(x/r))) * 180/np.pi
5 gama_der = derivative(gammafce,hs)
6 gamas = gammafce(hs)
7 plt.axhline(0,color ='k')
8 plt.plot(hs,gama_der)
9 plt.show()
Budeme hledat takovou hodnotu γ, kde platí dγ
dh = 0, tedy průnik derivace s osou y. Na obrázku 1.9
je znázorněn průběh derivace duhové funkce a zvýrazněna je osa y. V přiblížení vidíme, že funkce protne
nulu přibližně kolem hodnoty h = 0, 21, avšak získat přesné číslo je komplikovanější – musíme použít další
numerickou metodu, půlení intervalů (bisekce).
1.3.3 Hledání kořene - metoda půlení intervalů
Nejjednodušší metoda pro nalezení kořene rovnice se nazývá metoda půlení intervalů. Spočívá v jednoduchém
předpokladu – ve chvíli, kdy funkce překročí osu y, změní se znaménko její funkční hodnoty. Její
princip je velmi jednoduchý:
Kapitola 1 Zákon lomu a Duhová funkce 10
Obrázek 1.9: Průběh derivace duhové funkce. V obrázku je také přiblížení, které v kódu nenajdete - je jen
ilustrační.
£ Postup
• Vezmeme 2 funkční hodnoty f(x0) a f(x1). Pokud platí, že sgn (f(x0)) = −sgn (f(x1)), víme, že se
kořen nachází někde v intervalu mezi x1 a x0. Pokud to neplatí, musíme interval rozšířit.
• Určíme funkční hodnotu v bodě x2, který leží přesně mezi body x1 a x0, x2 = x0+x1
2 .
• Zjistíme, zda má funkční hodnota f(x2) stejné znaménko, jako f(x0) nebo jako f(x1). Podle toho
určíme nové dva okrajové body tak, aby měly funkční hodnoty opačné znaménko.
• Pokračujeme stejným způsobem, dokud nenajdeme takové xn, kde platí f(xn) = 0 (nebo alespoň
f(xn) = 0 ± , kde 1).
£
Vyzkoušíme si to tedy na naší derivaci duhové funkce. Využijeme opět lambda funkci a vytvoříme
funkci der, která počítá přímo derivaci duhové funkce
1 #Puleni intervalu
2 der = lambda x: derivative(gammafce,x)
3 x0 = hs[0]
4 x1 = hs[-5]
5 xs = np.array([x0,x1])
6 epsilon = 0.001
7 vysledek = False
8 plt.axhline(0,color ='k')
9 plt.plot(hs,der(hs))
10 while (vysledek == False):
Kapitola 1 Zákon lomu a Duhová funkce 11
11 x2 = (x0 + x1)/2.
12 xs = np.array([x0,x1])
13 plt.plot(xs,der(xs),'-or')
14 if((der(x2) < np.abs(0.0+epsilon))):
15 vysledek = x2
16 break
17 plt.plot(x2,0,'-ob')
18 if (np.sign(der(x0)*der(x2)) == 1): x0 = x2
19 elif (np.sign(der(x1)*der(x2)) == 1): x1 = x2
20 plt.show()
21 print(vysledek, gammafce(vysledek))
Obrázek 1.10: Hledání kořene pomocí metody půlení intervalů.
Výsledek, který nám kód vypočítá, je hodnota γ = 41.45◦.
C Úkoly
1. Spočítejte, jaké jsou maxima pro ostatní barvy podle tabulky 1.2.
2. Spočítejte minimum duhové funkce sekundární duhy pro různé barvy. Určete úhlovou šířku Alexandrova
tmavého pásu.
3.  Existuje také duha třetího a čtvrtého řádu, kterou ale pozorujeme na opačné straně, tedy kolem
Slunce. Jsou už opravdu velmi vzácné. Zkuste odvodit, jakým způsobem se v takovém případě paprsky
v kapkách odrážejí a jakou mají duhovou funkci.
C
Kapitola 1 Zákon lomu a Duhová funkce 12
Souhrn kapitoly
V této kapitole jsme nastínili základní principy numerického řešení fyzikálních problémů a základní vykreslování
grafů pomocí knihovny matplotlib v Pythonu. Naučili jsme se jednoduché metody numerické
derivace a hledání kořenů, které budou dále rozvedeny také v následujících kapitolách.
? Otázky
1. Proč je v Alexandrův pás tmavý?
2. Jaké jsou nevýhody metody bisekce?
3. Proč jsou barvy sekundární duhy otočené oproti primární?
4.  Jaké pořadí barev bude v duze 3. a 4. řádu?
?
Další informace
• Encyklopedie fyziky, Duha, J. Reichl
• Optické úkazy v atmosféře na www.ukazy.astro.cz
•  Duha 3. a 4. řádu (anglicky)
Kapitola 2
Modelování dat
Kateřina Klimovičová, Gabriel Török
Rychlý náhled: V kapitole si povíme něco o modelování dat. Především se zaměříme na metodu
nejmenších čtverců.
Klíčová slova: Modelování dat, metoda nejmenších čtverců.
ºº Cíle studia: Po prostudování kapitoly by jste měli být schopni používat metodu nejmenších čtverců.
Mějme dáno N bodů
[x1; y1], [x2; y2], ....[xN ; yN ], (2.1)
kde x1, x2, ...xN jsou pevně dané hodnoty veličiny x. y1, y2, ...yN jsou hodnoty získané měřením veličiny y.
Někdy se setkáváme s označením nezávisle (veličina x) a závisle (veličina y) proměnná.
M Příklad
Veličinou y může být poloha vozidla jedoucího konstantní rychlostí po rovné silnici (označme x ). Veličina
x může být čas měření daného x (označme t). Máme N bodů
[t1; x1], [t2; x2], ..., [tN ; xN ]. (2.2)
Hodnota x1 je poloha vozidla měřená v čase t1 atd.
M
Předpokládejme, že y je známou funkcí x. Neboli, že platí
y = f(x). (2.3)
M Příklad
V našem příkladu máme
x = f(t) = v0t + x0, (2.4)
kde v0 a x0 jsou konstanty (rychlost a poloha vozidla v čase t = 0).
M
13
Kapitola 2 Modelování dat 14
Pokud jsme schopni veličinu y měřit přesně, platí
y1 = f(x1), y2 = f(x2), ..., yN = f(xN ). (2.5)
Jinými slovy, všechny body v rovině x − y leží na křivce dané vztahem y = f(x).
Běžně se však setkáváme se situací, kdy veličinu y nejsme schopni měřit přesně.
Funkce f(x) obecně závisí na P parametrech. Předpokládejme, že počet parametrů P je menší než
počet bodů N. Označme parametry jako a1, a2, ..., aP . Můžeme napsat y(x) = f(x; a1, a2, ..., aP ). Platí
yi = f(xi; a1, a2, ..., aP ) + ei. (2.6)
Veličiny e1, e2, ..eN jsou náhodné (chyby měření).
M Příklad
Vraťme se k našemu příkladu. Pokud polohu vozidla měříme přesně, graf závislosti polohy x na čase t může
vypadat jako ten umístěný na levém panelu obr. 2.1. Pokud polohu nejsme schopni naměřit zcela přesně,
graf bude vypadat jako ten na pravém panelu stejného obrázku. Pro jednotlivé xi platí
xi = f(ti; v0, x0) + ei, (2.7)
kde
f(ti; v0, x0) = tiv0 + x0. (2.8)
0
50
100
150
200
250
300
0 2 4 6 8 10 12 14 16
t [s]
0
50
100
150
200
250
300
0 2 4 6 8 10 12 14 16
x’[m]
t [s]
x’[m]
Obrázek 2.1: Poloha vozidla x jako funkce času měření t. Levý panel: Poloha vozidla měřená přesně.
Černé body odpovídají jednotlivým měřením. Červená křivka je dána rovnicí (2.7). Parametry jsme zvolili:
v0 = 15 m.s−1 a x0 = 50 m. Pravý panel: obdoba levého panelu pro případ kdy rychlost není měřena přesně.
M
Pokud bychom měli veličinu y změřenou přesně, pak proložit ji křivkou y = f(x; a1, a2, ...aP ) není
problém. Pokud nejsme schopni veličinu y měřit přesně pak musíme nejdříve odhadnout hodnoty zatím
neznámých parametrů a1, a2, ..., aP . Přitom požadujeme aby křivka co nejlépe přiléhala k bodům
[x1; y1], [x2; y2], ...., [xN ; yN ]. Jelikož nalezené hodnoty nebudou přesnými hodnotami parametrů, ale pouze
jejich odhady, je třeba si zavést nové označení. Například je můžeme značit velkým písmenem (A1, A2, ...AP ).
Dále je třeba zvolit kritérium určující jak moc křivka přiléhá k bodům.
Kapitola 2 Modelování dat 15
2.1 Metoda nejmenších čtverců
V případě metody nejmenších čtverců je kritériem přiléhavosti křivky k bodům minimálnost veličiny
S =
N
i=1
[yi − f(xi; A1, A2, ...Ap)]2
. (2.9)
Jedná se o součet kvadrátů (čtverců) rozdílu mezi yi a f(xi; A1, A2, ...Ap).
Odhady A1, A2, ..., Ap získáme z podmínky
N
i=1
[yi − f(xi; A1, A2, ...Ap)]2
= min (S) (2.10)
M Příklad
V našem příkladu jedoucího vozidla je
x = v0t + x0 (2.11)
Neboli
S =
N
i=1
[xi − v0t − x0]2
. (2.12)
Odhady na V0 a X0 získáme z podmínky
N
i=1
[xi − V0t − X0]2
= min (S). (2.13)
Nyní aplikujeme výše uvedený postup na body (viz pravý panel obr. 2.1):
[1; 66, 52], [3; 77, 14], [6; 128, 01], [8; 174, 53], [12; 230, 79], [15; 287, 52]. (2.14)
Jinými slovy najdeme takové V0 a X0, pro které veličina S nabývá minima. Tímto postupem dostaneme
V0 = 16, 3297 m.s−1 a X0 = 38, 27890 m . K výpočtu jsme použili program jehož přepis je uveden v kapitole
2.2. Odpovídající křivku je možné vidět na obrázku 2.2.
t [s]
x’[m]
0
50
100
150
200
250
300
0 2 4 6 8 10 12 14 16
Obrázek 2.2: Poloha vozidla x jako funkce času měření t proložená křivkou nejlépe přiléhající k bodům.
Křivka je dána rovnicí (2.7) s parametry V0 = 16, 3297 m.s−1 a X0 = 38, 27890 m.
M
Kapitola 2 Modelování dat 16
2.1.1 Lineární regresní funkce
Nyní se zabývejme případem, kdy y = f(x; a1, a2, ...aP ) je tvaru
y = a1φ1(x) + a2φ2(x) + ... + aP φP (x). (2.15)
φ1(x), φ2(x), ...φP (x) jsou funkce x, avšak neobsahují žádné další parametry. Jinými slovy, y „závisí“ na
parametrech a1, a2, ..., aP lineárně.
Pro yi obecně platí
yi = a1φ1(xi) + a2φ2(xi) + ... + apφp(xi) + ei, (2.16)
kde ei jsou náhodné veličiny. Odhad parametrů A1, A2, ..., AP určíme z podmínky
S =
N
i=1
[yi − a1φ1(xi) − a2φ2(xi) − ... − apφp(xi)]2
= min S. (2.17)
K výpočtu můžeme použít rovnic
∂S
∂A1
= 0,
∂S
∂A2
= 0, ... ,
∂S
∂AP
= 0.
Neboli
∂S
∂A1
= (−2)
N
i=1
[yi − A1φ1(xi) − A2φ2(xi) − ... − AP φP (xi)]φ1(xi) = 0,
∂S
∂A2
= (−2)
N
i=1
[yi − A1φ1(xi) − A2φ2(xi) − ... − AP φP (xi)]φ2(xi) = 0,
...
∂S
∂AP
= (−2)
N
i=1
[yi − A1φ1(xi) − A2φ2(xi) − ... − AP φP (xi)]φP (xi) = 0.
Zavedeme substituce
Bjh =
N
i=1
φh(xi)φj(xi), Bh =
N
i=1
yiφh(xi).
Dostáváme soustavu p lineárních rovnic o p neznámých:
B11A1 + B21A2 + ... + BP1AP = 0,
B12A1 + B22A2 + ... + BP2AP = 0,
...
B1P A1 + B2P A2 + ... + BPP AP = 0.
Řešením jsou odhady parametrů A1, A2, ..., AP .
Kapitola 2 Modelování dat 17
M Příklad
V příkladu, který nás provází celou kapitolou, je funkce
x = v0t + x0 (2.18)
skutečně lineární ve výše uvedeném smyslu.
Již výše bylo uvedeno, že
S =
6
i=1
xi − v0ti − x0
2
. (2.19)
Dostáváme
∂S
∂v0
= (−2)
6
i=1
xiti − v0
6
i=1
titi − x0
6
i=1
ti = 0, (2.20)
∂S
∂x0
= (−2)
6
i=1
xi − v0
6
i=1
ti − x0
6
i=1
= 0. (2.21)
(2.22)
Po úpravě máme dvojici rovnic
6
i=1
xiti − v0
6
i=1
titi − x0
6
i=1
ti = 0, (2.23)
6
i=1
xi − v0
6
i=1
ti − x0
6
i=1
= 0. (2.24)
Zaveďme substituce
B1 =
6
i=1
tixi, B2 =
6
i=1
xi, B11 =
6
i=1
titi, B12 =
6
i=1
ti, B22 =
6
i=1
= 6. (2.25)
Rovnice se zredukují do tvaru
B1 − v0B11 − x0B12 = 0, (2.26)
B2 − v0B12 − x0B22 = 0. (2.27)
Řešením rovnic dostáváme
x0 = 38.27890 m, v0 = 16.3297m.s−1
. (2.28)
K výpočtu byl použit program uvedený v kapitole 2.2.
M
2.1.2 Příklad z reálného světa - Určení stáří vesmíru
V roce 1929 Edwin Hubble uvedl (viz [3]), že čím je galaxie vzdálenější od Země, tím se od ní rychleji
vzdaluje1. Zjištění lze shrnout v podobně tzv. Hubbleova (či Hubbleova-Lemaîterova) zákona
v = Hr. (2.29)
1
Ponechme nyní stranou spory o skutečném prvenství této práce.
Kapitola 2 Modelování dat 18
Hubbleův zákon dává do souvislosti rychlost vzdalování galaxie v a vzdálenost galaxie r od země. Konstanta
úměrnosti, H, se nazývá Hubbleova konstanta. Může být interpretována jako relativní rychlost rozpínání
vesmíru.
Z hodnoty Hubbleovy konstanty je teoreticky možné určit stáří vesmíru. Potřebujeme však vědět nejen
současnou hodnotu H, ale i vývoj s časem. Předpokládejme nyní, že hodnota H se s časem nemění. Jinými
slovy, že vesmír se už od svého počátku rozpíná stejně rychle.
C Úkol
Na konci této podkapitoly jsou uvedeny rychlosti a vzdálenosti galaxií od Země. Vaším úkolem je určit
z těchto pozorování hodnotu Hubbleovy konstanty a z ní předpokládané stáří vesmíru. Pokud si s tím
nebudete vědět rady níže uvádíme podrobnější návod. Závislost rychlosti vzdalování na vzdálenosti od
Země je vykreslena na obr. 2.3.
C
0
5
10
15
20
25
30
35
40
0 100 200 300 400 500 600 700
r [Mpc]
r[km.s]x10
-1-3
Obrázek 2.3: Hubbleův diagram. Závislost rychlosti vzdalování galaxií v na jejich vzdálenosti od Země r.
Hondoty jsme převzali z [5].
Kapitola 2 Modelování dat 19
Hodnoty rychlostí v galaxií meřených v km.s−1
a vzdálenosti galaxií r v megaparsecích (založeno
na [5]):
ri = [35.101, 40.501, 42.864, 55.689, 58.39, 66.152, 68.515, 72.903, 76.615,
84.378, 84.716, 93.153, 98.891, 103.954, 111.041, 115.429, 117.117,
118.129, 127.917, 128.592, 134.667, 138.717, 145.805, 148.168, 150.193,
155.593, 163.693, 167.743, 169.769, 186.982, 188.332, 190.357, 202.845,
214.658, 227.821, 231.533, 235.583, 237.271, 241.996, 259.547, 280.473,
292.96, 303.761, 341.225, 347.975, 385.439, 420.878, 499.18, 635.535]
vi = [2535.6, 2704.8, 3380.5, 3429.7, 4057.1, 4130, 5408.8, 4781.9, 5940.1,
4927.4, 5385.8, 5917.2, 7389.2, 6448.7, 7824.3, 6956.1, 7631.7, 8307.3,
8549.2, 9755.5, 7801.8, 8718.9, 8212.7, 9443.3, 10722.1, 11542.7,
8841.2, 10723.3, 10361.5, 12847.6, 13740.4, 10893.7, 11956.1, 14996.7,
14635.7, 14949.6, 13550.6, 16132.2, 17387, 15602.9, 18258.2, 16690.9,
19055.9, 23666.5, 22412.4, 21160.5, 26374, 30215.3, 37293.5]
Poznámka
Rada: Předpokládejte, že rychlost vzdalování známe přesně (i když tomu tak není) a vzdálenost je měřena
s chybou. Jak spočítat stáří vesmíru? Při velkém třesku by vzdálenost mezi všemi objekty ve vesmíru měla
jít do nuly. Tedy můžeme říci, že pro vzdálenost galaxií v čase t platí
r = vt + r0, kde r0 = 0. (2.30)
Z výše uvedeného je patrné, že čas, za který se galaxie dostanou do současného stavu, vypočteme z rovnice
t =
r
v
=
1
H
. (2.31)
Předpokládané stáří vesmíru je tedy nepřímo úměrné hodnotě Hubbleovy konstanty. Dávejte si však pozor
na jednotky.
Řešení: Mě vyšlo stáří vesmíru přibližně 15.7 mld let .
2.2 Příklady v Pythonu
Níže uvádíme několik příkladů, jak lze řešit výše uvedené úkoly.
1. Hledání parametrů pro které křivka nejlépe přiléhá k bodům. Viz sekce 1.1.
1 from scipy.optimize import minimize
2
3 def square(a,b, data):
4 x, y = data
5 return ((y-(a+b*x))*(y-(a+b*x))).sum()
6
Kapitola 2 Modelování dat 20
7 d = numpy.array([1, 3, 6, 8, 12, 15])
8 v = numpy.array([66.5248, 77.1354, 128.014, 174.528, 230.793, 287.515])
9
10 data = (d,v)
11
12 def squarekon(p):
13 x,y=p
14 return square(x,y,data)
15
16 minimize(squarejedna,(10.0,10.0))
2. Obdobně jako 1., avšak tentokrát je výpočet proveden metodou popsanou v sekci 1.2.
1 import numpy as np
2 from scipy.optimize import minimize
3
4 d = numpy.array([1, 3, 6, 8, 12, 15])
5 v = numpy.array([66.5248, 77.1354, 128.014, 174.528, 230.793, 287.515])
6
7 data = (d,v)
8
9 A = np.array([[-(d*d).sum(), -(d).sum()], [-(d).sum(), -6]])
10 B = np.array([-(d*v).sum(), -(v).sum()])
11 X = np.linalg.solve(A,B)
12
13 print(X)
3. Příklad kódu o kterém jsme dosud nehovořili, ale který může sloužit ke kontrole výše uvedených
výsledků. Zde jsme použili knihovnu Kapteyn
1 import numpy
2 from kapteyn import kmpfit
3
4 def residuals(p, data):
5 x, y = data
6 a, b = p
7 return (y-(a+b*x))
8
9 d = numpy.array([1, 3, 6, 8, 12, 15])
10 v = numpy.array([66.5248, 77.1354, 128.014, 174.528, 230.793, 287.515])
11
12 paramsinitial = [1.0, 1.0]
13 fitobj = kmpfit.Fitter(residuals=residuals, data=(d,v))
14 fitobj.fit(params0=paramsinitial)
15 print('Parametry: ', fitobj.params)
Kapitola 2 Modelování dat 21
Další informace
• Kapteyn Astronomical Institute. Least squares fitting with kmpfit. Dostupné z
https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html
• PRESS William H., TEUKOLSKY Saul A., VETTERLING William T., FLANNERY Brian P. Numerical
recipes. The art of scientific computing. Third edition. Cambridge. 2007, ISBN-13 978-0-521-
88068-8
• REKTORYS, Karel a kolektiv. Přehled užité matematicy. SNTL. 1980, ISBN 04-003-81.
Souhrn kapitoly
V kapitole jsme se věnovali především využití metody nejmenších čtverců.
? Otázky
1. Vysvětlete důvody, proč jsme na počátku předpokládali, že počet bodů N je větší než počet
parametrů P.
?
Kapitola 3
Interpolace a extrapolace
Gabriela Urbancová
Rychlý náhled: V této kapitole se naučíme, co je interpolace a extrapolace, jak je využívat v praxi,
tedy jak spojit zadanou množinu bodů křivkou tak, aby tato křivka byla hladká, procházela všemi zadanými
body a co nejvíce se blížila původní, námi hledané, funkci. Uvedeme si zde příklady různých interpolací,
k jejichž výpočtu a k tvorbě grafů využijeme programovací jazyk Python a jeho vestavěné knihovny a
funkce. V těchto krátkých programech se naučíte porovnávat různé interpolační metody, jak numericky,
tak graficky. Taky si ukážeme i jiný způsob, jak ověřit kvalitu zvolené interpolace, konkrétně pomocí
Riemannova určitého integrálu.
Klíčová slova: Interpolace, extrapolace, interpolační uzel, interpolace polynomem, stupeň polynomu,
koeficient polynomu, lineární interpolace, Lagrangeova metoda, lineární splajn, kubický splajn, Python,
Scipy, quad, Riemannův určitý integrál, Rungeova funkce.
ºº Cíle studia: Po prostudování této kapitoly budete vědět, k čemu a jak se používá interpolace a
extrapolace. Budete znát několik druhů interpolací, v jakých případech se dají použít a kdy je lepší zvolit
jinou aproximační metodu. Také budete umět zvolit vhodné knihovny a funkce v programu Python, spočítat
a vykreslit jednotlivé interpolace a srovnat jejich výsledky. Nakonec si spočítáte integrály známých funkcí
a srovnáte přesné výsledky s jejich interpolovanými hodnotami.
3.1 Co je to interpolace a extrapolace
Jednou z aproximačních metod, která se hojně využívá v numerické matematice, je interpolace, resp. extrapolace.
Hlavním znakem interpolace je, že nově sestrojená aproximační funkce g(x) musí přesně procházet
všemi známými body. Tato množina bodů, tzv. interpolačních uzlů, je výsledkem původní funkce f(x). Ve
většině případů je to tak, že tato původní funkce f(x) je neznámá a není spočítána ve všech bodech. Mějme
tedy jednoduchou tabulku obsahující napočítané hodnoty yi náležící jednorozměrné funkci f(xi) v bodech
x0, x1, ..., xN , xi = xj pro i = j, (3.1)
ale my potřebujeme zjistit hodnotu y v jiném bodě x, a to tak, aby x0 < x < xN . Interpolací se tedy
snažíme nalézt hodnoty, které se nachází mezi dvěma interpolačními uzly uvnitř známého intervalu.
22
Kapitola 3 Interpolace a extrapolace 23
Obrázek 3.1: Extrapolace funkce f(x) = sin(x) na 20ti interpolačních uzlech.
Někdy se stane, že nám některé okrajové body v požadovaném intervalu chybí a potřebujeme je znát.
Takovému procesu se říká extrapolace a snažíme se nalézt přibližné hodnoty funkce mimo interval známých
hodnot. Funguje na stejném principu jako interpolace, tedy opět hledáme vhodnou náhradní funkci. Díky
této funkci, která odpovídá okrajovým hodnotám nám neznámé funkce, můžeme určit i hodnoty mimo
původní rozsah měření, ale pouze s jistou pravděpodobností. Čím dál od okrajových bodů extrapolujeme,
tím více nepřesné máme výsledky (roste odchylka od skutečných hodnot). Extrapolace je metoda méně
využívaná, protože není tak přesná a funguje spíše jako hrubý odhad toho, jak by se mohla funkce chovat,
když je za okrajovými uzly. Ukázku extrapolované funkce f(x) = sin(x) můžete vidět na obr. 3.1, použili
jsme lineárního a Lagrangeův interpolační polynom. Na námi uvedeném příkladě se odhad minulého ani
budoucího chování pomocí extrapolace moc nezdařil, i když jsem zvolili poměrně hodně interpolačních uzlů.
 Definice
Interpolace je tedy takový typ aproximace, u které požadujeme, aby hledaná funkce přesně procházela
zadanými body neboli hledáme takovou funkci g(x), která se shoduje s funkcí f(x) v uzlových bodech xi
g(xi) = f(xi) = yi i = 1, . . . , N. (3.2)

Existuje celá řada různých typů interpolací [6], u některých požadujeme i shodnost derivací v uzlových
bodech, tj. g (xi) = f (xi), u jiných i spojitost druhých derivací (např. Hermiteova interpolace). Každá z
existujících interpolací má své výhody i omezení.
V tomto textu si uvedeme některé základní interpolace [7], omezíme se na případ jedno dimenzionální,
kdy všechny body leží v rovině v kartézském souřadném systému a jsou seřazeny vzestupně podle první
souřadnice. Samozřejmě se dají interpolovat i funkce více proměnných, ale pro zachování jednoduchosti
Kapitola 3 Interpolace a extrapolace 24
této kapitoly se jimi zabývat nebudeme.
Na konci této kapitoly najdete části kódu psaného v Pythonu, pomocí kterého si můžete vyzkoušet
druhy interpolací, které zde uvedeme.
3.2 Polynomiální interpolace
Při interpolování se snažíme najít takovou třídu funkcí, která bude v nějakém smyslu výhodná. Nejvíce
využívané funkce v běžné praxi jsou ty ve tvaru polynomu
P(x) = A0 + A1x + A2x2
+ . . . + AN xN
=
N
i=1
Ai · xN−i
. (3.3)
Požadujeme, aby byla splněna interpolační podmínka (3.2)
PN (xi) = yi i = 1, . . . , N (3.4)
a dostáváme tedy předpis pro interpolační polynom nejvýše stupně N.
Jak můžeme vidět, tak matematický tvar polynomu vypadá jednoduše. Stačí, když určíme stupeň polynomu,
napočítáme příslušné koeficienty pomocí řešení soustavy lineárních rovnic a dostaneme jednoznačný
tvar hledaného interpolačního polynomu. Stupeň polynomu je definován jako nejvyšší mocnitel proměnné v
nenulovém členu. Výhodou tohoto tvaru funkce je také následné počítání derivace a integrace. Volbou různého
stupně polynomu jsme schopni dosáhnout libovolné přesnosti, ale pozor, u vysokých stupňů polynomu
už dochází k velkých zaokrouhlovacím chybám a k nežádoucím oscilacím.
3.2.1 Lineární interpolace
Polynomem prvního stupně (pro N = 1) interpolujeme lineárně. Dva sousední uzly se jednoduše proloží
úsečkou. Hodnotu interpolační funkce v bodě x spočítáme tak, že nejdříve zjistíme, mezi kterými uzly
námi hledaný bod leží. Je potřeba vědět, zda mám uzly rozmístěny rovnoměrně či nerovnoměrně. Na tento
výpočet existují různé algoritmy. Když nakonec znám polohu xi a xi+1, námi hledaná interpolační funkce
se spočítá ze vztahu
g(x) − f(xi)
x − xi
=
f(xi+1) − f(xi)
xi+1 − xi
⇒ g(x) = f(xi) +
f(xi+1) − f(xi)
xi+1 − xi
(x − xi). (3.5)
Takto vzniklá interpolační funkce je sice spojitá, ale není hladká, tedy její první derivace není spojitá, proto
je její použití značně omezené. Na obr. 3.2 vidíme srovnání lineární a polynomiální interpolace..
3.2.2 Lagrangeův interpolační polynom
Při výpočtu jednotlivých koeficientů polynomu může docházet ke vzniku velkých zaokrouhlovacích chyb,
hlavně v případech vysokého stupně polynomu. Lepším řešením je se tomuto výpočtu úplně vyhnout a
použít tzv. Lagrangeovu metodu. Ta spočívá ve složení jednotlivých elementárních polynomů
li(x) =
(x − x0) . . . (x − xi−1)(x − xi+1) . . . (x − xN )
(xi − x0) . . . (xi − xi−1)(xi − xi+1) . . . (xi − xN )
(3.6)
splňujících podmínku
li(xj) = δij, (3.7)
Kapitola 3 Interpolace a extrapolace 25
Obrázek 3.2: Srovnání lineární a polynomiální interpolace pro případ 4 interpolačních uzlů.
do výsledného Lagrangeova interpolačního polynomu
LN (xk) =
N
i=0
yili(xk) = yk, pro i, k = 0, . . . , N. (3.8)
Pokud máme tabulku obsahující vysoký počet interpolačních uzlů, kterými výsledný polynom musí
projít, tak Lagrangeova metoda není vhodná, její výpočet se stává složitý a zdlouhavý. Naopak je metoda
dobře využitelná, pokud interpolujeme různé funkce, ale ve stejných uzlových bodech.
3.3 Interpolace pomocí splajnu
Jak už jsme uvedli dříve, v případech vysokého počtu interpolačních uzlů musíme k interpolaci zvolit
polynom vysokého stupně. Výsledná funkce se sice v blízkosti uzlů chová celkem hezky, ale mezi uzly
dochází ke kmitům, viz obr. 3.3, které jsou samozřejmě nežádoucí, ale jsou přímým důsledkem vysoké
zaokrouhlené chyby ve výpočtech. Zde přichází na řadu metoda, které se česky říká splajn (z angl. spline).
Opět využíváme funkce ve tvaru polynomu, ale tentokrát tak, že hledaná interpolační funkce je po celém
intervalu rozdělená na malé části. Jednotlivé části mezi uzly jsou popsány polynomy nízkého stupně a v
uzlech na sebe hladce navazují.
Pro představu si zde uvedeme pouze dva druhy splajnu, konkrétně Lineární a Kubický splajn. Existuje
jich samozřejmě více a mají různé výhody i nevýhody, další informace najdete např. v doporučené literatuře.
3.3.1 Lineární interpolační splajn
Nejjednoduším druhem interpolace splajnem je lineární a neznamená nic jiného, než proložení dvou sousedních
uzlů úsečkou. Mějme tedy dva sousední uzly dané souřadnicemi [xi−1, yi−1] a [xi, yi], které proložíme
Kapitola 3 Interpolace a extrapolace 26
Obrázek 3.3: Interpolace Rungeovy funkce pomocí Lagrangeova polynomu 10. stupně.
úsečkou, tedy lineárním interpolačním polynomem ve tvaru [7]
Si(x) = yi−1 +
(yi − yi−1)
(xi − xi−1)
(x − xi−1) = yi−1 + sδi, (3.9)
kde S(x) je spojitá funkce, ale její první derivace ve vnitřních uzlech obecně spojitá není. Pokud známe
dostatečné množství uzlů, pak lze chybu interpolace ovlivňovat a zmenšovat.
C Úkol
Vysvětlete vlastními slovy, čím se od sebe liší lineární interpolace a lineární splajn.
C
3.3.2 Kubický splajn
Velmi používanou aproximací je interpolace pomocí kubického splajnu. Využíváme tak vlastností polynomu
3. stupně, který je spojitý a má také spojité první i druhé derivace a zároveň má nízké riziko vzniku
nežádoucích kmitů. Obr. 3.4 potvrzuje, že jsme pomocí kubického splajnu dostali naprosto vyhovující
aproximaci Rungeovy funkce i při zachování počtu interpolačních uzlů.
Podrobnou definici zde pro jednoduchost uvádět nebudeme, ale krásně vysvětlený kubický splajn
i s postupem řešení můžete najít na http://kfe.fjfi.cvut.cz/ psikal/nme/cv04/kubspline.pdf.
C Úkol
Mějte zadanou funkcí pomocí 12ti interpolačních uzlů. Která z intepolačních metod bude
nejvhodnější?
C
Kapitola 3 Interpolace a extrapolace 27
Obrázek 3.4: Interpolace Rungeovy funkce pomocí kubického splajnu.
3.4 Interpolace v programu Python
V programovacím jazyce Python existuje velká knihovna Scipy, která obsahuje velké množství různých
druhů interpolací, jejich výčet lze najít na SciPy.org. Pro naše účely si výběrem jen některé druhy interpolací
a jejich použití si ukážeme v následujícím příkladě.
M Příklad
Nechť máme funkci f(x0) = sin(x0) a známe její funkční hodnotu v 5ti interpolačních uzlech, jejichž
souřadnice jsou x0 = [0.0, 0.785, 1.571, 2.356, 3.142] na intervalu < 0, π >.
• Srovnejte, jakou výslednou hodnotu spočítají různé druhy interpolací dostupné v knihovně Scipy
(konkrétně použijeme interpolaci lineárním polynomem, Lagrangeovu interpolaci, lineární splajn a
také splajny 3. a 4. stupně), a to konkrétně v bodech x = [0.2, 0.5, 0.9, 2.4] s hodnotami vypočítanými
originální funkcí.
• Na základě získaných výsledků rozhodněte, která z interpolací je pro náš příklad nejvhodnější.
K výpočtu použijte níže uvedený kód v Pythonu.
Uvedeme zde řešení pouze pro x = 0.2 (ostatní jsou za úkol):
1 # Importujeme knihovny, které chceme použít
2 import numpy as np
3 from scipy import interpolate as inter
4 import math
5
6 # Poloha 5ti interpolačních uzlů a originální funkce, ze které dopočítáme hodnoty
7 x0 = np.array([0.0 , 0.785, 1.571, 2.356, 3.142])
Kapitola 3 Interpolace a extrapolace 28
8 y0 = [math.sin(i) for i in x0]
9
10 # Poloha bodů, ve kterých chceme znát interpolační hodnoty
11 x = np.array([0.2, 0.5, 0.9, 2.4])
12
13 # Naše interpolační funkce
14 funkce_linear = inter.interp1d(x0, y0, kind='linear') # lineární interpolace
15 funkce_lagrange = inter.lagrange(x0, y0) # Lagrangeův polynom
16 funkce_spline = inter.InterpolatedUnivariateSpline(x0, y0, k=1) # splajn 1. stupně
17 funkce_cubic_spline = inter.InterpolatedUnivariateSpline(x0, y0, k=3) # splajn 3. stupně
18 funkce_spline_4 = inter.InterpolatedUnivariateSpline(x0, y0, k=4) # splajn 4. stupně
19
20 print("Interpolovaná hodnota 'linear' v bodě", x[0], '=', funkce_linear(x[0]))
21 print("Interpolovaná hodnota 'lagrange' v bodě", x[0], '=', funkce_lagrange(x[0]))
22 print("Interpolovaná hodnota 'spline' v bodě", x[0], '=', funkce_spline(x[0]))
23 print("Interpolovaná hodnota 'cubic spline' v bodě", x[0], '=', funkce_cubic_spline(x[0]))
24 print("Interpolovaná hodnota 'cubic spline 4' v bodě", x[0], '=', funkce_spline_4(x[0]))
25 print("Originální hodnota v bodě", 0.2, '=', math.sin(0.2))
Program na obrazovku vypíše:
Interpolovaná hodnota ’linear’ v bodě 0.2 = 0.18008284868926522
Interpolovaná hodnota ’lagrange’ v bodě 0.2 = 0.19689743442955368
Interpolovaná hodnota ’spline’ v bodě 0.2 = 0.18008284868926525
Interpolovaná hodnota ’cubic spline’ v bodě 0.2 = 0.20516280705344686
Interpolovaná hodnota ’cubic spline 4’ v bodě 0.2 = 0.19689743442955354
Originální hodnota v bodě 0.2 = 0.19866933079506122
M
Zkuste si najít nápovědu k funkci interp1d z knihovny scipy.interpolate, tato funkce je totiž
hodně specifická a má velmi široký záběr. V našem příkladě ji používáme pouze k výpočtu lineární interpolace,
ale v nápovědě najdete, že existují i volby ‘nearest’, ‘nearest-up’, ‘zero’, ‘slinear’,
‘quadratic’, ‘cubic’, ‘previous’ a ‘next’ a každá z nich počítá interpolaci odlišným způsobem.
M Příklad
Nechť stále máme funkci f(x0) = sin(x0). Na intervalu < 0, π > si pomocí příkazu v Pythonu
1 x0 = np.linspace(0, math.pi, num=10, endpoint=True)
napočítáme různý počet interpolačních uzlů (budeme měnit číslo num) v ekvidistantních vzdálenostech od
sebe. Pomocí podobného příkazu
1 x = np.linspace(0, math.pi, num=100, endpoint=True)
si navolíme počet nových bodů, ve kterých chceme znát hodnotu interpolace. Vysledujte, jakým způsobem
se mění předchozí výsledky, pokud zvyšuji nebo snižuji počet interpolačních uzlů. Diskutujte výsledky a
použitelnost jednotlivých druhů interpolací.
M
Kapitola 3 Interpolace a extrapolace 29
C Úkol
Existuje i jiný způsob, jak ověřit na kolik vybraná interpolace odpovídá původní funkci?
C
M Příklad
Mějme Rungeovu funkci g(x0) = 1/(1 + 25x2
0). Na intervalu < −1, 1 > si pomocí příkazu v Pythonu
spočítejte její integrační hodnotu. Poté funkci interpolujte všemi již uvedenými způsoby a spočítejte jejich
integrační hodnoty. Výsledné hodnoty porovnejte a rozhodněte, která z interpolací je nejvhodnější.
1 # importované části knihoven
2 from math import sin, pi
3 from scipy.integrate import quad
4
5 # definice funkce, kterou chceme numericky integrovat
6 def g(x):
7 return 1./(1.+25.*x**2.)
8
9 # použijeme vestavěnou funkci quad pro výpočet určitého integrálu od 0 do Pi
10 vysledek, chyba = quad(g, -1., 1.)
11 print("Numerický výsledek funkce je {:f} (+-{:g})" .format(vysledek, chyba))
Na obrazovku se nám vypíše:
Numerický výsledek funkce je 0.549360 (+-2.86683e-09).
Nyní si spočítáme, jaké integrační hodnoty dostaneme použitím různých interpolací, jinými slovy spočítáme
obsahy ploch ohraničené interpolovanou křivkou. Ze znalostí o počítání určitých integrálů zvolíme
Riemannův integrál. Na obr. 3.5 vidíme grafické srovnání levé, středové a pravé Riemannovy sumy spočítané
pro případ Rungeovy funkce. Zvolili jsme dělení intervalu na 11 oddílů. Hned na první pohled vidíme
jasné rozdíly ve způsobu výpočtu obsahů ploch pod křivkou.
Obrázek 3.5: Graficky znázorněný postup při počítání Riemannova integrálu Rungeovy funkce.
Lze samozřejmě použít i jiný druh integrálu, např. Simpsonovo pravidlo (více informací v knihovně
scipy.integrate.simpson), které k výpočtu používá kvadratický polynom na každém oddílu křivky. Také
Kapitola 3 Interpolace a extrapolace 30
existuje lichoběžníkového pravidlo (scipy.integrate.trapezoid), které aproximuje každý oddíl křivky
přímkou.
Z výše uvedeného kódu víme, jaká je přesná hodnota integrálu Rungeovy funkce na intervalu < −1, 1 >.
Použitím Riemannovy sumy dostáváme hodnoty (pro případ dělení na 11 oddílů):
Levá Riemannova suma: 0.5477044454139318,
Středová Riemannova suma: 0.4956374418136541,
Pravá Riemannova suma: 0.33904407478069454.
Pro znázornění výpočtu jsme použili interpolační funkci nazvanou funkce_spline, která nám interpoluje
Rungeovu funkci pomocí splajnu 1. stupně. Výpočty ostatních druhů interpolací jsou za úkol.
1 # Importujeme knihovny, které chceme použít
2 import numpy as np
3 import matplotlib.pyplot as plt
4 from scipy import interpolate as inter
5
6 # Poloha 11ti interpolačních uzlů a originální funkce, ze které dopočítáme hodnoty
7 x0 = np.linspace(-1, 1, num=11, endpoint=True)
8 y0 = [1/(1.+(25.*(i)**2)) for i in x0]
9
10 # Poloha bodů, ve kterých chceme znát interpolační hodnoty
11 x = np.linspace(-1, 1, num=100, endpoint=True)
12
13 # Naše interpolační funkce
14 funkce_linear = inter.interp1d(x0, y0, kind='linear', bounds_error=False, fill_value='extrapolate')
15 funkce_lagrange = inter.lagrange(x0, y0) # Lagrangeův polynom
16 funkce_spline = inter.InterpolatedUnivariateSpline(x0, y0, k=1) # splajn 1. stupně
17 funkce_cubic_spline = inter.InterpolatedUnivariateSpline(x0, y0, k=3) # splajn 3. stupně
18 funkce_spline_5 = inter.InterpolatedUnivariateSpline(x0, y0, k=5) # splajn 5. stupně
19
20 # Do grafu vykreslíme modrou gridovací mřížku
21 plt.plot(x0, y0, 'o', label='Data')
22 plt.grid(linestyle="-", color=(0.7, 0.8, 1.0))
23
24 # Vykreslíme si všechny použité interpolace do jednoho grafu
25 plt.plot(x, funkce_linear(x), '-', label='lineární polynom')
26 plt.plot(x, funkce_lagrange(x), '-', label='Lagrangeův polynom')
27 plt.plot(x, funkce_spline(x), '-', label='splajn 1. stupně')
28 plt.plot(x, funkce_cubic_spline(x), '-', label='kubický splajn')
29 plt.plot(x, funkce_spline_5(x), '-', label='splajn 5. stupně')
30
31 # Označíme osy v grafu
32 plt.xlabel('x')
33 plt.ylabel('f(x)')
34 # Legendu v grafu
35 plt.legend()
36
37 plt.show()
38
39 # Předpisy pro výpočet Riemannovy sumy ve 3 různých polohách
40 a = -1
Kapitola 3 Interpolace a extrapolace 31
41 b = 1
42 N = 11
43 dx = (b-a)/N
44 x_vlevo = np.linspace(a, b - dx, N)
45 x_uprostred = np.linspace(dx/2,b - dx/2,N)
46 x_vpravo = np.linspace(dx,b,N)
47
48 # Výpočet Riemannova integrálů pro různé polohy
49 Riemannova_suma_vlevo = np.sum(funkce_spline(x_vlevo) * dx)
50 Riemannova_suma_uprostred = np.sum(funkce_spline(x_uprostred) * dx)
51 Riemannova_suma_vpravo = np.sum(funkce_spline(x_vpravo) * dx)
52
53 # Tisk výsledků
54 print("Levý Riemannův integrál:", Riemannova_suma_vlevo)
55 print("Středový Riemannův integrál:", Riemannova_suma_uprostred)
56 print("Pravý Riemannův integrál:", Riemannova_suma_vpravo)
Obrázek 3.6: Různé druhy interpolací aplikovaných na Rungeovu funkci.
Spuštěním výše uvedeného kódu obdržíme obr. 3.6 a text:
Levý Riemannův integrál: 0.5473018959649976
Středový Riemannův integrál: 0.5027747653416103
Pravý Riemannův integrál: 0.35668449197860963
Porovnáme nyní získané hodnoty integrálů Rungeovy funkce:
quad funkce: 0.549360(+ − 2.86683e − 09)
Riemannova levá suma aplikovaná na originální funkci: 0.5477044454139318
Levý Riemannův integrál z funkce_spline : 0.5473018959649976
Kapitola 3 Interpolace a extrapolace 32
Všechny získané výsledky takto porovnejte a diskutujte.
M
C Úkoly
1. Jak můžeme na obr. 3.6 vidět, tak např. interpolace splajnem 5. stupně už nebude ideální, zkuste si
spočítat její Riemannův intergrál a opět diskutujte výsledky.
2. Zkuste změnit počet oddílů, na které se plocha pod křivkou dělí, jak moc to ovlivní výsledky?
3. Vyměňte Rungeovu funkci např. za funkci f(x0) = sin(x0) a zkuste spočítat výše uvedené Riemannovy
integrály.
i) Který z Riemannových integrálů bude nejblíže odpovídat hodnotě originální funkce?
ii) Zkuste změnit i počet oddílů, na které si plochu pod křivkou rozdělíte. Diskutujte výsledky.
C
Souhrn kapitoly
V této kapitole jste se dozvěděli, k čemu a jak se používá interpolace a extrapolace. Vysvětlili jsme Vám
různé druhy interpolací, popsali jsme, kde a jak se dají využít a pro jaké případy jsou či nejsou vhodné.
Pomocí programů psaných v programovacím jazyce Python, využitím hlavně knihovny Scipy a jejich
vestavěných funkcí, jste měli za úkol spočítat interpolační hodnoty v různých bodech křivky za použití
různých interpolačních metod. Naučili jste se, jakým způsobem ověřit kvalitu Vámi zvolené interpolační
metody i za použití Riemannova integrálu.
Další informace
• PRESS William H., TEUKOLSKY Saul A., VETTERLING William T., FLANNERY Brian P. Numerical
recipes. The art of scientific computing. Third edition. Cambridge. 2007, ISBN-13 978-0-521-
88068-8
• ČERMÁK, Libor a Rudolf HLAVIČKA. Numerické metody I. In: VUT, FSI, Ústav matematiky
[online]. Brno, 1.11.2006. Dostupné z:
http://mathonline.fme.vutbr.cz/Numericke-metody-I/sc-8-sr-1-a-11/default.aspx
Kapitola 4
Numerická řešení rovnic, hledání kořenů
Rychlý náhled: Tato kapitola se zabývá hledáním kořenů funkcí. Ukážeme si dvě základní metody a
to metodu půlení intervalu (bisekci) a Newtonovu Raphsonovu metodu (metodu tečen). Metody budeme
nejprve aplikovat na hledání kořenů kvadratické funkce. V druhé části si ukážeme aplikaci na hledání
rovnovážných bodů gravitačně vázaných systémů.
ºº Cíle studia:
• Pochopit matematické pozadí základních numerických metod pro hledání kořenů funkcí.
• Umět používat metody hledání kořenů, které jsou součástí knihovny scipy.
• Seznámit se s Rocheovým potenciálem a jeho vlastnostmi.
Všechny rovnice lze převézt na tvar f(x) = 0 a proto se budeme v této kapitole zabývat pouze rovnicemi
tohoto typu. Úloha tedy spočívá v nalezení x0 pro které platí f(x0) = 0. Hledání řešení této rovnice se
často označuje jako hledání kořenů (root finding). Pro řadu případů může být takových řešení více. Např.
kvadratická rovnice x2 − 4 = 0 má tato řešení dvě, a to x1 = 2 a x2 = −2. Některé úlohy naproti tomu
řešení v oboru reálných čísel nemají, např. rovnice x2 + 4 = 0.
Při řešení rovnice numericky je proto třeba znát o řešení nějaké informace (např. řešení existuje přibližně
zde, nebo řešení existuje na intervalu (a,b), který známe). Toto se může zdát jako značné omezení nicméně
při řešení reálných problémů v přírodních vědách takové informace často máme k dispozici.
My si nejdříve ukážeme numerické hledání kořenů kvadratické funkce1 a posléze se zaměří na složitější
funkci, které je ovšem zajímavá z hlediska výzkumu Vesmíru, konkrétně dvojhvězd.
4.1 Hledání kořenů kvadratické funkce
Začněme zopakováním základních vlastností kvadratické rovnice. Mějme kvadratickou rovnici
x2
+ ax + b = 0, (4.1)
1
Z terminologického hlediska je hledání kořenů funkce f(x) ekvivalentní s řešením rovnice f(x) = 0. Vždy ale hovoříme o
řešení (řešeních) rovnice, případně o kořenu (kořenech) funkce.
33
Kapitola 4 Numerická řešení rovnic, hledání kořenů 34
3 2 1 0 1 2 3 4 5
5
0
5
10
15
20 x2 2x+4
x2 2x+1
x2 2x 3
0 1 2 3 4 5 6 7 8
10
0
10
20
30
40
x0 x2 x1
x2 2x 5
p1
0 1 2 3 4 5 6 7 8
10
0
10
20
30
40
x2
x1
x2 2x 5
p1
Obrázek 4.1: Vlevo: Znázornění vlivu koeficientu b na kvadratickou funkci. Znázorněny jsou situace, kdy
existují dva kořeny (modře), jeden kořen (zeleně) a situace kdy kořen neexistuje (červeně). Uprostřed:
Znázornění metody půlení intervalů. Vpravo: Znázornění Newtonovy Raphsonovy metody.
kde uvažujeme, že člen u x2 je roven jedné. Pokud by byl člen u x2 různý od jedničky, můžeme celou rovnici
vydělit tímto koeficientem a převézt ji na tvar (4.1).
Z hlediska hledání řešení kvadratické rovnice mohou nastat tři situace v závislosti na hodnotě tzv.
diskriminantu, který v případě rovnice ve tvaru (4.1) vypadá D = a2 − 4b
• D < 0 rovnice (4.1) řešení nemá.
• D = 0 tj. b = a2/4 má rovnice jedno řešení x0 = −a/2 a rovnici lze přepsat na tvar (x + a/2)2 = 0.
• D > 0 má rovnice dvě řešení x1 = −a/2 +
√
D/2 a x2 = −a/2 −
√
D/2.
Podíváme-li se podrobně na výše uvedené vidíme, že pokud kořeny kvadratické funkce existují leží jeden z
nich vlevo od −a/2 a druhý leží vpravo od a/2. K podobnému závěru bychom došli, pokud bychom využili
známých vztahů koeficientů a a b ke kořenů rovnice x1 a x2
a = −(x1 + x2), (4.2)
b = x1x2. (4.3)
I z těchto vlastností vyplývá, že koeficient −a/2 je možné brát jako průměr koeficientů x1 a x2 a že intervaly
na kterých se nachází právě jeden kořen jsou tedy (−∞, −a/2) a (−a/2, ∞). Ke stejnému závěru jsme došli
již dříve ze znalosti kořenů.
4.1.1 Hledání kořenů metodou půlení intervalů
Metoda půlení intervalů je jednou ze základních metod hledání kořenů. Postup vychází z počáteční znalosti
intervalu na kterém se kořen nachází. V ukázce budeme hledat větší kořen tj. kořen který se nachází v
intervalu (−a/2, ∞). Nekonečno je z pochopitelných důvodů jako okrajový bod nevhodné, proto zkusíme
najít vhodnější bod. Uděláme odhad, že tímto bodem by mohl být bod x = a/2 a spočítáme hodnotu levé
strany rovnice v tomto bodě. Vzhledem k naší volbě koeficientu u x2 víme, že hodnota levé strany rovnice v
x0 = −a/2 musí být záporná, aby kořeny existovaly2. Pokud je hodnota v x = a/2 kladná, máme vyhráno
a víme, že kořen leží na intervalu (x, x ). Pokud je hodnota levé strany rovnice v bodě x záporná stejně
jako v bodě x0 proložíme těmito body přímku a nalezneme průsečík x1 této přímky z osou x. Vzhledem
k tomu, že přímka leží na intervalu (x , ∞) pod parabolou musí i v bodě, kde přímka protíná osu x ležet
2
Je-li koeficient u x2
> 0 je parabola obrácená vrcholem dolu (tj. v bodě x0 = −a/2 leží minimum kvadratické funkce).
Nachází-li se toto minimum na osou x, průsečíky s osou x neexistují a tím pádem nelze kořeny nalézt.
Kapitola 4 Numerická řešení rovnic, hledání kořenů 35
parabola nad ní a hodnota musí levé strany rovnice musí být v bode x1 kladná. Nyní tedy již bezpečně
známe interval na kterém kořen leží, označme ho (x0, x1).
Při zmenšování intervalu postupujeme následovně. Levou stranu rovnice si označíme jako f(x) a levé
strana rovnice je rovna nula. Zvolíme x2 = (x1 − x0)/2. Spočítáme f(x2). Je-li f(x2) > 0 víme, že kořen
leží na intervalu (x0, x2). Je-li f(x2) < 0 víme, že kořen leží na intervalu (x2, x1). Nový interval označíme
(x0, x1) spočítáme nové x2 a postupuje do doby než je splněna podmínka x1 − x0/(x1 + x0) < η, kde η
značí námi požadovanou přesnost. Řešení potom můžeme psát jako x = (x0 + x1)/2.
M Příklad
Následující část ukazuje využití metody půlení intervalů v pythonu.
1 from scipy import optimize as opt
2
3 def f(x):
4 return x**2-2*x - 5
5
6 koren=opt.bisect(f,1,7)
7
8 print(koren)
Program vypíše hodnotu 3.4494897427834985. Na řádku 1 importujeme z knihovny scipy rutinu optimize.
Na řádku 3 definujeme kvadratickou funkci. Hledaný koren potom najdeme metodou optimize.bisect,
kde prvním argumentem je funkce jejíž kořen hledáme, druhým a třetí argument určují interval na kterém
se kořen hledá – v našem případě se jedná o interval (1, 7).
M
C Úkol
Vyřešte následující úkoly.
1. Výše uvedený postup jsme aplikovali na hledání kořenu vpravo od minima kvadratické funkce. Jak
by se postup změnil pro hledání kořenu vlevo od minima? Jak byste poznali zda hledaný kořen leží
v intervalu (x0, x2) či (x2, x1) v obecném případě, kdy nevíte zda funkce f(x) na intervalu (x0, x1)
roste či klesá.
2. Pročtěte si dokumentaci scipy.optimize.bisect
3. Jakým způsobem můžeme nastavit přesnost s jakou se kořen hledá? Nastavte přesnost na různé
hodnoty a porovnejte výsledky.
4. Nastavte přesnost na poměrně vysokou hodnotu (např. xtol=0.01 a zkuste změnit pravý okraj intervalu
na kterém kořen hledáme. Má toto nastavení vliv na výsledek? Je výsledek vždy v požadované
toleranci od přesného řešení, které získáte přímým výpočtem?
5. Upravte jednoduchý skript tak, aby našel druhý kořen kvadratické rovnice.
C
Kapitola 4 Numerická řešení rovnic, hledání kořenů 36
4.1.2 Hledání kořenů Newtonovou – Raphsonovou metodou
Metoda půlení intervalů je vhodná pro úlohy kde známe interval na kterém se kořen nachází. Někdy můžeme
hledat kořeny funkce u které víme, kde přibližně kořen leží. Toho poté můžeme využít pro najití kořenu s
libovolnou přesností.
Metoda je znázorněna vlevo na Obr.4.1. Vyjdeme z prvotního odhadu místa, kde se kořen nachází.
V případě kvadratické funkce f(x) = x2 − 2x − 5 zvolíme bod vpravo od vrcholu paraboly (x0 = 1). V
ukázce jsme zvolili x1 = 7. V bodě x1 nalezneme funkční hodnotu (f(x1)) a numericky spočítáme derivaci
v tomto bodě. Z znalosti derivace spočítáme tečnou přímku v bodě x1. Na obrázku je znázorněna červeně
a označena jako p1. Spočítáme průsečík p1 s osou x a průsečík použijeme jako nový odhad kořenu x2.
Pokračujeme dále iterativně dokud nenalezneme kořen s požadovanou přesností.
M Příklad
Následující část ukazuje využití Newtonovy Raphsonovy metody v pythonu.
1 from scipy import optimize as opt
2
3 def f(x):
4 return x**2-2*x - 5
5
6 koren=opt.newton(f,7)
7
8 print(koren)
Program vypíše hodnotu 3.449489742783178. Na řádku 1 importujeme z knihovny scipy rutinu optimize.
Na řádku 3 definujeme kvadratickou funkci. Hledaný koren potom najdeme metodou optimize.newton,
kde prvním argumentem je funkce jejíž kořen hledáme, druhým argumentem je prvotní odhad na kořen, v
našem případě jsme zvolili x1 = 7. .
M
C Úkol
Vyřešte následující úkoly.
1. Pročtěte si dokumentaci scipy.optimize.newton
2. Jakým způsobem můžeme nastavit přesnost s jakou se kořen hledá? Nastavte přesnost na různé
hodnoty a porovnejte výsledky.
3. Upravte jednoduchý skript tak, aby našel druhý kořen kvadratické rovnice.
C
Další informace
• Podívejte se na další metody hledání kořenů funkcí. Při řešení vlastních problémů vždy dbejte o to,
abyste měli odhad jaké řešení hledáte. Např. si můžete vykreslete funkci, jejíž kořen hledáte apod.
Pokud používáte jiné programovací jazyky než python můžete využít např. numerical recipes nebo
GSL
Kapitola 4 Numerická řešení rovnic, hledání kořenů 37
• Dokumentaci k probíraným knihovnám naleznete online pro půlení intervalů i pro Newtonovu Raphsonovu
metodu
4.2 Pohyb částic v poli dvou kolem sebe obíhajících objektů
Zabývejme se situací kdy máme dvě gravitačně vázaná tělesa, která obíhají po kruhových trajektoriích
kolem vzájemného těžiště. Nás zajímá pohyb mnohem lehčího tělesa v jejich blízkosti. Taková situace je
například soustava Slunce a Jupiter, které obíhají kolem společného těžiště a nás zajímá pohyb asteroidů,
které jsou gravitačně ovlivňovány oběma tělesy.
Z hlediska výpočetní náročnosti je vhodné zvolit souřadnou soustavu ve které je Slunce i Jupiter v
klidu a výpočet provádět v ní. Asteroid, který je v klidu v této soustavě ve skutečnosti obíhá kolem středu
soustavy stejnou úhlovou rychlostí jako Jupiter.
4.2.1 Teorie
Uvažujeme dvě tělesa o hmotnostech M1 a M2 ve vzdálenosti a, kde jsme hmotnější těleso označili indexem
1 a platí tedy M1 > M2. My se omezíme na jednoduchý příklad, kdy nás budou zajímat drobná tělesa,
která se pohybují ve stejné rovině jako vzájemně se obíhající tělesa. V případe soustavy Slunce – Země by
se jednalo o tělesa která se pohybují v rovině ekliptiky. My si tuto rovinu označíme jako rovinu x, y. Tělesa
umístíme na osu x a do počátku soustavy souřadnic umístíme hmotný střed.
Je-li hmotný střed v bodě [0, 0] pak pro souřadnice těles musí platit
M1r1 = M2r2, (4.4)
kde ri je vzdálenost daného tělesa od hmotného středu. Zároveň známe vzdálenost obou těles a, tedy platí
r1 + r2 = a. Hmotnější těleso má tedy souřadnice [x1, 0], těleso méně hmotné má souřadnice [x2, 0], kde
x1,2 jsou dány následovně
x1 = −r1 = −
M2a
M1 + M2
,
x2 = r2 =
M1a
M1 + M2
. (4.5)
Úhlová rychlost, kterou tělesa obíhají kolem společného hmotného středu je dána Keplerovým zákonem
ω2
= G
M1 + M2
a3
, (4.6)
kde G = 6.6724 × 10−11 je gravitační konstanta3.
Pohybové rovnice v námi uvažované soustavě je možné zapsat ve tvaru
¨x = −
∂UR
∂x
+ 2ωvy,
¨y = −
∂UR
∂y
− 2ωvz, (4.7)
kde UR je Rocheův potenciál definovaný následujícím způsobem
UR(x, y) = −
GM1
d1
−
GM2
d2
−
1
2
ω2
(x2
+ y2
), (4.8)
3
Hodnotu gravitační konstanty uvádíme dle aktuálních měření ze stránek NIST.
Kapitola 4 Numerická řešení rovnic, hledání kořenů 38
kde d1,2 jsou vzdálenosti bodu o souřadnicích (x, y) od těles 1, 2 a můžeme je psát následujícím způsobem
d2
1 = (x − x1)2
+ y2
, (4.9)
d2
2 = (x − x2)2
+ y2
. (4.10)
V Rocheově potenciálu (4.8) mají první dva členy význam gravitačních potenciálů polí generovaných tělesy
o hmotnostech M1, M2. Třetí člen odpovídá potenciálu odstředivé síly, která působí v rotující soustavě i
na částici která je v této soustavě v klidu.
V rovnicích (4.7) odpovídá druhý člen příslušné složce Coriolisova zrychlení
aCor = −2ω × v. (4.11)
V našem případě je ω = (0, 0, ω) a v = (vx, vy, 0). Coriolisovo zrychlení nepůsobí na částice, které jsou v
soustavě v klidu, ale působí na pohybujících se částice.
Rocheův potenciál graficky znázorňujeme na Obr.4.2. Potenciál má pět lokálních extrémů, které značíme
Li a nazýváme je Lagrangeovými body. V těchto bodech je ∂UR/∂x = ∂UR/∂y = 0. Bod L1 leží mezi
tělesy, body L2 a L3 leží na ose x vně soustavy. Body L4 a L5 leží ve vrcholech rovnostranných trojúhelníků
se základnou danou tělesy M1 a M2.
Body L1, L2 a L3 jsou nestabilní a částice v nich umístěná by vlivem malé perturbace uletěla nakonec
pryč. Body L4 a L5 jsou stabilní pokud je4
M1 > 25M2
1 + 1 − 4/625
2
. (4.12)
K tomu dochází např. v soustavě Slunce – Jupiter, kde se v bodech L4 a L5 pohybují tzv. Trójské asteroidy.
Naproti tomu v binárních systémech dvou hvězd tato situace většinou nenastává.
4.2.2 Implementace v Pythonu
V repozitáři na githubu si můžete stáhnout skript, který řeší pohyb asteroidu v blízkosti bodu L4 soustavy
Slunce – Jupiter. My se zde zaměříme na hledání bodů L1, L2 a L3 pomocí Newtonovy - Raphsonovy
metody.
Začneme definováním základních veličin
1 DistanceJ = 778570000000. # m JUPITER FROM SUN
2 G = 6.6724*10**-11
3 Jupiter_mass = 1.8982*10**27 # kg
4 Sun_mass = 1.989*10**30 # kg
5
6
7 M1=Sun_mass
8 M2=Jupiter_mass
9 a=DistanceJ
10 Ang_vel=math.sqrt(G*(M1+M2)/(a**3))
11 r2=M1*a/(M1+M2)
12 r1=M2*a/(M1+M2)
Nadefinujeme Rocheův potenciál UR(x, y, vx, vy)
4
Zájemci o odvození se mohou podívat např. zde.
Kapitola 4 Numerická řešení rovnic, hledání kořenů 39
2.01.51.00.50.00.51.01.52.0
2.0
1.5
1.0
0.5
0.0
0.5
1.0
1.5
2.0
-3.01-2.73-2.45-2.18-1.90-1.62-1.34-1.07-0.79-0.51
Obrázek 4.2: Znázornění Rocheova potenciálu. Souřadnice na osách x a y jsou v jednotkách a.Nahoře:
Hodnoty Rocheova potenciálu (v jednotkách ω2a2 pro M2 = 0.8 M1. Vlevo dole: Ekvipotenciály pro
hodnoty, které odpovídají Lagrangeovým bodům. Ty jsou označeny znázorněny obrázku jako L1...L5. I zde
je M2 = 0.8 M1. Vpravo dole: Stejné jako vlevo dole s tím, že zde je M2 = 0.2 M1
Kapitola 4 Numerická řešení rovnic, hledání kořenů 40
1 def pot(x,y):
2 x1=-r1
3 x2=r2
4 d1=math.sqrt((x-x1)**2 + y**2 )
5 d2=math.sqrt((x-x2)**2 + y**2 )
6 result = -G*(M1/d1 + M2/d2) -1.*Ang_vel*Ang_vel*(x**2 + y**2)/2.
7 return result
Zrychlení ve směru osy x ax = ∂UR(x, y, vx, vy)/∂x je definováno ve funkci acx, zatímco acxax je speciální
případ zrychlení působící na částici v klidu (vx = vy = 0) na ose x (y = 0). My zde implementujeme vlastní
derivaci, ale mohli bychom využít i předdefinované (např. derivative z knihovny scipy.misc).
1 def acx(x,y,vx,vy):
2 dx=a/1000.
3 result=-(pot(x+dx,y) - pot(x-dx,y))/(2.*dx) + 2.* Ang_vel*vy
4 # result=-(-pot(x+2.*dx,y) + 8.*pot(x+dx,y) - 8.*pot(x-dx,y) + pot(x-2.*dx,y))/(12.*dx)
5 # + 2.* Ang_vel*vy
6 return result
7
8 def acx_ax(x):
9 return acx(x,0,0,0)
Hledání extrémů potenciálu UR na ose x můžeme přepsat na hledání kořenů funkce ax(x). Tedy hledáme
místa, kde na umístěnou částici nepůsobí žádná síla.
Definujme parametr α = M2/(M1 + M2). Je-li tento parametr malý α 1 potom je možné vyjádřit
pozici bodů Li[xi, 0] vztahy
x1 = a 1 −
α
3
1/3
(4.13)
x2 = a 1 +
α
3
1/3
(4.14)
x3 = −a 1 + a
5
12
α (4.15)
. Tyto body využijeme jako odhad kořenu pro Newtonovu Raphsonovu metodu.
1 alpha=M2/(M1+M2)
2 P0=a*a*Ang_vel*Ang_vel
3
4 L1X_0=a*(1.-(alpha/3.)**(1./3.))
5 L1X=opt.newton(acx_ax,L1X_0,tol=a*10**-5)
6 L1Y=0.
7 P1=pot(L1X,L1Y)/P0
8
9 print(L1X_0/a,L1X/a)
V uvedené funkci je P0 jednotka, která se používá pro měření potenciálu. L1X_0 značí původní odhad a
L1X přesnější řešení. Jako požadovanou přesnost jsme uvedli 10−5a. Je důležité vědět s jak velké hodnoty
kořen má a nastavit přesnost na rozumnou hodnotu. Pro soustavu Slunce – Jupiter dostaneme
Kapitola 4 Numerická řešení rovnic, hledání kořenů 41
L1X_0=0.931757139893148 a L1X=0.9323711678782919. Zkusíme si nyní nastavit hodnotu M2 = 0.1M0.
Dostáváme L1X_0=0.688234046118128 a L1X=0.6266027208867219. Vidíme, že nyní je přesná hodnota
od původní dále, což jsme očekávali, neboť poměr hmotností byl nyní menší.
C Úkol
Vyřešte následující úkoly.
1. Najděte polohy dalších Lagrangeových bodů.
2. Prozkoumejte chování pro další poměry hmotností.
3. Upravte kód pro soustavu Slunce – Země a pro soustavu Země – Měsíc.
C
Další informace
• Stáhněte si kód z Github a zkuste si vykreslit trajektorie pro různé kombinace hmotností a různé
počáteční polohy asteroidů.
Souhrn kapitoly
V této kapitole jsme se zabývali metodami hledání kořenu funkcí. Newtonovu Raphsonovu metodu jsme
poté aplikovali na hledání rovnovážných poloh v poli dvou kolem sebe obíhajících objektů.
Kapitola 5
Statistický popis dat, distribuční funkce
Adam Hofer
Rychlý náhled: Tato kapitola se zabývá statistikou, jež nám slouží k popisu libovolného datového
souboru. Jedná se jak o popis dat za pomoci charakteristik, což je popsáno základní deskriptivní statistikou.
Dále jsou zde popsány nejčastěji používané distribuční funkce
Klíčová slova: statistika, popis dat, distribuční funkce, normální rozdělení, PDF, CDF, variabilita,
průměr
ºº Cíle studia: V této kapitole je cílem studenta naučit pracovat se základními statistickými metodami
pro popis dat. Po přečtení tak bude student schopen říct základní informace o libovolných datech a popsat
je za pomocí základních charakteristik. Dále budete schopní v Pythonu naprogramovat jednoduché funkce
pro zpracování dat a jejich základní statistický popis.
Pokud chceme dobře porozumět světu okolo nás, je vhodné znát alespoň základní statistické metody.
Proto si v této kapitole vybrané metody popíšeme. Jelikož naše okolí lze často popisovat za pomoci různých
měření a následně tak naměřenými daty. Ty však mohou být poměrně rozsáhlé, tudíž je vhodné začít
s první částí statistiky a tou je statistika popisná neboli deskriptivní. Ta se nám naměřený soubor dat
snaží popisovat pomocí několika číselných hodnot, jež o daných datech něco vypovídají. Příkladem může
být měření teploty v průběhu dne. Dejme tomu, že budeme měřit teplota každou čtvrthodinu. Za jeden den
tak získáme 96 hodnot. Z těchto hodnot máme velmi podrobnou představu o počasí v daném dni. Nicméně
když budeme chtít někomu dalšímu popsat jaký bylo ten daný den počasí sdělovat mu všech 96 hodnot by
bylo jistě nadbytečné. Proto si vystačíme s několika charakteristikami jež nám ten soubor teplot popíšou.
Může to být například minimum, maximum a průměr. K získaní představy o počasí v daném dni mi tak
stačí tyto tři hodnoty a to je v základu princip popisné statistiky. Samozřejmě charakteristik existuje vícero
než v tomto našem jednoduchém příkladě. Proto se s nimi pojďme seznámit:
5.1 Základní charakteristiky popisné statistiky
Tyto charakteristiky jsou nám schopny dobře popsat datový soubor a říct nám základní informace o
obsažených datech.
42
Kapitola 5 Statistický popis dat, distribuční funkce 43
5.1.1 Minimum
nejnižší hodnota v datovém souboru. Pro její nalezení nám v Pythonu slouží funkce min(list)
5.1.2 Maximum
nejvyšší hodnota v datovém souboru, kterou můžeme nalézt použitím funkce max(list)
5.1.3 Počet hodnot
Neboli celkový počet vzorků v datovém souboru. Tento údaj o počtu hodnot v datovém souboru je nezbytný
pro většinou dalších výpočtů. V Pythonu jej můžeme najít funkcí len(list), kterou aplikujeme na datový
soubor.
5.1.4 Aritmetický průměr
Základní charakteristika popisující střed dat. Anglicky jej nazýváme mean a označujeme jej ¯x. Spočteme
jej tak, že sečteme veškeré hodnoty v datovém souboru a tento součet vydělíme počtem hodnot. Tato
charakteristika je nicméně náchylná na odlehlé hodnoty, které nám mohou výsledný průměr zkreslovat.
Příkladem budiž platy. Například máme 9 dělníku pracujících za minimální mzdu a jednoto top manažera
s milionovým příjmem. Pokud spočteme aritmetický průměr jejich platů získáme částku která rozhodně
nepopisuje reálný příjem většiny zaměstnanců v dané firmě, ale je zkreslena odlehlou hodnotou platu
manažera. V Pythonu můžeme (díky statistic module) aritmetický průměr spočíst funkcí mean(list).
Případně můžeme využít knihovnu Numpy a funkci numpy.mean().
¯x =
1
n
n
i=1
xi (5.1)
C Úkol
Napište vlastní funkci pro výpočet aritmetického průměru. Takže nám ale nic nebrání napsat si vlastní
funkci za pomoci sum() a length(). Funkcí sum prvně sečteme položky dle předchozího vzorce a následně
vydělíme počtem prvků o což se nám postará funkce length.
C
5.1.5 Medián
Je charakteristika středu a je to hodnota u níž platí, že 50 % dat je větších nebo rovno než medián a 50 %
dat je menších nebo rovno. Při sudém počtu vzorků pak spočítáme aritmetický průměr dvou prostředních
čísel. V Pythonu jej spočteme funkcí median(list). Označujeme jej ˜x.
5.1.6 Modus
Další z charakteristik, kdy spočítáme počty výskytů jednotlivých prvků a modus je ten nejčastější. V Pythonu
se funkce nazývá mode(list). Označujeme jej ˆx.
Kapitola 5 Statistický popis dat, distribuční funkce 44
5.1.7 Kvantily
Jde o hodnoty, které nám rozdělují setříděná data podobně jako medián. Medián je tedy speciální případ
kvantilů. Tedy kvantil x0.33 rozděluje data tak, že 33 % dat je menších nebo rovno než x0.33a 66 % je větších
nebo rovno než x0.33. V Pythonu jej díky knihovně numpy můžeme spočítat funkcí numpy.quantile().
5.1.8 Rozpětí
Základní a nejjednodušší míra variability. Získáme jej tak, že stačí odečíst minimum a maximum. Na data
bez extrému (outliers) to stačí, ale co třeba na (1, 2, 3, 4, 500)? Rozpětí vychází 500 − 1 = 499, přitom
většina dat je z rozsahu 1 až 4.
5.1.9 Mezikvartilové rozpětí a odchylka
Řeší problémy s extrémy. Budeme pracovat s rozdíly kvartilů na rozdíl od klasického rozpětí. Pokud budeme
chtít popsat, jak se data odchylují od mediánu, vydělíme jej dvěma.
5.1.10 Rozptyl
Nejpoužívanější míra rozptýlenosti. Spočteme jej tak, že sečteme druhé mocniny všech odchylek od průměru.
Jde o celkovou míru rozptýlenosti ta nám ale moc neříká, jak jsou data průměrné rozptýlená. Označujeme
jej σ2 či s2 či V ar(X). Lze jej spočíst za pomoci následujícího vzorce
s2
=
1
n − 1
n
i=1
(xi − ¯x)2
(5.2)
Pro výpočet za pomoci Pythonu můžeme použít funkci numpy.var().
5.1.11 Směrodatná odchylka
Určuje průměrnou odchylku od středu. Stačí odmocnit rozptyl. Bývá označována řeckým písmenem σ či s.
Jedná se o často používanou míru variability, která nám vypovídá o tom, nakolik se od sebe navzájem liší
jednotlivé prvky v souboru zkoumaných dat. V praxi lze říct je-li směrodatná odchylka malá, jsou si prvky
souboru dat povětšinou navzájem podobné, naopak je-li velká jsou prvky odlišné.
σ = s =
1
n − 1
n
i=1
(xi − ¯x)2 (5.3)
Směrodatnou odchylku lze vypočíst pomocí funkce numpy.std() v Pythonu.
5.1.12 Variační koeficient
Slouží ke studiu, zda není s daty něco podivného. Vydělíme směrodatnou odchylku průměrem a počítáme
jej v procentech. Hodnoty větší než 15 % až 30 % (záleží jakého charakteru mám data) svědčí o nějakém
problému (třeba kdybychom porovnávali stáří babičky a vesmíru).
vx =
sx
¯x
× 100 (5.4)
Kapitola 5 Statistický popis dat, distribuční funkce 45
5.2 Distribuce pravděpodobnosti
Pro popis dat využíváme obvykle číselných hodnot, proto je pro popis některých dějů v přírodě nutno
stanovit pravidla a každému jevu přiřadit hodnotu či pravděpodobnost s jakou ten jev nastává.
5.2.1 Náhodná veličina
Náhodná veličina X je reálná funkce definovaná na množině všech elementárních jevů, která každému jevu
přiřadí reálné číslo. Můžeme ji rozdělit na spojitou a diskrétní. V případě spojité náhodné veličiny je obor
hodnot M otevřený nebo uzavřený interval. Pro diskrétní náhodnou veličinu je obor hodnot M konečná
nebo nekonečná posloupnost.
Náhodnou veličinu lze také popsat jako libovolná veličinu, kterou je možné opakovaně měřit u různých
objektů, v různých místech nebo v různém čase ale její hodnoty závisí na náhodě. Tyto výsledky lze pak
následně zpracovat.
Příkladem může být házení mincí, kde jevu když padne rub přiřadíme číselnou hodnotu 0, a jevu, kdy
padne líc přiřadíme 1.
5.2.2 Pravděpodobnostní funkce
Pravděpodobnostní funkce (anglicky probability mass function, zkráceně PMF) je funkce, která udává
pravděpodobnost, že diskrétní náhodná veličina se přesně rovná nějaké hodnotě [1]. Pravděpodobnost
jednotlivých hodnot x označujeme P(x) a obvykle je zobrazujeme pomocí tabulky a grafu.
Obrázek 5.1: Pravděpodobnostní funkce (vlevo) a pravděpodobnostní funkce pro poctivou hrací
kostku(vpravo).
5.2.3 Distribuční funkce diskrétní veličiny
Pomocí pravděpodobnostní funkce lze zavést distribuční funkci vztahem
F(x) = P[X ≤ x] (5.5)
Distribuční funkce je vždy neklesající, je spojitá zprava a její hodnoty leží v intervalu od 0 do 1.
Průběhy distribučních funkcí a jejich vlastnosti jsou dány typem jejich rozdělení, základní seznam
diskrétních rozdělení je zde:
• Alternativní rozdělení
• Binomické rozdělení
Kapitola 5 Statistický popis dat, distribuční funkce 46
• Poissonovo rozdělení
• Negativně binomické rozdělení
• Pascalovo rozdělení
• Geometrické rozdělení
• Hypergeometrické rozdělení
• Logaritmické rozdělení
Více informací o jednotlivých rozděleních naleznete v libovolné literatuře o teorii pravděpodobnosti či
statistice.
5.2.4 Distribuční funkce spojité veličiny
Distribuční funkce anglicky označovaná jako cumulative distribution function zkráceně CDF je definována
vztahem:
FX(x) = P[X ≤ x] (5.6)
V případě, že X je spojitá náhodná proměnná s hustotou f, potom platí
FX(x) =
x
−∞
f(t)dt (5.7)
Průběhy distribučních funkcí a jejich vlastnosti jsou dány typem jejich rozdělení, základní seznam
spojitých rozdělení je zde:
• Rovnoměrné rozdělení
• Normální (Gaussovo) rozdělení
• Logaritmicko-normální rozdělení
• Exponenciální rozdělení
• Cauchyho rozdělení
• Maxwellovo-Boltzmannovo rozdělení
• Studentovo rozdělení
• Fisherovo-Snedecorovo rozdělení
• χ2 rozdělení (Chí kvadrát)
Více informací o jednotlivých rozděleních naleznete v libovolné literatuře o teorii pravděpodobnosti či
statistice.
Pro vykreslení funkce za pomoci Pythonu využijeme knihovnu SciPy a její funkci cdf(x, loc=0,
scale=1).
Kapitola 5 Statistický popis dat, distribuční funkce 47
Obrázek 5.2: Příklady distribučních funkcí pro normální rozdělení.
5.2.5 Hustota pravděpodobnosti
Hustota pravděpodobnosti se anglicky nazývá Probability Density Function proto bývá označována jako
PDF. Tuto funkci, když zintegrujeme na kterémkoli jejím vzorku, vyjde nám relativní pravděpodobnost,
že hodnota náhodné proměnné by se rovnala tomuto vzorku.[2] Zobrazuje nám hustotu rozložení pravdě-
podobnosti.
Ω
ρ(x)dx = 1 (5.8)
Pro vykreslení funkce za pomoci Pythonu využijeme knihovnu SciPy a její funkci pdf(x, loc=0,
scale=1).
Obrázek 5.3: Příklady hustot pravděpodobností pro normální rozdělení.
Kapitola 5 Statistický popis dat, distribuční funkce 48
5.2.6 Centrální limitní věta
jsou-li Xi nezávislé náhodné veličiny s konečným rozptylem, pak výběrový průměr má při dostatečně velkém
počtu pozorování přibližně normální rozdělení, ať už Xi pocházejí z libovolného rozdělení.
5.2.7 Normální rozdělení
Normální rozdělení neboli Gaussovo rozdělení je jedno z nejvýznamnějších rozdělení a to díky Centrální
limitní větě. Toto rozdělení tak popisuje velké množství jevů a pro většinu aplikací si s ním vystačíme. Normální
rozdělení je definováno 2 parametry (neboli charakteristikami) a to rozptylem σ2 a střední hodnotouµ.
Hustota pravděpodobnosti tohoto rozdělení má tvar
f(x) =
1
√
2πσ2
e−
(x−µ)2
2σ2 (5.9)
Typickým příkladem normálního rozdělení je například výška dospělých lidí či velikost jejich inteligence.
Jelikož nejvíce lidí má průměrnou výšku či inteligenci je křivka v tomto bodě nejvyšší a do obou stran
postupně klesá směrem k nule. To vypovídá o tom, že extrémně malých či velkých lidí je velmi málo a to
stejné platí i o inteligenci.
C Úkoly
1. Budeme losovat z 5 bílých a 9 černých míčků. Náhodná veličina X nám představuje počet bílých
míčků mezi 5 vybranými. Vytvořte pravděpodobnostní a distribuční funkci této náhodné veličiny.
2. Zkusme si následující pokus, kde budeme házet virtuální poctivou mincí. Hodíme třicetkrát mincí a
budeme počítat kolikrát padne rub (1) či líc (0). Tento následně pokus budeme opakovat 100 krát.
K tomu využijeme Python a vykreslete si hustotu pravděpodobnosti.
C
Souhrn kapitoly
Statistické zpracování dat je obsáhlá disciplína, jež by s klidem zabrala samostatný předmět, tato kapitola
je pouze shrnutím základních a nejpoužívanějších pojmů používaných při základním zpracování dat.
? Otázky
1. Jaké známe charakteristiky středu a jak se liší?
2. K čemu nám slouží směrodatná odchylka?
3. Proč je zrovna tak důležité Normální rozdělení?
?
Další informace
• PRESS William H., TEUKOLSKY Saul A., VETTERLING William T., FLANNERY Brian P. Numerical
recipes. The art of scientific computing. Third edition. Cambridge. 2007, ISBN-13 978-0-521-
88068-8
Kapitola 5 Statistický popis dat, distribuční funkce 49
• HUFF, Darrell. Jak lhát se statistikou. Nakladatelství BRÁNA A.S.. Praha,2013. ISBN-978-80-7243-
623-1
• FEYNMAN, Richard Phillips, Robert B. LEIGHTON a Matthew SANDS. Feynmanovy přednášky
z fyziky: revidované vydání s řešenými příklady. 2. vydání. Přeložil Ivan ŠTOLL. Praha: Fragment,
2013. ISBN 978-80-253-1642-9.
Kapitola 6
Řešení okrajové úlohy relaxační metodou
Jiří Horák
Rychlý náhled: Tato kapitola se zabývá numerickým řešením systému obyčejných diferenciálních rovnic
1. řádu relaxační metodou, konkrétně řešení okrajové úlohy. Studenti se podrobně seznámí se způsobem
implementace této metody v jazyce Python. Součástí je také vzorové řešení Gorillas problému.
Klíčová slova: numerické řešení diferenciálních rovnic, obyčejné diferenciální rovnice, relaxační metoda,
okrajová úloha
ºº Cíle studia: Po prostudování kapitoly budete schopní v Pythonu řešit systémy diferenciálních rovnic
se zadanými okrajovými podmínkami, například pohyb těles v gravitačním poli země. Kapitola je dále
základem pro následující kapitolu – Bondiho akrece, kde bude tato metoda využita na složitější systém
rovnic.
6.1 Okrajová úloha
Relaxační metoda je užitečný nástroj pro řešení jedno-dimenzionálních okrajových úloh. Dále uvažujeme
set obyčejných diferenciálních rovnic 1. řádu ve tvaru:
dY
dx
= F (x, Y) , xa ≤ x ≤ xb, (6.1)
kde x je nezávislá proměnná, Y(x) = [Y 0(x), . . . , Y N−1(x)] je vektor závislých proměných (proměnné
jednotlivých řešení Y α(x) jsou značeny řeckými indexy) a F : RN+1 → RN jsou nelineární funkce popisující
pravé strany rovnic.
Set N diferenciálních rovnic by měl být ohraničen N okrajovými podmínkami (omezeními). Obecně
mohou být okrajové podmínky z obou stran oblasti řešení (ve speciálním případě, kdy jsou všechny jen na
jedné straně, se úloha stává počáteční úlohou). Obecně mohou být okrajové podmínky psány
A [Y (xa)] = 0, B [Y (xb)] = 0, (6.2)
kde A a B jsou dvě nelineární funkce, A : RN → RNa a B : RN → RNb , kde Na a Nb je počet omezení
aplikován na xa a xb, a Na + Nb = N.
50
Kapitola 6 Řešení okrajové úlohy relaxační metodou 51
6.1.1 Diskretizace
Soustavu obyčejných diferenciálních rovnic budeme řešit na nepravidelné síti tvořené M body xa =
x0 . . . xi, . . . xM−1 = xb. Rovnice vyhodnotíme na středových bodech (xi+1/2) každého segmentu. Derivace
jsou jen přibližné, vypočítané pomocí metody centrální diference, hodnoty řešení Yi+1 jsou aritmetickým
průměrem Yi a Yi+1 na dvou sousedních bodech. To je
dY
dx i+1/2
≈
Yi+1 − Yi
xi+1 − xi
, Fi+1/2 ≈ F
1
2
(xi + xi+1) ,
1
2
(Yi+1 + Yi) . (6.3)
Takto obdržíme N × (M − 1) algebraických rovnic popisujících obyčejné diferenciální rocnice v každém
segmentu (0 ≤ i ≤ M − 2),
Ei (Yi, Yi+1) ≡ Yi+1 − Yi − Fi+1/2 (xi+1 − xi) = 0, (6.4)
plus v celkem N okrajových oblastech na koncích domény,
A (Y0) = 0, B (YM−1) = 0 (6.5)
pro N × M hodnot všech nezávislých proměnných na všech bodech sítě Y α
i . Formálně to můžeme zapsat:
G(Y) ≡









A (Y0)
E0 (Y0, Y1)
...
EM−2 (YM−2, YM−1)
B (YM−1)









= 0. (6.6)
kde
Y ≡ [Y0, . . . , YM−1] = Y 0
0 , Y 1
0 , . . . , Y α
i , . . . Y N−1
M−1 . (6.7)
6.1.2 Iterativní řešení
Rovnice G(Y) = 0 může být řešena iterativně použitím Newton-Ralphsonovy metody. Předpokládejme, že
máme počáteční odhad řešení Y(n). Hledáme vylepšené řešení Y(n+1) = Y(n) +δY takové, že G(Y(n+1)) = 0.
Expanzí v lineárním řádu dostaneme
G Y(n+1)
= G Y(n)
+
i,α
∂G
∂Y α
i
δY α
i = 0. (6.8)
To je lineární rovnice, která může být řešena použitím Gaussovy eliminace.
Kapitola 6 Řešení okrajové úlohy relaxační metodou 52
Všimněte si, že struktura matice L ≡ ∂G/∂Y α
i je následující:
L =


























∂Aβ
∂Y α
∂Eβ
0
∂Y α
0
∂Eβ
0
∂Y α
1
∂Eβ
1
∂Y α
1
∂Eβ
1
∂Y α
2
. . . . . .
∂Eβ
M−2
∂Y α
M−2
∂Eβ
M−2
∂Y α
M−1
∂Bβ
∂Y α


























. (6.9)
Je tomu tak proto, že Ei závisí pouze na hodnotách řešení ve dvou sousedních bodech sítě. Máme
∂Eβ
i
∂Y α
j
= δβ
α δj
i+1 − δj
i −
1
2
∂Fβ
∂Y α
i+1/2
(xi + xi+1) δj
i + δj
i+1 (6.10)
6.2 Implementace
Nejprve importujeme nezbytné moduly numpy pro manipulaci s vektory a maticemi, a matplotlib pro
tvorbu grafů.
1 import numpy as np
2 import math
3 import matplotlib
4 import matplotlib.pyplot as plt
5
6 # some settings of matplotlib:
7 from matplotlib import rc
8 rc('text', usetex=True) # ... in order to use LaTex fonts in the graphs.
9 rc('legend', frameon=False) # ... as a default, legends without the boundary frame.
Nyní implementujeme jeden krok relaxace. Funkce relaxation_step bere jako vstup:
• Array hodnot nezávislé proměnné, X = [x0, . . . , xM−1].
• Vektor řešení všech hodnot závislé proměnné, Y = [Y 0
0 , Y 1
0 , . . . , Y α
i , . . . Y N−1
M−1 ] popisující počáteční
odhad.
• Pravou stranu rovnice F(x, Y) = [F1(x, Y), . . . FN−1(x, Y)].
• Jakobián pravé strany rovnice jacF = ∂(F0, . . . , FN−1)/∂(Y 0, . . . , Y N−1)
• Dvě funkce A a B, udávající okrajové podmínky na obou stranách oblasti řešení, A(Y) = 0 v x = xa
a B(Y) = 0 v x = xb, a jejich jakobiány ∂(A0, . . . , ANa−1)/∂(Y 0, . . . , Y N−1)
a ∂(B0, . . . , BNb−1)/∂(Y 0, . . . , Y N−1).
Kapitola 6 Řešení okrajové úlohy relaxační metodou 53
Funkce vyhodnotí vektor G a jakobián L. Nakonec řeší rovnici
LδY = −G (6.11)
Výstup funkce je vektor δY.
1 def relaxation_step(X, Y, F, jacF, A, jacA, B, jacB):
2
3 # We find the number of mesh points M and number of the solution variables N
4 # from the sizes of the solution vectors X and Y:
5 M = len(X)
6 N = len(Y)//M
7 if len(Y)%M != 0:
8 raise ValueError('Incompatible sizes of X and Y arrays.')
9
10 # Create the empty vector G and empty matrix L:
11 G = np.zeros(N*M)
12 L = np.zeros((N*M, N*M))
13
14 # Filling G and L:
15 #------------------
16
17 # We start with the constrains at left border of the domain (the block dA/dY):
18 # We assign the number of constraints Na as a size of A(Y_0)
19
20 # Evaluating A(Y_0).
21 # Note that vector Y_0 corresponds to first N elements of the vector Y.
22 Aval = A(Y[0:N])
23 Na = len(Aval)
24 if Na > 0:
25 G[0:Na] = Aval
26 L[0:Na, 0:N] = jacA(Y[0:N])
27
28 # Now we calculate and place contributions of the individual segments:
29 # (the index i numbers the segments)
30 for i in range(M-1):
31 # Find the midpoint (x) and the size of the segment (dx):
32 x = (X[i+1] + X[i])/2
33 dx = X[i+1] - X[i]
34 # Find the average value of the solutions y = (Y_i + Y_{i+1})/2
35 # and the difference dy = Y_{i+1} - Y_i.
36 # Note that the vector Y_i is located in the solution vector Y[] at
37 # indices i*N, i*N+1, ..., i*N+N-1.
38 y = (Y[i*N:(i+1)*N] + Y[(i+1)*N:(i+2)*N])/2
39 dy = Y[(i+1)*N:(i+2)*N] - Y[i*N:(i+1)*N]
40 # Fill the vector G: Note that the contributions has to be placed AFTER the
41 # boundary conditions (first Na elements)
42 G[Na+i*N:Na+(i+1)*N] = dy - F(x, y)*dx
43 # Fill the matrix L:
44 # First calculate the gradient of F with respect to the solution variables,
45 jF = jacF(x, y)
46 delta = np.identity(N)
47 L[Na+i*N:Na+(i+1)*N, i*N:(i+1)*N] = -dx*jF/2 - delta
Kapitola 6 Řešení okrajové úlohy relaxační metodou 54
48 L[Na+i*N:Na+(i+1)*N, (i+1)*N:(i+2)*N] = -dx*jF/2 + delta
49
50 # Finally the constrains at right end of the domain:
51 if Na < N:
52 # There should be N-Na constraints.
53 # Evaluate the boundary function B(Y_{M-1}) & Jacobian:
54 Bval = B(Y[N*(M-1):])
55 Nb = len(Bval)
56 if Na + Nb != N:
57 raise
ValueError('The number of boundary constraints ({}) does not match the number of variables
N))
→
→
58 G[Na+(M-1)*N:] = Bval
59 L[Na+(M-1)*N:,N*(M-1):] = jacB(Y[N*(M-1):])
60
61 # The matrix L and the vector G is now assembled.
62 # Solve L.dY = -G for the improvement dY:
63 dY = np.linalg.solve(L, -G)
64
65 return dY
Tato funkce je implementována v modulu relaxation.py.
MM Gorilla problem
Kdysi dávno existovala velmi populární počítačová hra napsaná v MS-DOS q-basic, jmenovala se Gorillas
nebo GORILLAS.BAS, podle jména souboru se zdrojovým kódem. Hra byla velmi jednoduchá. Byli v ní dvě
gorily, které po sobě házely vybuchující banány nad obrysem města. Hra byla určená pro dva hráče, každý
musel zadat počáteční rychlost a úhel hodu tak, aby trefil druhého hráče.
Pro ty z vás, kteří se zajímají o počítačovou prehistorii, tady je odkaz na online verzi:
https://www.retrogames.cz/play654 − DOS.php
Chceme řešit stejný problém použitím relaxační metody. Podívejme se na problém takto: Vezmeme
střelce (první gorilu) a cíl (druhou gorilu) ve vzdálenosti d. Pro jednoduchost, předpokládejme že střelec
i cíl jsou ve stejné výšce a zanedbejme jakoukoli interakci objektů (banánů) s okolním vzduchem. Jaká má
být počáteční rychlost a úhel tak, abychom trefili cíl? Protože problém vede na nekonečný počet řešení,
musíme specifikovat čas letu objektu tb.
Ze znalostí základní fyziky víme, že pohyb bude určen tímto systémem diferenciálních rovnic:
¨x = 0, ¨y = −g, (6.12)
kde x(t) a y(t) popisují trajektorii banánu a g je tíhové zrychlení. Naše situace je popsána čtyřmi okrajovými
podmínkami
x(0) = 0, y(0) = 0, x(tb) = d, y(tb) = 0. (6.13)
Dvě z nich jsou určeny v t = 0, delší dvě v t = T. Tyto představují okrajové podmínky úlohy.
Chceme úlohy řešit s použitím relaxační metody. Proto musíme transformovat systém dvou obyčejných
diferenciálních rovnic 2. řádu na systém čtyř obyčejných diferenciálních rovnic 1. řádu. Nejjednodušší
Kapitola 6 Řešení okrajové úlohy relaxační metodou 55
způsob jak to provést, je zavést rychlosti u = ˙x a v = ˙y jako nové závislé proměnné. Pak dostaneme
d
dt






x
u
y
v






=






u
0
v
−g






. (6.14)
Dále zavedeme vzorkování časového intervalu [0, tb]. Použijeme rovnoměrně rozloženou síť bodů, takže
vektor hodnot nezávislé proměnné je T = [0, . . . , ti, . . . , tM−1 = tb]. Vektor lokálního řešení je jednoduše
Y = [x, u, y, v]T . Vektor celkového řešení má tvar
Y = [x0, u0, y0, v0, . . . , xi, ui, yi, vi, . . . xM−1, uM−1, yM−1, vM−1] .
Funkce F, A a B popisující pravé strany rovnic a okrajové podmínky jsou
F (t, [x, u, y, v]) = [u, 0, v, −g], A ([x0, u0, y0, v0]) = [x0, y0] = 0,
B ([xM−1, uM−1, yM−1, vM−1]) = [xM−1, yM−1] = 0,
a odpovídající Jakobiány:
∂J
∂Y
=






0 1 0 0
0 0 0 0
0 0 0 1
0 0 0 0






,
∂A
∂Y
=
∂B
∂Y
=
1 0 0 0
0 0 1 0
Pro použití relaxační metody musíme poskytnout počáteční odhad, který potenciálně povede ke skutečnému
řešení. Počáteční odhad není pravým řešením diferenciálních rovnic, ale obvykle stačí, když splňuje základní
topologické vlastnosti předpokládaného řešení (například, že spojitě spojuje počáteční a koncové okrajové
podmínky). My použijeme obzvláště jednoduchý první odhad – přímku. Pak předpokládáme, že se banán
hýbe s konstantní rychlostí horizontálním směrem a dosáhne cíle v čase tb. Tato situace odpovídá stavu
bez gravitačního pole. Horizontální rychlost je v tomto případě dána jako u = d/tb, vertikální rychlost je
nulová a trajektorie je popsána x(t) = td/tb a y = 0. Řešení je implementováno dále.
1 class GorillasProblem:
2
3 def __init__(self, d, g, tb, M):
4 """
5 Initialization.
6 `d` is the distance of the target,
7 `g` is the gravitational acceleration
8 `tb` is the time, when the should reach the target,
9 `M` is the number of mesh points in time domain.
10 """
11 self._d = d
12 self._g = g
13 # create evenly spaces sample for the independent variable (time):
14 self._T = np.linspace(0, tb, M)
15 # put as the solution vector the initial guess:
16 self._Y = self.initial_guess()
Kapitola 6 Řešení okrajové úlohy relaxační metodou 56
17 self._norm_dY = None
18
19 def initial_guess(self):
20 """ Return the global solution vector Y corresponding
21 to the initial guess. """
22 u0 = self._d/self._T[-1] # ...constant horizontal speed
23 Y = [0, u0, 0, 0]
24 for t in self._T[1:]:
25 Y += [u0*t, u0, 0, 0]
26 return np.array(Y)
27
28 def F(self, x, Y):
29 """ The right-hand sides of the differential equations. """
30 [x, u, y, v] = Y
31 return np.array([u, 0, v, -self._g])
32
33 def jacF(self, x, Y):
34 """ The Jacobian of the right-hand side. """
35 return np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]])
36
37 def A(self, Y):
38 """ The boundary constrains at t = 0: x=0 and y=0 """
39 [x, u, y, v] = Y
40 # We return vector of just 2 items as there are only
41 # two constraints: x = 0 and y = 0.
42 return np.array([x, y])
43
44 def jacA(self, Y):
45 """ Jacobian of the boundary constrains at t=0. """
46 return np.array([[1, 0, 0, 0], [0, 0, 1, 0]])
47
48 def B(self, Y):
49 """ The boundary constrains at t = tb: z=d and y=0 """
50 [x, u, y, v] = Y
51 return np.array([x - self._d, y])
52
53 def execute_one_step(self):
54 """ Executes the relaxation method. """
55 dY = relaxation_step(self._T, self._Y, self.F, self.jacF, self.A, self.jacA, self.B, self.jacA)
56 self._norm_dY = np.linalg.norm(dY)
57 self._Y += dY
58
59 def norm_dY(self):
60 """ Returns the L2-norm of the last update. """
61 return self._norm_dY
62
63 def initial_velocity(self):
64 """ Returns the initial velocity of the current solution and the angle. """
65 u = self._Y[1]
66 v = self._Y[3]
67 return math.hypot(u, v), math.atan2(v, u)*180/math.pi
68
69 def plot_solution_xy(self, ax, **kwargs):
Kapitola 6 Řešení okrajové úlohy relaxační metodou 57
70 """ Plots the solution into the axes object `ax` """
71 ax.plot(self._Y[0::4], self._Y[2::4], **kwargs)
A jeden příklad pro parametry d = 10 m, g = 9.8 m · s−2 a tb = 5 s.
1 gorillas = GorillasProblem(10, 9.8, 5, 20)
2 fig, ax = plt.subplots(1)
3 gorillas.plot_solution_xy(ax, ls='-', marker='x', label='initial guess')
4
5 gorillas.execute_one_step()
6 gorillas.plot_solution_xy(ax, ls='-', marker='x', label='1-st step')
7 print('1-st step norm = {}'.format(gorillas.norm_dY()))
8
9 gorillas.execute_one_step()
10 gorillas.plot_solution_xy(ax, ls='-', marker='x', label='2-nd step')
11 print('2-nd step norm = {}'.format(gorillas.norm_dY()))
12
13 ax.set_xlabel(r'$x [\mathrm{m}]$')
14 ax.set_ylabel(r'$y [\mathrm{m}]$')
15 ax.legend()
16 print('initial velocity = {} m/s, angle = {} deg'.format(*gorillas.initial_velocity()))
V tomto případě konvergovala relaxační metoda hned po prvním kroku. Je to díky tomu, že systém
obyčejných diferenciálních rovnic je homogenní a lineární. Newton-Raphsonova metoda aproximující funkci
G(Y) tak dává přímo přesné řešení. To však není obvyklý případ, jak uvidíte při řešení následujícího
problému.
6.3 Numerické nebo automatické derivace
Vypočítat Jakobián je jednoduché, avšak někdy to může být zdlouhavé. Místo přímého zapsání vzorce
můžeme použít numerickou diferenciaci:
Kapitola 6 Řešení okrajové úlohy relaxační metodou 58
Zavedeme numerickou derivaci funkce jedné proměnné (numder) a gradient ∂kF pro funkci více proměnných.
F může být také vektor nebo tenzor.
1 def numder(f, x, dx=1e-10):
2 return (f(x+dx) - f(x-dx))/(2*dx)
3
4 def grad(F, X, eps=1e-5):
5 n = len(X)
6 result = []
7 for i in range(n):
8 FF = lambda t: F(np.concatenate((X[:i], [t], X[i+1:]), axis=None))
9 result.append(numder(FF, X[i], dx=eps))
10 return np.array(result)
C Úkol
Zkuste znova vyřešit Gorilla problém, avšak bez zanedbání odporu vzduchu. Navíc uvažujte nenulovou
rychlost větru a různou výšku střelce a terče.
C
Souhrn kapitoly
V této kapitole jsme se seznámili s metodou řešení soustavy obyčejných diferenciálních rovnic pomocí
metody relaxace, také jsme si zopakovali použití Newton-Raphsonovy metody. Dále poznatky využijeme
v následující kapitole.
? Otázky
• Kolik okrajových podmínek musíme znát, abychom mohli řešit systém N obyčejných diferenciálních
rovnic 1. řádu?
• Jak se nazývá úloha, která má všechny okrajové podmínky na jedné straně oblasti řešení?
• Jakým způsobem je možné vyřešit problém, který je zadán soustavou obyčejných diferenciálních
rovnic vyššího řádu?
?
Kapitola 7
Bondiho akrece
Jiří Horák
Rychlý náhled: V této kapitole budeme vytvářet numerický model Bondiho akrece v jazyce Python.
Využijeme k tomu relaxační metodu popsanou v předcházející kapitole. Řešení budeme implementovat pro
síť souřadnic s logaritmickým rozložením.
Klíčová slova: Bondiho akrece, numerické modelování, numerické řešení diferenciálních rovnic, relaxační
metoda.
ºº Cíle studia: V této kapitole najdete konkrétní implementaci a použití relaxační metody na řešení
reálného fyzikálního problému, konkrétně Bondiho akrece. Také se naučíte definovat a používat různé
rozložení sítě souřadnic.
7.1 Akreční procesy
Akrece je proces, při němž se hmota akumuluje na povrchu hmotného tělesa díky gravitační přitažlivé síle.
K akreci dochází nejčastěji ve vesmíru v okolí velmi hmotných těles, jako jsou hvězdy, případně i hmotné
planety nebo extrémní objekty (černé díry, neutronové hvězdy). V případě hvězd nebo planet dochází
k akreci prachových částic, v případě kompaktních objektů se jedná o plasma, díky extrémním tlakům
a teplotám, které v jejich okolí nalezneme.
Specifický případ, který je zároveň nejjednodušší, je případ stacionární sféricky symetrické (Bondiho)
akrece.
7.2 Sférická akrece
Uvažujeme nerelativistickou, sféricky symetrickou akreci polytropické kapaliny na kompaktní objekt.
 Definice
Polytropická kapalina je popsána polytropickou stavovou rovnicí
p = Kργ
, (7.1)
kde ρ a p je hustota a tlak plynu, K je polytropická konstanta a γ je polytropický index.
59
Kapitola 7 Bondiho akrece 60

Polytropická kapalina je jednoduchý model kapaliny, pro který různé hodnoty polytropického indexu odpovídají
různým druhům kapalin. Rychlost a, kterou se v kapalině šíří zvuk, je dána
a2
=
dp
dρ
= γKργ−1
=
γp
ρ
. (7.2)
 Definice
Gravitační pole kompaktního objektu je popsáno Newtonovským gravitačním potenciálem
Φ(r) = −
GM
r
, (7.3)
kde M je hmotnost kompaktního objektu, G je gravitační konstanta a r je radiální souřadnice.

Předpokládáme, že ve velké vzdálenosti od centra je tok v klidu, takže radiální rychlost ur vymizí. Asymptotické
hodnoty hustoty a tlaku označíme ρ∞ a p∞.
Dynamika kapaliny je dána rovnicí kontinuity popisující zachování hmoty, a radiální Eulerovou rovnicí
popisující zachování hybnosti v radiálním směru. Rovnice kontinuity má tvar
· (ρv) =
1
r2
d
dr
r2
ρur
= 0 (7.4)
V tomto případě může být dále integrována
r2
ρur
= const ≡ −
˙M
4π
. (7.5)
Význam konstanty ˙M je jasný: udává, jak rychle je hmota dodávána centrálnímu objektu. Rovnice (7.5)
udává pouze celkové množství hmoty, které překročí sféru o poloměru r za jednotku času, ˙M = 4πr2ρ(−ur)
je konstantní bez ohledu na velikost sféry. Tato podmínka musí být splněna, jinak by došlo k akumulaci
hmoty v některých místech toku a proces by nebyl stacionární.
Radiální Eulerova rovnice popisuje pohyb elementu kapaliny pod vlivem tlakové a gravitační síly
ur dur
dr
= −
1
ρ
dp
dr
−
GM
r2
= −a2 d ln ρ
dr
−
GM
r2
. (7.6)
Použijeme–li logaritmus radiální derivace rovnice kontinuity, dostaneme
d ln ρ
dr
= −
2
r
−
1
ur
dur
dr
(7.7)
a radiální Eulerova rovnice má pak tvar
d(ur)2
dr
=
4a2
r − 2GM
r2
1 − a2
(ur)2
. (7.8)
7.2.1 Bezrozměrná verze
Definujeme-li následující bezrozměrné proměnné:
y ≡
−ur
a∞
2
, ˜a ≡
a
a∞
, x ≡
r
r∗
, ˙m ≡
˙M
4πr2
∗ρ∞a∞
, (7.9)
Kapitola 7 Bondiho akrece 61
kde
r∗ ≡
GM
a2
∞
(7.10)
je přirozená délková, je možné celou situaci popsat jednou obyčejnou diferenciální rovnicí (Bondiho rovnice)
a jednou algebraickou rovnicí
dy
dx
=
4˜a2
x − 2
x2
1 − ˜a2
y
, ˜a2
=
˙m
x2y1/2
γ−1
. (7.11)
Okrajové podmínky při x → ∞ vyplývají z požadavku nulové radiální rychlosti ve velké vzdálenosti,
(y → 0), rychlost zvuku pak nabývá asymptotické hodnoty ˜a → 1. Druhá podmínka dává asymptotické
chování
y ≈
˙m
x2
2
(x → ∞). (7.12)
7.2.2 Kritický bod a podmínka regularity
Bezrozměrná Bondiho rovnice má kritický bod, ve kterém rychlost proudění nabývá lokální rychlosti zvuku,
y = ˜a2. V tom případě je jmenovatel na pravé straně rovnice (7.11) nulový a řešení rovnice je může být
regulární pouze pokud je ve stejném bodě nulový také čitatel. Takové místo nazýváme jako sonický bod a
příslušnou radiální souřadnici označíme xs. Aby bylo řešení regulární, musí tedy platit
4˜a2
xs
−
2
x2
s
= 0 ⇔ y = ˜a2
. (7.13)
Z toho vyplývá další podmínka na řešení Bondiho rovnice. Obecně platí, že integrační křivka y(x), která
prochází spojitě sonickým bodem (a splňuje podmínku regularity určenou výše), nemusí být v souladu
s požadovaným asymptotickým chováním řešením při x → ∞. Obě tyto podmínky jsou tak splněny jen při
určitých hodnotách akreční rychlosti ˙m, která hraje roli vlastní hodnoty problému. Stacionární akrece je
možná jen pro tuto hodnotu akrečního toku, při jiných hodnotách nemusí být akrece nutně stacionární.
7.3 Analytické řešení
Další podmínka regularity představuje zajímavou výzvu při numerickém řešení. Dříve než se o něj pokusíme,
nastíníme si analytické řešení, abychom měli něco k porovnání. Hlavní výhoda sférické akrece je, že je možné
integrovat také rovnici pro hybnost, takže se jedná o algebraický problém (přestože je dán transcendentní
rovnicí).
Začněme s bezrozměrnou rovnicí pro hybnost
y1/2 dy1/2
dx
+
1
ρa2
∞
dp
dx
+
1
x2
= 0. (7.14)
Je snadné integrovat první a třetí člen podle x. Druhý člen lze integrovat následovně:
1
a2
∞
1
ρ
dp
dx
dx =
1
a2
∞
dp
ρ
=
γK
a2
∞
ργ−2
dρ =
γKργ−1
(γ − 1)a2
∞
=
˜a2
γ − 1
. (7.15)
Pak dostaneme
y
2
+
˜a2
γ − 1
−
1
x
= const ≡ Ψ, (7.16)
Kapitola 7 Bondiho akrece 62
což je bezrozměrná Bernoulliho rovnice. Integrační konstanta Ψ označuje odlišné integrační křivky, které
jsou řešením y(x) Bernoulliho rovnice (7.11). Ta z nich, která odpovídá naší okrajové podmínce, musí
splňovat y → 0 a ˜a → 1 pro x → ∞, je tedy dána hodnotou
Ψ ≡ Ψ0 = 1/(γ − 1). (7.17)
Které hodnotě Ψ odpovídají řešení, které spojitě prochází kritickým bodem? Z podmínky regularity víme,
že
ys = ˜a2
s , xs =
1
2˜a2
s
. (7.18)
Rovnice pro zachování hmoty nám dále dovoluje vyjádřit rychlost zvuku v sonickém bodě pomocí akreční
rychlosti jako
˜as = (4 ˙m)
γ−1
5−3γ . (7.19)
Po dosazení do Bernoulliho rovnice dostáváme
Ψ ≡ Ψs =
˜a2
s
2
+
˜a2
s
γ − 1
− 2˜a2
s =
5 − 3γ
2γ − 2
˜a2
s =
5 − 3γ
2γ − 2
(4 ˙m)
2γ−2
5−3γ . (7.20)
Jelikož požadujeme, aby řešení procházelo spojitě sonickým bodem, a zároveň dávalo správné chování na
velkých poloměrech, musí platit Ψs = Ψ0. To nám dává jednoduché analytické vyjádření polohy sonického
bodu xs a lokální rychlost zvuku ˜as v sonickém bodě,
xs =
5 − 3γ
4
, ˜a2
s =
2
5 − 3γ
. (7.21)
Tato rychlost zvuku odpovídá akreční rychlosti ˙m0
˙m0 =
1
4
2
5 − 3γ
5−3γ
2γ−2
, (7.22)
což je vlastní hodnota problému.
7.3.1 Asymptotické chování řešení:
Z předchozích odstavců už víme, jak se má řešení chovat ve velkých vzdálenostech od kompaktního objektu
(když x → ∞). Jaké chování lze očekávat, když x → 0? K nalezení odpovědi se skvěle hodí Frobeniova teorie.
Předpokládejme, že řešení blízko bodu x = 0 vystihuje závislost y ≈ A/xα. Z Bernoulliho rovnice vybereme
nejrychleji rostoucí členy, které musí být v rovnováze a jejich porovnáním určíme α a A. Z rovnice kontinuity
je vidět, že ˜a2 ∼ x(α−2)(γ−1). První člen Bernoulliho rovnice (člen odpovídající kinetické energii) roste jako
x−α, druhý člen (odpovídající vnitřní energii) roste jako x(α−2)(γ−1) a poslední třetí člen (odpovídající
potenciální energii) roste jako x−1. Porovnáním exponentů za předpokladu že 1 ≤ γ ≤ 5/3 zjistíme,
že nejvyšší hodnota α odpovídá rovnováze mezi kinetickým a potencionálním členem a odpovídá α = 1
a A = 2, takže
y ≈
2
x
(x → 0). (7.23)
Odsud vidíme, že velmi blízko kompaktního objektu má největší vliv gravitace a tlak způsobuje jen zanedbatelné
perturbace, zatímco ve velkých vzdálenostech je nejvýraznější rovnováha mezi gravitační potenciální
a vnitřní energií. V obou případech má navíc zanedbávaný člen tendenci snižovat hodnoty rychlosti. To se
bude hodit při numerickém řešení Bernoulliho rovnice pomocí metody půlení intervalů.
Kapitola 7 Bondiho akrece 63
7.3.2 Grafické znázornění analytických výsledků
Na obrázku 7.1 je znázorněno analytické řešení popsáno dříve. Abychom graf mohli vytvořit, začneme
nahráním potřebných modulů:
1 import numpy as np
2 import math
3 import scipy.optimize as opt
4 import matplotlib
5 import matplotlib.pyplot as plt
6
7 # in order to use LaTex fonts in the graphs:
8 from matplotlib import rc
9 rc('text', usetex=True)
10 rc('legend', frameon=False)
11
12 # introduce "nice" colors:
13 purple = '#9467bd'
14 blue = '#1f77b4'
V Pythonu dále definujeme funkce pro různé fyzikální veličiny popisující akreční tok. V obecném případě
je akreční tok popsán dvěma parametry – polytropickým indexem γ a bezrozměrnou akreční rychlostí ˙m.
Regulární řešení však existuje jen pro jednu konkrétní hodnotu ˙m0, která je implementována jako funkce
mdot0. Bernoulliho konstanta Ψ, jako funkce bezrozměrného poloměru x a kvadratické rychlosti y, je
implementována pomocí funkce Psi. Její dvě speciální hodnoty Ψs a Ψ0, odpovídající integrální křivce
procházející sonickým bodem a řešení se správným chováním na velkých vzdálenostech, jsou dány funkcemi
Psis a Psi0. Asymptotická chování řešení blízko kompaktního objektu a v nekonečnu popisují funkce
yasy_0 a yasy_inf, poloha kritického bodu [xs(γ, ˙m), ys(γ, ˙m)] je dána funkcemi xs a ys.
1 def mdot0(gamma):
2 return (2/(5-3*gamma))**((5-3*gamma)/(2*gamma-2))/4.
3 def Psi(x, y, gamma, mdot):
4 return y/2. - 1./x + (mdot/(x*x*math.sqrt(y)))**(gamma-1)/(gamma-1)
5 def Psis(gamma, mdot):
6 return (5-3*gamma)/(2*gamma-2)*(4.*mdot)**((2*gamma-2)/(5-3*gamma))
7 def Psi0(gamma):
8 return 1/(gamma-1)
9 def yasy_inf(x, mdot):
10 return (mdot/(x*x))**2
11 def yasy_0(x):
12 return 2/x
13 def xs(gamma, mdot):
14 return 1/(2*ys(gamma, mdot))
15 def ys(gamma, mdot):
16 return (4*mdot)**((2*gamma-2)/(5-3*gamma))
Nyní vykreslíme grafické řešení pro různé hodnoty akrečního toku. Místo řešení Bernoulliho rovnice,
abychom získali specifické funkce y = y(x; Ψ), vykreslíme čáry konstantního Ψ v rovině (x, y). Ukážeme si
tři případy. První dva jsou nefyzikální, s ˙m < ˙m0 odpovídající Ψ0 < Ψs a ˙m > ˙m0, což odpovídá Ψ0 < Ψs.
Nakonec třetí případ, kdy ˙m = ˙m0, dává fyzikální smysl.
Kapitola 7 Bondiho akrece 64
1 def flow_plot(ax, gamma, mdot):
2 xvals = np.logspace(-2, 2, 300)
3 yvals = np.logspace(-5, 3, 300)
4 n = len(xvals)
5 X, Y = np.meshgrid(xvals, yvals)
6 Z = np.array([[Psi(x, y, gamma, mdot) for x in xvals] for y in yvals])
7 ax.contour(X, Y, Z, levels=[Psis(gamma, mdot)], colors=[blue])
8 ax.contour(X, Y, Z, levels=[Psi0(gamma)], colors=[purple])
9 ax.plot(xvals[2*n//3:], [yasy_inf(x, mdot) for x in xvals[2*n//3:]] ls=':', color=purple)
10 ax.plot(xvals[:n//5], [yasy_0(x) for x in xvals[:n//5]], ls=':', color=purple)
11 ax.set_xlabel(r'Radial coordinate $x$')
12 ax.set_ylabel(r'Squared velocity $y$')
13 ax.axhline(ys(gamma, mdot), ls='--', color='gray')
14 ax.axvline(xs(gamma, mdot), ls='--', color='gray')
15 ax.set_ylim(1e-5,1e+3)
16 ax.set_xscale('log')
17 ax.set_yscale('log')
18
19 fig, panels = plt.subplots(1, 3, figsize=(14,4))
20 flow_plot(panels[0], 4./3., 0.5*mdot0(4/3))
21 flow_plot(panels[1], 4./3., 1.5*mdot0(4/3))
22 flow_plot(panels[2], 4./3., mdot0(4/3))
23 panels[0].annotate(r'$\dot{m}<\dot{m}_0$', xy=(0.98, 0.98), xycoords="axes fraction", va="top", ha='right')
24 panels[1].annotate(r'$\dot{m}>\dot{m}_0$', xy=(0.98, 0.98), xycoords="axes fraction", va="top", ha='right')
25 panels[2].annotate(r'$\dot{m}=\dot{m}_0$', xy=(0.98, 0.98), xycoords="axes fraction", va="top", ha='right')
26
27 fig.tight_layout()
Obrázek 7.1: Tři případy pro různé hodnoty akrečního toku.
Je vhodné zavést funkci, která vrací transonické řešení y v daném bodě x – tady je: Hodnota kvadrátu
rychlosti je nalezena jako řešení Bernoulliho rovnice pro dané x. Použijeme metodu bisekce. Jako okrajové
hodnoty použijeme rychlost zvuku v sonickém bodě a dvě asymptotické formy (pro x → 0, ∞) vynásobené
„bezpečnostním“ faktorem
1 def exact_solution(x, gamma):
2 mdot = mdot0(gamma)
3 x0 = xs(gamma, mdot)
Kapitola 7 Bondiho akrece 65
4 y0 = ys(gamma, mdot)
5 if x < x0:
6 ya, yb = y0, max(yasy_0(x), 1e2*y0)
7 return opt.bisect(lambda y: Psi(x, y, gamma, mdot)-Psi0(gamma), ya, yb)
8 else:
9 ya, yb = min(1e-2*yasy_inf(x, mdot), 1e-2*y0), y0
10 return opt.bisect(lambda y: Psi(x, y, gamma, mdot)-Psi0(gamma), ya, yb)
11
12 fig, ax = plt.subplots(1)
13 flow_plot(ax, 4/3, mdot0(4/3))
14 xvals = np.logspace(-2, 2)
15 ax.plot(xvals, [exact_solution(x, 4/3) for x in xvals], ls='--', color=blue);
Obrázek 7.2: Přesné řešení pro kvadrát rychlosti.
Kapitola 7 Bondiho akrece 66
7.4 Řešení s použitím relaxační metody, řešení problému vlastních hodnot
na neznámé oblasti řešení
Chceme řešit Bondiho úlohu numericky. Nejprve si pojďme shrnout základní vlastnosti úlohy.
• Bondiho akrece je „okrajová úloha“. Řešení je vymezeno na velkých radiálních souřadnicích a navíc je
další omezení dáno podmínkou regularity v sonickém bodě. Numerické řešení musí splňovat omezení
na obou stranách oblasti řešení. Jak jsme viděli, relaxační metoda může být v tomto případě velmi
nápomocná.
• V případě Bondiho akrece známe předem pozici sonického bodu, ale v obecném případě nemusí být
jeho pozice známá. Abychom mohli lépe demonstrovat numerické řešení, budeme dále předstírat,
že nevíme kde se sonický bod nachází. To však dále komplikuje situaci. Pokud bychom považovali
bezrozměrnou radiální souřadnici za nezávislou proměnou, neznali bychom předem oblast řešení (nevěděli
bychom, kde se nachází vnitřní okraj). Musíme tedy použít nějaký trik a tento problém vyřešit
dříve než začneme přímo implementovat relaxační metodu.
• Protože je řešení dáno jednoduchou diferenciální rovnicí prvního řádu, je problém očividně příliš
omezený a řešení existuje pouze pro určité hodnoty akrečního toku ˙m. To znamená, že Bondiho úloha
je také problém vlastních hodnot. Hledáme konkrétní hodnotu parametru (bezrozměrného akrečního
toku ˙m) tak, že existuje řešení kompatibilní s oběma okrajovými podmínkami. Jak můžeme použít
relaxační metodu pro nalezení vlastních hodnot problému?
7.4.1 Shrnutí problému
Pojďme si nejprve shrnout problém v jeho matematické formě. Hledáme řešení y(x) a hodnotu parametru ˙m
problému zadaného jednou obyčejnou diferenciální rovnicí,
dy
dx
=
4˜a2
x − 2
x2
1 − ˜a2
y
, ˜a2
=
˙m
x2y1/2
γ−1
, (7.24)
s dvěma okrajovými omezeními
y ≈
˙m
x2
2
(x → ∞) (7.25)
a
4˜a2
xs
−
2
x2
s
= 0 ⇔ y = ˜a2
. (7.26)
7.4.1.1 Úprava problému vlastních hodnot
Máme dvě okrajové podmínky pro jednu obyčejnou diferenciální rovnici, které závisí na jednom dodatečném
parametru ( ˙m). Hledáme řešení y(x) a odpovídající hodnotu ˙m takovou, že obě okrajové podmínky jsou
splněny. Nejjednodušší způsob, jak najít řešení, je považovat parametr ˙m za další závislou proměnou,
˙m = ˙m(x), danou další obyčejnou diferenciální rovnicí d ˙m/dx = 0 zajišťující, že se bude jednat o konstantu.
Pak dostaneme soustavu dvou obyčejných diferenciálních rovnic pro y(x) a ˙m(x) s dvěma okrajovými
omezeními,
Kapitola 7 Bondiho akrece 67
d
dx
˙m
y
=


0
4˜a2
x − 2
x2 1 − ˜a2
y
−1

 , A ( ˙m, y) =
4y
xa
−
2
x2
a
= 0, B ( ˙m, y) = y −
˙m2
x2
b
= 0. (7.27)
Pozice vnějšího okraje xb je ve své podstatě volná, v praxi bychom měli volit rádius dostatečně daleko,
tak aby skutečné řešení souhlasilo s asymptotickým tvarem, který jsme odvodili v předchozí části.
Podíváme-li se na analytické řešení, vidíme, že už na xb ≈ 102 odpovídá skutečné řešení asymptotickému
tvaru velmi dobře (alternativně můžeme zkusit najít další člen v asymptotickém rozvoji, který indikuje
chybu předešlého řádu, a podle toho se rozhodnout, kam umístíme vnější okrajovou podmínku). Na druhou
stranu, poloha vnitřní okrajové podmínky je stále neznámá – místo toho máme vztah
y(xa) − ˜a2
( ˙m, xa, y) = 0. (7.28)
To je mnohem lepší, ale stále ne perfektní. Pozice vnitřní okrajové podmínky xa není předem známá, místo
toho závisí na lokální hodnotě závislých proměnných y a ˙m.
7.4.1.2 Co udělat s vnitřní okrajovou podmínkou?
Vnitřní okrajová podmínka je očividně také neznámá, ale máme pro ní jedno omezení navíc. Co kdybychom
do závislých proměnných zahrnuli také radiální souřadnici? Jinak řečeno, zavedeme novou závislou proměnnou
q, tak že x = x(q), y(x) ≡ y[x(q)] = y(q), a obdobně pro ˙m. Pak lze derivaci podle x rozložit na
derivaci složené funkce
d
dx
≡
dq
dx
d
dq
(7.29)
Pojďme se podívat, co se stane s naším systémem rovnic. Nyní máme tři neznámé funkce ˙m(q), x(q)
a y(q) a systém obyčejných diferenciálních rovnici vypadá takto:
d
dq



˙m
x
y


 =




0
ψ
ψ 4˜a2
x − 2
x2 1 − ˜a2
y
−1



 . (7.30)
Objevuje se zde nová funkce ψ(q) = dx/dq, která vychází ze změny nezávislé proměnné x na q. Tato funkce
je víceméně volitelná, pokud zvolíme ψ ≡ 1, dostaneme původní problém. Pokud jí zvolíme rovnou nějaké
jiné konstantě (např. α), dostaneme jen jiné škálování radiální souřadnice (q = x/α). Na první pohled se
zdá, že jsme nic nevylepšili: měřit radiální souřadnici v x nebo v násobcích x – jaký je v tom rozdíl? Avšak to
je právě smysl této práce. Všimněte si, že stále nevíme kde se nachází vnitřní okrajová podmínka xa, avšak,
jakákoli je její hodnota, vždy bude existovat taková hodnota ψ, která bude transformovat (neznámý) interval
x, [xa, xb] do daného předepsaného intervalu q, třeba [0, 1]. Zaveďme x(0) = xa a x(1) = xb. Hodnota ψ,
která definuje tuto transformaci je neznámá, ale můžeme jí považovat za další vlastní podmínku problému,
úplně stejně jako jsme to udělali s ˙m. Takže, v konečném důsledku řešíme problém pro čtyři závislé proměnné
˙m, x, y a ψ danými rovnicemi
d
dq






˙m
x
y
ψ






=






0
ψ
ψ 4˜a2
x − 2
x2 1 − ˜a2
y
−1
0






, ˜a2
( ˙m, x, y) =
˙m
x2y1/2
γ−1
. (7.31)
Kapitola 7 Bondiho akrece 68
s dvěma omezeními řešení na obou stranách,
A [ ˙m, x, y, ψ]T
=
4˜a2
x − 2
x2
y − ˜a2
= 0, B [ ˙m, x, y, ψ]T
=
x − xb
y − ˙m2
x4
= 0 (7.32)
V této formě je již možné problém řešit použitím relaxační metody. Avšak než to uděláme, musíme
přidat několik poznámek.
7.4.1.3 Obecná poznámka k síti souřadnic
Transformace z neznámého intervalu [xa, xb] na pevný interval v q:
V předchozím textu jsme předpokládali, že [xa, xb] → [0, 1], takže x(q = 0) = xa a x(q = 1) = xb. Ve
skutečnosti je mnohem obvyklejší transformace [xa, xb] → [0, M − 1], kde M je počet bodů sítě. Pak
jednotlivé body sítě odpovídají celočíselným hodnotám nezávislé proměnné q (např. xi ≡ x(q = i)). Dosud
jsme předpokládali lineární transformaci q → x takovou, že x(0) = xa a x(M −1) = xb. Takže transformace
má tvar
x = xa +
q
M − 1
(xb − xa) , ψ =
dx
dq
=
xb − xa
M − 1
(7.33)
Můžeme ale také zvolit nelineární funkci. Například pokud vezmeme místo dx/dq = ψ,
dx
dq
= ψθ(x), (7.34)
kde ψ je konstanta a θ(x) je funkcí x. Po integraci pak dostaneme
dx
θ(x)
= ψq + const (7.35)
Funkce θ očividně ovlivní transformaci q → x. Jelikož celočíslené hodnoty q označují jednotlivé body sítě,
θ(x) určuje umístění těchto bodů. Například volba θ(x) = x definuje logaritmické rozložení bodů sítě, které
bude využito při implementaci dále.
 Můžeme použít i pokročilejší rozložení bodů, například metodu adaptivní sítě. Například
θ(x, Y ) = ∆ + βY
dY
dx
−1
(7.36)
nám dává uniformní síť, kde je vzdálenost mezi body rovná ∆ právě tehdy, když β = 0. Pro nenulové β
je rozložení bodů ovlivněno také magnitudou funkční derivace: síť je hustší v místech, kde se řešení mění
rychleji.
7.5 Implementace
Při implementaci využijeme relaxační metodu pro řešení Bondiho akrece. Použijeme 4 proměnné [ ˙m, x, y, ψ]T .
Dvě z nich, ˙m a ψ jsou vlastní hodnoty problému. Pro transformaci q → x použijeme dx/dq = xψ, tedy
logaritmické rozložení bodů sítě. Hodnota q = 0 odpovídá vnějšímu okraji umístěnému na x = 102. Vnitřní
okraj je dán funkcí q = M − 1, kde M je počet bodů sítě. Voláme relaxační metodu implementovanou
v relaxation.py (viz předchozí kapitola). Jako počáteční odhad použijeme dvě vlastní hodnoty ˙m a ψ.
Poté integrujeme obyčejné diferenciální rovnice použitím metody Runge-Kutha až do bodu, kdy rychlost
dosáhne 0.8 lokální rychlosti zvuku. S tímto řešením můžeme začít následný relaxační proces.
Kapitola 7 Bondiho akrece 69
Jelikož relaxační metoda vyžaduje implementaci jakobiánu pravé strany rovnic a hraniční omezení,
musíme buďto najít přesný vztah analyticky, nebo jej vypočítat numericky – tak, jak to uděláme dále. Dvě
funkce numder a numjac berou funkci jedné nebo více proměnných a vrací funkci, která počítá numerickou
derivaci metodou centrální diference. Elegantnější řešení by však bylo využít automatickou diferenciaci.
1 def numder(f, dx=1e-10):
2 return lambda x: (f(x+dx) - f(x-dx))/(2*dx)
3
4 def numjac(F, eps=1e-5):
5 def jac(X):
6 n = len(X)
7 result = []
8 for i in range(n):
9 FF = lambda t: F(np.concatenate((X[:i], [t], X[i+1:]), axis=None))
10 result.append((FF(X[i]+eps) - FF(X[i]-eps))/(2*eps))
11 return np.transpose(np.array(result))
12 return jac
13
14 # example: der is derivative of the function x --> x^2
15 der = numder(lambda x: x**2)
16 print(der(1))
17
18 # example: jac is jacobian of the function [x, y] --> [x^2, x+y]
19 jac = numjac(lambda X: np.array([X[0]**2, X[0]+X[1]]))
20 print(jac([1,1]))
Zde je implementace:
1 import scipy.integrate as integrate
2 import relaxation
3
4 class BondiFlow:
5 """
6 Solution class for finding regular solutions of the Bondi equation
7
8 Initialization:
9 ---------------
10 gamma : polytropic index of the flow
11 mdot : initial guess of the accretion rate
12 xinf : position of the outer boundary
13 M : rough number of the grid points (logarithmically spaced)
14 """
15
16 def __init__(self, gamma, mdot, xinf, M):
17 """ Initialization. """
18 # set the private variables:
19 self._gamma = gamma
20 self._mdot = mdot
21 self._xinf = xinf
22 # set the Jacobians of RHS and boundary constraints:
23 # Note since F depends on both dependent and independent variables (q and X),
Kapitola 7 Bondiho akrece 70
24 # and Jac is related only to X, we have to do the following trick:
25 self.jacF = lambda q, X: numjac(lambda Y: self.F(0, Y))(X)
26 self.jacA = numjac(self.A)
27 self.jacB = numjac(self.B)
28 # finally, calculate the initial guess:
29 self.make_initial_guess(M)
30
31 def make_initial_guess(self, M):
32 """ Calculating the initial guess using direct RK integration. """
33 # initial values of the solution variables = first 4 entry of the solution vector _Ys
34 # (here psi0 is calculated by integrating dx/dq equation and assuming that the
35 # sonic point is approximately at xs ~ 1/2)
36 xs = 0.5
37 psi0 = math.log(xs/self._xinf)/(M-1)
38 self._Ys = np.array([self._mdot, self._xinf, math.log(self._mdot**2/self._xinf**4), psi0])
39 self._qs = [0]
40
41 # introduce the scipy.integrate ode solver:
42 solver = integrate.ode(self.F).set_integrator('vode', method='bdf', rtol=1e-10, nsteps=10000)
43 solver.set_initial_value(self._Ys, self._qs[0])
44
45 # integrate until we the velocity is greater than sqrt(0.8) sound speed or number of
46 # mesh point exceeds 1.5 times the initially setted. The local values are storred in
47 # the local solution vectors _qs and _Ys:
48 q = 1
49 finish = False
50 while not(finish):
51 Y = solver.integrate(q)
52 self._qs.append(q)
53 self._Ys = np.concatenate((self._Ys, Y))
54 q += 1
55 [mdot, x, y, psi] = Y
56 a2 = (mdot/(x*x*math.exp(y/2)))**(self._gamma - 1)
57 finish = (q > 1.5*M) or (math.exp(y) > 0.8*a2)
58
59 def F(self, q, Y):
60 """ The right-hand sides of the ODEs. """
61 [mdot, x, y, psi] = Y
62 a2 = (mdot/(x*x*math.exp(y/2)))**(self._gamma - 1)
63 f = (4*a2/x - 2/x**2)/(math.exp(y) - a2)
64 return np.array([0, x*psi, f*x*psi, 0], float)
65
66 def A(self, Y):
67 """ Constraints at x_inf. (corresponding to q=0) """
68 [mdot, x, y, psi] = Y
69 xinf = 1e2
70 return np.array([x - xinf, math.exp(y) - mdot**2/x**4], float)
71
72 def B(self, Y):
73 """ Constraints at the sonic point (corresponding to q=M-1). """
74 [mdot, x, y, psi] = Y
75 a2 = (mdot/(x*x*math.exp(y/2)))**(self._gamma - 1)
76 return np.array([math.exp(y) - a2, 4*a2/x - 2/x**2], float)
Kapitola 7 Bondiho akrece 71
77
78 def calculate(self, k=1., atol=1e-3, maxiter=50):
79 """
80 Iterative use of the relaxation method.
81 Returns (success, i), where success is True if the iteration converged and
82 i is number of the iteration steps.
83 """
84 success = False
85 i = 1
86 # The iterative process continues until number of iteration exceeds `maxiter`
87 # or the L2-norm of the improvements < 'atol'
88 while not(success) and (i < maxiter):
89 dY = k*relaxation.relaxation_step(self._qs, self._Ys,
90 self.F, self.jacF, self.A, self.jacA, self.B, self.jacB)
91 self._Ys += dY
92 self._mdot += dY[0]
93 i += 1
94 success = (np.linalg.norm(self.B(self._Ys[-4:])) < atol)
95 return success, i
96
97 def modify_gamma(self, gamma):
98 """ Change the value of polytropic index. """
99 self._gamma = gamma
100
101 @property
102 def mdot(self):
103 """ Actual accretion rate. """
104 return self._mdot
105
106 @property
107 def gamma(self):
108 """ Actual polytropic index. """
109 return self._gamma
110
111 @property
112 def xs(self):
113 """ Actual position of the sonic point. """
114 return self._Ys[-3]
115
116 def plot_solution_xy(self, ax, **kwargs):
117 """ Place the solution into the axes object ax. """
118 ax.plot(self._Ys[1::4], np.exp(self._Ys[2::4]), **kwargs)
... a příklad použití. Zkoušíme najít řešení pro γ = 1.3 s počátečním odhadem překračujícím správnou
hodnotu o 10%.
1 # set the params and initialize the solver:
2 gamma = 1.3
3 mdot_guess = 1.1*mdot0(gamma)
4 bondi = BondiFlow(gamma=gamma, mdot=mdot_guess, xinf=100., M=100)
5
6 # plot the initial guess (the relaxation is not called yet):
Kapitola 7 Bondiho akrece 72
7 fig, ax = plt.subplots(1)
8 flow_plot(ax, bondi.gamma, mdot0(bondi.gamma))
9 bondi.plot_solution_xy(ax, ls='-', color='gray', marker='', label='initial guess')
10
11 # make the relaxation process:
12 success, i = bondi.calculate(atol=1e-3, maxiter=100, k=0.2)
13
14 # plot the solution:
15 bondi.plot_solution_xy(ax, ls='-', color=blue, marker='', label='relaxed solution')
16 ax.legend()
17
18 # print the report:
19 print('success: {}, number of iterations: {}'.format(success, i))
20 print('initial mdot = {}'.format(mdot_guess))
21 print('mdot found = {}'.format(bondi.mdot))
22 print('correct mdot = {}'.format(mdot0(gamma)))
Obrázek 7.3: Numerické řešení
C Úkol
 Zkuste implementovat řešení pomocí adaptivní sítě bodů.
C
7.5.1  Zkoumání prostoru parametrů
Relaxační metoda je úžasný nástroj pro průzkum prostoru parametrů. Jediný volný parametr je polytropický
index γ. Pokud jiné řešení začneme s počáteční hodnotou γ0, zjistíme, že řešení stejně jako předtím
začíná s počátečním odhadem. Poté změníme poltropický index z γ0 na γ0 + dγ. A použijeme předchozí
řešení jako počáteční odhad pro nové řešení. Tato procedura je implementována níže:
1 # various values of the polytropic index:
2 gammas = np.linspace(1.1, 1.6, 30)
3
4 # initial guess for the first solution, and relaxed solution:
5 mdot_guess = 1.1*mdot0(gammas[0])
Kapitola 7 Bondiho akrece 73
6 bondi = BondiFlow(gamma=gammas[0], mdot=mdot_guess, xinf=100, M=100)
7 sucess, i = bondi.calculate(atol=1e-3, maxiter=100, k=0.2)
8
9 # we store the accretion rate, position of the sonic point and number of iterations:
10 mdots = [bondi.mdot]
11 xss = [bondi.xs]
12 iters = [i]
13
14 # now we go through all values of the polytropic index, modify the solver
15 # and upadate & store the solutions:
16 for gamma in gammas[1:]:
17 bondi.modify_gamma(gamma)
18 sucess, i = bondi.calculate(atol=1e-3, maxiter=100, k=0.2)
19 mdots.append(bondi.mdot)
20 xss.append(bondi.xs)
21 iters.append(i)
22
23 # finally, we plot the resuts. Black solid line is the solution of the analytic theory:
24 fig, panels = plt.subplots(1,3, figsize=(16,4))
25 panels[0].plot(gammas, mdots, ls='', marker='x', label='numeric')
26 panels[0].plot(gammas, [mdot0(gamma) for gamma in gammas], ls='-', color='black', label='analytic')
27 panels[0].set_xlabel(r'Polytropic index $\gamma$')
28 panels[0].set_ylabel(r'Accretion rate $\dot{m}$')
29 panels[0].legend()
30
31 panels[1].plot(gammas, xss, ls='', marker='x', label='numeric')
32 panels[1].plot(gammas, [xs(gamma, mdot0(gamma)) for gamma in gammas], ls='-', color='black', label='analytic')
33 panels[1].set_xlabel(r'Polytropic index $\gamma$')
34 panels[1].set_ylabel(r'Sonic radius $x_\mathrm{s}$')
35 panels[1].legend()
36
37 panels[2].plot(gammas, iters, ls='-', marker='x')
38 panels[2].set_xlabel(r'Polytropic index $\gamma$')
39 panels[2].set_ylabel(r'Number of iterations');
Kapitola 7 Bondiho akrece 74
Souhrn kapitoly
V této kapitole jsme implementovali problém Bondiho akrece a vyřešili jej pomocí relaxační metody. Také
jsme zavedli různé sítě bodů, od lineární, přes logaritmickou až po adaptivní.
? Otázky
1. Jaký má význam rovnice kontinuity?
2. Co reprezentuje konstanta ˙m?
3. Jak se změní řešení při použití jiného polytropického indexu?
4. Jakým způsobem lze určit vnější okrajovou podmínku pro radiální souřadnici?
5. Jaká je výhoda použití logaritmické sítě oproti lineární?
?
Kapitola 8
Některé aproximační metody řešení diferenciálních
rovnic
Martin Blaschke
Rychlý náhled: V této kapitole si představíme aproximační řešení diferenciálních rovnic, konkretně
Poissonovu a Linstedtovu metodu a především Krylovu-Bogoliubovu metodu. Implementace je nastíněna
v programu Wolfram Mathematica. Kapitola je určená pro pokročilé studenty.
Klíčová slova: numerické řešení diferenciálních rovnic, aproximační metody, Poissonova metoda, Linstedtova
metoda, Krylova-Bogoliubova metoda, Wolfram Mathematica
ºº Cíle studia: Cílém této kapitoly je ukázat aproximativní řešení diferenciálních rovnic, která jsou
užitečná v případě, kdy nelineární členy jsou malé (zanedbatelné) vzhledem k lineární části diferenciální
rovnice. Po prostudování kapitoly bude student schopen použít tuto metodu na řešení diferenciálních rovnic.
8.1 Řešení v podobě Talorovy řady generované v Mathematice
Než přejdeme k perturbační metodě, ukážeme si, že program Mathematica je schopen generovat řešení ve
tvaru Taylorovy řady, které bude validní na určitém rozpětí parametrů.
Abychom zjistili, jak lze program Mathematica použít k vygenerování Taylorovy řady, uvažujeme o konkrétním
příkladu jednoduchého (to znamená g = l) matematického kyvadla: θ + sin θ = 0 s počátečním
úhlem θ(0) = π/3 a počáteční rychlostí θ (0) = 0.
MM Řešený příklad
Odvozte řešení Taylorovy řady jednoduchého matematického kyvadla podle dané počáteční podmínky, do
osmého řádu (t8) přesnosti.
Řešení
Za¡neme tím, že rozvineme θ[t] v Taylorovu řadu kolem času t = 0,
p i¢emž řády vyšší než t8
budeme ignorovat a vložíme počáteční podmínky do ř—£¤ F
75
Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 76
Clear["Global`*"]
s1[t_] = Series[θ[t], {t, 0, 8}] /. {θ[0]  π / 3, θ′
[0]  0}
π
3
+
1
2
θ′′
[0] t2
+
1
6
θ(3)
[0] t3
+
1
24
θ(4)
[0] t4
+
1
120
θ(5)
[0] t5
+
1
720
θ(6)
[0] t6
+
θ(7)
[0] t7
5040
+
θ(8)
[0] t8
40 320
+ O[t]9
Toto následně vložíme do diferenciální rovnice a automaticky seřadíme jednotlivé členy podle řády
v proměnné t. Toto seřazení zajistí příkaz Simplify[], nebo v některých případech i příkaz FullSimplify[]
.
ode = Simplify[D[s1[t], {t, 2}] + Sin[s1[t]]  0]
3
2
+ θ′′
[0] + θ(3)
[0] t +
1
4
θ′′
[0] + 2 θ(4)
[0] t2
+
1
12
θ(3)
[0] + 2 θ(5)
[0] t3
+
1
48
-3 3 θ′′
[0]2
+ θ(4)
[0] + 2 θ(6)
[0] t4
+
1
240
-10 3 θ′′
[0] θ(3)
[0] + θ(5)
[0] + 2 θ(7)
[0] t5
+
-10 3 θ(3)
[0]2
- 15 θ′′
[0]3
+ 3 θ′′
[0] θ(4)
[0] + θ(6)
[0] + 2 θ(8)
[0] t6
1440
+ O[t]7
 0
Jelikož proměnná t (čas) může nabývat libovolných hodnot, je jediným způsobem jak zajisti obecné
řešení to, že všechny koeficienty ve všech řádech t jsou nulové. Tato podmínka může být vyjádřena pomocí
příkazu LogicalExpand[].
eq1 = LogicalExpand[ode]
3
2
+ θ′′
[0]  0 && θ(3)
[0]  0 &&
1
4
θ′′
[0] + 2 θ(4)
[0]  0 &&
1
12
θ(3)
[0] + 2 θ(5)
[0]  0 &&
1
48
-3 3 θ′′
[0]2
+ θ(4)
[0] + 2 θ(6)
[0]  0 &&
1
240
-10 3 θ′′
[0] θ(3)
[0] + θ(5)
[0] + 2 θ(7)
[0]  0 &&
-10 3 θ(3)
[0]2
- 15 θ′′
[0]3
+ 3 θ′′
[0] θ(4)
[0] + θ(6)
[0] + 2 θ(8)
[0]
1440
 0
Symbol && objevující se mezi rovnicemi reprezentuje logický součet. Tento systém obyčejných rovnic
můžeme nyní přímo vyřešit.
Solve[eq1]
θ′′
[0]  -
3
2
, θ(3)
[0]  0, θ(4)
[0] 
3
4
,
θ(5)
[0]  0, θ(6)
[0]  3 , θ(7)
[0]  0, θ(8)
[0]  -
49 3
8

Řešení ve formě řady pak nalezneme přímím použitím jednotlivých výsledků.
sol = s1[t] /. %[[1]]
π
3
-
3 t2
4
+
t4
32 3
+
t6
240 3
-
7 t8
15 360 3
+ O[t]9
Povšimněme si, že vyjádření neobsahuje liché členy řady. Ke konečnému vyjádření výsledku můžeme
použít příkazu Normal[] k odstranění symbolu reprezentujícího vyšší řády přesnosti.
Normal[sol]
Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 77
π
3
-
3 t2
4
+
t4
32 3
+
t6
240 3
-
7 t8
15 360 3
Toto řešení ve tvaru řady bude přesné (lépe řečeno užitečné) pouze pro omezené množství času t,
o čemž se v tomto případě lze přesvědčit porovnáním s analytickým řešením.
C Úkoly
1. Spočítejte řešení příkladu až do O(t15) a prodiskutujete možnost a praktičnost nalezení tohoto řešení
metodou papír - tužka.
2. Vykreslete (Plot[]) řešení pomocí řady pro různé řády přesnosti vzhledem k analytickému řešení.
3. Nahraďte sin (x) členem (x2 +x3) cos (x) a proveďte celý výpočet znovu. Opět uvažujte do jaké míry
by tento postup byl užiteční pokud by byl prováděn metodou papír - tužka.
C
8.2 Poissonova metoda
Perturbační metody jsou použitelné, když je nultý řád analyticky řešitelný a vyšší řády jsou malé. Tuto
malost zajišťuje parametr , který je do problému vložen ručně. Uvažujme například Van der Polovu rovnici
D[x[t], {t, 2}] - ϵ (1 - x[t]^2) D[x[t], t] + x[t]  0
x[t] - ϵ 1 - x[t]2
 x′
[t] + x′′
[t]  0
Parametr byl vložen před závorku tak, aby nultý řád ( = 0) byla jednoduchá lineární diferenciální
rovnice.
Základní myšlenka perturbačních metod spočívá v tom, že i řešení se dá rozepsat pomocí mocninné
řady v (ne nutně Taylorovi).
Jestliže je skutečně malé, obyčejně potřebujeme pouze několik členů řady k nalezení přibližného řešení
dostatečné přesnosti.
Mocninné řady obvykle divergují, a proto je občas potřeba užít dalších metod, které extrahují užitečnou
informaci z divergující řady (Padé aproximace).
MM Řešený příklad
Nalezněte řešení do 6-tého řádu v pomocí Poissonovy perturbační metody pro
x[τ] + ϵ x[τ]2
+ x′
[τ]  0
Řešení
Clear["Global`*"]
n = 6;
x[τ] = Sum[xi[τ] ϵ^i, {i, 0, n}]
x0[τ] + ϵ x1[τ] + ϵ2
x2[τ] + ϵ3
x3[τ] + ϵ4
x4[τ] + ϵ5
x5[τ] + ϵ6
x6[τ]
Pozn. Index se vytvo¡í kombinací kláves Ctrl + - . Nelineární rovnice je poté rozvedena
eq = Expandx[τ] + ϵ x[τ]2
+ D[x[τ], τ]
Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 78
x0[τ] + ϵ x0[τ]2
+ ϵ x1[τ] + 2 ϵ2
x0[τ] x1[τ] + ϵ3
x1[τ]2
+ ϵ2
x2[τ] +
2 ϵ3
x0[τ] x2[τ] + 2 ϵ4
x1[τ] x2[τ] + ϵ5
x2[τ]2
+ ϵ3
x3[τ] + 2 ϵ4
x0[τ] x3[τ] +
2 ϵ5
x1[τ] x3[τ] + 2 ϵ6
x2[τ] x3[τ] + ϵ7
x3[τ]2
+ ϵ4
x4[τ] + 2 ϵ5
x0[τ] x4[τ] +
2 ϵ6
x1[τ] x4[τ] + 2 ϵ7
x2[τ] x4[τ] + 2 ϵ8
x3[τ] x4[τ] + ϵ9
x4[τ]2
+ ϵ5
x5[τ] +
2 ϵ6
x0[τ] x5[τ] + 2 ϵ7
x1[τ] x5[τ] + 2 ϵ8
x2[τ] x5[τ] + 2 ϵ9
x3[τ] x5[τ] +
2 ϵ10
x4[τ] x5[τ] + ϵ11
x5[τ]2
+ ϵ6
x6[τ] + 2 ϵ7
x0[τ] x6[τ] + 2 ϵ8
x1[τ] x6[τ] +
2 ϵ9
x2[τ] x6[τ] + 2 ϵ10
x3[τ] x6[τ] + 2 ϵ11
x4[τ] x6[τ] + 2 ϵ12
x5[τ] x6[τ] +
ϵ13
x6[τ]2
+ x0
′
[τ] + ϵ x1
′
[τ] + ϵ2
x2
′
[τ] + ϵ3
x3
′
[τ] + ϵ4
x4
′
[τ] + ϵ5
x5
′
[τ] + ϵ6
x6
′
[τ]
A jednotlivé koeficienty jsou dány dohromady do eq2.
eq2 = Collect[eq, ϵ]
x0[τ] + 2 ϵ12
x5[τ] x6[τ] + ϵ13
x6[τ]2
+
ϵ7
x3[τ]2
+ 2 x2[τ] x4[τ] + 2 x1[τ] x5[τ] + 2 x0[τ] x6[τ] +
ϵ8
(2 x3[τ] x4[τ] + 2 x2[τ] x5[τ] + 2 x1[τ] x6[τ]) +
ϵ9
x4[τ]2
+ 2 x3[τ] x5[τ] + 2 x2[τ] x6[τ] + ϵ10
(2 x4[τ] x5[τ] + 2 x3[τ] x6[τ]) +
ϵ11
x5[τ]2
+ 2 x4[τ] x6[τ] + x0
′
[τ] + ϵ x0[τ]2
+ x1[τ] + x1
′
[τ] +
ϵ2
(2 x0[τ] x1[τ] + x2[τ] + x2
′
[τ]) + ϵ3
x1[τ]2
+ 2 x0[τ] x2[τ] + x3[τ] + x3
′
[τ] +
ϵ4
(2 x1[τ] x2[τ] + 2 x0[τ] x3[τ] + x4[τ] + x4
′
[τ]) +
ϵ5
x2[τ]2
+ 2 x1[τ] x3[τ] + 2 x0[τ] x4[τ] + x5[τ] + x5
′
[τ] +
ϵ6
(2 x2[τ] x3[τ] + 2 x1[τ] x4[τ] + 2 x0[τ] x5[τ] + x6[τ] + x6
′
[τ])
Protože je libovolný parametr musí být v řešení každý koeficient roven nule. Příkaz CoefficientList[]
nám udává seznam koeficientů, které je vhodné zobrazit v tabulce.
Table[CoefficientList[eq2, ϵ][[i]]  0, {i, 1, n + 1}] // TableForm
x0[τ] + x0
′
[τ]  0
x0[τ]2
+ x1[τ] + x1
′
[τ]  0
2 x0[τ] x1[τ] + x2[τ] + x2
′
[τ]  0
x1[τ]2
+ 2 x0[τ] x2[τ] + x3[τ] + x3
′
[τ]  0
2 x1[τ] x2[τ] + 2 x0[τ] x3[τ] + x4[τ] + x4
′
[τ]  0
x2[τ]2
+ 2 x1[τ] x3[τ] + 2 x0[τ] x4[τ] + x5[τ] + x5
′
[τ]  0
2 x2[τ] x3[τ] + 2 x1[τ] x4[τ] + 2 x0[τ] x5[τ] + x6[τ] + x6
′
[τ]  0
sol0 = Flatten[DSolve[{CoefficientList[eq2, ϵ][[1]]  0, x0[0] == 1}, x0[τ], τ]]
{x0[τ]  -τ
}
sol1 =
Flatten[DSolve[{CoefficientList[eq2, ϵ][[2]]  0 /. sol0, x1[0]  0}, x1[τ], τ]]
sol2 = Flatten[
DSolve[{CoefficientList[eq2, ϵ][[3]]  0 /. sol0 /. sol1, x2[0]  0}, x2[τ], τ]]
sol3 = Flatten[DSolve[{CoefficientList[eq2, ϵ][[4]]  0 /. sol0 /. sol1 /. sol2,
x3[0]  0}, x3[τ], τ]]
sol4 = Flatten[DSolve[{CoefficientList[eq2, ϵ][[5]]  0 /. sol0 /. sol1 /. sol2 /.
sol3, x4[0]  0}, x4[τ], τ]]
sol5 = Flatten[DSolve[{CoefficientList[eq2, ϵ][[6]]  0 /. sol0 /. sol1 /. sol2 /.
sol3 /. sol4, x5[0]  0}, x5[τ], τ]]
sol6 = Flatten[DSolve[{CoefficientList[eq2, ϵ][[7]]  0 /. sol0 /. sol1 /. sol2 /.
sol3 /. sol4 /. sol5, x6[0]  0}, x6[τ], τ]]
x1[τ]  --2 τ
(-1 + τ
)
x2[τ]  -3 τ
(-1 + τ
)2

x3[τ]  --4 τ
(-1 + τ
)3

x4[τ]  -5 τ
(-1 + τ
)4

Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 79
x5[τ]  --6 τ
(-1 + τ
)5

x6[τ]  -7 τ
(-1 + τ
)6

Simplify[x[τ] /. sol0 /. sol1 /. sol2 /. sol3 /. sol4 /. sol5 /. sol6]
-7 τ
6 τ
- 5 τ
(-1 + τ
) ϵ + 4 τ
(-1 + τ
)2
ϵ2
3
τ
(-1 + τ
)3
ϵ3
+ 2 τ
(-1 + τ
)4
ϵ4
- τ
(-1 + τ
)5
ϵ5
+ (-1 + τ
)6
ϵ6

GeneralSol = FullSimplify Sum(-ϵ)^k -(k+1) τ
(-1 + τ
)^(k), {k, 0, Infinity}
1
-ϵ + τ
(1 + ϵ)
V konečném kroku, můžeme jednotlivé řády perturbativního řešení porovnat s analytickým výsledkem.
Zde se obvykle pokládá hodnotě 1.
Plot
1
-ϵ + τ
(1 + ϵ)
/. ϵ  1, Sum(-ϵ)^k -(k+1) τ
(-1 + τ
)^(k), {k, 0, 1} /. ϵ  1,
Sum(-ϵ)^k -(k+1) τ
(-1 + τ
)^(k), {k, 0, 2} /. ϵ  1,
Sum(-ϵ)^k -(k+1) τ
(-1 + τ
)^(k), {k, 0, 3} /. ϵ  1, {τ, 0, 4}
1 2 3 4
0.2
0.4
0.6
0.8
1.0
C Úkoly
1. Nalezněte perturbační řešení pro nelineární tvrdou strunu danou rovnicí:
x′′
[τ] + ω2
x[τ] + x3
[τ]  0
2. Nalezněte Poissonovu expanzi do desátého řádu v pro ODE:
x′
[τ] + x[τ] + x4
[τ]  0
x′
[τ] - x[τ] + x5
[τ]  0
x′
[τ] - x[τ] - x4
[τ]  0
C
8.3 Lindstedtova metoda (Poincarého-Lindstedtova metoda)
V poruchové teorii je Lindstedtova metoda (technika) užitá pro uniformní aproximaci periodických řešení
obyčejných dif. rovnic (ODE), když ostatní metody selhávají.
Demonstrujme metodu na jednoduchém příkladu tvrdé nelineární struny. Abychom mohli vyřešit rov-
nici
x′′
[t] + ω2
x[t] + ϵ x3
[t]  0
pro malé hodnoty , předpokládáme, že frekvence se změní příspěvkem nelineárního členu, a že může
být vyjádřena jako poruchová řada. Nechť Ω je neznámá frekvence periodického řešení strunové rovnice,
pak provedeme substituci τ = Ωt tak, aby řešení mělo periodu 2π vzhledem k nové proměnné τ.
Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 80
x[τ + 2 π] = x[τ]
Takže nelineární strunová rovnice přejde od tvaru
Ω2
x′′
[τ] + ω2
x[τ] + ϵ x3
[t]  0
Lindstedtova metoda předpokládá, že
x (τ) = x0 (τ) + ϵ x1 (τ) + ...
a
Ω = Ω0 + ϵ Ω1 + ...
Pro = 0, musí zřejmě platit Ω0 = ω
MM Řešený příklad
Nalezněte Lindstedtovo řešení pro tvrdou strunu v prvním řádu .
x′′
[t] + ω2
x[t] + ϵ x3
[t]  0
Řešení:
Clear["Global`*"]
n = 1;
x[τ] = Sum[xi[τ] ϵ^i, {i, 0, n}]
Ω = Sumωi ϵi
, {i, 0, n}
x0[τ] + ϵ x1[τ]
ω0 + ϵ ω1
Vložíme rovnici pro tvrdou strunu a řešení rozvedeme.
eq1 = Expand[Ω^2 D[x[τ], {τ, 2}] + ω0 ^2 x[τ] + ϵ (x[τ])^3]
ω0
2
x0[τ] + ϵ x0[τ]3
+ ϵ ω0
2
x1[τ] + 3 ϵ2
x0[τ]2
x1[τ] + 3 ϵ3
x0[τ] x1[τ]2
+ ϵ4
x1[τ]3
+
ω0
2
x0
′′
[τ] + 2 ϵ ω0 ω1 x0
′′
[τ] + ϵ2
ω1
2
x0
′′
[τ] + ϵ ω0
2
x1
′′
[τ] + 2 ϵ2
ω0 ω1 x1
′′
[τ] + ϵ3
ω1
2
x1
′′
[τ]
Seřadíme jednotlivé členy podle řádu v a vygenerujeme potřebný počet dif. rovnic.
eq2 = Collect[eq1, ϵ]
odsystem =
Table[Expand[CoefficientList[eq2, ϵ][[i]] / ω0 ^2]  0, {i, 1, n + 1}] // TableForm
ω0
2
x0[τ] + ϵ4
x1[τ]3
+ ω0
2
x0
′′
[τ] + ϵ x0[τ]3
+ ω0
2
x1[τ] + 2 ω0 ω1 x0
′′
[τ] + ω0
2
x1
′′
[τ] +
ϵ2
3 x0[τ]2
x1[τ] + ω1
2
x0
′′
[τ] + 2 ω0 ω1 x1
′′
[τ] + ϵ3
3 x0[τ] x1[τ]2
+ ω1
2
x1
′′
[τ]
x0[τ] + x0
′′
[τ]  0
x0[τ]3
ω0
2
+ x1[τ] +
2 ω1 x0
′′
[τ]
ω0
+ x1
′′
[τ]  0
Nyní jsi tyto ODE separátně zobrazíme
od0 = odsystem[[1, 1]]
od1 = odsystem[[1, 2]]
x0[τ] + x0
′′
[τ]  0
x0[τ]3
ω0
2
+ x1[τ] +
2 ω1 x0
′′
[τ]
ω0
+ x1
′′
[τ]  0
První rovnici od0 analyticky vyřešíme za určitých okrajových podmínek, např.
Sol0 = DSolve[{od0, x0[0]  A, x0
′
[0]  0}, x0[τ], τ]
x0 = x0[τ] /. Sol0[[1]]
{{x0[τ]  A Cos[τ]}}
A Cos[τ]
Tento nulový řád řešení je následně vložen do rovnice od1 jako funkce:
Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 81
od1 = od1 /. x0  Function[τ, A Cos[τ]]
A3
Cos[τ]3
ω0
2
-
2 A Cos[τ] ω1
ω0
+ x1[τ] + x1
′′
[τ]  0
Poté je i od1 vyřešeno za použitím nulových okrajových podmínek, neboť informaci o okrajových
podmínkách systému již obsahuje nultý řád řešení.
Sol1 = FullSimplify[DSolve[{od1, x1[0]  0, x1
′
[0]  0}, x1[τ], τ]]
x1[τ]  A
Sin[τ] A2
(6 τ + Sin[2 τ]) - 16 τ ω0 ω1
16 ω0
2

Členy obsahující τ jsou seskupeny
x1 = Collect[x1[τ] /. Sol1[[1]], τ]
-
A3
Sin[τ] Sin[2 τ]
16 ω0
2
A
τ Sin[τ] 6 A2
- 16 ω0 ω1
16 ω0
2
Koeficient stojící před τ je položen nule čímž nalézáme první řád opravy k frekvenci.
freq = Solve[Coefficient[x1, τ]  0, ω1]
ω1 
3 A2
8 ω0

Frekvence je do prvního řádu v dána vztahem:
Ω = Ω /. freq[[1]]
3 A2
ϵ
8 ω0
+ ω0
K vyjádření prvního řádu řešení můžeme například položit τ sin τ → 0 a zjednodušit výraz pomocí
příkazu TrigReduce[], nebo příkazem Simplify[]:
x1 = TrigReduce[x1 /. τ Sin[τ]  0]
x1p = Simplify[x1 /. τ Sin[τ]  0]
-A3
Cos[τ] + A3
Cos[3 τ]
32 ω0
2
-
A3
Cos[τ] Sin[τ]2
8 ω0
2
Celkové řešení do patřičného řádu je:
x = x0 + ϵ Collect[x1, A^3]
xp = x0 + ϵ Collect[x1p, A^3]
A Cos[τ] +
A3
ϵ (-Cos[τ] + Cos[3 τ])
32 ω0
2
A Cos[τ] -
A3
ϵ Cos[τ] Sin[τ]2
8 ω0
2
C Úkoly
1. Analyticky (pomocí programu Mathematica) dokažte následující trigonometrickou identitu:
Cos[5 τ] == 16 (Cos[τ])^5 - 20 (Cos[τ])^3 + 5 Cos[τ]
2. Nalezněte řešení tvrdé struny do třetího řádu a ten porovnejte s analytickým řešení. Diskutujte, proč
metoda selhává!
C
Kapitola 8 Některé aproximační metody řešení diferenciálních rovnic 82
Souhrn kapitoly
V této kapitole jsme představili některé aproximační metody řešení diferenciálních rovnic a ukázali si jejich
použití při řešení fyzikálních problému v programu Wolfram Mathematica.
? Otázky
1. Kdy můžeme použít aproximativní metody? Jaký typ diferenciálních rovnic tak můžeme řešit?
2. Pro jaké případy se používá Lindstedtova metoda?
3. Jak se k řešení používá Taylorova řada?
?
Kapitola 9
Profil spektrální čáry od rotujícího prstence hmoty
Vladimír Karas
Rychlý náhled: Tato kapitola popisuje vznik spektrálních čar vznikajících v prstenci hmoty rotujícím
kolem kompaktního tělesa. Tato kapitola je pouze zadání úlohy, kterou student může vypracovat v rámci
závěrečné seminární práce.
Klíčová slova: Spektrální čáry, kompaktní objekty, numerická integrace
ºº Cíle studia: Cílem kapitoly je představit studentům princip vzniku spektrálních čar z rotujícího
prstence v okolí kompaktního objektu tak, aby byli studenti schopni samostatně navrhnout a implementovat
řešení.
MM Spektrální profil emisní čáry vznikající v rotujícím prstenci
Spektrální čáry jsou zpravidla vyzařovány určitým zářivým přechodem v atomech zahřátého plynu. Vlnová
délka λ (nebo ekvivalentně energie resp. frekvence vyzařovaných fotonů ν) je přesně dána. Následkem pohybu
jednotlivých atomů a také vlivem gravitačního pole v místě zdroje se však může výsledná pozorovaná
spektrální čára jevit vzdálenému pozorovateli rozšířená, s určitým deformovaným profilem. Také může dojít
k posunu její vlnové délky (či energie E = hν) vlivem Dopplerova jevu a gravitačního červeného posunu.
Pozorovaný zářivý tok v laboratorní soustavě spojené s teleskopem lze symbolicky zapsat jako integraci
F(ν) =
(Over observer’s plane)
I(ν)|θ=θobs
dS. (9.1)
M Příklad
Převeďte výše uvedený obecný výraz do vhodného tvaru sumy příspěvků od jednotlivých elementů zdroje
záření, který pro účely tohoto cvičení považujte za kruhový prstenec o poloměru R. Numerickou integrací
přes azimutální úhel podél prstence určete pozorovaný spektrální profil F(ν) od spektrální čáry emitované
na jednotkové energii prstencem s danou konstantní rotační rychlosti v(R) v ekvatoriální rovině. Objekt je
pozorován vzdáleným pozorovatelem pod daným úhlem (tzv. inklinací) v intervalu 0 (podél osy prstence)
až 90 (v rovině prstence) stupňů.
83
Kapitola 9 Spektrální profily 84
Pozn.: světelné paprsky považujeme z přímky. Neuvažujeme na tomto místě zakřivení paprsků
v silné gravitaci vlivem efektů obecné teorie relativity. Poloha středu spektrální čáry (tzv. centroidu)
zjevně nezávisí na inklinaci disku. Při úplném, self-konzistentním relativistickém řešení tohoto problému
není možné oddělovat gravitační posuv od Dopplerova. Vlnová délka centroidu bude na inklinaci záviset,
uvážíme-li úplný vztah posuv vlnové délky, anizotropní emisivitu prostředí, závislou na úhlu vyzařování,
i gravitační fokusaci světelných paprsků.
££ Řešení: Vztah mezi pozorovatelem (obs) zaznamenanou a zdrojem emitovanou (em) vlnovou délkou
lze zapsat v přibližném tvaru (Gerbal & Pelat 1981; Karas et al. 1995)
λobs = λem
(1 − βy)
1 − β2
≈ λem 1 +
1
2
β2
− βy , (9.2)
kde β = v/c, y = sin φ sin θobs (φ je azimutální úhel v rovině prstence a θobs je inklinace pozorovatele, c
značí rychlost světla).
Obdobně lze odvodit vztah pro gravitační červený posuv
λobs =
λem
1 − β2
. (9.3)
Komentář: Záření plynu, obíhajícího u vnitřního okraje disku, podléhá silnějšímu gravitačnímu
a příčnému Dopplerovu jevu, než záření z odlehlejších oblastí. Výsledkem je nesouměrnost pozorovaných čar,
které jsou zároveň celkově posunuty k červenému konci spektra. K tomu přistupuje dopplerovsky zesílené
záření přibližujícího se materiálu, jež zvýrazňuje vrchol na modré straně čáry. Graf 9.1 ilustruje spektrální
profil, kde elementy prstence vydávají záření o frekvenci ν = ν0, představující v lokální klidové soustavě
úzkou spektrální čáru s charakterem Diracovy funkce δ. Profil se potom vytváří překrytím jednotlivých
příspěvků ovlivněných Dopplerovým jevem a gravitačním červeným posuvem. Pozorovaný profil závisí na
pozorovatelově inklinaci θobs.
M
C Úkoly
1. Diskutujte význam největší a nejnižší frekvence pozorované čáry. Za předpokladu keplerovské rotace
kolem centrálního tělesa o hmotnosti M splňuje orbitální rychlost látky, odpovídající vrcholům
v profilech čar, vztah
vKα|R=Rin
vKα|R=Rout
=
Rout
Rin
nezávisle na úhlu inklinace. Tento poměr poskytuje první odhad rozměru disku.
2. Integrací přes R zjistěte spektrální profil čáry od radiálně rozlehlého disku. V případě disku se profil
čáry získá jeho rozdělením do jednotlivých prstenců o radiální šířce dr, z nichž každý vyzařuje se svou
vlastní lokální frekvencí a intenzitou, a následným složením výsledného signálu. Integrace probíhá od
určitého vnitřního okraje Rin do vnějšího poloměru Rout. Prostudujte polohu centroidu energie Ec
spektální čáry rozlehlého disku.
C
Kapitola 9 Spektrální profily 85
Obrázek 9.1: Spektrální čáry plynného prstence na poloměru 6 Rg v Newtonovské aproximaci normalizovaný:
vlevo – tak, že Fνdν = 1, vpravo – k velikosti peaku modrého posunu na dané spektrální čáře.
Černá tečkovaná čára ukazuje vlastní frekvenci.
Obrázek 9.2: Spektrální čáry plynného prstence na poloměru 6 Rg ve speciální relativitě, normalizovaný:
vlevo – tak, že Fνdg = 1, vpravo – k velikosti peaku modrého posunu na dané spektrální čáře. Černá
tečkovaná čára ukazuje vlastní frekvenci.
Obrázek 9.3: Spektrální čáry plynného prstence na poloměru 6 Rg ve obecné relativitě, normalizovaný:
vlevo – tak, že Fνdg = 1, vpravo – k velikosti peaku modrého posunu na dané spektrální čáře. Černá
tečkovaná čára ukazuje vlastní frekvenci.
Kapitola 9 Spektrální profily 86
Další informace
1 Gerbal, D., & Pelat, D. 1981, “Profile of a line emitted by an accretion disk. Influence of the geometry
upon its shape parameters”, Astronomy & Astrophysics 95, 18;
2 Karas, V., Lanza, A., & Vokrouhlický, D. 1995, “Emission-line profiles from self-gravitating thin
disks”, The Astrophysical Journal 440, 108;
3 Kojima, Y. 1991, “The effects of black hole rotation on line profiles from accretion discs”, Monthly
Notices of the Royal Astronomical Society 250, 629
4 Laor, A. 1991, “Line profiles from a disk around a rotating black hole”, The Astrophysical Journal
376, 90
5 Stolc, M. 2019, “Accretion discs in the context of tidal disruption of stars in nuclei of galaxies”, Master
Thesis, Prague
Souhrn kapitoly
V této kapitole jsme uvedli studenty do problematiky spektrálních čar, které jsou základem pro samostatnou
práci.
? Otázky
1. Co je zdrojem záření vytvářející spektrální čáry?
2. Jak se mění tvar čar, používáme-li stále přesnější teorii gravitace? Čím je to způsobeno?
?
Literatura
[1] Stewart, William J. Probability, Markov Chains, Queues and Simulation: The Mathematical Basis of
Performance Modeling. [s.l.]: Princeton University Press, 2011. Dostupné online. ISBN 978-1-4008-
3281-1
[2] GRINSTEAD, Charles M.; SNELL, J. Laurie. Grinstead & Snell, Úvod do pravděpodobnosti. [s.l.]:
Orange Grove Texts, 2009. ISBN 161610046X.
[3] HUBBLE, Edwin. A Relation between Distance and Radial Velocity among Extra-Galactic Nebulae.
Proceedings of the National Academy of Sciences of the United States of America. 1929, 15, 168-173.
[4] Kapteyn Astronomical Institute. Least squares fitting with kmpfit. Dostupné z
https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html
[5] KIRSHNER, Robert P. Hubble’s diagram and cosmic expansion. Proceedings of the National Academy
of Science. 2004, 101, 8-13.
[6] PRESS William H., TEUKOLSKY Saul A., VETTERLING William T., FLANNERY Brian P. Numerical
recipes. The art of scientific computing. Third edition. Cambridge. 2007, ISBN-13 978-0-521-
88068-8
[7] ČERMÁK, Libor a Rudolf HLAVIČKA. Numerické metody I. In: VUT, FSI, Ústav matematiky [online].
Brno, 1.11.2016. Dostupné z: http://mathonline.fme.vutbr.cz/Numericke-metody-I/sc-8-sr-1-a-
11/default.aspx
87