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 2016 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 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 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 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 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á. 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 termodynamický zákon PA081: Programování numerických výpočtů ► 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č. 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 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 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 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 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 Dynamika pevného tělesa Komplexní příklad 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ů float k = 1.3, m = 0.2, p dt = 0.001; = 1, y = 0, t = 0, while (...) { float pi = yl = } p - k*y*dt, y + p/m*dt; t += dt; p = pl; y = yl; Dynamické systémy Diferenciální rovnice Příklady Dopředná metoda Zpětná metoda Semi implicitní metoda Leap frog Metody vyšších řádů Výpočet integrálu Dynamika pevného tělesa Komplexní příklad <□► « ď ► < 1 ► 4 1 ► 9/46 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 float u = 0, U = 14, R = 22e3, C = le-6, dt = 0.00001; whi1 e (...) { float dq = (U-u)/R * dt; u += dq/C; } ► stejně i pro proměnlivý vstup, např. U = sin(27T ■ 10~3 ■ t 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 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 5 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/46 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 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 ► řešení lze vyjádřit analyticky a naprogramovat Komplexní příklad 15/46 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/46 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/46 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/46 Zábo, skoč ► 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 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 20/46 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 ■ v 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 ► metoda 4. řádu ► posloupnost výpočtu fci = hf(xn,yn) k2 = 1 i hf(xn + ^yn + 2fel) 1 1 fe4 = hf(xn + h,3/n + fc3) !/ !, !, 1, o 3 3 o 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 23/46 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/46 Pokročilé metody PA081: Programovaní 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 5 ^^(V 24/46 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/46 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/46 Pohybové rovnice pevného tělesa ► stavový vektor y = (x, q, P, L)T ► x - poloha ► q - orientace ► P - lineární hybnost (vyjadruje rychlost) ► L - úhlová hybnost (vyjadruje úhlovou rychlost) ► rovnice / —P \ m dy F V úhlovou rychlost co vypočteme co = q(M-1q-1(L)) 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 ► coq je vnoření co do prostoru kvaternionů ► integrace dopřednou nebo semiimplicitní metodou ► lineární rovnice o 13 proměnných 28/46 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/46 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 PA081: Programování numerických výpočtů ► 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 ► 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/46 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/46 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/46 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 ligandu 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/46 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 řádů Výpočet integrálu Dynamika pevného tělesa Komplexní příklad 40/46 Komprese vektorového pole ► série snímků, každý 3D mřížka 2563 vektorů ► výrazná časová i prostorová spojitost ► algoritmus komprese vycházející z audia a videa ► sekvence snímků IBBBPBBBP apod. ► I - komprimovaný samostatně ► P - závisí na předchozím I nebo P ► B - interpolace mezi I-P resp. P-P 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 Komprese vektorového pole ► I-snímek ► analogie komprese audia ► každá složka zvlášť ► logaritmická transformace ► kvantizace na 16 bitů ► průchod 3D mřížkou definovaným způsobem ► lineární prediktor ► předpokládáme, že následující bod je stejný ► zaznamenáváme odchylku ► Riceovo kódování odchylek ► í = n div 2fc, j = n mod 2k ► unární zakódování í, následované 0 a fc bity j ► vhodné pro geometrické rozložení 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 Komprese vektorového pole ► B,P snímky ► opět Riceovo kódování odchylek ► adaptivní mřížka ► strukturu fixujeme pro frameset ► změna mezi framesety může způsobit nespojitosti ► zpětné sjednocení mřížek ( A ■ / + 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 Online dekomprese ► časově náročná komprese příliš nevadí (off-line) ► rychlá dekomprese je kritická ► ligand je malý - není třeba celý snímek ► aktualizace 1-10 Hz ► paralelní implementace ► sady procesů pro I a P/B snímky ► celý I s předstihem (závislosti mezi bloky) ► B/P jen potřebné bloky (podle polohy ligandu) 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 On-line dekomprese Snímky v čase Průběh dekomprese 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 Použiti dat ve snímku: Data pro další snímek: Proměnlivý ligand ► změny v tvaru důsledkem interakce s receptorem ► nelze napočítat dopředu ► řádově menší molekula ► popsat změny tvaru dalšími volnými proměnnými ► zahrnout do on-line simulace ► vypsaná diplomová práce 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