Univerzita Karlova v Praze Matematicko-fyzikální fakulta

Univerzita Karlova v Praze Matematicko-fyzik´aln´ı fakulta

´ RSK ˇ ´ PRACE ´ BAKALA A

Vojtˇech Miloˇs Vizualizace vypoˇ cten´ eho proudˇ en´ı pomoc´ı metody sledov´ an´ı ˇ c´ astic Katedra numerick´e matematiky

Vedouc´ı bakal´aˇrsk´e pr´ace: Studijn´ı program: Studijn´ı obor:

RNDr. Ing. Jaroslav Hron, Ph.D. Matematika Obecn´a matematika

Praha 2013

Na tomto m´ıstˇe bych velmi r´ad vˇrele podˇekoval ochotn´emu vedouc´ımu m´e pr´ace, RNDr. Ing. Jaroslav Hron, Ph.D., za veˇsker´ y ˇcas, rady a trpˇelivost, kterou pro mˇe naˇsel.

Prohlaˇsuji, ˇze jsem tuto bakal´aˇrskou pr´aci vypracoval samostatnˇe a v´ yhradnˇe s pouˇzit´ım citovan´ ych pramen˚ u, literatury a dalˇs´ıch odborn´ ych zdroj˚ u. Beru na vˇedom´ı, ˇze se na moji pr´aci vztahuj´ı pr´ava a povinnosti vypl´ yvaj´ıc´ı ze z´akona ˇc. 121/2000 Sb., autorsk´eho z´akona v platn´em znˇen´ı, zejm´ena skuteˇcnost, ˇze Univerzita Karlova v Praze m´a pr´avo na uzavˇren´ı licenˇcn´ı smlouvy o uˇzit´ı t´eto pr´ace jako ˇskoln´ıho d´ıla podle §60 odst. 1 autorsk´eho z´akona.

V . . . . . . . . . . . . . . . . . . . . . . dne . . . . . . . . . . . . .

Podpis autora

N´azev pr´ace: Vizualizace vypoˇcten´eho proudˇen´ı pomoc´ı metody sledov´an´ı ˇc´astic Autor: Vojtˇech Miloˇs Katedra: Katedra numerick´e matematiky Vedouc´ı bakal´aˇrsk´e pr´ace: RNDr. Ing. Jaroslav Hron, Ph.D., Katedra numerick´e matematiky Abstrakt: Na re´aln´ ych datech naskenovan´e c´evy pomoc´ı CT (poˇc´ıtaˇcov´e tomografie) ˇci MRI (magnetick´e rezonance) je vytvoˇrena s´ıt’ ze ˇctyˇrstˇen˚ u a ve vrcholech jsou vypoˇcteny vektory rychlosti pomoc´ı metody koneˇcn´ ych prvk˚ u. C´ılem pr´ace je navrhnout a implementovat algoritmus stopov´an´ı ˇc´astice v t´eto s´ıti s d˚ urazem na maxim´aln´ı pˇresnost. Za t´ım u ´ˇcelem uvedeme metodu analytick´e integrace, kter´a je moˇzn´a d´ıky po ˇca´stech line´arn´ımu vektorov´emu poli. Prezentujeme vlastn´ı modifikace t´eto metody a v´ ysledky porovn´ame s tradiˇcn´ım zp˚ usobem numerick´e integrace. D´ale je v pr´aci zm´ınˇena modifikace obecn´ ych krokov´ ych metod pro urˇcen´ı buˇ nky, v n´ıˇz ˇc´astice nach´az´ı, bez nutnosti prohled´an´ı cel´e s´ıtˇe po kaˇzd´em kroku. Kl´ıˇcov´a slova: sledov´an´ı ˇca´stice, lok´alnˇe exaktn´ı metoda, vizualizace proudˇen´ı tekutiny

Title: Visualization of computed flows using particle tracing methods Author: Vojtˇech Miloˇs Department: Department of Numerical Mathematics Supervisor: RNDr. Ing. Jaroslav Hron, Ph.D., Department of Numerical Mathematics Abstract: A mesh of tetrahedrons is build from a real data that were obtained by CT (computer tomography) or MRI (magnetic resonance imaging) scan of a blood vessel and then velocity vectors are computed in every vertex using the finite element method. The objective of the work is to design and implement a particle tracing algorithm in such a mesh with an emphasis on a maximum accuracy of the result. For this purpose, an analytic integration method is presented, which can be used thanks to piecewise linear vector field. We propose our own modification of this method and compare its results with the conventional numerical integration. Furthermore, the work describes modifications of common step methods that are used to determine the cell that contains a particular particle without the need to search the whole mesh in each step. Keywords: particle tracing, local exact method, visualization of fluid flow

Obsah ´ Uvod

2

1 Z´ akladn´ı myˇ slenka sledov´ an´ı ˇ c´ astice 1.1 Analyticky zadan´e vektorov´e pole . . . . . . . . . . . . . . . . . . 1.2 Vektorov´e pole vznikl´e interpolac´ı . . . . . . . . . . . . . . . . . .

3 5 6

2 Pouˇ zit´ e metody a jejich srovn´ an´ı 2.1 Optimalizovan´a Eulerova metoda . . . . . . . . . . . 2.1.1 Z´akladn´ı varianta stopov´an´ı ˇca´stice . . . . . . 2.1.2 Modifikovan´a varianta stopov´an´ı ˇca´stice . . . 2.2 Lok´alnˇe exaktn´ı metoda . . . . . . . . . . . . . . . . 2.2.1 Explicitn´ı vyj´adˇren´ı trajektorie v r´amci buˇ nky 2.2.2 V´ ypoˇcet u ´nikov´eho bodu z buˇ nky . . . . . . . 2.2.3 Moˇzn´e typy bunˇek . . . . . . . . . . . . . . . 2.3 Srovn´an´ı metod . . . . . . . . . . . . . . . . . . . . . 2.3.1 Pˇresnost metody . . . . . . . . . . . . . . . . 2.3.2 Ostatn´ı krit´eria . . . . . . . . . . . . . . . . . 2.4 Program . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

8 8 8 10 11 11 17 21 22 22 23 26

Z´ avˇ er

27

Apendix 2.5 Struktura programu a pouˇzit´e algoritmy 2.5.1 TTetBasicVol . . . . . . . . . . . 2.5.2 TracingEuler . . . . . . . . . . . 2.5.3 TracingLocalExact . . . . . . . . 2.5.4 TiskCastice . . . . . . . . . . . . 2.5.5 Main . . . . . . . . . . . . . . . . 2.6 Uˇzivatelsk´a pˇr´ıruˇcka . . . . . . . . . . .

28 28 28 30 30 31 31 31

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

Literatura

32

Seznam obr´ azk˚ u

33

1

´ Uvod C´ılem t´eto pr´ace je pˇresn´ y v´ ypoˇcet a vizualizace proudˇen´ı tekutiny v jist´e oblasti z´ajmu. Potˇrebu zpr˚ uhlednit chov´an´ı tekutiny maj´ı napˇr´ıklad leteˇct´ı inˇzen´ yˇri pˇri konstrukci Boeingu ˇci designeˇri v automobilk´ach. Zde jsou aplikov´any klasick´e simulaˇcn´ı algoritmy, nebot’ je nesrovnatelnˇe levnˇejˇs´ı simulovat proudˇen´ı vzduchu kolem modelu nov´eho auta na poˇc´ıtaˇci neˇz vytvoˇrit skuteˇcnou maketu vozidla, uv´est do chodu cel´ y aerodynamick´ y tunel a zamˇestnat spoustu lid´ı. Nakonec se mus´ı pˇrikroˇcit k n´akladn´emu testu v tunelu pro potvrzen´ı spr´avnosti simulac´ı, avˇsak pˇredchoz´ı zkouˇsky mohou prob´ıhat jen virtu´alnˇe, viz [12]. Nicm´enˇe z´ajem o grafick´e zn´azornˇen´ı proudu tekutiny sah´a mnohem d´al, napˇr´ıklad k l´ekaˇrsk´e praxi pˇri zkoum´an´ı cirkulace krve v krevn´ım ˇreˇciˇsti. M´ali pacient aneurysma (v´ ydut’ dut´eho org´anu, zpravidla tepny), jeˇstˇe nen´ı jist´e, zda by se mˇelo pˇristoupit k operaci. Rozhodnut´ı velmi citlivˇe z´aleˇz´ı na charakteru proudˇen´ı v dutinˇe, coˇz pˇredstavuje pr´avˇe motivaci ke vzniku t´eto pr´ace, jej´ıˇz n´apln´ı je navrhnout a implementovat algoritmus pro simulaci s d˚ urazem na maxim´aln´ı pˇresnost. V prvn´ı kapitole je nast´ınˇen proces z´ısk´av´an´ı dat ze skuteˇcn´eho svˇeta a n´asledn´ y zp˚ usob uchov´an´ı pro dalˇs´ı pr´aci. Jedn´a se o sezn´amen´ı s ˇsirˇs´ı t´ematikou, kter´a poskytne z´aklad pro pochopen´ı myˇslenky a z´ajmu sledov´an´ı ˇc´astic. D´ale se uˇz vˇenujeme pouze zobrazen´ı jiˇz zpracovan´ ych dat bez ohledu na jejich p˚ uvod. Druh´a kapitola nab´ız´ı konkr´etn´ı koncepty algoritm˚ u pro k´ yˇzenou simulaci a jejich srovn´an´ı pˇredevˇs´ım ve velikosti chyby.

Obr´azek 1: Uk´azka simulace obt´ek´an´ı formule vzduchem. 2

Kapitola 1 Z´ akladn´ı myˇ slenka sledov´ an´ı ˇ c´ astice Tato kapitola nab´ız´ı u ´vod do cel´e problematiky sledov´an´ı ˇc´astice. Za t´ım u ´ˇcelem zm´ın´ıme nˇekolik uˇz´ıvan´ ych pojm˚ u, kter´e ˇcerp´ame z knihy Mechanika kontinua [3]. Dom´ ena Ω je oblast z´ajmu vyplnˇen´a tekutinou. M˚ uˇze se jednat napˇr´ıklad o vnitˇrek potrub´ı ˇci prostor okolo leteck´eho profilu pˇri zkoum´an´ı vztlaku kˇr´ıdla. Proudˇen´ı tekutiny zcela popisuje vektorov´ e pole ~v = ~v (x,t), kter´e v kaˇzd´em bodˇe x ∈ Ω ud´av´a vektor rychlosti pro kaˇzd´ y ˇcas t. Existuje nˇekolik zp˚ usob˚ u - kˇrivek, jak proudˇen´ı zviditelnit. • Proudnice (Streamline) je kˇrivka v kaˇzd´em bodˇe teˇcn´a ke smˇeru rychlosti. • Emisn´ı ˇ c´ ara (Streakline) pˇredstavuje stopu barviva vpouˇstˇen´eho do dom´eny v jednom bodˇe. • Trajektorie (Pathline) je stopa nehmotn´e ˇc´astice un´aˇsen´e proudem tekutiny. R´adi bychom poznamenali, ˇze proudnice z´avis´ı pouze na aktu´aln´ım vektorov´eho pole, zat´ımco emisn´ı ˇca´ry i trajektorie jsou ovlivnˇeny jeho ˇcasov´ ym pr˚ ubˇehem. Pro intuitivn´ı pˇredstavu jsou tyto charakteristiky kˇrivek jasn´e, ovˇsem pro u ´plnost uvedeme jeˇstˇe jejich pˇresn´e definice.

3

Obr´azek 1.1: Demonstrace rozd´ıl˚ u mezi jednotliv´ ymi kˇrivkami popisuj´ıc´ımi vektorov´e pole v pr˚ ubˇehu ˇcasu. Kr´atk´e ˇc´arky pˇredstavuj´ı smˇer vektorov´eho pole. ˇ Cerven´ a kˇrivka je trajektorie, modr´a kˇrivka je emisn´ı ˇca´ra a pˇreruˇsovan´e ˇca´ry ukazuj´ı proudnice.

4

Definice 1. Necht’ Ω je otevˇren´a souvisl´a podmnoˇzina R3 , 0 < T ∈ R a ~v ∈ R3 je spojit´e vektorov´e pole na Ω × (0, T ). Pak parametrick´a kˇrivka ~xp : [a,b] → Ω pro a < b ∈ R se naz´yv´a proudnice v ˇcase t ∈ (0,T ), jestliˇze plat´ı ∀s ∈ [a,b]

d~xp (s) × ~v (~xp (s),t) = 0, ds

kde ~u × ~v je vektorov´y souˇcin vektor˚ u ~u a ~v . Definice 2. Necht’ Ω je otevˇren´a souvisl´a podmnoˇzina R3 a x0 ∈ Ω. Necht’ 0 < T ∈ R a ~v ∈ R3 je spojit´e vektorov´e pole na Ω × (0, T ). Necht’ ~xp : (0,T ) → Ω je parametrick´a kˇrivka, kter´a splˇ nuje d~xp (t) = ~v (~xp (t),t). dt Parametrick´a kˇrivka ~xp (t) se naz´yv´a trajektorie proch´azej´ıc´ı bodem x0 , jestliˇze existuje t0 ∈ (0,T ), ˇze je splnˇena poˇc´ateˇcn´ı podm´ınka ∀t ∈ (0, T )

