PA081: Programování numerických
výpočtů
Dynamické systémy
PA081: Programování numerických
výpočtů
6. Simulace dynamických dějů
jaro 2019
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Něco, co se hýbe
► objekty reálného světa kolem nás
► velká část počítačových her
► vesmírná tělesa
► plyn a kapaliny
► předpověď počasí
► simulace požárů, povodní, ...
► elektrický proud
► interakce atomů a molekul
Programování numerických výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Jak se to hýbe
PA081: Programování numerických
výpočtů
► fyzikální zákony
► postihují principy chování reálného světa
► jejich dodržení v simulaci dodává důvěryhodnost
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
3/47
Jak se to hýbe
► fyzikální zákony
► postihují principy chování reálného světa
► jejich dodržení v simulaci dodává důvěryhodnost
► Newtonovy
1. setrvačnosti
2. síly - F = ma
3. akce a reakce
► termodynamické
1. Celková energie izolovaného systému zůstává konstantní.
2. Entropie izolovaného systému nikdy neklesá.
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
2. termodynamický zákon
► energie má snahu rovnoměrně se rozptýlit.
► psychologické vnímání času je souhrn zkušeností s jeho projevy.
► pozpátku pustený film poznáme na první pohled.
► nelze zcela snadno vyjádřit proč.
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
2. termodynamický zákon
► energie má snahu rovnoměrně se rozptýlit.
► psychologické vnímání času je souhrn zkušeností s jeho projevy.
► pozpátku pustený film poznáme na první pohled.
► nelze zcela snadno vyjádřit proč.
► dovoluje uchovávat energii (palivo/potrava + kyslík)
► de-facto podmínka existence života
► ale také příčina smrti
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
2. termodynamický zákon
► energie má snahu rovnoměrně se rozptýlit.
► psychologické vnímání času je souhrn zkušeností s jeho projevy.
► pozpátku pustený film poznáme na první pohled.
► nelze zcela snadno vyjádřit proč.
► dovoluje uchovávat energii (palivo/potrava + kyslík)
► de-facto podmínka existence života
► ale také příčina smrti
► příčina Murphyho zákonů
► více viz http: //secondl aw. oxy. edu
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Stavba simulace
► stavové proměnné
► popisují, v jakém stavu se systém nachází
► poloha, rychlost, ...
► jejich změny v čase
► jak na sebe stavové proměnné působí vzájemně (rychlost změna polohy)
► jak se projevují vnější impulsy (síla mění rychlost)
► diferenciální rovnice (obyčejné)
► vztah stavových proměnných a jejich časových derivací
► postupným řešením (integrací) získáme vývoj systému v čase
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Řešení diferenciálních rovnic
► analytické
► jen speciální případy
► i tak komplikované metody
► dopředná (explicitní) Eulerova metoda
► použita v následujících příkladech
► v praxi téměř nepoužitelná (viz dále)
► zpětná (implicitní) Eulerova metoda
► semiimplicitní metody
► komplexnější metody
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Padající kulička
► stavové proměnné: výška z, rychlost pádu v
► změna výšky odpovídá rychlosti: dz = -vdt
► změna rychlosti - gravitační síla: dv = gdt
► korektní zápis systému rovnic
dz/dt = -v     dv/dt = g
► v kódu přímočará náhrada dt konečným krokem
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Pružinový oscilátor
PA081: Programování numerických
výpočtů
► hmotný bod na pružině, pohyb jen v jednom směru
► vše ostatní zanedbáno
► změna polohy y odpovídá rychlosti, tj. hybnost/hmotnost
d y = —d t m
► změna hybnosti odpovídá působící síle, ta je úměrná vychýlení pružiny
dP = F = -kydt
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Pružinový oscilátor
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
9/47
Průlet komety
► poloha x, nenulová počáteční rychlost v
► pro zjednodušení vektory ve 2D
► změna polohy odpovídá rychlosti: dx = \dt
► změna rychlosti dána gravitačním působením Slunce:
PA081: Programování numerických
výpočtů
X
X
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Jednoduchý elektrický obvod
► vstupním napětím U nabíjíme kondenzátor přes odpor
R
u
C u
► výstupní napětí je úměrné kapacitě a náboji: u = q/C
► změna náboje je velikost proudu: dq/dt = í
► proud podle Ohmová zákona: í = (U - u)/R
oat u = 0, U = 14, R = 22e3, C = le-6, dt = 0.00001
i 1 e (...) {
float dq = (U-u)/R * dt; u += dq/C;
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
► stejně i pro proměnlivý vstup, např. U = sin(27r -10 3 ■ t
Problémy dopředné metody
► nejvíce se projeví na pružinovém oscilátoru
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
10 15 20 25 30 35 40 45 50
Ař = 0.0001 (zelená), Ař = 0.01 (červená)
Problémy dopředné metody
► řešíme rovnici (resp. systém rovnic)
dy dt
= f(t,y(t))
► počítáme posloupnost
y-n+i = yn + At
dy dt
► výsledná sin x je ideálně nešikovná funkce
► pro x e [0,7t] je sinx > 0 a je konkávni, pro x e [tt, 2tt] je sinx < 0 a je konvexní
► použitá derivace vždy nadhodnocuje
► výpočet do systému "přidává energii"
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
13/47
Zpětná Eulerova metoda
► použijeme derivaci až v bodě tn+i
y-n+i = yn + At
dy_ dt
tn+l
► dosazením vede na řešení rovnice
y-n+i = yn + Atf(tn+1,yn+1))
► implicitní metoda
► pro výpočet dalšího kroku musíme řešit rovnici nebo systém rovnic
► obecně nemusí být triviální
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Zpětná Eulerova metoda
PA081: Programování numerických
výpočtů
► pro pružinový oscilátor
yn+i — y n +
m
Pn+1
Pn+1 = Pn~ Atkyn+1
► jednoduchý systém lineárních rovnic
► řešení lze vyjádřit analyticky a naprogramovat
Dynamické systémy
Diferenciální rovnice Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
15/47
Zpětná Eulerova metoda
► výstup pro Ař = 0.01
10 15 20 25 30 35 40 45 50
PA081: Programování numerických
výpočtů
► metoda také není numericky stabilní
► má sklon systém tlumit
► pro simulace v principu vyhovuje
► nutné řešení systému rovnic v každém kroku může být problém
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
16/47
Semiimplicitní metoda
► kompromis mezi dopřednou (explicitní) a zpětnou (implicitní)
► nevyžaduje řešení komplikovaného systému rovnic
► nahrazuje ho vhodnou aproximací
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Semiimplicitní metoda
► kompromis mezi dopřednou (explicitní) a zpětnou (implicitní)
► nevyžaduje řešení komplikovaného systému rovnic
► nahrazuje ho vhodnou aproximací
► časovou derivaci dy/dt chápeme jako funkci f (y)
► dosadíme ji do zpětné metody
Yn+l = Yn + Atf (y
►
použijeme Taylorův rozvoj
f(Yn+l) = f(Yn) +
3y
(Yn+l ~Yn) + * * *
►
zanedbáme vyšší členy
Kometa semiimplicitní metodou
► stavový vektor y = (x, v)r
► pohybová rovnice
f(y) =
dy
dx/dt dv/dt
kde g = -Gmkfnsx/r3 a r = ||x ► Jakobián pro Taylorův rozvoj
df_
dy
dv/dx dv/dv dg/dx dg/d\
v g
0 1
dg/dx 0
PA081: Programování numerických
výpočtů
► výpočet dg/dx podle pravidel o derivování složené funkce
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
18/47
Kometa semiimplicitní metodou
PA081: Programování numerických
výpočtů
► v čase n máme přímo k dispozici hodnoty stavových proměnných y = (x, v)r
► na základě pohybové rovnice dy/dt = ... (zahrnuje i externí silové působení) spočteme jejich derivace, tj. f (yn)
► podle předem spočteného vzorce získáme dg/dx v čase n
► dosadíme do rovnice
- Ař
3y
y
)
(Yn+i ~ Yn) = Atf (yn)
► řešením je Ay = yM+1 - yn
► narozdíl od zpětné metody je to jen systém 4 lineárních rovnic
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
19/47
Zábo, skoč
PA081: Programování numerických
výpočtů
► metoda „leap frog"
► specializovaná pro systémy tvaru
dx/dt = v     d\/dt = F(x)
► v sudém kroku integrujeme pozici, v lichém rychlost
X2n = X2n-2 + AtV2n-l V2n+1 = V2n-1 + AtF(x2n)
► pro harmonické systémy je numericky stabilní
► nepřidává energii
► chyba dopředně metody v sudém kroku je kompenzována chybou v lichém kroku
► v praxi velmi rozšířená metoda
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
20/47
Metody vyšších řádů
► obecná formule metody 1. řádu
y-n+i = yn + hf(xn,yn)
pro f(xn,yn) = je to Eulerova dopředná
► lze ukázat, že její chyba je 0(h2)
► Taylorovým rozovjem hypoteticky známeho výsledku
► o 1 vyšší exponent kroku, proto 1. řád
► metody vyšších řádů mají chybu 0(hk+1)
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Ilustrativní metoda 2. řádu
► začneme stejné
fci = hf(xn,yn)
► pokusný mezikrok do poloviny intervalu h
1 i
k2 = hf(xn + 7>h,yn + -fci)
► definitivní krok s použitím výsledků mezikroku
► lze ukázat, že chyba je 0(h?)
► k posunu o h potřebujeme 2x vyhodnotit /
► lepší než metoda 1. řádu s polovičním krokem
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Runge-Kutta
PA081: Programování numerických
výpočtů
► metoda 4. řádu
► posloupnost výpočtu
	= hf(xn	Jíl)
k2	= hf(xn	+ -^Kyn + -ki)
	= hf(xn	+             + 2fe2)
k.4	= hf(xn	+ fi,yn + k3)
yn+\	= yn + l	111 fei + -k2 + -k3 + 77 k4 D            D D
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádu
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
23/47
Pokročilé metody
► Která metoda je nejlepší?
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
24/47
Pokročilé metody
PA081: Programování numerických
výpočtů
► Která metoda je nejlepší?
► jak kdy, záleží na konkrétním problému
► neplatí, že čím vyšší řád, tím delší krok si můžeme dovolit
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
24/47
Pokročilé metody
► Která metoda je nejlepší?
► jak kdy, záleží na konkrétním problému
► neplatí, že čím vyšší řád, tím delší krok si můžeme dovolit
► adaptivní délka kroku
► dlouhé kroky v plochých oblastech, krátké v bohatších
► zdvojnásobení kroku metody Runge-Kutta
► Fehlberguv algoritmus
► zpětnovazební metody
►
► Bulirsch-Stoer
► dlouhý krok s vnitřní extrapolací na nekonečně krátké mezikroky
► více viz Numerical Recipes
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
24/47
Výpočet integrálu
► neurčitý integrál
f (t) =
g(t)dt
odpovídá nalezení libovolného řešení / rovnice
df(t)/dt = g(t)
► určitý integrál
Ch
y =
g(t)dt = f(t1)-f(t0)
to
odpovídá konkrétnímu řešení téže rovnice s počáteční podmínkou f (to) = 0
► položíme yo = 0 a vhodnou numerickou metodou integrujeme od to do t\
Dynamika pevného tělesa
► hmotný bod je příliš zjednodušující abstrakce
► pevné těleso
► nenulový objem, nepodléhá žádným deformacím
► k popsání stavu systému jsou nutné další veličiny
► aktuální orientace (natočení)
► matice 3x3 nebo kvaternion
► úhlová rychlost
► vektor co
► velikost - rychlost rotace
► směr - osa rotace
► matice (tenzor) setrvačnosti
► analogie hmotnosti
► prostorové rozložení hmotnosti
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Matice setrvačnosti a úhlová hybnost
► lineární hybnost P = mv je relativně intuitivní veličina
► úhlová hybnost L = Mco už ne
► jednoduché vyjádření pohybové rovnice
dL(t)/dt = T(t)
► analogie dP(t)/dt = F(t)
► pro volně rotující těleso zůstává L konstantní (co ne)
► matice setrvačnosti
► vychází z T = F x r
► pro množinu diskrétních bodů
/ m(rl + r})
y
-mryrx
\ -mrzr:
X
-mrxry m{rl + rj) -mrzry
-mrxrz
\
-mrvrz m{rl +r|) )
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
► ve spojitém případě analogicky integrováním
► závisí na aktuální orientaci
► více viz http://www.cs.cmu.edu/-baraff/sigcourse/
27/47
Pohybové rovnice pevného tělesa
► stavový vektor y = (x, q, P, L)1
► x - poloha
► q - orientace
► P - lineární hybnost (vyjadřuje rychlost)
► L - úhlová hybnost (vyjadřuje úhlovou rychlost)
► rovnice
dy
dt
V    T /
► úhlovou rychlost co vypočteme
co = qOVrV^L))
► cog je vnoření co do prostoru kvaternionů
► integrace dopřednou nebo semiimplicitní metodou
► lineární rovnice o 13 proměnných
Model interakce molekul
PA081: Programování numerických
výpočtů
► dvě molekuly
► pevná, velká - receptor
► menší, pohyblivá - ligand
► pohyblivou molekulu ovládáme
► myš a klávesnice,
► zařízení silové zpětné vazby
► molekuly silově interagují
► přitažlivé i odpudivé síly
► každý atom s každým
► význam
► nalezení místa, kde si ligand „sedne"
► pochopení vlastností interakce v okolí
► simulace musí působit důvěryhodně
► co nej dokonalejší iluze manipulace s reálnými objekty
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
29/47
Zařízení silové zpětné vazby
► vyvinula se z teleoperátorů
► manipulace na dálku, např. v nebezpečném prostředí
► manipulace s virtuálním prostředím
► impedanční (odporové)
► zařízení snímá polohu, případně rychlost
► působí vypočtenou silou
► modelovaný systém je vyjádřen jako „odpor" - silová reakce systému na změnu polohy
► admitanční (vodivostní)
► zařízení snímá silové působení
► přesune se do vypočtené polohy
► modelovaný systém je vyjádřen jako „vodivost" - ochota systému reagovat na vnější silové působení
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Zařízení silové zpětné vazby
► v ideálním případě lze modelovat vše obojím
► vstupují nedokonalosti zařízení
► prakticky jsou impedanční zařízení lepší v modelu, kde převládá aktivita operátora
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Zařízení silové zpětné vazby
► v ideálním případě lze modelovat vše obojím
► vstupují nedokonalosti zařízení
► prakticky jsou impedanční zařízení lepší v modelu, kde převládá aktivita operátora
► hmat má výrazně nižší latenci než zrak
► nutná obnovovací frekvence 1 kHz
► vysoká náročnost na numerické vyhodnocení modelu
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Přímé mapování síly
► model je de facto statický
► vstup (souřadnice zařízení) jsou přímo promítnuty do modelu
► na jejich základě je vypočtena silová interakce
► síla je opět přímo aplikována na zařízení
► např. kolize dvou koulí, modelujeme tvrdou pružinou
F =
0
r > #1 + R
k(Rľ +R2-r)   r < Rľ + R2
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Přímé mapování síly
► model je de facto statický
► vstup (souřadnice zařízení) jsou přímo promítnuty do modelu
► na jejich základě je vypočtena silová interakce
► síla je opět přímo aplikována na zařízení
► např. kolize dvou koulí, modelujeme tvrdou pružinou
PA081: Programování numerických
výpočtů
F =
0
r > #1 + R
k(Rľ +R2-r)   r < Rľ + R2
► jednoduché, přehledné, ale takhle reálný svět nefunguje
► veškeré dynamické chovaní (tj. jaký je efekt aplikované sily) jsme z modelu vyloučili
► zejména nedokážeme modelovat různé hmotnosti kolidujících objektů
► dynamicky se chová fyzické zařízení atd., ale to nestačí
► nezatracujeme úplně, na některé aplikace stačí
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Virtuální spojení
► „virtual coupling"
► manipulovaný objekt jek zařízení připojen viskoelastickým spojem
► paralelní pružina a tlumič
F = feAx + bA\
► silová zpětná vazba
► jakou silou má zařízení působit na operátora
► vstup pro model
► jak působí operátor na manipulovaný virtuální objekt
► tj. 3. Newtonův zákon v praxi
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
33/47
Silové interakce molekul
► molekulu Ugandu modelujeme jako pevné těleso
► vzájemná poloha atomů se nemění
► hmotnosti atomů koncentrovány v jádrech — matice setrvačnosti
► silová interakce atomu vždy působí v místě jádra
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Silové interakce molekul
► molekulu Ugandu modelujeme jako pevné těleso
► vzájemná poloha atomů se nemění
► hmotnosti atomů koncentrovány v jádrech — matice setrvačnosti
► silová interakce atomu vždy působí v místě jádra
► van der Waalsova - LJ potenciál
12 A 6B
FLJ =--7^ +
r
13
r
7
► elektrostatická
Fei = e
r2
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Silové interakce molekul
► molekulu Ugandu modelujeme jako pevné těleso
► vzájemná poloha atomů se nemění
► hmotnosti atomů koncentrovány v jádrech — matice setrvačnosti
► silová interakce atomu vždy působí v místě jádra
► van der Waalsova - LJ potenciál
12 A 6B
FLJ =--7^ +
r
13
r
7
► elektrostatická
Fei = e
r2
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
► n ^ 104 atomů receptoru, m « 100 atomů Ugandu
► m x n interakcí nespočítáme 1000x za vteřinu
34/47
van der Waalsova interakce
► předvýpočet potenciálu na mřížce
► běžně používáno v chemickém software
► interpolace příliš vyhladí strmé stěny
5 4			
	-	1          i          1 i	
