Numerické metody II Jan Schee April 17, 2018 Contents 1 Parciální diferenciální rovnice 1 1.1 Numerické metody řešení................ 1 1.1.1 Problém rovnováhy................ 3 1.1.2 Problém šíření.................. 8 1.1.3 Numerické schéma Lax-Wendroff: o(Ax2) a o(At2) 28 2 Eigen-value problém 33 2.1 Schrôdingerova rovnice.................. 33 2.2 Vlastní frekvence schwarzchildovy černé díry ..... 41 3 Generátor náhodných čísel a rozdělovači funkce 42 3.1 Číslo 7T........................... 44 3.2 Náhodná procházka.................... 46 4 Vícerozměrné integrály 49 4.1 Opakování jednorozměrné integrace........... 49 4.2 Monte-Carlo........................ 51 1 1 Parciální diferenciální rovnice Rovnice určující pohyb nebo stacionární stav fyzikálnych systémů lze shrnout do generického, kvazilieárího tvaru d2f d2f d2f / p df df\ + 2fli2^-^- + «22^- + F [x,y,j,—, — \ = 0 (1) oxz oxoy oyz \ ox oy J kde ciij jsou konstanty vedoucí ke klasifikaci PDR prostřednictvím parametru A = auaii — Je-li A < 0 mluvíme o hyperbolické rovnici (vlnová rovnice), pro A > 0 dostáváme eliptickou rovici (Poissonova rovice), a pro A = 0 obdržíme parabolickou rovnici (rovnice difúze). 1.1 Numerické metody řešení Fyzikální problémy, které popisují parciální diferenciální rovnice, můžeme rozdělit do tří kategorií: • problém rovnováhy - eliptické PDR, např. Laplaceova nebo Poissonova rovnice, • problém šíření - parabolické a hyperbolické rovnice, např. rovnice difúze nebo vlnová rovnice, • Eigenvalue problém - příkladem je Schrodingerova rovnice a její stacionární stavy. Obecně mají PDE nekonečné monožství řešení. Když stanovíme počáteční a okrajové podmínky tak z toho nekonečného množství vybíráme právě jedno. Při řešení problému rovnováhy stanvujeme pouze okrajové podmínky, protože zde nehledáme vývoj v čase. Naproti tomu problém šíření 2 reprezuntuje úlohy kdy se fyzikální systém vy vyjí v čase a tedy potřeba specifikovat počáteční podmínky. Fyzikální jevy jsou zpravidla řízeny PDE druhého řádu. Při formulaci matematického modelu fyzikálního problému jsou uvažovány tři typy okrajových podmínek • Dirichletova okrajová podmínka - na hranicích integrační oblasti je specifikována hodnota hledané funkce f(xb,yb), kde body (xb,Vb) představují hranici. • Neumanova okrajová podmínka - je specifikována hodnota derivace funkce df/dň ve směru normály ň k hranici. • Smíšené okrajové podmínky - jedná se o lineární kombinaci hodnot a derivací ve směru normály okraje domény integrace, tj. af{xb, Vb) + bdf/dň(xb, yb). Zastavme se na okamžik u Neumanovy okrajvé podmínky. Derivaci ve směru ň určíme ze vztahu d f _ d f d f d f - = 71 ■ - = t).™- -+- 71.-.- (2) Dále, necht je hranice určená implicitně rovnicí F(x,y) = 0. (3) Normálový vektor k tého hranici je potom určen vztahem ň = ± 1 (4) 2 ex ^ kde je dF/dx dF/dy' (5) 3 1.1.1 Problém rovnováhy Numerické schéma které řeší problém rovnováhy vysvětlíme na příkladu rozdělení teploty ve vodiči. Uvažujme vodič se čtvercovým průřezem 10 x 10 cm a tlouštce dz « 1. Horní hrana vodiče je udržována na stálé teplotě T = 100° C a zbylé tři hrany udržujeme na teplotě 0° C. Rozdělení teploty ve vodiči určuje Laplaceova rovnice — — - 0 (6) dx2 dy2 K vyřešení této úlohy využijeme metodu konečných diferencí. Doménu D(x, y) nad kterou hledáme rozdělení teploty T(x, y) rozdělíme na knečnou mřížku NI x N J bodů, přičemž souřadnice x a, y budou x = xmin + idx ay = ymin + jdy, (7) kde je Ó.X = [Xmax — Xmin) / NI ä dy = {y max ~ ymin) / N J. (8) Rovnici (6) převedeme na diferenční tvar takto. Rozvineme ře šení T v okolí bodu (x, y) do Taylorovy řady dT 1 dT T(x + Ax,y) = T(x,y) + —Ax + -—Ax2 + o(Ax)3, (9) ox 2 ox dT 1 dT T(x - Ax, y) = T(x, y) - —Ax + -—-Ax2 - o(Ax)3. (10) ox 2 ox Poslední dvě rovnice sečteme a vyjádříme d2T/dx2, obdržíme vztah d2T ^T(x + Ax,y) + T {x - Ax, y) - 2T(x, y) dx2 ~ Ax2 ■ K ) 4 Analogicky, pro d2T/dy2 dostaneme vztah d2T T(x,y + Ay) + T(x,y-Ay)-2T(x,y) (12) Funkci T(x,y) budeme vyčislovat pouze v uzlových bodech, tj. T{x^ ž/) = T{xmin + iáx,ymin+ jdy), (13) přejdeme k novému značení T(x,y)^Tij. (14) Dosazením vztahů (11) a (12) do (6) a po osamostatnění Tíj dostaneme rovnici 1 ThJ 2(1 + /c2) Ti+ljj + Tí-xj + k (TiJ+1 + Ti,-{) (15) kde jsme už využili přeznačení (14) a zavedli k = Ax/Ay. K nalezení řešení Tíj na celé mřížce NI x N J použijeme iterační schéma. Štandartne se diskutují tři iterační schémata: • Jacobiho, • Gaussovo-Seidelovo, • SOR (Succesive-Over-Relaxation). Zde se omezíme na poslední dvě iterační schémata. Gaussovo-Seidelovo (GS) schéma Zavedeme nový index m kterým budeme číslovat jednotlivé iterace. GS schéma pro rovnici (15) bude mít tvar (k rovnici (15) přičteme ±Tíj) Trn+1 = Trn ^ ATrn+l 5 kde je a rpm+1 2(1 + k2) ZTříj + + *2PJj+i + ^í!) - 2(1 + k2)T, (17) ^ V tomto schématu, slouží k určení hodnoty T[^+1 jak hodnoty předchozí iterace m tak i aktuální iterace m + 1. To znamená, že k výpočtu vystačíme s jedním dvourozměrným polem ve kterém jsou jeho prvky postupně přepisovány, v našem případě zleva-doprava a sdola-nahoru (viz obr.l) m 6 3 X 0 "D -*—> °3 1 á jm r i P f i g P f l i P f | i P f T-A77+1 ' / J P f P^ ^ l g P^ ^ l á P^ ^ | á P^ ^ jm+1 '/J-1 | i P f | i P^ ^ l é P^ ^ l 4 P^ ^ | á P^ ^ | i P^ ^ | á P^ ^ k á P^ ^ P^ ^ p^ ^ P^ ^ É á růst indexu i —> Figúre 1: Schéma rozdělení integrační domény D(x,y) a směr generování řešení GS metodou. Aktuálně určovaná teplota v bodě i, j (červený bod) je určená z b... Okrajové podmínky naši úlohy lze reprezentovat následovmě T0j=TNI-ltj = 0 V j = 0,...,NJ-2 (18) 7 a rij0 = O arijJVj-i = 100 V i = 0,...,iVJ-l. (19) Iniciální stav iteračního procesu nastavíme na 0, tj., Tíj = 0 V i = 1,... ,NI-2 a j = 1,.. .,NJ-2. (20) Iterační proces zastavíme tak, že v každém iteračním kroku určíme A = MAX(AT™+1) (21) a srovnáme se zvolenou přesností e, tzn. m + 1 iterace je poslední pokud je splněna podmínka A < e. (22) SOR schéma SOR iterační schéma vychází ze schématu Gausse-Saidela přidáním over-ralaxačního parametru oj ^+1=?íj+wAIIj+1- (23) Je-li uj = 1, dostaneme se zpět ke GS schématu. Maximální rychlost konvergence je dosažena vhodnou volbou parametru ujopt G [1.0, 2.0]. SOR schéma a smíšené okrajové podmínky V předchozích dvou úlohách jsme při hledání rozdělení teploty použili Dirichletovy okrajové podmínky, tj. na hranici integrační oblasti jsme určili hodnoty funkce T. To odpovídalo situaci, kdy jsme včechny čtyři hrany destičky udržovali na konstantních teplotách. Nyní náš experiment poněkud změníme. 8 1.1.2 Problém šíření Parciální diferenciální rovnice spadající do této kategorie jsou parabolické a hyperbolické. Typickým představitelem parabolických PDR je rovnicd vedení tepla a typickým představitelem hyperbolických PDR je vlnová rovnice. Rovnice vedení tepla Uvažujme následující fyzikální pokus. Máme měděnou tyč délky L (třeba L = 2m) a průměru d « L (třeba d = 2cm). Tyč tepelně izolujeme od okolí , vyjma jejích konců (teplo vstupuje/vystupuje pouze na koncích tyče). Tuto tyč vložíme do prostředí s teplotou Tq a necháme ji tam tak dlouho, dokud nedojde k tepelné rovnováze s okolím (tyč bude mít konstantní teplotu Tq). Po vyjmutí z prostředí, v čase t = 0 přiložíme ke koncům tyče dva tepelné elementy s teplotami Tx a T2 (např. Tx = 0°C a T2 = 50°C). Monitorujeme jak se vy vyjí rozložení teploty v tyči v čase t. Abychom mohli kvantitativně předovědět výsledek pokus musíme vytvořit jeho matematický model. Tento model se skládá ze tří složek: • PDR - řídí vývoj teploty v tyči, • BCs - okrajové podmínky popisují fyzikální problém na krajích tyče, • ICs - počáteční podmínky popisují rozložení teploty v tyči na začátku pokusu. 9 PDR Rovnici vedení tepla lze odvodit následující úvahou. Uvažujme tenkou tyč o příčném průřezu A. Určíme celkovou změnu tepla v tyči áQi/át v části od x do x + Ax. Zřejmě bude platit dQ - = Qt(x) + Qt(x + Ax) (24) dt kde je Qt(x) = -kA^^ (25) a ^ , * \ , . dT(x + Ax) Qt(x + Ax) = kA—-J- 26 dx Celkové vnitřní teplo daného úseku materiálu s hustotou p bude dQ/ = cTdra = cApTds Qi = í + cATds. (27) Pro jeho celkovou časovou zenu nakonec obdržíme vztah cpAp fX+AX Ttds = k A [Tx{x + Ax) - Tx{x)]. (28) J X Integrál na levé straně rozřešíme pomocí věty o střední hodnotě, podle které existuje číslo x < £ < x + Aa; takové, že platí / Tt(s,t)ds = Tx(t,t)Ax. (29) J x Dostáváme tak rovnici rr(c+\ uTx(x + Ax-Tx(x)) cpTt(£,t) = k--r-. (30) Ax 10 &(x)«V A a(x+Ax)>0 x x+Ax Figuře 2: Schéma elementu tyče, oblast mezi x a x + Ax, které analyzujeme tok tepla. V limitě Aa: —> 0 je zřejmě £ —> x a dostaneme rovnici vedení tepla v homogenním materiálu s koeficientem vodivosti a ve tvaru du d2u a (31) dt ~ 5a;2 kde je 0 < a: < L, 0 < t < oo a • u(x,t) - teplota v tyči, v místě x a čase t. • du/dt - rychlost změny teploty v místě x a čase t. • d2u/dx2 - konkávnost teplotního profilu u(x,t). Význam poslední zmíněné veličiny snadno nahlédneme, když rozvineme u(x,t) do Taylorovy řady při fixovaném čase t, tj. c)u 1 cféu u(x + Ax,t) = u(x,t) + —Aa: + --^ (Aa:)2 + o(Aa:)3, (32) u(x — Ax, t) c)u 1 cféu u(x,t) - ^-Ax + - — (As)2 - o(Aa:)3. (33) dx ' 2 dx2 Sečtením obou rovnic a osamostatněním d2u/dx2 dostaneme vztah d2u _ 2 dx2 Ax2 u(x + Ax, t) + u(x — Ax, t) u(x, t) (34) 11 Z této rovnice odvodíme následující výsledky: • pokud bude v místě x a case t teplota u(x, t) 0 a celkový tepelný tok bude kladný, • pokud bude v místě x a case t teplota u(x, t) =průměrná teplota v sousedních bodech, pak bude d2u/dx2 = 0 a tepelný tok bude nulový, • pokud bude v místě x a case t teplota u(x, t) >průměrná teplota v sousedních bodech, pak bude d2u/dx2 < 0 a celkový tepelný tok bude záporný. Jinak řečeno, když bude teplota v bodě x a čase t větší než je průměrná teplota dvou sousedních bodů x + Ax a x — Ax pak teplota v bodě x bude klesat a to s rychlostí a2d2u/dx2. BCs Ze zkušenosti víme, že fyzikální problémy jsou nějakým způsobem ohraničeny. Hranicemi našeho experimentu jsou konce tyče x\ = 0 a x2 = L kde udržujeme stálé teploty T\ = 0°C a T2 = 50°C, tj. BCs=f"íľ'?r? <35> [ u(L,t) = T2 pro 0 < t < 00. ICs Dále je zřejmé, že fyzikální systém se musí vyvíjet z výchozího známého stavu. Tento stav je určen, zpravidla, pro t = 0. V našem případě představuje počáteční stav rozložení teploty v tyči v čase t = 0, tj. ICs: u(x, 0) = T0 pro 0 < x < L. (36) 12 Numerické řešení K nalezení numerického řešení problému, který matematicky popsán systémem du 2 d2u dt-aä^' (37) BCs=f"íľ'?r? <38> [ u(L,t) = T2 a ICs: u(x, 0) = T0 pro 0 < x < L, (39) použijeme opět metodu konečných diferencí Tyč délky L rozdělíme na N stejných dílů, každý o velikosti Ax = L/N (40) a časový integrační krok bude At. Budeme určovat řešení v uzlových bodech u(x, t) = u(xo + iAx, kAt) (41) a zavedeme označení u(x0 + iAx, kAt) u\. (42) Dále nahradíme parciální derivaci aproximativními vztahy duí uÍ+l - uí dt At a (43) Ôx2 ~ (Ax)2 " V ' 13 Ze vztahů (43) a (44) dosadíme do rovnice (37) a osamostatníme uf+1, obdržíme vztah u^+1 = u\ + 2a2 At u i+1 + Utl k 2 1 (45) Právě odvozené integrační schéma se nazývá FTCS (Forward Time Centered Space). Figuře 3: Snímky časového vývoje rozdělení teploty v tenké měděné tyči. Vlnová rovnice Určeme, jak se bude šířit rozruch v jednorozměrné struně délky L a na vodní hladině. Rovnicí, která popisuje časový a prostorový vývoj iniciální perturbace, je vlnová rovnice. Pro veličinu u = u(x, t), která zde má roli amplitudu výchylky, vlnovou rovnici píšeme ve tvaru d2u 2^2"u dť2=cd^2 ( } kde parametr c je rychlost šíření vlny. 14 Figuře 4: Struna upevněná v bodech A a B. Ukážeme si jednoduché odvození rovnice pro šíření příčného vlnění struny. Strunu napneme mezi body A a, B (viz obr. 4). Vychýlíme strunu z rovnovážné polohy. Ve struně se vyvolané vnitřní napětí se bude snažit uvést strunu opět do rovnovážné polohy. Analyzujme detailně okolí bodu u(x). 15 Figure 5: Detail struny. Soustředíme hmotnost struny délky AI do bodu m, budou na tento hmotný bod síly T a —T, takže výsledná působící síla bude (obr. 5) Tnet = T(cos ol<2 — cos ol\)ex + T(sin ol 0 přejde na rovnici otz oxz kde jsme zavedli rychlost šíření příčné vlny ve struně c = yjT/cr. V našem experimentu aplikujeme na konce struny dva typy okraji vých podmínek: • von Neumanova okrajová podmínka - struna je pevně ukotvena na koncích, u(0,t) = u(L,t) = 0. (53) • Dirichletova okrajová podmínka - rozruch (vlna) se odráží od levého (pravého) okraje nádrže, du , du \x=0 o \x=L ox ox Iniciální rozruch bude popsán Gaussovou funkcí 0. (54) u(x,0) = exp ( ———2~^~ ) ' (55) kde bod xq = L/2. 17 Diskretizace úlohy Definiční obor nad kterým hledáme řešení opět rozdělíme na konečnou mříž a vlnovou rovnici převedeme na diferenční tvar. Druhé derivace podle x a, t aproximujeme vztahy a d2u dx2 d2u /±x2 \Uj+l + Uj-1 dt2 - At2 O ^U3 Dosazením do vlnové rovnice (46) dostaneme rovnici pro odhad hod-noty Uj ve tvaru 2u: 2ut (56) (57) u k+l 2ukAl-(32) .fc-1 (58) kde jsme zavedli parametr (3 = cAt/Ax. Pozornou analýzou tohoto vztahu zjistíme, že k určení hodnoty uk+1 potřebujeme znát hodnoty uk i uk i. Toto schéma je dvouúrovňové (viz obr. 6). .fc-i k+1 k-1 J-1 J Figuře 6: Schématický diagram integračního schématu, které je dáno vztahem (58). 18 Počáteční podmínky Protože je integrační schéma (58) dvouúrovňové, musíme specifikova počáteční podmínky ve dvou časových úrovních. Necht má první časová úroveň index k = 0 a druhá časová urveň index k = 1. Veličinu u v uzlových bidech mřížky nastavíme v první časové úrovni na Gaussův impulz, tj. Zbývá vyřešit, jak určit hodnoty u ve druhé časové úrovni, tj. itj. Pro fixovanou pozici x rozvineme u v okolí první časové úrovně do prvního řádu v At, tj. d \ u)~u°+(-^j Ač, pro j = {1,2,...,iV-2}. (60) "3 — 3 Z vlnové rovnice (46) plyne platnost vztahu d (du\ nd2u du nd2u . . _ dt (m) =cd^^m^ cVAí- (61) Dosazením posledního vztahu do (60) obdržíme výraz pro inicializaci veličiny u ve druhé časové úrovni, tj. dostaneme 1 =^ + /32(^+1+^_1-2íx°), pro j = {1,2,..., iV-2}. (62) Zde jsme za d2udx2 dosadili z diferenčního vztahu (56). Ještě poznamenejme, že počáteční podmínky jsou stejné pro oba naše příklady šíření vlny, tj. pro strunu a vodní hladinu. 19 Okrajové podmínky Diskutujme nyní nastavení okrajových podmínek v případě struny upevněné na koncích, tj. v bodech x = 0a,x = Lve kterých je u(0,t) = u(L,t) = 0. V diskretizačním formalizmu mají okrajové podmínky tvar uo = uN-i = 0? Pr0 všechny hodnoty k. (63) Dostáváme se k druhé úloze a to k šíření vlny na vodní hladině. Hladina je ohraničena břehem na a; = 0 a a; = L. V těchto bodech platí podmínky (ŕto)o (dajjy-i °" ^ Rozvineme řešení u v okolí bodů j = 0 a j = N — 1 do Taylorovy řady do prvního řádu v Aa; a obdržíme u\ = ul+(^)k &x (65) o a .k „.k J-i ~ \ TTľ I jv-l dx, u%_2 = u%_1-[ — ) Ax. (66) Uplatněním podmínek (64) tak dostáváme hodnoty u v krajních bodech intervalu Uq = Ui a u%_1 = u%_2, pro všechny hodnoty fc. (67) 20 Figuře 7: Snímky šíření poruchy v napnuté struně při použití von Neumanovy okrajové podmínky. 21 Figuře 8: Snímky šíření vlny na napnuté struně při použití Dirichle-tovy okrajové podmínky. Vlnová rovnice ve 2D Předchozí analýzu lze snadno rozšířit na 2D případ šíření vln (vodní hladina, blána bubnu). Vlnová rovnice má v tomto případě tvar (68) 22 Při konstrukci odpovídající diferenční rovnice rozvineme, opět, řešení v okolí bodu (£, x, y) do Taylorovy řady s přesností do 2. řádu v Ax. Diferenční rovnice je potom dána výrazem , n i n _ on At2 v «J ÍJ ÍJ/ ~ V Ax2 + 2"^) (69) odkud vyjádříme u™^1 ui,~j — ^Ui,j Ui,j + Px {^i+lj + ui-i,j ^ui,j) + " 2 -1. (87) Parametry Ax, At a v jsou potom svázány vztahem \vAt/Ax\ < 1 a v < 0. (88) 25 Numerické schéma FTCS(Forward Time Centered Space): o(Ax2), o(At) První derivaci v prostorové souřadnici získáme z Taylorova rozvoje v okolí bodu t, x ve směru ±Ax, tj. (Ji I <+i = u7l + —\7lAx + o(Ax2), (89) J J ox J un3_x = u]-^\]Ax + o(Ax2), (90) odkud plyne aproximativní vztah pro první parciální derivaci n _ n uj-i dxlj 2Ax y J Po dosazení do advektivní rovnice (??) dostaneme FTCS integrační schéma Aplikujeme na toto schéma von Neumanovu analýzu stabitity, tj. dosadíme do tohoto schématu poruchu (76) a určíme \G\. Nejprve pro G obdržíme vztah G = l-iasiYi(kAx) (93) a následně pro \G\2 obdržíme výsledek \G\2 = \GG*\ = l + a2sm2(kAx) > 1, (94) který implikuje nestabilitu FTCS schématu. Tuto nestabilitu řeší následující schéma. 26 Numerické schéma Lax-Friedrichs: o(Ax2), o(At2) Nestabulitu FTCS eliminujeme, když v rovnici (92) nahradíme první člen na pravé straně prostorovým průměrem sousedních bodů, t j J 2 Nové schéma pro naši advektikví rovnici, Lax-Friedrich, má potom tvar + i vAt / \ n+l „ J + l j-l _ fn _ un\ (g6) i ~ 2 2Ax V^'+1 Opět podrobíme toto numerické schéma analýze stabitity a pro zesilovací parametr G obdržíme vztah G = cos(kAx) — iasin(kAx) (97) a pro GG* dostaneme výraz \G\2 = \GG*\ = 1 - sm2(kAx)(l - a2). (98) Pokud bude \a\ < 1 tak je \G\2 < 1 a schéma je stabilní. Modifikací členu u™ jsme do systému vložili numerickou disipaci, jak následně ukážeme. Rovnici (96) stačí upravit na tvar - "i - «ľ + 3+ 2 3 ~ 2ÄÍ " """O (") resp., lépe, na tvar ~ ^ ~ f ^+1 7 ^ + ^+1 + Ul~l ~ 2UlAt (100) At V 2Aa; / 2At2 v ' což je v limitě A —> 0 převedeno na rovnici du du 1 d2u . di + Vd-X = 2WAt (101) 27 Výraz na pravé straně představuje numerickou disipaci. Lax-Friedrichs je tedy schéma pro poslední uvedenou rovnici. Pro At « 1 jej lze s opatrností aplikovat na rovnici (72). Numerické schéma Leap frog: o(At2), o(Ax2) Druhý řád v At i v Ax vede na další numerické schéma - Leap frog. Advektivní rovnice přejde na tvar _J-^— ~ _v J+1 (102) 2At 2Ax v ' který vede na schéma uTl - «r1 - K1 - ■ (103) Figure 9: Leap-frog schéma je v časové oblasti dvouúrovňové. Když chceme určit funkci u v čase n +1 (šedé body), musíme znát hodnutu funkce v časových úrovních n (červené body) a n — 1 (červené body). Navíc touto metodou získáváme dvě na sobě nezávislá řešení (bílé a černé bloky.) 28 Schéma leap-frog taky podrobíme analýze stability. Pro parametr G dostaneme vztah a G = ±yj 1 — a2 sm2(kAx) — ia sm(kAx) GG* = l-a2 sm2(kAx) + a2 sm2(kAx) = 1. (104) (105) Schéma Leap-frog je mezně stabilní. K určení hodnoty it- musíme znát hodnoty u v časových úrovních n a n — 1. Jedná se tedy o dvouúrovňové integrační schéma. Hlavní nevýhodou je však existence dvou nezávislých řešení, jak je vidět na obrázku (obr. 9). 1.1.3 Numerické schéma Lax-Wendroff: o(Ax2) a o(At2) Zachovat vysokou míru stability a přitom mít k dispozici časově jedno úrovňové schéma nám umožní Lax-Wendroff schéma. - n+1 n+1/2 j-1 j-1/2 j+1/2 n j+1 Figuře 10: Schéma numerické metody integrace Lax-Wendroff. Odd-vozené pomocí " 1/2"-Lax-Friedrich schématu a "l/2"-Leap frog schématu. 29 Toto schéma odvodíme pomocí "l/2"-Lax-Friedrichs a "1/2"-Leapfrog schémat: "1/2"'-Leapfrog: Standardní Leapfrog má tvar = »r1 - v£ (u1+1 - uU) . (106) Pro Ať —>• Ač/2 a Ax —> Ax/2 přejde toto schéma na tvar 7/n+l _ vn 7;AV2 ( n+l/2 n+1/2 A^/2 0-+V2 "j-i/i Nyní stačí určit u™*^ a u^j^- K tomu použijeme "1/2"-Lax-Friedrichs. " 1 /2"-Lax-Friedrichs: Nejprve obaz výše zmíněné body určíme pomocí "l/2"-FTCS. Obdržíme ^+1=^-^(^-^) (108) a pro bod (n + 1/2, j — 1/2) dostaneme a pro bod [n + 1/2, j + 1/2) dostaneme - nj+l/2 U Az/2 W + 1 J / " 1 } 30 Aplikací Lax-Friedrichs na poslední dva vztahy dostaneme konečné výrazy n+l/2 _ «" + Ar ,.n a n+l/2 Uj+1+Uj At / \ ^+1/2 = 2 " ^ K+1 " M"J • (U2) Poslední vztahy dosadíme do (107) a obdržíme finální schéma Lax-Wendroff ve tvaru un+l = un _ ^ X i J 2A (113) Opět je ihned vidět, že Lax-Firedrichs vnáší do nového schématu numerickou disipaci, kterou reprezentuje podlední člen na pravé straně. Podrobme nyní toto schéma analýze stability. Pro chybu e bude zesilovací faktor G dán výrazem G = 1 — a2(l — cos kAx) — žen sin kAx. (H4) a GG* výrazem \GG*\ = 1 - a2(l - a2)(l - cos2(kAx)). (115) Schéma Lax-Wendroff bude stabilní pokud bude a2 < 1. Vlnová rovnice jako systém lineárních PDR prvního řádu Vrátíme se zpátky k ID vlnové rovnici pro veličinu u = u(t, x) ve tvaru d2u 2^2u dt2 dx2 ťlT3 ("6) 31 a převedeme ji na soustavu lineárních PDR prvního řádu. Rovnici (116) převedeme na tvar (d d \ ( d d\ n nebo lépe, na tvar d d \ (du du\ n Z posledního výrazu je vidět, že bude výhodné zavést substituci du du /-i -i r\\ r=m a s = cYx- (119) Když tuto substituci dosadíme do původní vlnové rovnice (116) tak obdržíme rovnici dr ds M=Cä-x (120) Rozvineme nyní rovnici (118). Dostaneme se tak k rovnici ds dr ds dr (121) dt dx dx dt Z rovnice (120) však dále vyplývá platnost vztahu ds dr ds dr (122) dt dx dt dx Soustava lineárních PDR prvního řádu, představujících vlnovou rovnici (116), je du ä* = r- (123) 32 Tuto soustavu nyní můžeme řešit jedním z diskutovaných schémat. Naší volbou v tomto textu bude schéma Lax-Wendroff. Nejprve převedeme poslední dvě rovnice a obdržíme r n+l _ _n r -i— 4-i + «(r?+i-2r? + r?-i) (126) a 4 + r 3 + i_rn+a{sn_2snJrSn_i) První rovnici ze soustavy pak snadno určíme ze vztahu u n+l _ n ~3 3 Parametr a je definován vztahem a = cAt ~Ax (127) (128) (129) 33 2 Eigen-value problém 2.1 Schrôdingerova rovnice Schrô dingerova rovnice je pohybovou rovnicí pro nerelativistický elektron. Budeme se zabývat jak časovou rovnicí a hledat řešení vyvíjející se v čase tak i nečasovou Schrôdingerovou rovnicí popisující stacionární stavu vázaných elektronů v potenciálu V. Jednorozměrná rovnice má tvar ih^EA = Hi,{x,y). (130) s Hamiltonovým operátorem hH, který nezávisí na čase. Řešení pak hledáme v separovaném tvaru ý(x,t) = il){x)(j>{t). (131) Toto řešení dosadíme do rovnice (130) a obdržíme rovnici Pravá i levá strana poslední rovnice se musejí rovnat jedné a téže konstantě rozměru energie, řekněme E. Tímto způsobem jsme obdrželi dvě diferenciální rovnice Hif;{x) = Eif;{x) (133) a ih^- = Eftt). (134) Rovnice (133) je hledaná bezčasová Schrôdingerova rovnice. 34 Bezčasová Schrodingerova rovnice Nyní se zabývejme hledáním vázaných stavů nerelativistické částice v potenciálu V(x), které jsou řešením bezčasové Schródingerovy rovnice Ťi2 d2 \ + V{x) )iIj = EiIj, (135) 2m dx2 ^ * 2 j 2 kde zřejmě je H = — 2mďř2 + V(x) a Ľ je energie částice. Částice v nekonečné potenciálové jámě Nekonečnou potenciálovou jámu charakterizuje potenciál ve tvaru T// x í oo pro x < -1 a x > 1 , , F(X) = { 0 pro s e [-1,1]. (136) Částice nemůže existovat v oblasti x < — 1 a, x > 1, proto je zde ijj{x) = 0. (137) Vlnová funkce je spojitá, v bodech x = —1 a x = 1 musí také platit ^(-1) = ^(+1) = 0. (138) Uvnitř potenciálové jámy má Schrodingerova rovnice tvar d2é J = -* ^ a její obecné řešení je t/j(x) = acos(kx) + bsin(kx). (140) 35 Z okrajových podmínek dostaneme = acos(—k) + 6sin(—k) = acos(k) — fesin(fc) = 0 (141) a ^(1) = acos(fc) + bsm(k) = 0. (142) Když sečteme poslední dvě rovnice tak dostaneme rovnici Tí 7T 2a cos k = 0 k = —, pro n = 1,3,.... (143) Pokud odečteme poslední dvě rovnice tak dostáváme rovnici Tí 7T 26sin/c = 0 k = —, pro n = 2,4,.... (144) Existují dvě větve řešení: • (1) - z požadavku ^,'(0) = 0^6 = 0 (145) obdržíme řešení i\){x) = acoskx (146) kde k splňuje rovnici (143). • (s) - z požadavku ^(0) = 0 => a = 0 (147) obdržíme řešení ifj(x) = bsiukx (148) kde k splňuje rovnici (144). 36 Po dosazení řešení (147) nebo (148) do rovnice (139) zjistíme, že energie částice bude kvantována a bude platit e=(— , pro n = 1,2,3,.... (149) Při hledání numerického řešení využijeme symetrii úlohy a budeme integrovat rovnici (139) od x = 0 směrem k x = 1 pro danou energii e. Hledáme takové e pro které bude = 0. Nejprve provedeme experiment, kdy numericky určíme, jak se chová i/j(x = l;e) pro energie e (obrázek 11 ilustruje, vývoj i/j(x = — 1; e)). Je vidět, že kvantované hodnoty energie částice v potenciálové jámě leží na intervalech [ei,e2J, kde platí ^(l, €i)i/j(1, 62) < 0. 0.5 0 -0.5 0 10 20 30 40 50 60 70 80 90100 Figure 11: Ilustrace vývoje funkcí if;(8\l,e) a if;(l\l,e). K určení jedné kvantované hodnoty en budeme postupovat podle následujícího algoritmu: 37 (150) • Pro takto nalezený interval použijeme jednu z metod hledání kořene na intervalu monotónnosti(např. půlení intervalu). Lineární harmonický oscilátor Potenciál lineárního harmonického oscilátoru je dán rovnicí V(x) = ÍcjV. (151) Rovnice (135) přejde na tvar a po převedení do bezrozměrného tvaru, máme rovnici de (A-a^ = 0. (153) Nejprve se analyzujme asymtotické chování řešení rovnice (153), tj. v oblastech £ —> ±oo. Rovnice (153) přejde na tvar " ^ = °- (154) 38 Snadno se přesvědčíme, že tuto rovnici řeší funkce ijj = aexp(£2/2) + 6exp(-£2/2). (155) První sčítanec v této funkci zřejmě způsobí, že i\) pro £ —> ±00 bude divergovat, což je v rozporu s požadavkem aby vlnová funkce byla kvadraticky integrabilní (musí platit —> ±00) -)• 0). V tomto případě bude mít vlnová funkce tvar ^(0 = &exp(-£2/2). (156) Obecné řešení rovnice (153) budeme tedy hledat ve tvaru ^(0 = ^(Oexp(-f2/2). (157) Dosazením tohoto řešení do rovnice (153) obdržíme diferenciální rovnici pro nenámou funkci v(£) ve tvaru v"(0 - 2£i/(0 + (A - 1M0 = 0. (158) Řešení budeme hledat pomocí mocninné řady 00 v(0 = J2a^k- (159) k=0 První a druhá derivace jsou 00 At) = E akK"-1 (160) fc=l a 00 ""(í) = E»t^-l)íM. (161) k=2 39 Po dosazení do (158) a přeindexaci dostaneme pro koeficienty ak rekurentní vztah 2fc-(A-l) °*+* = {k + 2)(k + l)ak- (162) Existují dvě větve řešení: • sudá - ao = 1 a a\ = 0, • lichá - ao = 0 a a\ = 1. Pro velká platí pro ak+2/ak přibližně afc+2 2/c 2 (163) afc (fc + 2)(fc+l) fc + 2" Ovšem pro funkci exp(£2), kterou rozvedeme do mocninné řady dostaneme oo exp(£2) = £Wfe> (164) k=0 kde platí ^ = (165) bk k + 2 v 7 Odtud vidíme, že pokud funkci i>(£) rorvineme do mocninné řady pak se asymptoticky bude chovat jako funkce exp£2 a způsobí divergenci i\) v ±oo. Funkce v(£) musí mít formu polynomu, tj. existuje takové n, od kterého je an+2 = 0. Takto se dostavme ke kvantování dosud libovolné hodnoty veličity A. Aby byla funkce ip kvadraticky konvergentní musí A splňovat podmínku 2n - (Xn - 1) = 0 \n = 2n + 1. (166) Při numerickém řešení budeme, díky symetrii potenciálu V, postupovat analogicky jako v případě nekonečné potenciálové jámy. Rozdíl 40 oo. je v hodnotě £ kdy vlnová funkce vymizí, zde je to pro £ —> Nekonečno zde nahradíme konečným ale "velkým" číslem 100. Když sudé a liché řešení vlnové rovnice pro lineární harmonický oscilátor v bodě £ = 100 vyneseme do grafu pro různé hodnoty parametru A G [0,10], obdržíme obrázek 12, ze kterého je ihned vidět, že nulové body i\) nastávají pro diskrétní hodnoty parametru A = 1, 3, 5, ... . Což přesně odpovídá výsledku (166). 0.5 0 -0.5 i—i—i—r j_i_i_L _L J_L 01 23456789 10 Figuře 12: Hodnoty vlnových funkcí v bodě £ = 100 jako funkce parametru A. 41 2.2 Vlastní frekvence schwarzchildovy černé díry Zabývejme se určením vlastních modů Schwarzchildovy černé díry. Konkrétně, budeme hledat řešení Regge-Wheelerovy rovnice dx2 + a2 - V(x) i\) = 0 (167) pro —oo < x < oo. Potenciál V(x) je na celém definičním oboru kladná a konečná, tj. /oo V(a;)da; = A (168) -oo přičemž A je konečné. V případě Schwarzchildovy geometrie má potenciál V tvar V{r) = 2(1- 2/r) ^ + ^ ~ 3 (169) kde je r = r(a^) a je dáno rovnicí z = r + 21og|r/2-l|. (170) Pro numerický výpočet řešení rovnice (167) je výhodné zavést novou funkci ó vztahem ip = exp (i J (jxíx^j . (171) Dosazením (171) do (167) obdržíme odpovídající rovnici pro funkci 4> ve tvaru i^+a2-2-V = 0. (172) dx Řešení, představující kvazinormální mody černé díry odpovídají okrajovým podmínkám 4> —> —a pro x —> oo (173) 42 a 4> —> + — oo. (174) Z praktických důvodů rozdělíme rovnici (172) na reálnou a imaginární část. Nejprve rozdělíme frekvenci g a řešení 0, tj. 0 = 01 + Ž02 a '2exp (-g) d,'. (186) 48 Figure 14: Náhodná procházka v M-B plynu. Vratme se k problému náhodné procházky. Mezi dvěmi srážkami se pohybuje částice plynu rovnoměrně přímočaře rychlostí v. Střední volná dráha pro částici v plynu s MB rozdělením rychlostí je i = (\/W) (187) kde <7 je plocha jedné částice (účinný průřez), n je hustota počtu částic v plynu. Simulace náhodné procházky v MB plynu je dána následujícím algoritmem: • generátor náhodných čísel určí hodnoty £v a £a, 49 • řešením rovnice (185) určíme v , rychlost odpovídající £v podle MB rozdělení, a jednotlivé složky rychlosti vx = v cos(27r£a) a vx = vsm(27ľ^a). • novou polohu určíme ze vztahu ř=vAt + r0, (188) kde je f = (x,y), v = (vx,vy) a časový interval po který se částice pohybuje mezi srážkami je At =< v > /£ (189) přičemž střední kvadratická rychlost < v > MB rozdělení určuje rovnice < v >= (190) 4 Vícerozměrné integrály 4.1 Opakování jednorozměrné integrace V této kapitole se omezíme na vícerozměrné integrály které budou maximálně 3D. Uvažujme následující určitý integrál funkce f(x,y,z) fX2 fV2{x) rz2(x,y) 1= dx dy dzf(x,y,z). (191) J xx Jyi{x) J zx{x,y) Tento integrál řešíme tak, že postupně řešíme ID integrály přes x, y and z. Tj. rozdělíme integraci přes 3D oblast na integraci tří ID integrálů, které mají následující tvar G(x,y) = / dzf(x,y,z), (192) Jzx (x,y) 50 n/2 O) H (x) = / G(x,y)dy, (193) Jyi(x) I=j H(x)dx. (194) Následující pseudokód ukazuje realizaci integrace přes 3D oblast: I = 0 FOR x = xl TO x2 H = 0 FOR y = yl TO y2 G = 0 FOR z = zl TO z2 G = G + f (x,y , z)*dz END FOR H = H + G*dy END FOR 1 = 1+ H*dx END FOR Přitom každá ze smyček může realizovat soŕiiovanoulD metodu integrace (Simpsonova pravidlo, Rombergova integrace, metoda Gaussovy kvadratury). Příklad: Necht je dána funkce f(Xly) = x2y + b. (195) Máme určit objem pod touto plochou nad kruhem o poloměru r = 2 v počátku souřadnic, tj. řešíme dvojný integrál rx2 ry2{x) I dx dyx y. (196) J xx Jyi{x) Meze zřejmě budou x\ = —r, x = -££/(*<)> (/2) = 7?E/2fe)- (198) 1 N N i=i i=i Uvažujme funkci jedné proměnné f (x) na intervalu x G [a, b] a funkční hodnoty sspadají do intervalu y G [2/1,2/2]- Určitý integrál funkce / je obsah plochy pod /, tj. ŕ A= f(x)dx. (199) J a Dále, necht platí f (x) > 0. Celková plocha obdélníku Atot = (b — o) * f max a označme Ntot celkový počet náhodných bodů rovnměrně 52 rozdělených na ploše odelníku a Na je počet náhodných bodů, které padnou pod křivku /. Za těchto předpokladu platí ľ f(x)dx~Atot^. (200) _A_ „ Na . ľb „, N , . , NA Atot Ntot 53