~xp (t0 ) = x0 . Definice 3. Necht’ Ω je otevˇren´a souvisl´a podmnoˇzina R3 , 0 < T ∈ R a ~v ∈ R3 je spojit´e vektorov´e pole na Ω×(0, T ). Necht’ φ(t,x0 ,t0 ) je trajektorie parametrizovan´ a parametrem t ∈ (0,T ) proch´azej´ıc´ı bodem x0 ∈ Ω v ˇcase t0 ∈ (0,T ). Necht’ 0 < a < b < T . Pak parametrick´a kˇrivka ψ(s) : (a,b) → Ω se naz´yv´a emisn´ı ˇc´ara v ˇcase t, pokud ∀s ∈ (a,b) ψ(s) = φ(t,x0 ,s). V pˇr´ıpadˇe stacion´arn´ıho (nemˇenn´eho v ˇcase) proudˇen´ı tyto tˇri kˇrivky spl´ yvaj´ı v jednu. To je uˇziteˇcn´e pozorov´an´ı z toho d˚ uvodu, ˇze se trajektorie daj´ı experiment´alnˇe zjiˇst’ovat nalezen´ım emisn´ıch ˇcar, kter´e jsou empiricky mnohem l´epe pozorovateln´e. D´ale se budeme jiˇz zab´ yvat pouze zkoum´an´ım trajektori´ı ˇca´stic.

1.1

Analyticky zadan´ e vektorov´ e pole

Jak je vidˇet pˇr´ımo z definice, zjiˇst’ovat trajektorie ˇca´stic znamen´a ˇreˇsit soustavu obyˇcejn´ ych diferenci´aln´ıch rovnic. Proto se trajektorie naz´ yvaj´ı integr´ aln´ı ’ kˇ rivky. Nejjednoduˇsˇs´ı zp˚ usob numerick´e integrace, zde myˇsleno jako zjiˇst ov´an´ı trajektorie ˇca´stice pˇri znalosti jej´ı poˇca´teˇcn´ı pozice a rychlosti v kaˇzd´em bodˇe, pˇredstavuje Eulerova metoda. Spoˇc´ıv´a v rozdˇelen´ım ˇcasu na mal´e u ´seky ∆t (tzv. ˇcasov´e kroky) a nach´az´ı-li se ˇca´stice v ˇcase t0 v bodˇe x0 ∈ Ω, nov´a pozice ˇca´stice se vypoˇcte takto xn+1 = xn + ~v (xn ,tn )∆t.

(1.1)

Toto se naz´ yv´a krok metody. V pr´avˇe vypoˇcten´e pozici ˇc´astice x1 se opˇet zjist´ı rychlost (nyn´ı jiˇz v ˇcas t0 + ∆t) a stejn´ ym zp˚ usobem se vypoˇcte nov´a pozice x2 . Takto postupnˇe vznikne posloupnost bod˚ u aproximuj´ıc´ı integr´aln´ı kˇrivku. Chyba zp˚ usoben´a jedn´ım krokem odpov´ıd´a velikosti O(∆t2 ) a d´ıky diskretizaci T krok˚ u je chyba metody O(∆t). Volbou menˇs´ıho O(∆t) pak zajist´ıme na n = ∆t pˇresnˇejˇs´ı v´ ypoˇcet simulace (menˇs´ı odchylku aproximaˇcn´ı metody od analytick´eho 5

ˇreˇsen´ı soustavy obyˇcejn´ ych diferenci´aln´ıch rovnic). Bohuˇzel je tato metoda pˇr´ıliˇs pomal´a v tom smyslu, ˇze pro dosaˇzen´ı pˇrijateln´e pˇresnosti vyˇzaduje vˇetˇsinou nepˇrijateln´ y v´ ypoˇcetn´ı ˇcas. Pˇri dostateˇcn´e hladkosti vektorov´eho pole ~v lze pouˇz´ıt v´ıce metod pro numerickou integraci tˇechto obyˇcejn´ ych diferenci´aln´ıch rovnic, viz [4], napˇr´ıklad metody Runge-Kutta r˚ uzn´ ych ˇra´d˚ u nebo v´ıcekrokov´e metody.

1.2

Vektorov´ e pole vznikl´ e interpolac´ı

Pˇri skuteˇcn´ ych datech z experimentu nikdy nejsme schopni dostat u ´plnou informaci o vektorov´em poli rychlost´ı v dan´e dom´enˇe. Nast´ın´ıme, na jak´ ych datech pracuje algoritmus sledov´an´ı ˇc´astice, a heuristicky pˇribl´ıˇz´ıme proces vzniku takov´ ych dat. Konkr´etnˇe se budeme vˇenovat pˇr´ıpadu simulace proudˇen´ı krve v c´evˇe ˇclovˇeka, viz [1]. Prvn´ı v´ ysledky experimentu jsou sn´ımky napˇr´ıklad z CT (poˇc´ıtaˇcov´a tomografie), kter´e poskytnou 3D model c´evy. Existuje zp˚ usob sledov´an´ı ˇc´astic kontrastn´ı l´atky pr´avˇe pomoc´ı CT, dokonce i jejich rychlost´ı, nicm´enˇe u ´skal´ı t´eto metody spoˇc´ıv´a v samotn´em procesu z´ısk´av´an´ı sn´ımku. Poˇc´ıtaˇcov´a tomografie skenuje k´ yˇzen´ y objekt pˇr´ıliˇs pomalu, tedy nen´ı schopna zachytit dynamick´e jevy dostateˇcnˇe pˇresnˇe. Umoˇzn´ı uspokojivˇe zobrazit tvar c´evy, avˇsak charakter proudˇen´ı krve nikoliv. Na tomto m´ıstˇe je moˇzn´e pouˇz´ıt metodu koneˇ cn´ ych prvk˚ u (FEM - finite element method), kter´a vyˇzaduje rozdˇelen´ı objemu c´evy na mnoho ˇctyˇrstˇen˚ u1 , ˇc´ımˇz vznikne tzv. s´ıt’ (mesh). T´ımto je pˇripravena p˚ uda pro FEM, podrobnˇe popsanou v knize [2], pro v´ ypoˇcet rychlost´ı ve vrcholech s´ıtˇe. V´ ypoˇcet proudˇen´ı metodou FEM je v souˇcasn´e dobˇe rozs´ahl´ y obor tzv. CFD (computational fluid dynamics), kter´ y se neust´ale zkoum´a. Jedn´a se o diskutovan´ y probl´em bez jednoznaˇcn´eho ˇreˇsen´ı2 . Nyn´ı pˇredpokl´adejme, ˇze jiˇz um´ıme urˇcit hodnoty ~v ve vrcholech s´ıtˇe. Nejpˇrirozenˇejˇs´ı cestou rozˇs´ıˇren´ı informace o rychlosti je line´arn´ı interpolace, ˇc´ımˇz z´ısk´ame po ˇc´astech line´arn´ı vektorov´e pole ~v v cel´e dom´enˇe Ω. Abychom dostali co nejpˇresnˇejˇs´ı informaci o proudˇen´ı tekutiny v jist´em kritick´em m´ıstˇe (typicky u ostr´ ych hran, pˇredpokl´adan´ ych vysok´ ych rychlost´ı, turbulenc´ı apod.), je nutn´e v jeho okol´ı m´ıt ˇ jemnou s´ıt’ ˇctyˇrstˇen˚ u. Casto se st´av´a, ˇze se dodateˇcnˇe dle vypoˇcten´ ych v´ ysledk˚ u ’ mˇen´ı s´ıt podle lok´aln´ıho z´ajmu a znovu se poˇc´ıtaj´ı rychlosti kv˚ uli novˇe vznikl´ ym vrchol˚ um.

Pro zm´ınˇen´e metody numerick´e integrace je nutn´e volit d´elku ˇcasov´eho kroku. V´ıcekrokov´e metody maj´ı nev´ yhodu, ˇze mus´ı m´ıt d´elku kroku konstantn´ı. Naproti tomu jednokrokov´e metody (napˇr. Runge-Kutta) mohou flexibilnˇe reagovat na to, zda-li dan´ y krok prov´adˇej´ı v numericky zaj´ımav´e oblasti, ˇci nikoliv. Metody, kter´e zkracuj´ı nebo prodluˇzuj´ı sv˚ uj ˇcasov´ y krok dle potˇreby, se naz´ yvaj´ı adap1 Metoda koneˇcn´ ych prvk˚ u v principu nevyˇzaduje, aby z´akladn´ı buˇ nkou s´ıtˇe byl pr´avˇe ˇctyˇrstˇen. Pro v´ ypoˇcet a dalˇs´ı pr´aci je to efektivn´ı tvar. Pˇri pouˇzit´ı libovoln´eho mnohostˇenu napˇr´ıklad jiˇz nelze obecnˇe uvnitˇr interpolovat vektorov´e pole line´arnˇe (za pˇredpokladu zadan´ ych rychlost´ı ve vrcholech). 2 Pot´ıˇz je v hled´ an´ı jednoznaˇcn´eho ˇreˇsen´ı Navierov´ ych-Stokesov´ ych rovnic.

6

tivn´ı. Jen nen´ı zcela jasn´e, jak takov´e metody konstruovat. Intuice napov´ıd´a volit d´elku kroku podle velikosti buˇ nky ˇc´astic´ı pr´avˇe okupovanou. Ovˇsem tak´e z´aleˇz´ı ’ na velikosti rychlosti, nebot m˚ uˇze b´ yt pˇr´ıliˇs vysok´a a n´asleduj´ıc´ı vypoˇcten´a pozice ˇca´stice m˚ uˇze vyj´ıt mimo dom´enu. Dokonce by mˇel b´ yt br´an zˇretel na velikost zmˇeny vektorov´eho pole v dan´em bodˇe, ponˇevadˇz u t´emˇeˇr konstantn´ıho pole nen´ı tˇreba volit pˇr´ıliˇs mal´e ˇcasov´e kroky. Velkou v´ yhodou lok´alnˇe exaktn´ı metody, kterou uvedeme v druh´e kapitole, je adaptivita vypl´ yvaj´ıc´ı pˇr´ımo z definice, aniˇz by uˇzivatel musel volit d´elku kroku.

7

Kapitola 2 Pouˇ zit´ e metody a jejich srovn´ an´ı Hlavn´ım objektem naˇseho z´ajmu je lok´alnˇe exaktn´ı metoda uveden´a v ˇcl´anku [10]. Jedn´a se o ne pˇr´ıliˇs rozˇs´ıˇrenou metodu sledov´an´ı ˇca´stic, aˇc apriori nab´ız´ı nˇekolik v´ yhod oproti klasick´ ym krokov´ ym metod´am. C´ılem kapitoly je implementovat tuto metodu s nˇekolika modifikacemi podle ˇcl´anku [7] a srovnat v´ ysledky s klasickou jednokrokovou metodou. Zde pouˇzijeme Eulerovu metodu (viz pˇredchoz´ı kapitola) s optimalizac´ı lokalizace ˇca´stice podle ˇcl´anku [9], kter´a je detailnˇe probr´ana v prvn´ı ˇca´sti t´eto kapitoly. V jej´ım z´avˇeru uv´ad´ıme z´akladn´ı strukturu programu pro vizualizaci. Vˇenujeme se pouze z´asadn´ım rys˚ um, implementaˇcn´ı detaily jsou diskutov´any v apendixu.

2.1

Optimalizovan´ a Eulerova metoda

Klasick´a Eulerova metoda se implementuje n´asledovnˇe. Pro danou pozici ˇc´astice x~0 je nalezena buˇ nka (v naˇsem pˇr´ıpadˇe se jedn´a o ˇctyˇrstˇen) s´ıtˇe, ve kter´e se ˇca´stice nach´az´ı. Pot´e je zjiˇstˇena hodnota rychlostn´ıho vektorov´eho pole v m´ıstˇe ˇca´stice ~v (x~0 ). Spoˇcte se nov´a pozice ˇca´stice ze vztahu 1.1, kde ∆t je zvolen´ y ˇcasov´ y krok. V´ yhodou metody je jej´ı snadn´a implementace a ˇz´adn´e pˇredpoˇc´ıt´av´an´ı dat. Ovˇsem z´asadn´ım nedostatkem je ˇcasov´a n´aroˇcnost, a to obzvl´aˇst’ kv˚ uli nutnosti prohled´an´ı cel´e s´ıtˇe pro urˇcen´ı okupovan´e buˇ nky. T´ımto probl´emem se zab´ yv´a ˇcl´anek [9] a navrhuje algoritmus vyˇzaduj´ıc´ı glob´aln´ı vyhled´av´an´ı pouze jednou na zaˇca´tku simulace. D´ale je d´ıky pˇredpoˇc´ıtan´e konektivitˇe (sousednosti) bunˇek zjiˇst’ov´ana c´ılov´a buˇ nka pˇr´ımo, tento jev naz´ yv´ame stopov´an´ı ˇc´astice“. ”