3 S 2 >	-\		
1		-	
0			
-e			
	i    l               i          , i		
		s  r -4»        1.5                   2 2.5 m          ' tg _	
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
35/47
van der Waalsova interakce
► díky vysokým exponentům relativně krátký dosah
► můžeme si dovolit ořez
► uvažujeme pouze dvojice atomů o vzájemné vzdálenosti menší než daná konstanta
► typicky 4-8 Á
► vypočítat m x n vzdáleností dvojic atomů je pořád příliš
van der Waalsova interakce
PA081: Programování numerických
výpočtů
► prostor rozdělíme na pravidelné buňky
► tak, aby se „největší" atom receptoru vešel do jedné buňky
► při náhodné pozici zasahuje nanejvýš do 8 buněk
► hashovací funkce souřadnic buňky
► struktura připravená před simulací
► projdeme všechny atomy receptoru
► atomy, které padnou (i částečně) do buňky zařadíme do příslušného řádku
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
van der Waalsova interakce
► prostor rozdělíme na pravidelné buňky
► tak, aby se „největší" atom receptoru vešel do jedné buňky
► při náhodné pozici zasahuje nanejvýš do 8 buněk
► hashovací funkce souřadnic buňky
► struktura připravená před simulací
► projdeme všechny atomy receptoru
► atomy, které padnou (i částečně) do buňky zařadíme do příslušného řádku
► atom Ugandu uměle zvětšíme o vzdálenost ořezu
► překryje « 10-500 buněk
► uvažujeme pouze atomy receptoru v příslušných řádcích
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
van der Waalsova interakce
PA081: Programování numerických
výpočtů
► i po ořezu zůstává příliš mnoho atomu receptoru
► silové působení (síla i točivý moment) jsou aditivní
► první derivace jsou také aditivní
► lze vše počítat paralelně a na závěr sečíst
► round-robin rozdělení atomů receptoru mezi jádra CPU
► statisticky rovnoměrné
► technicky používáme MPI
► musí běžet na jednom stroji nebo na velmi rychlé síti
► vysoké nároky na latenci
► rozeslání dat, výpočet, sesbírání výsledku za méně než lms
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Elektrostatická interakce
PA081: Programování numerických
výpočtů
► nízký exponent (2) - dlouhý dosah
► nelze efektivně použít ořez
► vypočítáme silové pole na mřížce
► působení na jednotkový náboj
► lze spočítat dopředu
► při běhu simulace interpolace pro všechny atomy Ugandu
► interpolace ořezává vysoké hodnoty
► v kombinaci s přesnou VdW interakcí to není problém
► extrémy jsou ve stejných místech
► VdW nedovolí tohoto stavu dosáhnout
► mřížka je příliš velká
► 2563 ■ 3 ■ 4B « 200MB
► adaptivní rozlišení mřížky
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Celková struktura simulace
► off-line výpočet
► elektrostatického pole
► hashovací tabulka pro VdW interakce
► řídící simulační vlákno (1 kHz)
► elektrostatické působení na všechny atomy ligandu
► součet dílčího silového působení VdW sil
► působení zařízení silové zpětné vazby na ligand (3 nebo 6 DoF)
► integrace pohybové rovnice (semiimplicitní, systém 13 lineárních rovnic)
► Nx výpočty dílčích VdW sil
► synchronně s řídícím vláknem
► pevně mezi sebe rozdělené atomy receptoru
► každý pro všechny atomy ligandu
► řízení zařízení silové zpětné vazby (1 kHz)
► asynchronní
► implementace viskoelastického spojení s ligandem
► grafika (25 Hz)
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
39/47
Pohyblivý receptor
► model tuhé molekuly neodpovídá realitě
► při nenulové teplotě mění svůj tvar spontánně
► při vzájemné interakci se ovlivňují a mění tvar
► spontánní pohyby receptoru umíme napočítat dopředu
► ~ 1000 a více snímků
► hash tabulka pro všechny snímky
► 1000 x 200 kB není problém
► neřešíme, i když umíme přepočítávat změny dostatečně rychle
► elektrostatické pole
► 1000 x 50 MB, kdo neplýtvá diskem s námi...
► komprese proměnlivého vektorového pole
PA081: Programování numerických
výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádu
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Domácí úkol
Scénář
► do blízkosti planety Jupiter vletí kometa a exploduje
► simulujeme trajektorii úlomků v gravitačním poli planety a čtyř měsíců
► detekujeme kolize úlomků s planetou a měsíci
► zjednodušující předpoklady
► vše jen ve 2D
► stabilní kruhový oběh všech měsíců
► úlomky na sebe vzájemně nepůsobí
► pro výpočet gravitační síly nahrazujeme tělesa bodem
► zanedbáváme vliv Slunce
► žádné relativistické efekty
Programování numerických výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Domácí úkol
Naivní implementace
► k dispozici v materiálech k předmětu
► pouze Eulerova dopředná metoda
► nepoužitelná, není dostatečně stabilní pro potřebný počet kroků
► parametry programu
► časový krok (10 s)
► počet úlomků (generovány deterministicky)
► maximální délka simulace
► metoda integrace
► možnost vypsat polohy všech těles (vizuální kontrola)
Programování numerických výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Domácí úkol
Grafické znázornění
PA081: Programování numerických
výpočtů
le+09
5e+08
-5e+Q8
-le+09
-1.5e+09
-2e+09 ■
-2.5e+09
'sim2' u 2:3 '■imZ' u 10:11 'sim2' u 12:13
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda Zpětná metoda
Semiimplicitní
metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
-le+09      -5e+08 0 5e+08       le+09      1.5e+09      2e+09 2,5e+0S
43/47
Domácí úkol
Zadání
► doimplementujte do programu stabilnější metody integrace
► hodnotí se i rychlost implementace
► základní optimalizace výpočtu
► vhodně volená metoda integrace, délka kroku apod.
► parametry spuštění
► počet úlomků, délka simulace
► tolerance na přesnost výsledku
Programování numerických výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Domácí úkol
Hodnocení
PA081: Programování numerických
výpočtů
► 1 bod za funkční implementaci (v toleranci přesnosti)
► 0.5 bodu za komentované optimalizační kroky vedoucí k využití alespoň 10 % teoretického výkonu jádra
► 0.5 za paralelní implementaci s efektivitou na 16 jádrech alespoň 70 %
► připočteno k výsledku zkoušky (písemka 12, na A třeba 9)
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Domácí úkol
Optimalizace
► základní
► použijte AVX instrukce (-mavx, případně -mfma)
► umožněte kompilátoru vektorizovat
► dbejte na využití cache
► double vs. float
► vyhněte se zbytečným výpočtům náročných funkcí (sqrt, sin, cos, ...)
► paralelní implementace
► pthreads, OpenMP, MPI, ...
► vhodné rozdělení úlomků mezi vlákna
► nevyužitá jádra po kolizích s měsíci
Programování numerických výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad
Domácí úkol
Odhad výkonu
► spočítat přesně operace je problematické
► netriviální úsilí, snadné zanesení chyby — falešný výsledek
► čitače v procesorech
► v Linuxu perf z balíku 1 i nux-tool s-generi c perf stat -e r5307c6 program argumenty ...
► magické r5307c6 funguje pro Haswell, více viz např. http://www.bni kol i c.co.uk/bl og/
hpc-howto-measure-flops.html
► ověřte si vztah hodnoty čitače ke skutečnému počtu provedených operací
► instrukce AVX - 8 operací ve float zároveň
► i při velmi dobrém využití cache se skutečně provede cca. 70 % instrukcí
Programování numerických výpočtů
Dynamické systémy
Diferenciální rovnice
Příklady
Dopředná metoda
Zpětná metoda
Semiimplicitní metoda
Leap frog
Metody vyšších řádů
Výpočet integrálu
Dynamika pevného tělesa
Komplexní příklad