2.1.1

Z´ akladn´ı varianta stopov´ an´ı ˇ c´ astice

Necht’ A je v´ ychoz´ı a B je c´ılov´a pozice ˇca´stice. Necht’ fi , i = 1,2,3,4, jsou ´ stˇeny pr´avˇe okupovan´e buˇ nky. Ukolem stopovac´ıho algoritmu je urˇcit, kterou stˇenou ˇca´stice opust´ı buˇ nku. Necht’ u ´seˇcka AB prot´ın´a stˇenu fi pro nˇejak´e i v bodˇe P . Pak P = A + λA (B − A), kde λA ∈ [0,1]. K´eˇz C je stˇred stˇeny fi a ~n jej´ı norm´alov´ y vektor. Pak vektor (P − C) n´aleˇz´ı do roviny urˇcen´e norm´alou ~n, tedy plat´ı 8

Obr´azek 2.1: Zn´azornˇen´ı situace pro 2D pˇr´ıpad s troj´ uheln´ıkovou s´ıt´ı. V´ ychoz´ı pozice ˇca´stice A se nach´az´ı v buˇ nce 1 a c´ılem je zjistit, ve kter´e buˇ nce se nach´az´ı c´ılov´a pozice B. Pro stˇeny f1 a f2 λA ∈ [0,1], kdeˇzto pro stˇenu f3 je λA < 0. Menˇs´ı λA n´aleˇz´ı stˇenˇe f1 , a proto je ˇca´stice posunuta do bodu P1 a aktu´aln´ı buˇ nka zmˇenˇena na buˇ nku 2. Analogick´ ym postupem se ˇc´astice posune do bodu P2 s okupovanou buˇ nkou 5. Zde ˇza´dn´a hodnota λA nen´aleˇz´ı intervalu [0,1], ˇc´ımˇz je algoritmus u konce.

(P − C) · ~n = 0, coˇz po dosazen´ı za P d´av´a λA =

(C − A) · ~n . (B − A) · ~n

(2.1)

Pozn´amka. Pro urˇcen´ı roviny stˇeny nen´ı nutn´e br´at pr´avˇe stˇred. V algoritmu staˇc´ı vz´ıt libovoln´ y vrchol stˇeny. Takto urˇc´ıme λA pro kaˇzdou stˇenu fi a pokud n´aleˇz´ı v intervalu [0,1], znamen´a to, ˇze u ´seˇcka AB prot´ın´a rovinu urˇcenou stˇenou fi . To pro urˇcen´ı u ´nikov´e stˇeny jeˇstˇe nestaˇc´ı, z aspiruj´ıc´ıch stˇen vybereme tu s nejmenˇs´ım λA , coˇz geometricky odpov´ıd´a rovinˇe, jeˇz u ´seˇcka AB protne v nejniˇzˇs´ım ˇcase. Pot´e d´ıky konektivitˇe mˇr´ıˇzky zjist´ıme, do kter´e buˇ nky se sledovan´a ˇc´astice pˇrem´ıstila. Za pozici A dosad´ıme pr˚ useˇc´ık trajektorie s pˇr´ısluˇsnou stˇenou a algoritmus opakujeme pro aktualizovanou v´ ychoz´ı pozici a okupovanou buˇ nku. Je tˇreba si d´at pozor na to, ˇze pr´avˇe dosaˇzen´a stˇena jiˇz nesm´ı b´ yt uvaˇzov´ana. Jestliˇze j´ı u ´seˇcka AB jiˇz jednou protnula, nem˚ uˇze ji protnout znovu. Z d˚ uvodu zaokrouhlovac´ı chyby m˚ uˇze nastat pˇr´ıpad, kdy algoritmus identifikuje v´ ychoz´ı stˇenu jako c´ılovou. To znamen´a, ˇze λA je t´emˇeˇr rovn´e nule, coˇz je bohuˇzel validn´ı v´ ysledek v intervalu [0,1]. Nazvˇeme tento probl´em opakovan´e prostupov´an´ı stˇen“. ” 9

Pokud pro ˇz´adnou stˇenu fi hodnota λA nen´aleˇz´ı do intervalu [0,1], mus´ı nutnˇe c´ılov´ y bod B n´aleˇzet stejn´e buˇ nce jako bod A. T´ım p´adem je algoritmus u konce. Probl´emem z´akladn´ı metody stopov´an´ı nen´ı ani tak nutnost vynech´av´an´ı jiˇz navˇst´ıven´e stˇeny, jako je to pot´ıˇz se zaokrouhlovac´ımi chybami a kompatibilitou pozice ˇca´stice s okupovanou buˇ nkou. Konkr´etnˇe mohou nastat dva pˇr´ıpady, kdy ˇca´stice fyzicky n´aleˇz´ı jin´e buˇ nce, neˇz je j´ı okupovan´a buˇ nka uloˇzen´a v programu. Uvaˇzujme soused´ıc´ı buˇ nky 1 a 2 a jako v´ yˇse body A a B, pˇriˇcemˇz bod A m´a pˇriˇrazenou buˇ nku 1. • A i B n´aleˇz´ı 2. Algoritmus tedy spoˇc´ıt´a pro kaˇzdou stˇenu buˇ nky 1 λA , jenˇze vˇzdy mimo interval [0,1]. Tedy prohl´as´ı, ˇze B tak´e n´aleˇz´ı buˇ nce 1. • A n´aleˇz´ı 2 a B n´aleˇz´ı 1. Pro spoleˇcnou stˇenu bunˇek vyjde λA v intervalu [0,1] a tedy je ˇca´stice pˇresunuta na pr˚ useˇc´ık stˇeny a u ´seˇcky AB za zmˇeny okupovan´e buˇ nky na 2. Snadno tedy m˚ uˇze doj´ıt ke ztr´atˇe informace o okupovan´e buˇ nce a bylo by nutn´e prov´est prohled´an´ı cel´e s´ıtˇe pro obnoven´ı kompatibility polohy ˇca´stice a buˇ nky. Bohuˇzel by se jednalo o glob´aln´ı pr˚ uzkum po kaˇzd´em kroku pro opravu pˇr´ıpadn´e chyby, coˇz je pˇr´ıliˇs n´aroˇcn´ y proces. Proto je tˇreba algoritmus modifikovat.

2.1.2

Modifikovan´ a varianta stopov´ an´ı ˇ c´ astice

Zm´ınˇen´e probl´emy z´akladn´ıho stopovac´ıho algoritmu vyˇreˇs´ı v´ ypoˇcet dalˇs´ıho charakteristick´eho ˇc´ısla pro kaˇzdou stˇenu λStred . Jedn´a se o analogickou definici ke vztahu 2.1, pouze nam´ısto v´ ychoz´ı pozice A je pouˇzit stˇred buˇ nky S. λStred =

(C − S) · ~n . (B − S) · ~n

N´asleduj´ıc´ı tvrzen´ı hovoˇr´ı o relevanci z´amˇeny S za A. Tvrzen´ı 1. Necht’ n ∈ Z a M je konvexn´ı n-stˇen v R3 . Jeho stˇeny oznaˇcme F1 , ...,Fn . Necht’ S je geometrick´y stˇred M a bod A ∈ Mo , B ∈ (M o )C 1 . Pak plat´ı ∀i ∈ 1, ...,n

(AB ∩ Fi 6= ∅) ⇔ (SB ∩ Fi 6= ∅).

D˚ ukaz. Necht’ i ∈ 1, ...,n je libovoln´e. Pak nadrovina urˇcen´a stˇenou Fi rozdˇeluje prostor R3 na dva otevˇren´e poloprostory. Oznaˇcme H1 ten, ve kter´em leˇz´ı cel´ y vnitˇrek n-stˇenu M (to lze d´ıky konvexitˇe M). Druh´ y poloprostor necht’ je H2 . Plat´ı tedy, ˇze A ∈ H1 a S ∈ H1 . Pak jiˇz trivi´alnˇe plyne AB ∩ Fi 6= ∅

⇔

B ∈ H2

⇔

SB ∩ Fi 6= ∅.

k 1

M o znaˇc´ı vnitˇrek mnoˇziny M , M C znaˇc´ı doplnˇek mnoˇziny M a M znaˇc´ı uz´avˇer mnoˇziny

M.

10

Nebot’ pˇr´ım´ ym d˚ usledkem tvrzen´ı je λStred ∈ [0,1] pr´avˇe tehdy, kdyˇz λA ∈ [0,1], pro urˇcen´ı podezˇrel´ ych stˇen u ´niku je moˇzn´e pro kaˇzdou stˇenu ˇctyˇrstˇenu spoˇc´ıtat pouze λStred . Nakonec definujme λm = min(1, max(0,λA )). Nyn´ı uk´aˇzeme, co je motivac´ı definice λm a v´ ypoˇctu λStred . Opˇet tedy uvaˇzujme soused´ıc´ı buˇ nky 1 a 2 a body A a B, pˇriˇcemˇz bod A m´a pˇriˇrazenou buˇ nku 1. • A i B n´aleˇz´ı 2. Algoritmus tedy pro kaˇzdou stˇenu buˇ nky 1 spoˇc´ıt´a λStred , a pro spoleˇcnou stˇenu obou bunˇek vyjde λStred ∈ [0,1], coˇz pˇredznamen´a, ˇze u ´seˇcka AB mus´ı protnout spoleˇcnou stˇenu. Nyn´ı se tedy spoˇc´ıt´a λA uˇz jen pro tuto stˇenu a to vyjde mimo interval [0,1]. Nicm´enˇe pokud vezmeme m´ısto λA hodnotu λm , pˇresune se ˇca´stice z bodu A pˇr´ımo do bodu B. Buˇ nka okupace je zmˇenˇena na 2. • A n´aleˇz´ı 2 a B n´aleˇz´ı 1. Pro vˇsechny stˇeny buˇ nky 1 vyjde λStred mimo interval [0,1], a tedy algoritmus konˇc´ı s v´ ysledkem, ˇze B n´aleˇz´ı buˇ nce 1. V obou pˇr´ıpadech modifikovan´ y algoritmus d´a spr´avn´ y v´ ysledek okupovan´e buˇ nky. Pro u ´plnost dod´av´ame, ˇze s modifikovanou metodou tak´e miz´ı povinnost hl´ıdat probl´em opakovan´e prostupov´an´ı stˇen“ zm´ınˇen´ y v´ yˇse. ” Na z´avˇer jeˇstˇe poznamen´ame, ˇze modifikovan´ y algoritmus striktnˇe vyˇzaduje konvexn´ı mnohostˇeny. Pokud se v s´ıti vyskytuj´ı nekonvexn´ı mnohostˇeny, je tˇreba je rozdˇelit na nˇekolik menˇs´ıch konvexn´ıch.

2.2

Lok´ alnˇ e exaktn´ı metoda

Probl´em zobrazov´an´ı trajektorie ˇca´stice v rychlostn´ım vektorov´em poli ~v je ekvivalentn´ı probl´emu ˇreˇsen´ı soustavy obyˇcejn´ ych diferenci´aln´ıch rovnic (d´ale jen soustavy ODR) d~x(t) = ~v (~x(t)), ~x(t0 ) = ~x0 . (2.2) dt Lok´alnˇe exaktn´ı metoda se liˇs´ı od klasick´ ych metod t´ım, ˇze navrhuje ˇreˇsit 2.2 analyticky a nal´ezt tak explicitn´ı vyj´adˇren´ı trajektorie s nulovou akumulovanou chybou.

2.2.1

Explicitn´ı vyj´ adˇ ren´ı trajektorie v r´ amci buˇ nky

Uvaˇzujme dom´enu Ω a v n´ı s´ıt’ ˇctyˇrstˇen˚ u s vektory rychlosti ve vrcholech bunˇek. Pak v kaˇzd´e buˇ nce m˚ uˇzeme rychlosti line´arnˇe interpolovat, ˇc´ımˇz z´ısk´ame po ˇca´stech line´arn´ı spojit´e vektorov´e pole ~v v Ω. Uvaˇzujme ˇca´stici nach´azej´ıc´ı se v urˇcit´e buˇ nce. Existuje interpolaˇcn´ı matice2 A ∈ R3×3 a vektor posunut´ı ~o ∈ R3 , ˇze pro vˇsechny ~x n´aleˇz´ıc´ı t´eto buˇ nce plat´ı ~v (~x) = A~x + ~o, 2

Pˇredpokl´ ad´ ame, ˇze matice A je regul´arn´ı.

11

coˇz spolu s 2.2 d´av´a line´arn´ı soustavu ODR, pro kterou existuje explicitn´ı ˇreˇsen´ı vyj´adˇren´e pomoc´ı maticov´e exponenci´aly. Definice 4. Necht’ A ∈ R3×3 a s ∈ R. Pak definujeme maticovou exponenci´alu jako ∞ X sk k sA A . e = k! k=0 Pozn´amka. Tato ˇrada je absolutnˇe konvergentn´ı pro kaˇzdou matici A ∈ R3×3 . Podle teorie ODR [8] m´a line´arn´ı soustava ODR d~x(t) = A~x + ~o, dt

~x(t0 ) = ~x0

(2.3)

ˇreˇsen´ı ve tvaru ~x(t) = e(t−t0 )A~x0 + (e(t−t0 )A − Id)A−1~o.

(2.4)

Nyn´ı je ot´azkou, jak spoˇc´ıtat maticovou exponenci´alu. Jeden zp˚ usob nab´ız´ı pouˇzit´ı Cayleyovy-Hamiltonovy vˇety, ˇze matice A je koˇrenem sv´eho charakteristick´eho mnohoˇclenu. My ale zvol´ıme zp˚ usob v´ ypoˇctu pˇres Jordan˚ uv tvar matice. Pˇredpokl´adejme, ˇze vlastn´ı ˇc´ısla interpolaˇcn´ı matice A jsou navz´ajem r˚ uzn´a a nenulov´a. Rozliˇsujeme dva moˇzn´e pˇr´ıpady vlastn´ıch ˇc´ısel re´aln´e matice. Re´ aln´ a vlastn´ı ˇ c´ısla M´a-li matice A tˇri navz´ajem r˚ uzn´a re´aln´a vlastn´ı ˇc´ısla, lze ji zapsat ve tvaru A = SDS −1 , kde



 λ1 0 0 D =  0 λ2 0  0 0 λ3

a matice S obsahuje ve sloupc´ıch vlastn´ı vektory pˇr´ısluˇsn´e postupnˇe vlastn´ım ˇc´ısl˚ um λ1 , λ2 a λ3 . Tento tvar je vhodn´ y pro v´ ypoˇcet maticov´e exponenci´aly, protoˇze esA =

∞ X sk k=0

k!

(SDS −1 )k =

∞ X sk k=0

Oznaˇcme

∞ X 1 S(D)k S −1 = S( (sD)k )S −1 . k! k! k=0

 eλ1 t˜ 0 0 Λ(t) =  0 eλ2 t˜ 0  , 0 0 eλ3 t˜ 

kde t˜ je symbol pro (t − t0 ), pak maticov´a exponenci´ala m´a tvar ˜

etA = SΛ(t)S −1 .

Nyn´ı upravme druhou ˇc´ast v´ yrazu 2.4. ˜

(etA − Id)A−1 = S(Λ(t) − Id)S −1 A−1 , 12

pˇriˇcemˇz A−1 = SD−1 S −1 a inverze diagon´aln´ı matice se spoˇcte snadno, je opˇet diagon´aln´ı s pˇrevr´acen´ ymi hodnotami v˚ uˇci p˚ uvodn´ı matici. Oznaˇcme tedy   eλ1 t˜−1 0 0 λ1   eλ2 t˜−1 ∆(t) =  0 0 . λ2 0

0

eλ3 t˜−1 λ3

Explicitn´ı vyj´adˇren´ı trajektorie ˇca´stice m´a tedy tvar ~x(t) = SΛ(t)S −1~x0 + S∆(t)S −1~o.

(2.5)

Komplexn´ı p´ ar vlastn´ıch ˇ c´ısel Pokud matice A m´a jen jedno re´aln´e a p´ar komplexnˇe sdruˇzen´ ych vlastn´ıch ˇc´ısel, tak´e existuje diagon´aln´ı matice D a transformaˇcn´ı matice S, ˇze A = SDS −1 . Nicm´enˇe obˇe jsou komplexn´ı. N´asleduj´ıc´ı tvrzen´ı hovoˇr´ı o tzv. re´aln´em Jordanovˇe tvaru matice. Vˇ eta 2. Necht’ A ∈ R3×3 a λ1 ,λ2 a λ3 jsou jej´ı vlastn´ı ˇc´ısla, pˇriˇcemˇz λ1 = α + βi a λ2 = α − βi pro α,β ∈ R, β 6= 0. Necht’ ~u a ~v jsou postupnˇe re´aln´a a imagin´arn´ı ˇc´ast vlastn´ıho vektoru pˇr´ısluˇsn´emu vlastn´ımu ˇc´ıslu λ1 a w ~ je vlastn´ı vektor pˇr´ısluˇsn´y vlastn´ımu ˇc´ıslu λ3 . Pak plat´ı A = SBS −1 , kde

a



 | | | ~  S =  ~u ~v w | | |   α β 0 B =  −β α 0  . 0 0 λ3

D˚ ukaz. Protoˇze ~u + i~v je vlastn´ı vektor, plat´ı A(~u + i~v ) = (α + iβ)(~u + i~v ) a srovn´an´ım re´aln´ ych a imagin´arn´ıch ˇc´ast´ı dost´av´ame A~u = α~u − β~v

(2.6)

A~v = α~v + β~u

(2.7)

Nyn´ı uk´aˇzeme, ˇze ~u a ~v jsou line´arnˇe nez´avisl´e. Pro spor pˇredpokl´adejme c~u + d~v = 0

(2.8)

tak, ˇze bud’ c 6= 0 nebo d 6= 0. Aplikov´an´ım A na 2.8 dostaneme (cα + dβ)~u + (dα − cβ)~v = 0. Po dosazen´ı d~v = −c~u tedy dost´av´ame βd~u = βc~v . To pˇr´ımo implikuje d~u = c~v a tedy d2~u = dc~v . Z´aroveˇ n z 2.8 m´ame cd~v = −c2~u. Koneˇcnˇe tedy m´ame d2~u = −c2~u. 13

Nebot’ c a d jsou re´aln´a ˇc´ısla a alespoˇ n jedno z nich nen´ı nulov´e, nem˚ uˇze platit 2 2 d = −c , ˇcili mus´ı platit ~u = 0. Nyn´ı z 2.7 obdrˇz´ıme A~v = α~v , protoˇze α nen´ı vlastn´ı ˇc´ıslo matice A, nutnˇe ~v = 0. Coˇz je spor, ponˇevadˇz ~u + i~v = 0, ale z´aroveˇ n je to vlastn´ı vektor matice A.

Jeˇstˇe je tˇreba uk´azat, ˇze tak´e w ~ je line´arnˇe nez´avisl´ y s rovinou definovanou vektory ~u, ~v . Pˇredpokl´adejme pro spor, ˇze existuj´ı re´aln´a ˇc´ısla c a d tak, ˇze c~u + d~v = w ~ a c 6= 0 nebo d 6= 0. Nyn´ı na tuto rovnost aplikujeme oper´ator A a dost´av´ame A(c~u + d~v ) = λ3 (w) ~ = λ3 (c~u + d~v ). Ted’ podle 2.6 a 2.7 dosad´ıme za levou stranu a dostaneme (cα + dβ − cλ3 )~u + (dα − cβ − dλ3 )~v = 0. Vzhledem k tomu, ˇze jsou vektory ~u a ~v line´arnˇe nez´avisl´e, nutnˇe mus´ı platit      α − λ3 β d 0 = . −β α − λ3 c 0 Tato soustava rovnic m´a netrivi´aln´ı ˇreˇsen´ı pouze v pˇr´ıpadˇe, ˇze λ3 je vlastn´ı ˇc´ıslo matice   α − λ3 β . −β α − λ3 Jenˇze jej´ı vlastn´ı ˇc´ısla jsou α+iβ a α−iβ, coˇz znamen´a, ˇze re´aln´e ˇc´ıslo λ3 vlastn´ım ˇc´ıslem b´ yt nem˚ uˇze. Z toho d˚ uvodu nutnˇe c = 0 a d = 0, coˇz je spor.

Necht’ nyn´ı m´ame matici S, jej´ıˇz sloupce tvoˇr´ı postupnˇe vektory ~u,~v a w. ~ M´ame tedy AS = A[~u|~v |w] ~ = [A~u|A~v |Aw] ~ = [α~u − β~v |α~v + β~u|λ3 w] ~ = [~u|~v |w]B. ~ Vzhledem k tomu, ˇze jsme jiˇz ovˇeˇrili regularitu matice S, je d˚ ukaz hotov.

k

Lemma 3. Necht’ k ∈ N a A,B ∈ Rk×k . Pokud AB = BA, pak eA+B = eA eB = eB eA . ˇ D˚ ukaz. Rady eA ,eB i eA+B jsou absolutnˇe konvergentn´ı, proto je moˇzn´e napsat souˇcin dvou ˇrad jako jednu ˇradu, tzv. Cauchy˚ uv souˇcin ˇrad, v tomto tvaru: eA eB =

∞ ∞ X An X B m n! m=0 m! n=0

14

=

=

∞ X X

An B p−n n!(p − n)! p=0 0≤n≤p

∞ X 1 X p! ( An B p−n ). p! n!(p − n)! p=0 0≤n≤p

Zde pouˇzijeme AB = BA a z binomick´e vˇety dostaneme ∞ X 1 = (A + B)p . p! p=0

k Vˇ eta 4. Necht’ α , β, t˜ ∈ R a  B= Pak t˜B

e

=e

αt˜

α β −β α

 .

cos(β t˜) sin(β t˜) − sin(β t˜) cos(β t˜)



D˚ ukaz. Necht’



1 0 0 1



0 1 −1 0



I= a

 J=

 .

.

Plat´ı J 2 = −I. Pak B = αI + βJ, tedy ˜

˜

˜

etB = etαI+tβJ , Coˇz d´ıky IJ = JI a pˇredchoz´ımu lemmatu znamen´a ˜

˜

˜

etB = etαI etβJ βJ (βJ)2 + + ···] 1! 2! β2 β4 β β3 ˜ = etα [(I − I + I + · · · ) + ( J − J + · · · )] 2! 4! 1! 3! ˜ = etα [(cos(β t˜))I + (sin(β t˜))J]   cos(β t˜) sin(β t˜) αt˜ =e . − sin(β t˜) cos(β t˜) ˜

= etα I[I +

k 15

Nyn´ı jsme si pˇrichystali apar´at pro explicitn´ı vyj´adˇren´ı trajektorie ˇca´stice. Analogicky jako v re´aln´em pˇr´ıpadˇe definujeme matice Λ(t) a ∆(t), pomoc´ı kter´ ych αt˜ ˜ pˇr´ımo vyj´adˇr´ıme v´ yraz 2.4. Oznaˇcme pro zjednoduˇsen´ı Ec (t) ≡ e cos(β t), Es (t) ≡ eαt˜ sin(β t˜) a El (t) ≡ eλ3 t˜. D´ale   Ec (t) Es (t) 0 0 . Λ(t) =  −Es (t) Ec (t) 0 0 El (t) Uvaˇzujme rozklad matice A = SBS −1 na re´aln´ y Jordan˚ uv tvar jako ve vˇetˇe 2. Vˇeta 4 n´am d´av´a n´avod, jak spoˇc´ıtat maticovou exponenci´alu. Necht’ s je re´aln´e ˇc´ıslo, pak

sA

e

=

∞ X sk k=0

k!

(SBS

−1 k

) =

∞ X sk k=0

k!

k

S(B) S

−1

∞ X 1 = S( (sB)k )S −1 . k! k=0

Dosazen´ım t˜ za s dost´av´ame ˜

etA = SΛ(t)S −1 .

Inverze matice B je 

B −1

α −β 1  β α = 2 α + β2 0 0

 0 . 0 1 2 2 (α + β ) λ3

T´ımto m´ıˇr´ıme k u ´pravˇe druh´e ˇc´asti v´ yrazu 2.4. ˜

(etA − Id)A−1 = S(Λ(t) − Id)B −1 S −1 , to motivuje definici

  ∆(t) = 

(Ec (t)−1)α+Es (t)β α2 +β 2 (Ec (t)−1)β−Es (t)α α2 +β 2

−(Ec (t)−1)β+Es (t)α α2 +β 2 (Ec (t)−1)α+Es (t)β α2 +β 2

0 0

0

0

(El (t)−1) λ3

Explicitn´ı vyj´adˇren´ı trajektorie ˇca´stice m´a opˇet tvar 2.5.

16

  .

2.2.2

V´ ypoˇ cet u ´ nikov´ eho bodu z buˇ nky

V tomto okamˇziku jsme schopni urˇcit pˇresnou trajektorii ˇca´stice v r´amci jedn´e buˇ nky. Zˇrejmˇe je nutn´e zjistit, kde ˇca´stice buˇ nku opust´ı. Jin´ ymi slovy je tˇreba urˇcit, kterou stˇenu protne trajektorie ˇca´stice nejdˇr´ıv a kde. Uvaˇzujme stˇenu s s vrcholy A,B,C, a norm´alou ~n. Definujme skal´arn´ı funkci f (t) = (~x(t) − A) · ~n. Pak bod ~x(t) pro nˇejak´e t n´aleˇz´ı rovinˇe s pr´avˇe tehdy, kdyˇz f (t) = 0. Po dosazen´ı z rovnice 2.5 dost´av´ame f (t) = (Λ(t)S −1 x~0 ) · (S T ~n) + (∆(t)S −1~o) · (S T ~n) − A · ~n. Po vyn´asoben´ı pˇr´ısluˇsn´ ych matic a porovn´an´ım produkt˚ u skal´arn´ıch souˇcin˚ u doˇ jdeme k neline´arn´ı funkci obsahuj´ıc´ı exponenci´aly. Reˇsen´ı naˇs´ı u ´lohy je ekvivalentn´ı nalezen´ı nejmenˇs´ıho nez´aporn´eho3 koˇrenu t´eto funkce. Neexistuje analytick´e ˇreˇsen´ı obecn´e rovnice obsahuj´ıc´ı souˇcet nˇekolika exponenci´al, byt’ jen s line´arn´ım exponentem. Proto je v ˇcl´anku [7] doporuˇceno pouˇz´ıt Newtonovu iteraˇcn´ı metodu. Probl´ emy klasick´ e Newtonovy metody. Uved’me klasick´ y algoritmus pro hled´an´ı koˇrenu funkce f (t). Necht’ t0 je zvolen´a poˇca´teˇcn´ı hodnota. N´asleduj´ıc´ı body aproximace koˇrenu jsou d´any rekurentn´ım vztahem tn+1 = tn + kn , kde kn = − ff0(t(tnn)) je korekce aproximace v kroku n a symbol f 0 (t) znaˇc´ı derivaci funkce f v bodˇe t. Newtonovu metodu nelze pouˇz´ıt v klasick´em smyslu uv´adˇen´em v ˇcl´anku [7]. Funkce f (t) totiˇz m˚ uˇze nab´ yvat lok´aln´ıch extr´em˚ u, kde metoda teˇcen selˇze, viz obr´azky 2.2 a 2.3 zachycuj´ıc´ı skuteˇcn´e funkce f (t) z´ıskan´e v pr˚ ubˇehu simulace. Dalˇs´ı velkou komplikac´ı je hled´an´ı pouze nez´aporn´ ych koˇren˚ u. Klasick´a metoda nerozliˇsuje kladn´e a z´aporn´e koˇreny a m˚ uˇze se tedy st´at, ˇze ze sv´e podstaty spr´avnˇe nalezne z´aporn´ y koˇren, aˇc existuje kladn´ y koˇren odpov´ıdaj´ıc´ı skuteˇcn´emu bodu u ´niku z buˇ nky. Probl´emem je tak´e poˇzadavek, aby nalezen´ y koˇren byl co nejmenˇs´ı. D´ıky v´ yskyt˚ um stacion´arn´ıch bod˚ u toto v obecn´em pˇr´ıpadˇe tak´e klasick´a metoda nezaruˇcuje, viz ob´azek 2.3. Modifikovan´ a Newtonova metoda Prezentujeme nˇekolik naˇsich modifikac´ı klasick´e metody za u ´ˇcelem vyhovˇet v´ yˇse zm´ınˇen´ ym poˇzadavk˚ um a vyhnout se zm´ınˇen´ ym probl´em˚ um. • Neuvaˇzov´an´ı z´aporn´ych koˇren˚ u. Jako poˇc´ateˇcn´ı bod zvol´ıme hodnotu t0 = 0 a pot´e neuvaˇzujeme klasickou korekci k0 , n´ ybrˇz jej´ı absolutn´ı hodnotu. Stejnˇe tak budeme br´at absolutn´ı hodnoty korekc´ı dokud funkce f (ts ) nezmˇen´ı znam´enko pro nˇejak´e s ∈ N. Pak opˇet budeme br´at hodnoty (bez absolutn´ı hodnoty) korekce kn pro vˇsechna n ≥ s. Jedn´a se o heuristiku, kter´a umoˇzn´ı metodˇe teˇcen pˇrekroˇcit lok´aln´ı extr´em v kladn´em smˇeru“. Pokud ” hodnota funkce f (ts ) zmˇen´ı znam´enko oproti hodnotˇe f (ts−1 ), ze spojitosti funkce f (t) vpl´ yv´a, ˇze v intervalu [f (ts−1 ),f (ts )] leˇz´ı hledan´ y koˇren. 3

Z´ aporn´ y koˇren by tak´e urˇcoval pr˚ unik trajektorie se stˇenou, ale v z´aporn´em ˇcase.

17

Obr´azek 2.2: Uk´azka pˇr´ıpadu, kdy Newtonova metoda teˇcen nekonverguje.

Obr´azek 2.3: Uk´azka pˇr´ıpadu, kdy Newtonova metoda teˇcen konverguje ke ˇspatn´emu koˇrenu vlivem pˇr´ıtomnosti stacion´arn´ıho bodu.

18

• Rozumn´a d´elka korekce. Jak je zn´azornˇeno na obr´azku 2.3, hodnota aproximace ts bl´ızk´a stacion´arn´ımu bodu m˚ uˇze zp˚ usobit nesmyslnˇe velkou korekci 0 ks vlivem velmi mal´e hodnoty derivace f (ts ). T´ım by byl proces hled´an´ı nejmenˇs´ıho kladn´eho koˇrene ne´ uspˇeˇsn´ y, i kdyby pot´e metoda k jist´emu koˇrenu konvergovala. Proto navrhujeme zvolit hodnotu kmax > 0 a uvaˇzovat umˇelou korekci τ , pokud algoritmus doporuˇc´ı korekci v absolutn´ı hodnotˇe vˇetˇs´ı neˇz je tato mez. V dalˇs´ım kroku uˇz m˚ uˇze algoritmus navrhnout smysluplnou korekci a pokraˇcovat tak d´ale v hled´an´ı koˇrenu. • Povolen´ı oscilace. Obr´azek 2.2 ukazuje, ˇze m˚ uˇze nastat nekoneˇcn´ y cyklus v klasick´em algoritmu. V tomto pˇr´ıpadˇe pˇr´ıpadˇe pˇredpokl´ad´ame, ˇze pokud by existoval koˇren v bodˇe ts , bude se jednat o pˇr´ıliˇs velk´ y ˇcas ts , kter´ y jiˇz z hlediska simulace nen´ı zaj´ımav´ y. Nakonec stejnˇe hled´ame minim´aln´ı tmin mezi vˇsemi stˇenami buˇ nky, tedy s velkou pravdˇepodobnost´ı bude platit tmin < ts . Navrhovan´ ym ˇreˇsen´ım je omezit poˇcet iterac´ı jistou konstantou imax . Pokud je tato mez pˇrekroˇcena, staˇc´ı pˇr´ıpad prohl´asit za nekonver” gentn´ı“ a standardn´ım postupem pokraˇcovat. • V´ıce povolen´ych iterac´ı. Klasick´a metoda zpravidla vyˇzaduje 5 aˇz 6 iterac´ı pro dosaˇzen´ı k´ yˇzen´e pˇresnosti. Nicm´enˇe v naˇsem pˇr´ıpadˇe je tˇreba br´at ohled na moˇznou velkou vzd´alenost od koˇrene, dokonce br´at v u ´vahu stacion´arn´ı bod mezi v´ ychoz´ı pozic´ı a k´ yˇzen´ ym koˇrenem. Proto je tˇreba nastavit maxim´aln´ı poˇcet iterac´ı imax na vyˇsˇs´ı hodnotu. Experiment´alnˇe zjiˇstˇen´a vhodn´a hodnota imax je okolo 15-ti, ovˇsem je tˇreba imax ch´apat jako parametr metody. • Konvergence zprava. Je tˇreba si uvˇedomit, ˇze pro simulaci je d˚ uleˇzit´e, aby se ˇca´stice v dalˇs´ım kroku nach´azela uvnitˇr dan´e buˇ nky. Pˇredpokl´adejme, ˇze ˇca´stice se nach´az´ı v buˇ nce 1 a metoda vr´at´ı jako koˇren ˇcas t − , kde  > 0 s t´ım, ˇze ˇc´astice proˇsla stˇenou s a dostala se tak do buˇ nky 2. V dalˇs´ım kroku metody opˇet pro stˇenu s obdrˇz´ıme nez´aporn´ y koˇren epsilon a ˇca´stice t´emˇeˇr nezmˇen´ı svou polohu. Stane se pouze, ˇze se zmˇen´ı okupace buˇ nky z 2 na 1, nebot’ dle algoritmu ˇca´stice proˇsla stˇenou s. Simulace se t´ımto dostane do nekoneˇcn´e smyˇcky nebo bude dlouhou dobu setrv´avat na jednou m´ıstˇe. Proto je nutn´e, aby Newtonova metoda vr´atila vˇzdy t + . ˇ astice mˇela Uk´azka selh´an´ı klasick´e Newtonovy metody je na obr´azku 2.4. C´ uniknout z buˇ nky v ˇcase odpov´ıdaj´ıc´ımu nejmenˇs´ımu kladn´emu koˇrenu funkce f (t) z obr´azku 2.6. Nicm´enˇe klasick´a metoda vr´atila z´aporn´ y koˇren a t´ım p´adem byla ztracena informace o spr´avn´em“ koˇrenu v okol´ı bodu t = 0,004. Pro jinou ” stˇenu pak vyˇsel koˇren pˇr´ısluˇsn´e funkce t = 0,06, mezi vˇsemi stˇenami byl nejmenˇs´ı, a simulace jej vybrala jako ˇcas, kdy ˇca´stice unikla z buˇ nky. Tento omyl byl napraven heuristikou o neuvaˇzov´an´ı z´aporn´ ych koˇren˚ u a v´ ysledkem je spr´avn´ y obr´azek 2.5.

Uveden´e modifikace jsou pouze heuristiky, jak zv´ yˇsit pravdˇepodobnost u ´spˇechu klasick´e metody. Ale je tˇreba m´ıt kontrolu, zda-li nalezen´ y koˇren je skuteˇcnˇe nejmenˇs´ı, ˇcili skuteˇcnˇe odpov´ıd´a spr´avn´emu bodu u ´niku z buˇ nky. Funkce f (t) odpov´ıd´a vzd´alenosti bodu trajektorie ˇca´stice od roviny pˇr´ısluˇsn´e stˇeny. Proto i 19

Obr´azek 2.4: V´ ypoˇcet ˇcasu u ´niku klasickou Newtonovou metodou.

Obr´azek 2.5: V´ ypoˇcet ˇcasu u ´niku modifikovanou Newtonovou metodou.

20

Obr´azek 2.6: Funkce f (t), kde selhala klasick´a Newtonova metoda. pro f (t0 ) = 0 je tˇreba otestovat, zda ˇca´stice pˇr´ısluˇsn´e stˇenˇe buˇ nky 4 . Dokonce ani v tomto pˇr´ıpadˇe si nem˚ uˇzeme b´ yt jist´ı, ˇze jsme naˇsli spr´avn´ y bod u ´niku ˇca´stice 5 z buˇ nky , ale z˚ ustane zachov´ana kompatibilita pozice ˇca´stice a informace o okupovan´e buˇ nce. Dalˇs´ı pr˚ ubˇeh simulace probˇehne spr´avnˇe a jedin´ y d˚ usledek bude mal´e poskoˇcen´ı“ ˇc´astice. Ovˇsem tento pˇr´ıpad nastane velmi zˇr´ıdka. ” M˚ uˇze se st´at, ˇze nen´ı nalezen ˇz´adn´ y vhodn´ y ˇcas u ´niku nebo nalezen´ y bod neleˇz´ı ve stˇenˇe viz pˇredchoz´ı odstavec. Algoritmus oznaˇc´ı buˇ nku za zvl´aˇstn´ı“ a ” pˇrepne na klasickou integraˇcn´ı metodu (zde na modifikovanou Eulerovu metodu), dokud se ˇca´stice neobjev´ı v nov´e buˇ nce.

2.2.3

Moˇ zn´ e typy bunˇ ek

Popsan´ y postup jsme uvaˇzovali v pˇr´ıpadˇe, ˇze interpolaˇcn´ı matice A je regul´arn´ı a m´a navz´ajem r˚ uzn´a vlastn´ı ˇc´ısla. Zvaˇzme nyn´ı pˇr´ıpad, ˇze vlastn´ı ˇc´ıslo je nulov´e. Bez u ´jmy na obecnosti uvaˇzujme λ3 = 0. Pak tˇret´ı diferenci´aln´ı rovnice 2.2 m´a tvar d(~x(t))3 = (S~oS −1 )3 , (~x(t0 ))3 = (x~0 )3 , dt kde (w(t)) ~ c´ı tˇret´ı souˇradnici vektoru w. ~ Odtud plyne 3 znaˇ (~x(t))3 = (x~0 )3 + (t − t0 )(S~oS −1 )3 . λ t˜

Tedy form´alnˇe staˇc´ı pouze napsat do matice ∆(t) m´ısto v´ yrazu e 3λ3−1 v´ yraz t˜, do ˜ matice Λ(t) napsat 1 nam´ısto eλ3 t a algoritmus bude fungovat spr´avnˇe. Posledn´ım speci´aln´ım pˇr´ıpadem jsou buˇ nky, jej´ıˇz vektory rychlost´ı ve vrcholech jsou rovnobˇeˇzn´e. Pak je jist´e, ˇze ˇc´astice se bude pohybovat v tomto smˇeru po pˇr´ımce, ˇcili u ´lohou z˚ ust´av´a pouze naj´ıt pr˚ useˇc´ık pˇr´ımky a roviny. 4

Zde je tˇreba poˇc´ıtat se zaokrouhlovac´ımi chybami, tedy nestaˇc´ı br´at v u ´vahu konvexn´ı obal vrchol˚ u stˇeny, je nutn´e pˇripustit jistou v´ ychylku ve smˇeru norm´aly stˇeny 5 M˚ uˇze nastat pˇr´ıpad, ˇze trajektorie tuto stˇenu protne dvakr´at v ˇcasech t0 a t1 , pˇriˇcemˇz t0 < t1 . T´ım p´ adem zm´ınˇen´ ym testem projde i ˇcas t1 , aˇc se nejedn´a o spr´avn´ y ˇcas u ´niku ˇc´astice z buˇ nky.

21

Uved’me moˇzn´e typy bunˇek. Typ Krit´ erium Re´aln´a Vlastn´ı ˇc´ısla jsou re´aln´a a navz´ajem r˚ uzn´a Komplexn´ı Vlastn´ımi ˇc´ısly jsou komplexnˇe sdruˇzen´a ˇc´ısla Paraleln´ı Vektory rychlost´ı jsou rovnobˇeˇzn´e Zvl´aˇstn´ı Vˇsechny ostatn´ı buˇ nky Pozn´amka. V ˇcl´anku [7] je uvedeno, ˇze za zvl´aˇstn´ı“ by se mˇely oznaˇcit i buˇ nky ” obsahuj´ıc´ı kritick´ y bod. Jedn´a se o bod, kde je vektor rychlosti nulov´ y, tedy pˇredurˇcuje zvl´aˇstn´ı chov´an´ı vektorov´eho pole uvnitˇr buˇ nky. Autoˇri ˇcl´anku se zˇrejmˇe domn´ıvaj´ı, ˇze pouze tyto buˇ nky mohou zp˚ usobit probl´emy s konvergenc´ı Newtonovy iteraˇcn´ı metody. Ovˇsem experiment´alnˇe ovˇeˇrenou pravdou je, ˇze i po jejich odstranˇen´ı se st´ale vyskytuj´ı buˇ nky, kde klasick´a Newtonova metoda selh´av´a. Proto jsme se rozhodli ignorovat polohu kritick´eho bodu pˇri rozdˇelov´an´ı bunˇek do jednotliv´ ych kategori´ı.

2.3

Srovn´ an´ı metod

Porovn´avat metody m˚ uˇzeme kvantitativnˇe v ˇcasov´e ˇci pamˇet’ov´e sloˇzitosti. Dalˇs´ım krit´eriem m˚ uˇze b´ yt implementaˇcn´ı n´aroˇcnost, stejnˇe jako v´ yˇse kvalifikace poˇzadovan´a po pˇr´ıpadn´em uˇzivateli. Nicm´enˇe hlavn´ı a obt´ıˇznˇe porovnateln´ y znak, ve kter´em by se mˇely naˇse dvˇe metody liˇsit, je pˇresnost.

2.3.1

Pˇ resnost metody

M˚ uˇzeme se na probl´em pod´ıvat z analytick´eho hlediska a zjist´ıme, ˇze modifikovan´a Eulerova metoda (d´ale jen MEM“) m´a diskretizaˇcn´ı chybu O(h2 ). Lok´alnˇe ” exaktn´ı metoda (d´ale jen LEM“) m´a diskretizaˇcn´ı chybu prakticky nulovou. ” V tom spoˇc´ıv´a jej´ı hlavn´ı pˇr´ınos do oboru numerick´e matematiky, protoˇze t´eto pˇresnosti MEM nedos´ahne nikdy. Jedin´ y zp˚ usob jak zpˇresnit MEM je zmenˇsov´an´ı ˇcasov´eho kroku, t´ım ovˇsem zase nar˚ ust´a vliv zaokrouhlovac´ıch chyb. Naproti tomu parametry LEM jsou maxim´aln´ı poˇcet iterac´ı Newtonovy metody imax , maxim´aln´ı povolen´a korekce kmax a umˇel´a korekce τ . Zmˇenou hodnot parametr˚ u (vyˇsˇs´ı i( max), niˇzˇs´ı τ a kmax ) dost´av´ame pomalejˇs´ı metodu, ale velkou v´ yhodou je, ˇze se nezvyˇsuje vliv zaokrouhlovac´ıch chyb. Zmˇena hodnot parametr˚ u totiˇz pouze zv´ yˇs´ı u ´spˇeˇsnost modifikovan´e Newtonovy metody v hled´an´ı koˇren˚ u a t´ım p´adem sn´ıˇz´ı poˇcet pˇr´ıpad˚ u, kdy je nutn´e pˇrepnout z lok´alnˇe exaktn´ı integrace na klasickou metodu s ˇcasov´ ym krokem. Dalˇs´ı moˇznost´ı zkoum´an´ı pˇresnosti metody je konstrukce dom´eny, v n´ıˇz zn´ame spr´avn´e trajektorie ˇca´stic, a zkoumat, jak moc se od nich liˇs´ı vypoˇcten´a data. Typick´ ym pˇr´ıkladem je v´alec s vektorov´ ym polem takov´ ym, ˇze by se ˇca´stice mˇely pohybovat po kruˇznic´ıch. Na obr´azku 2.7 je vlevo zn´azornˇena simulace pohybu ˇca´stice MEM s ˇcasov´ ym krokem ∆t = 0.001. Simulovan´a ˇc´astice se postupnˇe oddaluje od osy v´alce. Vpravo je zn´azornˇena stejn´a situace se shodnou poˇca´teˇcn´ı polohou ˇca´stice vypoˇc´ıtan´a pomoc´ı LEM, kter´a ani v jednom pˇr´ıpadˇe nepˇrepnula na klasickou integraci s ˇcasov´ ym krokem. Tato ryz´ı exaktn´ı integrace tedy siˇ simulace MEM ˇcinil mulovala pˇresnou trajektorii ˇca´stice po uzavˇren´e kˇrivce. Cas 22

9 sekund, kdeˇzto LEM skonˇcila po 12-ti sekund´ach. Jeˇstˇe dopln´ıme, ˇze pro z´ısk´an´ı podobn´e pˇresnosti MEM jsme potˇrebovali ˇcasov´ y krok ∆t = 0.00002 a v´ ypoˇcetn´ı ˇcas 93 sekund, viz obr´azek 2.8.

Obr´azek 2.7: Srovn´an´ı MEM ∆t = 0,001 (vlevo) a ryz´ı LEM bez pˇrepnut´ı na MEM (vpravo).

Obr´azek 2.8: MEM ∆t = 0,00002.

2.3.2

Ostatn´ı krit´ eria

Dalˇs´ım d˚ uleˇzit´ ym srovn´avac´ım krit´eriem je jiˇz zm´ınˇen´a ˇcasov´a a prostorov´a n´aroˇcnost. MEM vyˇzaduje pˇredv´ ypoˇcet norm´al vˇsech stˇen a stˇredy bunˇek s´ıtˇe. Naproti tomu LEM je nesrovnatelnˇe n´aroˇcnˇejˇs´ı na pamˇet’ i ˇcas pˇredv´ ypoˇctu. Potˇrebuje rozdˇelit buˇ nky na jednotliv´e typy, u paraleln´ıch bunˇek uloˇzit smˇer vektorov´eho pole a u komplexn´ıch a re´aln´ıch bunˇek se jedn´a o v´ ypoˇcet vlastn´ıch ˇc´ısel, 23

transformaˇcn´ı matice sest´avaj´ıc´ı z vlastn´ıch vektor˚ u a jej´ı inverze. Pˇredv´ ypoˇcet na s´ıti ˇc´ıtaj´ıc´ı okolo 20000 bunˇek trval 4 sekundy MEM a 37 sekund LEM. Pro porovn´an´ı ˇcasov´e sloˇzitosti samotn´e simulace pouˇzijeme graf na obr´azku 2.9 popisuj´ıc´ı ˇcas simulace pro 100 ˇca´stic. R˚ ust ˇcasov´e n´aroˇcnosti LEM je mnohem pozvolnˇejˇs´ı, protoˇze pouze mal´a ˇca´st bunˇek je oznaˇcena za zvl´aˇstn´ı“. Typicky se ” jedn´a o m´enˇe neˇz 0,1% vˇsech bunˇek.

Obr´azek 2.9: Srovn´an´ı ˇcasov´e sloˇzitosti simulace pro 100 ˇc´astic. Ve prospˇech LEM hovoˇr´ı poˇcet napoˇcten´ ych bod˚ u trajektorie. Pokud uvaˇzujeme 1000 ˇca´stic, krok metody 0.001 a c´ılov´ y ˇcas T = 1, v´ ysledkem MEM bude 106 bod˚ u. V praxi se m˚ uˇze prodlouˇzit c´ılov´ y ˇcas, zkr´atit d´elka kroku a velikost dat se stane z´avaˇzn´ ym probl´emem. Je moˇzn´e napˇr´ıklad po vypoˇcten´ı deseti nov´ ych bod˚ u trajektorie devˇet zapomenout a uchovat si pouze polohu posledn´ıho. Ot´azkou z˚ ust´av´a, jestli nevynech´ame ty d˚ uleˇzit´e a naopak v oblasti pˇr´ım´e trajektorie nebude zbyteˇcnˇe mnoho bod˚ u. LEM si pamatuje pouze body pr˚ uniku se stˇenami bunˇek. To znamen´a, ˇze v oblasti z´ajmu bude bod˚ u trajektorie v´ıc a naopak. LEM je tedy pˇrirozenˇe adaptivn´ı metodou.

24

Obr´azek 2.10: Nahoˇre MEM ∆t = 0.001, uprostˇred MEM ∆t = 0.0001, dole LEM bez zvl´aˇstn´ıch bunˇek. Stoj´ı za povˇsimnut´ı, ˇze na horn´ım obr´azku zcela chyb´ı nejvyˇsˇs´ı smyˇcky“ trajektorie. ”

25

Nakonec nab´ıdneme srovn´an´ı tvar˚ u trajektori´ı ˇc´astic v re´aln´e dom´enˇe, konkr´etnˇe se jedn´a o c´evu v mozku s patologickou vydut´ı, viz obr´azek 2.10. Na skuteˇcn´ ych datech se jen velmi tˇeˇzko pozn´av´a, zda-li je vypoˇcten´a trajektorie spr´avn´a ˇci nikoliv, nicm´enˇe pˇresvˇedˇciv´ ym argumentem je fakt, ˇze se zmenˇsov´an´ım ˇcasov´eho kroku MEM, se napoˇcten´e trajektorie bl´ıˇz´ı v´ ysledku LEM.

2.4

Program

Na z´avˇer kapitoly uvedeme z´akladn´ı strukturu implementovan´eho simulaˇcn´ıho programu. Podrobnosti je moˇzn´e nal´ezt v apendixu. Projekt je napsan´ y v jazyce C++ a sest´av´a se ze tˇr´ı soubor˚ u: • uTetVol.h je rozhran´ı pro pr´aci s tˇr´ıdou TTetBasicVol reprezentuj´ıc´ı dom´enu rozdˇelenou na ˇctyˇrstˇeny. Tato tˇr´ıda disponuje metodami pro naˇcten´ı dat o dom´enˇe z .vtk souboru6 , pro pˇredv´ ypoˇcet dat typu norm´aly stˇen ” ˇctyˇrstˇen˚ u“ ˇci sousednosti bunˇek7“. D´ale zahrnuje pomocn´e metody pro ” stopovac´ı algoritmy pro vˇetˇs´ı pˇrehlednost hlavn´ıho souboru. • uTetVol.cpp je tradiˇcnˇe pouze implementace metod deklarovan´ ych ve stejnojmenn´em hlaviˇckov´em souboru. • MainFile.cpp obsahuje dvˇe funkce TracingEuler a TracingLocalExact, kter´e parametrem pˇrij´ımaj´ı tˇr´ıdu TTetBasicVol, v´ ychoz´ı pozici ˇca´stice, celkov´ y ˇcas simulace a parametry jednotliv´ ych sledovac´ıch metod. Kaˇzd´a simulaˇcn´ı funkce vrac´ı vektor trajektorie datov´eho typu double se souˇradnicemi jednotliv´ ych bod˚ u. Simulaˇcn´ı j´adro je navrˇzen´e tak, ˇze vygeneruje poˇc´ateˇcn´ı polohy ˇc´astic a postupnˇe je pos´ıl´a do jedno ze dvou stopovac´ıch algoritm˚ u, aby nazpˇet dostalo vektor trajektorie. Takto vznikne vektor dataStore vˇsech trajektori´ı, jenˇz je pomoc´ı funkce TiskCastice konvertov´an ve spr´avn´em form´atu do .vtk souboru. Vygenerovan´ y soubor lze zobrazit v libovoln´em programu podporuj´ıc´ı form´at .vtk. My jej zobrazujeme pomoc´ı opensource softwaru ParaView [6].

6 7

Form´ at souboru .vtk je rozˇs´ıˇren´ ym mediem pro sd´ılen´ı dat mezi rozliˇcn´ ymi programy. Pojmem sousednost bunˇek m´ın´ıme informaci, kter´e buˇ nky soused´ı se kter´ ymi.

26

Z´ avˇ er Navrhli a implementovali jsme program pro sledov´an´ı ˇc´astic zaloˇzen´ y pˇredevˇs´ım na lok´alnˇe exaktn´ı metodˇe se z´aloˇzn´ı variantou pro zvl´aˇstn´ı buˇ nky ve formˇe modifikovan´e Eulerovy metody. Jedn´a se o algoritmus vykresluj´ıc´ı trajektorii ˇc´astice pomaleji neˇz klasick´e krokov´e algoritmy, zato s nesrovnatelnˇe vyˇsˇs´ı pˇresnost´ı. Tyto poˇzadavky pr´avˇe na pˇresnost kladou napˇr´ıklad nˇekter´e simulace proudˇen´ı krve, coˇz m˚ uˇze pˇredstavovat prostor pro vyuˇzit´ı lok´alnˇe exaktn´ı metody. Z´aklad metody vych´az´ı z ˇcl´anku [7], kde vˇsak nejsou ˇreˇseny vˇsechny aspekty navrˇzen´eho algoritmu. Naˇs´ım hlavn´ım pˇr´ınosem pro metodu je n´avrh modifikace Newtonovy iteraˇcn´ı metody, aby konvergovala ke koˇrenu funkce vyhovuj´ıc´ımu dan´ ym poˇzadavk˚ um, a to i za pˇredpokladu pˇr´ıtomnosti stacion´arn´ıho bodu funkce mezi poˇca´teˇcn´ı hodnotou a k´ yˇzen´ ym koˇrenem funkce. Jedn´a se o heuristiky, nezaruˇcuj´ı tedy vˇzdy spr´avn´ y v´ ysledek, nicm´enˇe na re´aln´ ych datech se osvˇedˇcily ´ echem je kupˇr´ıkladu stopov´an´ı ˇc´astic v s´ıti c´evy, viz obr´azek velmi dobˇre. Uspˇ 2.10, kter´e probˇehlo v ryze exaktn´ım m´odu bez pouˇzit´ı numerick´e integrace, ˇcili s nulovou chybou pod´el cel´e trajektorie. Moˇzn´ ym pokraˇcov´an´ım pr´ace je n´avrh sofistikovanˇejˇs´ıch heuristik pro hled´an´ı k´ yˇzen´eho koˇrene funkce v pˇrijateln´em ˇcase. Dalˇs´ım rozˇs´ıˇren´ım pr´ace je prozkoum´an´ı metod typu Runge-Kutta nam´ısto Eulerovy metody, pˇr´ıpadnˇe v´ıcekrokov´ ych metod pro numerickou integraci.

27

Apendix K pr´aci dod´av´ame podrobnou dokumentaci k programu s implementac´ı lok´alnˇe exaktn´ı metody a modifikovan´e Eulerovy metody.

2.5

Struktura programu a pouˇ zit´ e algoritmy

V u ´vodu pop´ıˇseme tˇr´ıdu TTetBasicVol reprezentuj´ıc´ı s´ıt’ ze ˇctyˇrstˇen˚ u, pot´e pˇrejdeme k popisu funkc´ı nach´azej´ıc´ıch se v souboru MainFile.cpp. Nebudeme se zab´ yvat konkr´etn´ımi n´azvy promˇenn´ ych ˇci pˇresn´ ym v´ yˇctem parametr˚ u kaˇzd´e funkce. Naˇs´ım c´ılem je srozumitelnˇe sdˇelit myˇslenku jednotliv´ ych komponent.

2.5.1

TTetBasicVol

Ve tˇr´ıdˇe TTetBasicVol jsou uloˇzena veˇsker´a data o s´ıti v jednorozmˇern´ ych ’ pol´ıch datov´eho typu double, at uˇz se jedn´a o skal´arn´ı veliˇciny, vektory ˇci matice8 . Je tak uˇcinˇeno za u ´ˇcelem sn´ıˇzit ˇcasov´e n´aroky pˇri pr´aci s daty. Tˇr´ıda obsahuje nˇekolik funkc´ı pro pˇr´ıstup k dat˚ um, kaˇzdou z nich pˇredch´az´ı kl´ıˇcov´e slovo inline urˇcuj´ıc´ı, ˇze pˇred kompilac´ı se tˇelo funkce m˚ uˇze pˇrepsat do m´ısta, kde je funkce vol´ana (viz [11]). T´ımto je podpoˇrena pˇrehlednost programu bez zv´ yˇsen´ı ˇcasov´e sloˇzitosti. N´asleduje v´ yˇcet hlavn´ıch metod tˇr´ıdy. Pro jednoduchost oznaˇcme S1 ,S2 ,S3 vrcholy stˇeny a V1 ,V2 ,V3 ,V4 vrcholy ˇctyˇrstˇenu. LoadFromVTKFile je metoda pro naˇcten´ı dat z .vtk souboru v ASCII form´atu. Jedn´a se o souˇradnice vrchol˚ u, seznam ˇctyˇrstˇen˚ u, vektory rychlost´ı ve vrcholech a pˇr´ıpadnˇe dalˇs´ı data jako tlak nebo teplota. BuildData z naˇcten´ ych dat pro kaˇzd´ y ˇctyˇrstˇen vypoˇc´ıt´a n´asleduj´ıc´ı. • Stˇ red buˇ nky odpov´ıd´a line´arn´ı kombinaci vrchol˚ u s koeficienty 41 . • Norm´ aly stˇ en poˇc´ıt´ame pomoc´ı normovan´eho vektorov´emu souˇcinu ~n =

(S2 − S1 ) × (S3 − S1 ) . ||(S2 − S1 ) × (S3 − S1 )||

Nebot’ vyˇzadujeme vnˇejˇs´ı norm´alu, provedeme jeˇstˇe skal´arn´ı souˇcin se ˇctvrt´ ym vrcholem ˇctyˇrstˇenu a pokud je kladn´ y, uloˇz´ıme si vektor −~n. 8

Matice jsou ukl´ ad´ any po ˇr´ adc´ıch kv˚ uli optimalizaci ˇcasov´e sloˇzitosti n´asoben´ı matice a vektoru

28

• Vrcholov´ a pˇ r´ısluˇ snost pro kaˇzd´ y vrchol urˇcuje, do kter´e buˇ nky patˇr´ı. V prvn´ı f´azi je pro kaˇzd´ y vrchol zjiˇstˇen pouze poˇcet takov´ ych bunˇek. D´ıky tomu m˚ uˇze b´ yt alokov´ano pole, do kter´eho jsou v druh´e f´azi vzestupnˇe zaps´any identifikaˇcn´ı ˇc´ısla pˇr´ısluˇsn´ ych bunˇek. Vrcholov´a pˇr´ısluˇsnost slouˇz´ı jen pro v´ ypoˇcet sousednosti, neukl´ad´a se. • Sousednost je informace o okol´ı buˇ nky. Pro kaˇzdou stˇenu je tˇreba urˇcit, do kter´e buˇ nky bychom se dostali jej´ım pˇrekroˇcen´ım. Provedeme tedy pr˚ unik vrcholov´e pˇr´ısluˇsnosti vˇsech tˇr´ı vrchol˚ u. D´ıky uspoˇra´d´an´ı prvk˚ u v seznamech projdeme kaˇzd´ y seznam maxim´alnˇe jednou. • Interpolaˇ cn´ı matice a translaˇ cn´ı vektor je takov´a matice A a vektor ~o, ˇze pro vˇsechna i ∈ 1,2,3,4 plat´ı ~v (Vi ) = AVi + ~o. Matice A a vektor ~o jsou urˇceny jednoznaˇcnˇe tˇemito ˇctyˇrmi podm´ınkami. • Vlastn´ı ˇ c´ısla a vlastn´ı vektory interpolaˇcn´ı matice A jsou poˇc´ıt´any za pouˇzit´ı knihovny Eigen [5]. • Typ buˇ nky je urˇcen podle vlastn´ıch ˇc´ısel matice A. Pokud jsou alespoˇ n dva koˇreny shodn´e, buˇ nka je oznaˇcna za zvl´aˇstn´ı“. Je-li mezi vlastn´ımi ” ˇc´ısly komplexnˇe sdruˇzen´ y p´ar, je buˇ nka komplexn´ı“. Pro pˇr´ıpad navz´ajem ” r˚ uzn´ ych re´aln´ ych ˇc´ısel je buˇ nka oznaˇcena re´aln´a“. Posledn´ım pˇr´ıpadem je ” buˇ nka paraleln´ı“, pro kterou si uchov´ame smˇer spoleˇcn´ y vˇsem vektor˚ um ” rychlosti ve vrcholech. JeTu critical je metoda, kter´a pˇrijme parametrem bod X a identifikaˇcn´ı ˇc´ıslo buˇ nky. Vr´at´ı true, pokud bod leˇz´ı v buˇ nce a false v opaˇcn´em pˇr´ıpadˇe. Bod n´aleˇz´ı buˇ nce pr´avˇe tehdy, kdyˇz n´aleˇz´ı konvexn´ımu obalu vrchol˚ u. Hled´ame tedy koeficienty line´arn´ı kombinace vektor˚ u a(V2 − V1 ) + b(V3 − V1 ) + c(V4 − V1 ) = X, coˇz je soustava line´arn´ıch rovnic. Bod X n´aleˇz´ı dan´e buˇ nce, pr´avˇe tehdy kdyˇz plat´ı a + b + c ≤ 1 a z´aroveˇ n a ≥ 0, b ≥ 0, c ≥ 0. Kv˚ uli zaokrouhlovac´ım chyb´am uvaˇzujeme ostr´e nerovnosti nam´ısto neostr´ ych. SpravnaStena slouˇz´ı k ovˇeˇren´ı, zda bod X n´aleˇz´ı dan´e stˇenˇe. Ve svˇetˇe s koneˇcnou pˇresnost´ı to ale neznamen´a ˇreˇsit ot´azku, zda-li bod n´aleˇz´ı do konvexn´ıho obalu vektor˚ u (S2 − S1 ) a (S3 − S1 ). Je nutn´e pˇridat jeˇstˇe norm´alu ~n stˇeny pro pokryt´ı pˇr´ıpadu mal´e odchylky ˇc´astice od stˇeny. My tuto funkci budeme pouˇz´ıvat za pˇredpokladu, ˇze ˇca´stice je dostateˇcnˇe bl´ızko roviny urˇcen´e body S1 , S2 a S3 , z ˇcehoˇz vypl´ yv´a, ˇze koeficient line´arn´ı kombinace u vektoru ~n vyjde zanedbatelnˇe mal´ y. Interpoluj pro danou buˇ nku a pozici X vr´at´ı pomoc´ı interpolaˇcn´ı matice A a translaˇcn´ıho vektoru ~o interpolovanou rychlost ~v (X) = AX + ~o. 29

NewtonIter a NewtonIter komplex jsou metody pro numerick´e nalezen´ı nejmenˇs´ıho kladn´eho koˇrenu funkce f (t) vypl´ yvaj´ıc´ı z lok´alnˇe exaktn´ı metody pˇri hled´an´ı bodu u ´niku ˇca´stice z buˇ nky. Jedn´a se o modifikovanou Newtonovu metodu navrˇzenou v kapitole 2. Pokud Newtonova metoda nekonverguje, vr´at´ı hodnotu −50. Pokud konverguje k z´aporn´emu koˇrenu, vr´at´ı hodnotu −10. Pˇri konvergenci ke spr´avn´emu koˇrenu vr´at´ı hodnotu o m´alo vˇetˇs´ı, aby ˇc´astice skuteˇcnˇe byla v sousedn´ı buˇ nce (pro pˇr´ıpad konvergence zleva).

2.5.2

TracingEuler

Funkce TracingEuler simulaˇcn´ı j´adro Eulerovy metody stopov´an´ı ˇc´astic. Funkce pˇrijme parametrem instanci tˇr´ıdy TTetBasicVol, poˇca´teˇcn´ı pozici ˇc´astice, celkov´ y ˇcas simulace a ˇcasov´ y krok. Posledn´ım pˇrij´ıman´ ym parametrem je ukazatel na vektor datov´eho typu double, do kter´eho bude postupnˇe zapisov´ana poˇc´ıtan´a trajektorie ˇc´astice. Motivac´ı pro pos´ıl´an´ı ukazatele na jiˇz existuj´ıc´ı pole je u ´spora ˇcasu. Pro velk´e mnoˇzstv´ı bod˚ u trajektorie by kop´ırov´an´ı vektoru znamenalo zbyteˇcnou ˇcasovou z´atˇeˇz. Na zaˇca´tku simulace je proveden glob´aln´ı pr˚ uzkum, ve kter´e buˇ nce se ˇca´stice nach´az´ı. Pokud v ˇza´dn´e, na standardn´ı v´ ystup je vyps´ana chybov´a hl´aˇska V´ ychoz´ı ” pozice ˇca´stice je mimo dom´enu“. Pak n´asleduje posouv´an´ı ˇca´stice ve smyslu modifikovan´e Eulerovy metody uveden´e v kapitole 2, dokud ˇc´astice neopust´ı dom´enu a aktu´aln´ı ˇcas bude menˇs´ı neˇz celkov´ y ˇcas simulace. Dalˇs´ı zastavovac´ı krit´erium je pˇr´ıliˇs mal´a zmˇena polohy ˇca´stice. Motivac´ı pro toto krit´erium je moˇznost, ˇze ˇca´stice skonˇc´ı v m´ıstˇe s nulovou rychlost´ı a pak simulace jiˇz nem´a smysl.

2.5.3

TracingLocalExact

Simulaˇcn´ım j´adrem lok´alnˇe exaktn´ı metody je funkce TracingLocalExact. Funkce pˇrijme parametrem instanci tˇr´ıdy TTetBasicVol, poˇca´teˇcn´ı pozici ˇc´astice, celkov´ y ˇcas simulace a parametry metody imax , kmax , τ . Jako u funkce TracingEuler, i zde je vektor pro zapisov´an´ı trajektorie pˇred´an v parametru funkce ukazatelem. Simulace zaˇc´ın´a nalezen´ım aktu´aln´ı buˇ nky ˇca´stice (v opaˇcn´em pˇr´ıpadˇe je vyps´ana chybov´a hl´aˇska V´ ychoz´ı pozice ˇc´astice je mimo dom´enu“ a simulace ” konˇc´ı). T´ımto se simulace dostane do cyklu, kde je dle typu buˇ nky vybr´ana metoda pro zjiˇstˇen´ı nov´e pozice. Pro zvl´aˇstn´ı buˇ nku se pouˇzije modifikovan´a Eulerova metoda viz kapitola 2. Pro paraleln´ı buˇ nky se jedn´a o nalezen´ı pr˚ useˇc´ıku kaˇzd´e stˇeny s pˇr´ımkou proch´azej´ıc´ı aktu´aln´ı pozic´ı ˇca´stice dan´ ym smˇerem. Z tˇechto kandid´at˚ u na bod u ´niku se vybere ten nejbl´ıˇz aktu´aln´ı pozice ˇc´astice. Re´aln´e a komplexn´ı buˇ nky v k´odu spl´ yvaj´ı v jeden pˇr´ıpad liˇs´ıc´ı se pouze ve dvou funkc´ıch. Jedn´a se o funkce NewtonIter ˇci NewtonIter komplex a funkce pro vyˇc´ıslen´ı bodu trajektorie za pˇredpokladu zn´am´eho ˇcasu u ´niku. Zp˚ usobem popsan´ ym v kapitole 2 se vypoˇcte pro kaˇzdou stˇenu potenci´aln´ı ˇcas u ´niku a z nich se vybere ten minim´aln´ı. Pokud pro ˇza´dnou stˇenu nedostaneme kladn´ y ˇcas, prohl´as´ıme buˇ nku za zvl´aˇstn´ı a po ovˇeˇren´ı kompatibility pozice ˇca´stice a aktu´aln´ı buˇ nky pokraˇcujeme dalˇs´ı iterac´ı hlavn´ıho cyklu. V pˇr´ıpadˇe, ˇze jsme alespoˇ n pro jednu stˇenu buˇ nky obdrˇzeli kladn´ y ˇcas u ´niku, jeˇstˇe mus´ıme ovˇeˇrit, ˇze pˇr´ısluˇsn´ y bod u ´niku skuteˇcnˇe n´aleˇz´ı stˇenˇe ˇctyˇrstˇenu a nejen rovinˇe stˇenou urˇcen´e. Nakonec posuneme ˇca´stici do spoˇcten´eho bodu u ´niku, zmˇen´ıme aktu´aln´ı buˇ nku na pˇr´ısluˇsnou sousedn´ı buˇ nku a

30

cyklus opakujeme. Zastavovac´ı krit´eria jsou shodn´e s krit´erii v jiˇz zm´ınˇen´e funkci TracingEuler.

2.5.4

TiskCastice

D˚ uleˇzitou funkc´ı je TiskCastice, kter´a dostane parametrem ukazatel na vektor dataStore datov´eho typu double obsahuj´ıc´ı souˇradnice bod˚ u a dalˇs´ı vektor urˇcuj´ıc´ı, kolik trajektori´ı vektor dataStore obsahuje a jakou maj´ı d´elku. Posledn´ım parametrem je jm´eno souboru, do kter´eho tato data vytiskne ve form´atu .vtk.

2.5.5

Main

Hlavn´ı tˇelo programu obsahuje ovl´adac´ı panel pro volbu poˇctu ˇc´astic, parametr˚ u metod ˇci n´azv˚ u vstupn´ıch ˇci v´ ystupn´ıch soubor˚ u. Pot´e je inicializov´ana instance tˇr´ıdy TTetBasicVol s n´azvem mesh, kter´a naˇcte data ze zvolen´eho souboru a provede pˇredv´ ypoˇcet pomoc´ı metody BuildData. N´asleduje cyklus pro simulov´an´ı trajektorie ˇca´stic jak lok´alnˇe exaktn´ı metodou, tak modifikovanou Eulerovou metodou. Nakonec jsou vygenerov´any dva v´ ystupn´ı soubory pomoc´ı funkce TiskCastice a vyps´any u ´daje o pr˚ ubˇehu simulace (poˇcet bod˚ u trajektori´ı, ˇcas simulace jednotliv´ ych metod ˇci pr˚ umˇern´ y poˇcet iterac´ı Newtonovy metody).

2.6

Uˇ zivatelsk´ a pˇ r´ıruˇ cka

Na pˇriloˇzen´em CD se nach´az´ı adres´aˇr Chluap“, ve kter´em je stejnojmenn´ y ” projekt spustiteln´ y v softwaru Visual Studio C++ 2010. Pro kompilaci projektu je nutn´e do Visual Studia pˇridat volnˇe dostupnou knihovnu Eigen [5]. Jedn´a se pouze o experiment´aln´ı verzi programu. Parametry se zad´avaj´ı pˇr´ımo v k´odu a to v hlavn´ım souboru MainFile.cpp v hlavn´ı funkci Main. Soubory, ze kter´ ych m´a program ˇc´ıst, mus´ı b´ yt v adres´aˇri, kde se tak´e nach´az´ı zdrojov´ y soubor MainFile.cpp. V t´eto sloˇzce jsou jiˇz pˇripraven´e pˇr´ıklady vstupn´ıch soubor˚ u 0-178.vtk, 0-182.vtk, 0-890.vtk a 0-930.vtk. V´ ysledek simulace se opˇet uloˇz´ı do sloˇzky, kde se nach´az´ı i zdrojov´e soubory, a je moˇzn´e je zobrazit napˇr´ıklad volnˇe dostupn´ ym softwarem ParaView [6]. Pˇr´ıklady soubor˚ u *.vtk vygenerovan´ ych programem jsou ve sloˇze Vysledky“ i s uk´azkou, jak vypadaj´ı po zobrazen´ım v ” ParaView. Pozn´amka. Je vhodn´e si v ParaView souˇcasnˇe s v´ ysledkem simulace zobrazit i dom´enu, ve kter´e byla simulace prov´adˇena.

31

Literatura [1] a Alexandru M. Morega, A. A. D. Blood Flow – Vessel Interaction in a Subclavian Aneurysm. URL http://www.ima.ro/conf_cci/e-book_11/ pdf_file/dobre_art_cci_11.pdf. [2] Bhavikatti, S. (2005). Finite Element Analysis. New Age International (P) Limited, Publishers. ISBN 9788122415896. ˇka, M. a Sopko, B. a. S. L. (2000). Mechanika kontinua. Druh´e [3] Brdic opraven´e vyd´an´ı. ACADEMIA, Praha. ISBN 80-200-0772-5. [4] Butcher, J. (2008). Numerical Methods for Ordinary Differential Equations. Wiley. ISBN 9780470753750. [5] Guennebaud, G., Jacob, http://eigen.tuxfamily.org.

B.

a

kol.

(2010).

Eigen

v3.

[6] Henderson, A. (2007). The ParaView Guide: A Parallel Visualization Application. [7] Kipfer, P., Reck, F. a Greiner, G. G.: Local exact particle tracing on unstructured grids. Computer Graphics Forum, 22, 133–142. ˇ , J. (2004). Obyˇcejn´e diferenci´aln´ı rovnice v re´aln´em oboru. Uˇcebn´ı [8] Kofron texty. Karolinum. ISBN 9788024609461. [9] Macpherson, G. B., N. N. a Weller, H. G. (2009). Particle tracking in unstructured, arbitrary polyhedral meshes for use in cfd and molecular dynamics. Commun. Numer. Meth. Engng., pages 263––273. [10] Nielson, G. M. a Jung, I.-H. (1999). Tools for computing tangent curves for linearly varying vector fields over tetrahedral domains. IEEE Transactions on Visualization and Computer Graphics, 5(4), 360–372. ISSN 10772626. doi: 10.1109/2945.817352. URL http://dx.doi.org/10.1109/2945. 817352. [11] Stroustrup, B. (2000). The C++ Programming Language. AddisonWesley Longman Publishing Co., Inc., Boston, MA, USA, 3rd edition. ISBN 0201700735. [12] Yongling, L. Wu Guangqiang, X. Numerical Simulation of the External Flow Field Around a Bluff Car. URL http://web.mscsoftware.com/ support/library/conf/auto00/p02300.pdf.

32

Seznam obr´ azk˚ u 1

Uk´azka simulace obt´ek´an´ı formule vzduchem.

. . . . . . . . . . .

2

1.1

Demonstrace rozd´ıl˚ u mezi jednotliv´ ymi kˇrivkami popisuj´ıc´ımi vektorov´e pole v pr˚ ubˇehu ˇcasu. Kr´atk´e ˇca´rky pˇredstavuj´ı smˇer vektoˇ rov´eho pole. Cerven´ a kˇrivka je trajektorie, modr´a kˇrivka je emisn´ı ˇca´ra a pˇreruˇsovan´e ˇca´ry ukazuj´ı proudnice. . . . . . . . . . . . . .

4

2.1

2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10

Zn´azornˇen´ı situace pro 2D pˇr´ıpad s troj´ uheln´ıkovou s´ıt´ı. V´ ychoz´ı pozice ˇca´stice A se nach´az´ı v buˇ nce 1 a c´ılem je zjistit, ve kter´e buˇ nce se nach´az´ı c´ılov´a pozice B. Pro stˇeny f1 a f2 λA ∈ [0,1], kdeˇzto pro stˇenu f3 je λA < 0. Menˇs´ı λA n´aleˇz´ı stˇenˇe f1 , a proto je ˇca´stice posunuta do bodu P1 a aktu´aln´ı buˇ nka zmˇenˇena na buˇ nku 2. Analogick´ ym postupem se ˇc´astice posune do bodu P2 s okupovanou buˇ nkou 5. Zde ˇza´dn´a hodnota λA nen´aleˇz´ı intervalu [0,1], ˇc´ımˇz je algoritmus u konce. . . . . . . . . . . . . . . . . . . . . . . 9 Uk´azka pˇr´ıpadu, kdy Newtonova metoda teˇcen nekonverguje. . . . 18 Uk´azka pˇr´ıpadu, kdy Newtonova metoda teˇcen konverguje ke ˇspatn´emu koˇrenu vlivem pˇr´ıtomnosti stacion´arn´ıho bodu. . . . . . . . . . . . 18 V´ ypoˇcet ˇcasu u ´niku klasickou Newtonovou metodou. . . . . . . . 20 V´ ypoˇcet ˇcasu u ´niku modifikovanou Newtonovou metodou. . . . . 20 Funkce f (t), kde selhala klasick´a Newtonova metoda. . . . . . . . 21 Srovn´an´ı MEM ∆t = 0,001 (vlevo) a ryz´ı LEM bez pˇrepnut´ı na MEM (vpravo). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 MEM ∆t = 0,00002. . . . . . . . . . . . . . . . . . . . . . . . . . 23 Srovn´an´ı ˇcasov´e sloˇzitosti simulace pro 100 ˇc´astic. . . . . . . . . . 24 Nahoˇre MEM ∆t = 0.001, uprostˇred MEM ∆t = 0.0001, dole LEM bez zvl´aˇstn´ıch bunˇek. Stoj´ı za povˇsimnut´ı, ˇze na horn´ım obr´azku zcela chyb´ı nejvyˇsˇs´ı smyˇcky“ trajektorie. . . . . . . . . . . . . . 25 ”

33