UČEBNÍ TEXTY VYSOKÝCH ŠKOL Vysoké učení technické v Brně Fakulta strojního inženýrství Doc. RNDr. Libor Čermák, CSc. Numerické metody pro řešení diferenciálních rovnic Ústav matematiky FSI VUT v Brně Obsah 1 Obyčejné diferenciální rovnice: počáteční úlohy 5 1.1 Formulace, základní pojmy .............................. 5 1.2 Eulerovy metody.................................... 7 1.3 Explicitní Rungovy-Kuttovy metody......................... 12 1.4 Lineární mnohokrokové metody............................ 19 1.4.1 Obecná lineární mnohokroková metoda.................... 19 1.4.2 Adamsovy metody................................ 20 1.4.3 Metody zpětného derivování.......................... 24 1.5 Tuhé problémy..................................... 26 2 Obyčejné diferenciální rovnice: okrajové úlohy 33 2.1 Metoda střelby..................................... 33 2.2 Diferenční metoda................................... 35 2.3 Metoda konečných objemů............................... 40 2.4 Metoda konečných prvků ............................... 41 3 Parciální diferenciální rovnice 48 3.1 Úloha eliptického typu................................. 49 3.1.1 Formulace úlohy................................. 49 3.1.2 Diferenční metoda................................ 50 3.1.3 Metoda konečných objemů........................... 56 3.1.4 Metoda konečných prvků............................ 60 3.2 Úloha parabolického typu............ ................... 67 3.3 Úloha hyperbolického typu .............................. 71 3.4 Hyperbolická rovnice prvního řádu.......................... 73 Literatura 80 3 Předmluva Tato skripta jsou určena pro studium předmětu Numerické metody II. Výklad navazuje na úvodní kurz Numerické metody a proto se předpokládá, že čtenář má základní znalosti o numerických metodách lineární algebry, řešení nelineárních rovnic, interpolaci, numerickém derivování a integrování. Předkládaný text lze použít také jako studijní literaturu pro doktorandské studium na FSI VUT. Skripta uvádějí celou řadu algoritmů a doprovodný text [10] k tomu přidává příklady a cvičení. Pro implementaci algoritmů a experimentování s nimi se výtečně hodí prostředí MATLABu. První kapitola se věnuje numerickému řešení počátečních úloh pro obyčejné diferenciální rovnice. Dílčí témata jsou tradiční, tj. Rungovy-Kuttovy metody, Adamsovy metody a metody zpětného derivování. Ve druhé kapitole jsou uvedeny tři základní metody řešení okrajových úloh pro obyčejné diferenciální rovnice druhého řádu, a sice diferenční metoda, metoda konečných objemů a metoda konečných prvků. Třetí kapitola je věnována numerickým metodám řešení parciálních diferenciálních rovnic. Pro eliptickou parciální diferenciální rovnici ve dvou prostorových proměnných je uvedena diskretizace diferenční metodou, metodou konečných objemů a metodou konečných prvků. Řešení parabolické a hyperbolické parciální diferenciální rovnice druhého řádu v jedné prostorové proměnné se provádí metodou přímek. Pro hyperbolickou rovnici prvního řádu v jedné prostorové proměnné je použita metoda charakteristik. V rámci klasických tematických okruhů jsem se snažil do skript zařadit takové numerické metody, které se v současnosti skutečně používají. Vycházel jsem přitom z osvědčených učebnic numerické matematiky, jakými jsou např. knihy [9], [17], [19], [26], a z vynikajících monografií, mezi nimi zejména [12], [22], [23], [14], [8], [2], [28]. Pokud jde o české zdroje, nejvíce podnětů jsem čerpal z knih [26], [27] a ze skript [16], [3] a [18]. Za chyby a přepisy, které se bohužel ve skriptech jistě vyskytnou, se dopředu omlouvám. Budu vděčný všem čtenářům, kteří mně na ně upozorní. Brno, květen 2013 Libor Čermák 4 1. Obyčejné diferenciální rovnice: počáteční úlohy V této kapitole se budeme zabývat problematikou numerického řešení počátečních úloh pro obyčejné diferenciální rovnice. 1.1. Formulace, základní pojmy Počáteční problém pro ODR1 spočívá v určení funkce y(t), která vyhovuje diferenciální rovnici y'(t) = f(t,y(t)) (1.1) a splňuje počáteční podmínku y(a) = V. (1.2) Je-li v nějakém okolí D bodu [a,rj\ funkce f(t,y) spojitá a splňuje-li v tomto okolí Líp-schitzovu podmínku s konstantou L vzhledem k proměnné y, tj. platí-li \f(t,u)-f(t,v)\ 1, neboli když r A leží v oblasti absolutní stability Ra- r\ e RA = {z E C I \z - 1| > 1}. Oblast absolutní stability IE metody je tedy obrovská, je to celý vnějšek \z — 1| > 1 jednotkového kruhu komplexní roviny C se středem v bodě [1, 0]. Interval absolutní stability IE metody Ia = (—oo,0). Podmínka stability (1.11) délku kroku IE metody zřejmě nijak neomezuje. Je to právě mimořádná stabilita, která je onou hledanou předností IE metody ve srovnání s EE metodou. Tento klad je však třeba vykoupit nutností řešit obecně nelineární rovnici. yn+± získáme jako přibližné řešení rovnice g (z) =0, kde g (z) = z-yn- rf(tn+1,z). Protože dobrou počáteční aproximaci lze získat extrapolací z hodnot yn,yn-i,..., dá se očekávat rychlá konvergence Newtonovy metody (pro řešení nelineárních rovnic). Praktická zkušenost potvrzuje, že tomu tak skutečně je. Lichoběžníkovou metodu dostaneme jako aritmetický průměr EE metody a IE metody: Vn+i =Vn + \r[f{tn,yn) + f(tn+1,yn+1)] . (1.14) Pro lokální diskretizační chybu lten = y{tn+1) - y{tn) - \r[f{tn,y{tn)) + f(tn+1,y(tn+1))] užitím Taylorovy věty odvodíme lten = -^rV"(ín)+0(r4). Lichoběžníková metoda (stručně TR metoda podle anglického trapezoídal rule) je tedy implicitní metoda řádu 2. Stabilitu TR metody zjistíme řešením testovací úlohy (1.10): 2 + tA 2 -t\ n+l Vn+i =yn + \T\[yn + yn+1], odtud yn+1 = \ + T\n = ■■■■■ Z — TA Není těžké ověřit, že podmínka stability (1.11) platí, právě když tA G Ra = {z e C I Re(z) < 0}. Oblast absolutní stability TR metody tedy obsahuje celou zápornou polorovinu komplexní rovniny C, interval absolutní stability Ia = (—oo, 0). TR metoda (s podporou IE metody) je v Matlabu implementována jako program ode23t. 11 1.3. Explicitní Rungovy-Kuttovy metody Obecný tvar s-stupňové explicitní Rungovy-Kuttovy metody je Vn+i =Vn + r(6ifci + b2k2 H-----h bsks), (1.15) kde koeficienty k,n i = 1,2,..., s, jsou určeny předpisem h = f(tn,yn), k2 = f(tn + rc2, yn + ra21k1), h = f(tn + rc3,yn + r(a31k1 + a32k2)), (1.16) fcs = f(tn + rcs,í/n + r(asiA;i + as2k2 H-----h aS)S_i/cs_i)), a kde 6«, q, a^- jsou konstanty definující konkrétní metodu. Rungova-Kuttova metoda (1.15), (1.16) je explicitní: nejdříve spočteme k\, pak k2 pomocí k\, pak k3 pomocí k\, k2 a tak dále, až nakonec spočteme ks pomocí k\,k2,..., ks_\. Vypočtené koeficienty kÍ7 i = 1, 2,..., s, dosadíme do (1.15) a dostaneme yn+±. V dalším budeme hovořit jen o Rungových-Kuttových metodách (stručně RK metodách), tj. slůvko „explicitní" vynecháme. Je však třeba připomenout, že existují také implicitní Rungovy-Kuttovy metody, těmi se však zabývat nebudeme. RK metody jsou zřejmě jednokrokové: k výpočtu yn+± potřebujeme znát jen yn, předchozí hodnoty yn-i, yn-2, • • • v kroku od tn do tn+1 nepoužijeme. Koeficient ki je směrnicí lokálního řešení procházejícího bodem [ť*,y*], kde [ťi,yi\ = [tn^n], t*=tn+TCi, y* = yn+T^ankx + a^k^-----ha^-ifci-i), i = 2,3,...,s. Do bodu [tn+i,yn+i] se tedy dostaneme z bodu [tn,yn] tak, že se posuneme po přímce, jejíž směrnice k* = b\k\ + b2k2 + • • • + bsks je lineární kombinací se směrnic k\,k2,... ,ks, pro metodu řádu alespoň 1 jde o vážený průměr, neboť J^*=1 = 1, viz (1.17). Abychom dostali konkrétní metodu, musíme určit stupeň s a dále konstanty q, 6«, a^. Konstanty RK metod je zvykem zapisovat do tabulky známé jako Butcherova tabulka: c2 0.21 c3 a3i cs asi as2 ... as,s-l bi b2 ... bs-i bs Jedním z kritérií při volbě konstant RK metody je dosažení dostatečné přesnosti. Tu měříme pomocí lokální diskretizační chyby s lten = y(tn+i) - y(tn) - ^ biki(tn), 12 kde h(tn) = f(tn,y(tn)), í-1 ki(tn) — f(tn + tq, yn-\- t *y ^ aijkj(tn)), i — 2, 3,..., s. i=i Lokální diskretizační chyba je chyba, které se dopustíme v jednom kroku za lokalizačního předpokladu yn = y(tn). RK metoda je řádu p, pokud lokální diskretizační chyba je řádu 0(tp+1). Pro p = 1, 2, 3 lze odvodit následující tzv. podmínky řádu: řád 1: 5> = 1> řád 2: s 5> = 1> 1=1 s / j u^ 2 ' i=2 řád 3: s 5> = 1> 1=1 s s s i—l = 2 , &jCj = - , bidijCj = g i=2 i=2 i=2 i=2 (1.17) Odvození podmínek řádu pro p = 1, 2, 3,4, 5 lze najít třeba v [22]. Protože všechny prakticky používané metody splňují podmínku Ci = ail + ai2-\-----haí,í-i, i = l,2,...,s, (1.18) budeme i my předpokládat, že podmínka (1.18) platí. RK metoda řádu p > 1 má globální diskretizační chybu řádu 0(rp). Předpokladem pro platnost tohoto tvrzení je dostatečná hladkost pravé strany /, konkrétně je třeba, aby funkce f(t, y) měla spojité derivace až do řádu p včetně. Pokud / má spojité derivace jen do řádu s < p, pak lze pro globální chybu dokázat pouze řád 0(rs), viz [22]. Označme p(s) maximální dosažitelný řád s-stupňové RK metody. Platí p(s) = s pro s = 1, 2, 3,4, p(8) = 6, p(5) = 4, p(9) = 7 p(6) = 5, p(s) < s — 2 pro s = 10,11,... p(7) = 6, Vidíme, že s-stupňové RK metody řádu s existují jen pro 1 < s < 4. Například metoda řádu 5 je nejméně 6-ti stupňová. Uveďme si několik nejznámějších metod. Metoda řádu 1. Pro s = p = 1 existuje jediná explicitní metoda a tou je nám již známá EE metoda yi+1 = y,t + rf(th Metody řádu 2. Pro s = p = 2 má explicitní RK metoda Butcherovu tabulku Podmínky (1.17) pro metodu řádu 2 stanoví c2 021 h b2 61 + 62 = 1, hC2 = \ 2 ' 13 a protože ve shodě s (1.18) předpokládáme a21 = c2? dostáváme tabulku kde ab = |. Parametry a, b jsou tedy svázány jednou podmínkou, takže zvolíme-li a ^ 0, je b = l/(2a). a a 1-6 6 Pro a = |je6 = la dostáváme metodu Vn+i =Un + rk2, kde k2 = f(tn + \r, yn + \rk{), kx = f(tn, yn), (1.19) známou pod názvem modifikovaná Eulerova metoda. Budeme ji značit EM1 jako první moáífikace Eulerovy metody. V anglicky psané literatuře je metoda (1.19) známa jako míápoínt Euler formula. Pro a = 1 je 6 | a dostáváme metodu yn+i=yn + \r{k1 + k2), kde h = f(tn,yn), k2 = f(tn + r,yn + rkt). (1.20) Budeme ji značit EM2 jako áruhou moáífikací Eulerovy metoáy. Metoda (1.20) se také často uvádí pod názvem Heunova metoáa. Pro a = |je6 = |a dostáváme metodu Vn+i =Un + \r{kx + 3k2), kde kx = f(tn,yn), k2 = f(tn + \r,yn + fr/ci), známou jako Ralstonova metoáa řádu 2 (stručně R2 metoda). Metody řádu 3. Pro s = p = 3 dostáváme Butcherovu tabulku c2 a 4 podmínky pro metodu řádu 3: 6i + 62 + 63 = 1, b2c2 + 63C3 c2 c3 C3 — «32 ^32 6l 1 2 ' b2c\ + 63ŕ| = i , 63a32c2 1 6 • Když zvolíme dva parametry 0 < c2 < C3, jsou tím všechny koeficienty metody jednoznačně určeny. Volba c2 = |, c3 = | vede na Ralstonovu metoáu řádu 3 (stručně R3 metodu): i Ž/n+i =Vn + \r (2fci + 3A;2 + 4A;3), 0 f k1 = f(tn,yn), k2 = f(tn + |r,í/n + ^t/ci), 3, 2 14 9 3 9 ^3 = /(*n + |r,í/n + frfc2). Ralstonova metoda řádu 3 je základem Runge-Kutta-Bogackí-Shampíne metody, viz [22], která je implementována do Matlabu jako funkce ode23 a jejíž popis uvedeme v této kapitole později. Metody řádu 4. Pro s = p = 4 je nejznámější klasická Rungova-Kuttova metoáa Vn+i =yn + \r (fci + 2fc2 + 2A;3 + fc4), T,yn + ±Tk2), kA = f(t 1 1 2 2 1 0 1 to 2 1 0 0 1 1 1 1 1 6 3 3 6 h = f(tn,yn) h = f(tn + \ 5r> Ž/n V \rkx) Tk3). 14 Klasická Rungova-Kuttova metoda (stručně cRK4) byla velmi populární v době, kdy se ještě nepoužívaly samočinné počítače a kdy proto velmi významným kritériem byla jednoduchost metody. Toto hledisko však v současné době ztratilo na významu a proto se používají jiné metody. Kvalitní dvojice metod řádu 4 a 5 jsou součástí metod Runge-Kutta-Fehlberg nebo Runge-Kutta-Dormand-Prínce, viz např. [12]. Posledně jmenovaná dvojice metod je použita v matlabovské funkci ode45, popis uvedeme v této kapitole později. Řízení délky kroku. V profesionálních programech uživatel zadá toleranci e a program délku kroku vybírá tak, aby velikost odhadu estn lokální chyby len nabývala pořád přibližně stejné hodnoty e. Krok od yn do yn+i je úspěšný, když |estn| rmax, položí se r = rmax, kde je rmax = 0,1(6—a). 3) V případě neúspěšného kroku navrženou délku r* redukujeme: r* = max(r*,gminr) , kde qmin = 0,5 resp. 0,1 v ode23 resp. ode45 při prvním neúspěchu v rámci jednoho kroku a r* = 0,5r při opakovaném neúspěchu v temže kroku. 15 4) V případě úspěšného kroku délku kroku r* redukujeme předpisem r* = min(r*,qmaxr), kde qmax = 5. 5) V kroku bezprostředně následujícím po neúspěšném kroku se délka kroku nesmí zvětšit. 6) Počáteční délka kroku r = 0ey^ [max(\V\,ea/er)] /\f(a,V)\, (1.23) přičemž pro r < rmin změníme r na rmin a pro r > rmax změníme r na rmax. 7) Při programování je třeba postupovat opatrně. Příkazy programu je nutné uspořádat tak, aby nemohlo dojít k dělení nulou, viz (1.22) a (1.23). Podrobnější informace týkající se řízení délky kroku lze najít v [22]. Odhad lokální chyby. Základní myšlenka je jednoduchá. Použijí se dvě metody, z nichž jedna je řádu p a druhá řádu p+1. Z výchozí hodnoty yn spočteme y™+1 přesnější metodou a méně přesnou metodou. Použitelný odhad lokální chyby je estn — yn+1 — yn+i- Obě metody používají tutéž množinu koeficientů {kj}j=l. V tom případě je totiž získání odhadu „laciné". Dvojice Rungových-Kuttových metod se popisují pomocí rozšířené But-cherovy tabulky, přičemž Vn+i = Vn + r ]CLi tfkj je metoda řádu p, c2 0.21 Cs asi bl h* ■ K bl* h** u2 ■ K* Ex E2 ■ Es y*n+i = yn + r J2Sj=i bTkj Je metoda řádu p + ^ estn = r ^2Sj=1 Ejkj je odhad lokální chyby, takže Ej = b** - 6*, j = 1, 2,..., s. Pokud ve výpočtu pokračujeme přesnější metodou, tj. když yn+± = y^+í = y*n+\ +estn, říkáme, že jsme použili dvojici metod s lokální extrapolací. Tento postup se v současných programech upřednostňuje. Druhou možností je pokračovat méně přesnou metodou, tj. položit yn+i = Vn+i- V tom případě se přesnější metoda použije jen pro získání odhadu chyby a říkáme, že jsme dvojici metod použili bez lokální extrapolace. Obě níže uvedené metody BS32 a DP54 se používají jako metody s lokální extrapolací. Příkladem metody, která se obvykle používá bez lokální extrapolace, je Runge-Kutta-Fehlbergova metoda řádu 4, označovaná stručně jako RKF45, viz např. [12]. Na vysvětlenou k použitým zkratkám uveďme, že první číslo značí řád metody a druhé číslo řád pomocné metody použité pro odhad lokální chyby. 16 Bogacki Shampine metoda, stručně BS32 metoda, je implementována v Matlabu jako funkce ode23. Rozšířená Butcherova tabulka BS32 metody je Přesnější z obou metod páru je Ralstonova R3 metoda 1 1 2 2 3 0 3 4 4 1 2 1 4 X 9 3 9 7 1 1 1 24 4 3 8 2 1 4 n 9 3 9 u 5 1 1 1 72 12 9 8 y n+l Ž/n [§*1 řádu 3. Pomocná metoda y n+l Vn 5*2 J*2 fc3] ^3 Ž/n+1 řádu 2 používá kromě koeficientů ki, k2 a /c3 navíc ještě koeficient £4 = f(xn+i,yn+i). V každém kroku se tedy počítají jen 3 nové hodnoty funkce /, neboť koeficient k\ ve stávajícím kroku je roven koeficientu /c4 z kroku předchozího, takže nově se počítají jen koeficienty k2, k3 a /c4. Výjimkou je případ neúspěšného kroku, kdy se hodnota yn+± neakceptuje a krok se krátí. Tyto případy však nebývají časté. Zařazení koeficientu k± do metody řádu 2 nás tedy téměř nic nestojí, umožní však zlepšit vlastnosti této metody a tím i celého páru. Metoda, ve které ks = f(tn+i,yn+i), bývá označována jako FSAL podle anglického First Same As Last. Zdůrazněme, že BS32 metoda se používá jako metoda s lokální extrapolací, tj. y n+i = y*n+i- Hodnoty přibližného řešení pro t G {tn,tn+i) spočteme dostatečně přesně pomocí kubického Hermitova polynomu Hs(t) určeného podmínkami = yn+i, H'^tn+1, ki kň Metoda BS32 je tedy skvělá: je řádu 3, v každém úspěšném kroku se pravá strana / počítá jen 3-krát, a to stačí jak na řízení délky kroku tak na výpočet řešení mezi uzly tn a tn+i. Dormand-Prince metoda, stručně DP54 metoda, je definována rozšířenou Butchero-vou tabulkou 1 5 1 5 3 10 3 40 9 40 4 5 44 45 56 15 32 9 8 9 19 372 6 561 25 360 2187 64 448 6 561 212 729 1 9017 3168 355 33 46 732 5 247 49 176 5103 18 656 1 35 384 0 500 1113 125 192 2187 6 784 11 84 5179 57 600 0 7 571 16 695 393 640 92 097 339 200 187 2100 1 40 35 384 0 500 1113 125 192 2187 6 784 11 84 0 71 0 71 71 17 253 22 1 57 600 16 695 1920 339 200 525 40 Metoda DP54 je typu FSAL, neboť k7 = f(tn+1,yn+1). Proto se v každém úspěšném kroku 17 metody počítají jen koeficienty k2,..., k7, koeficient k\ byl už vypočítán jako koeficient k7 v předchozím kroku. Zdůrazněme, že DP54 metoda se používá jako metoda s lokální extrapolací, tj. yn+1 = y**+1. Hodnoty přibližného řešení pro t G {tn,tn+i) spočteme dostatečně přesně pomocí in-terpolačního polynomu H^s) = yn + rkBq, kde k = (k1} k2, k3, /c4, k5, k$, k7), B = " 1 183 64 37 12 145 -i 128 0 0 0 0 0 1500 371 1000 159 1000 371 0 125 32 125 12 375 64 0 9 477 3 392 729 106 25 515 6 784 0 11 7 11 3 55 28 . 0 3 2 -4 5 2 J Q = (Q, Q2, 1 a e 7^ 0 je limn^oo \yn\ —► oo, což je nepřijatelné: pro e = 0 je yn = 0 přesné řešení rovnice y' = 0, avšak pro s ^ 0, \e\ 1, tj. již pro nepatrnou poruchu startovacích hodnot yj, j = 0,1,..., k — 1, dostaneme řešení zcela nevyhovující. Dá se ukázat, že vyloučit musíme také případ, kdy \r\ = 1 je kořen násobnosti větší než 1. Konvergence. Uvažujme D-stabilní LMM řádu p > 1. Jestliže startovací hodnoty zadáme s chybou řádu 0(rp), pak globální diskretizační chyba je rovněž řádu 0(rp). Precizní formulaci a důkaz této věty lze najít např. v [27], [12]. Předpokladem její platnosti je dostatečná hladkost pravé strany /, konkrétně je třeba, aby funkce f(t,y) měla spojité derivace až do řádu p včetně. Pokud / má spojité derivace jen do řádu s < p, pak lze pro globální chybu dokázat pouze řád 0(rs). Absolutní stabilita. Řešíme-li testovací úlohu (1.10) LMM na rovnoměrném dělení s krokem r, dostaneme k ^(a3--rAft)yn+w = 0. (1.25) j=0 Řešení hledejme ve tvaru yn = rn. Po dosazení do diferenční rovnice (1.25) obdržíme k k ^{a3 - T\Pj)rn+1-3 = rn+1-k ^(a3- - TX/3j)rk-j = rn+1-k7i(r,rX) = 0 , j=0 j=0 k kde Tr(Z,z) = Y,(aj-z/3j)Zk-1 j=0 je tzv. polynom stability LMM. Jestliže 7r(r, r\) = 0, pak yn = rn je řešením diferenční rovnice (1.25) a podmínka stability rn —► 0 pro n —► oo platí, právě když |r| < 1. Oblast R a absolutní stability LMM metody definujeme jako množinu bodů z komplexní roviny, pro které z podmínky 7r(£, z) = 0 plyne |£| < 1. Podmínka absolutní stability tedy platí, když tA g Ra = {z e C I tt(£, z) = 0 |^| < 1}. 1.4.2. Adamsovy metody Integrací diferenciální rovnice (1.1) od tn do tn+1 dostaneme Ptn+l y(tn+1) - y(tn) = / f(t,y(t))dt. 20 Funkci f (t, y (t)) aproximujeme pomocí interpolačního polynomu Pk-i(t) stupně k—l, tj. ŕtn+1 y(tn+l) ~ y(tn) +/ Pk-l(t)dt, kde Pk-l(tn+l-j) = f(tn+l-j,y(tn+l-j))- J tn Když přibližnou rovnost nahradíme rovností a přesné řešení nahradíme řešením přibližným, dostaneme předpis Ptn+l yn+i = yn+ Pk-i(t)dt, kde Pk-i{tn+i-j) = f(tn+i-j,yn+i-j)- (1.26) Adams-Bashforthovy metody dostaneme, když v (1.26) zvolíme j = 1, 2,..., k. Konstrukci polynomu Pk-i(t) lze přehledně vyjádřit tabulkou t tn—1 tn+l—k Pk-l{t) f n fn-1 fn+l-k Adams-Bashforthovu metodu lze zapsat ve tvaru k y-a+1 =y-a+r^2 Pk,j fn+l-j • Stručně ji budeme označovat jako ABk metodu. ABk metoda je explicitní, fc-kroková, řádu k, D-stabilní. Pro konstantní délku kroku, tj. když tn+i-j = tn+i — 3Ti j = 1) 2, • • •, k, jsou koeficienty ABk metod pro k = 1,2,..., 6, spolu s chybovými konstantami Ck+1 a dolními mezemi a*k intervalů absolutní stability («£,0), uvedeny v následující tabulce: k @k,2 #,3 01,4 01,5 fikfi Uk+1 a*k 1 1 1 2 -2 2 3 2 1 2 5 12 -1 3 23 12 16 12 5 12 3 8 8 11 4 55 24 59 24 37 24 9 24 251 720 3 10 5 1901 720 2 774 720 2 616 720 1274 720 251 720 95 288 90 551 6 4 277 1440 7923 1440 9 982 1440 7 298 1440 2 877 1440 475 1440 19 087 60 480 5 57 Všimněte si, že AB1 metoda je totožná s EE metodou, tj. AB1=EE. Koeficienty ABk metod lze určit z formule fc-i yn+i =yn + r^2 rf^Jfn, j=0 21 kde VJ je operátor zpětné diference, V fn = fn, V/n = fn ~ fn—li ■ ■ ■ ^ fn = V /n — V fn—li i = 2,3, . . . . a koeficienty 7* jsou definovány rekurentním předpisem Chybová konstanta C^+1 = 7^. Podrobnosti lze najít například v [12]. Adams-Moultonovy metody dostaneme, když v (1.26) zvolíme j = 0,1,..., k — 1. Konstrukci polynomu Pk-i(t) lze přehledně vyjádřit tabulkou Čn+1 tn+2-k Pfc-l(t) fn+1 f n fn+2-k Adams-Moultonovu metodu lze zapsat ve tvaru fc-i Vn+l = ÍJn + t fn+1-j ■ j=0 Stručně ji budeme označovat jako AMk metodu. AMk metoda je implicitní, pro k = 1 je jednokroková a pro k > 1 je {k — l)-kroková, je řádu k a D-stabilní. Koeficienty AMk metod pro konstantní délku kroku a pro k = 1,2,... ,6, spolu s chybovými konstantami Ck+i a dolními mezemi intervalů absolutní stability («&, 0), jsou uvedeny v následující tabulce: k Pk,0 h,i Pk,2 Pk,3 Pk,4 As, 5 Ck+1 Ctk 1 1 1 2 — OO 2 1 2 1 2 1 12 — OO 3 5 12 8 12 1 12 1 24 -6 4 9 24 19 24 5 24 1 24 19 720 -3 5 251 720 646 720 264 720 106 720 19 720 3 160 90 49 6 475 1427 798 482 173 27 863 45 1440 1440 1440 1440 1440 1440 60 480 38 Všimněte si, že AM1 metoda je totožná s IE metodou, tj. AM1=IE, a že AM2 metoda je totožná s TR metodou, tj. AM2=TR. Koeficienty AMk metod lze určit prostřednictvím formule fc-i Vn+l =Vn+r^2 7jV3/n+l> j=0 22 kde koeficienty 7^ jsou definovány předpisem 70 = 1, 7j = lj ~ i = 1,2,.... Chybová konstanta Ck+i = 7fc- Podrobnosti viz [12]. Metody prediktor-korektor. AM metody jsou přesnější a stabilnější než AB metody. Nevýhodou AM metod je jejich implicitnost. Abychom určili yn+±, musíme řešit rovnici fc-i Vn+i = fiVn+i), kde s = l,2,...,S, a nakonec položíme Vn+l = ÍJn+1 1 fn+1 = /(Wl)!/n+l) ) dostaneme metodu prediktor-korektor, kterou schématicky označujeme P(EC)5E. Přitom P značí předpověď počáteční aproximace AB metodou, C korekci AM metodou a E vyhodnocení pravé strany /(ín+i, ž/i+i)- Zvolíme-li jako prediktor metodu ABk a jako korektor metodu AMk, dostaneme metodu, kterou značíme ABk-AMk-P(EC)5E. Její chybová konstanta je rovna chybové konstantě Cp+i korektoru, oblast absolutní stability se blíží oblasti absolutní stability korektoru AMk až pro S —► 00. Nejčastěji se používá schéma PECE, kdy se korekce provede jen jednou a pravá strana se počítá dvakrát. V dalším se omezíme právě na schéma PECE. Abychom mohli řídit délku kroku, potřebujeme znát odhad estn lokální chyby. Jestliže y*l+l = Vn+i spočteme prediktorem ABk a y**+í = y^+i korektorem AMk, pak tzv. Milnův odhad lokální chyby dává C estn = T7*--tel ~ ^n+l)' (L27) odvození viz [12]. Nakonec položíme y-a+i = y**+i + estn- (1-28) Korekce y™+1 pomocí odhadu lokální chyby estn se nazývá lokální extrapolace. Celý krok označujeme zkratkou ABk-AMk-PECLE, přičemž písmeno L vyznačuje použití lokální extrapolace. Metoda ABk-AMk-PECLE je řádu k + 1, oblast absolutní stability je větší než u prediktoru ABk ale menší než u korektoru AMk, viz [12], [1]. Alternativní odhad lokální chyby lze získat také tak, že yn+± spočteme metodou AM(k+l) a položíme estn = yn+1 - y**+1. (1.29) 23 Výslednou metodu lze označit jako ABk-AM(k+l)-PECE. Odhady (1.28) a (1.29) nejsou sice totožné, pro malé r jsou však prakticky nerozlišitelné. Konkrétně pro metodu AB2-AM2-PECLE organizujeme výpočet následovně: AB2: y*n+1 = yn + ±r(3/n - fn^) fn+l = f(tn+l,yn+l) AM2: y*n*+1 = yn + ±r(/*+1 + /„) estn = g (ž/n+l ~~ Vn+l)i ÍJn+1 = ÍJn+1 + eS^n fn+l = f^n+liVn+l)- V případě AB2-AM3-PECE metody nahradíme řádek L řádkem C: AM3: yn+1 = yn + ^r(5/*+1 + 8fn - /n_i), estn = yn+1 - y*n*+1. Startovací hodnotu y\ získat třeba pomocí EM2 metody. Řízení délky kroku a řádu metody. Kvalitní programy založené na metodách predik-tor-korektor mění jak délku kroku tak řád metody. Tak například matlabovský program odell3 používá metody ABk-AM(k+l)-PECE pro k = 1, 2,..., 12. Změna řádu se provádí současně se změnou délky kroku. Algoritmy tohoto typu se označují jako VSVO algoritmy (varíable step varíable order). Základní myšlenka je jednoduchá. Předpokládejme, že jsme vypočetli yn+± metodou ABk-AM(k+l)-PECE s krokem délky r. Určíme odhad est^ lokální chyby podle (1.28) nebo (1.29). Jestliže |est^| > e, jde o neúspěch a výpočet yn+i je třeba opakovat, v opačném případě pokračujeme výpočtem přibližného řešení yn+2- V každém případě je však třeba stanovit novou délku kroku a nový řád. Rád se může změnit nejvýše o jednotku, tj. v metodě ABk-AM(k+l)-PECE místo k může být nově také k — 1 nebo k + 1. Odhady odpovídajících lokálních chyb est^1 a est^+1 lze získat snadno, viz [22]. Nové délky kroků r^_1? t£ a t£+1 stanovíme podobně jako pro jednokrokovou metodu, r* = frie/lestfflV+V pro l = k - 1, k, k + 1, a největší z čísel r^_1? t£, r^+1 určí jak novou délku kroku tak nový řád. Konkrétně, je-li největší t£, k se nemění a jako novou délku kroku vezmeme r* = t£, je-li největší t£+1, zvětšíme k o jedničku a pokračujeme s krokem délky r* = r^+1 a je-li největší r^_1? k o jedničku snížíme a pokračujeme s krokem délky r* = r^_x. Pro krok délky r* 7^ r je třeba vypočítat hodnoty f^+i-j Pro Ki+i-j = ^n+i — ÍT* ■ To lze snadno provést interpolací pomocí fn+i, fn, ■ ■ ■, fn+i-k- Podrobný popis strategie VSVO lze najít např. v [12], [22]. Start metody není žádný problém: začneme metodou AB1-AM2-PECE, počáteční délku kroku určíme např. podle (1.23) a algoritmus VSVO se už sám rychle vyladí na správné hodnoty jak délky kroku tak řádu metody. 1.4.3. Metody zpětného derivování Při řešení tzv. tuhých problémů je třeba pracovat s metodami, které se vyznačují neomezenou oblastí absolutní stability. Metody zpětného derivování (stručně BDF podle 24 backward differentiation formulas) tuto vlastnost mají. Pro fc-krokovou metodu zpětného derivování použijeme zkrácený zápis BDFk metoda. Dostaneme ji tak, že v rovnici y'(tn+1) = f(tn+1,y(tn+1)) nahradíme derivaci y'(tn+1) pomocí derivace P'k(tn+i) interpolačního polynomu Pk(t) stupně k procházejícího body [tn+1,y(tn+1)], [tn, y(tn)], [tn+i-fc, y(tn+1-k)]. Když pak nahradíme y(tn+1_j) přibližnými hodnotami yn+1_j, j = 0,1,..., k, dostaneme metodu -Pfc(Wl) = /(Wb!/n+l)' Přehledné vyjádření aproximujícího polynomu Pk(t) uvádí následující tabulka t tn+1 tri tn+1—k Pk(t) Vn+1 Vn Vn+l-k BDFk metodu zapíšeme ve tvaru k ak,jVn+l-j = TPkfifn+l ■ j=0 Metoda BDFk je implicitní, fc-kroková, řádu k, pro k < 6 je D-stabilní. Pro konstantní délku kroku jsou koeficienty akj, Pk,o a chybové konstanty Ck+i BDFk metod uvedeny v následující tabulce: k "fc,2 "fc,3 "fc,4 "fc,5 "fc,6 Pk,0 Ck+1 ctk 1 1 -1 1 1 2 90° 2 1 4 3 1 3 2 3 2 9 90° 3 1 18 11 9 11 2 11 6 11 3 22 88° 4 1 48 25 36 25 16 25 3 25 12 25 12 125 73° 5 1 300 137 300 137 200 137 75 137 12 137 60 137 10 137 52° 6 1 360 147 450 147 400 147 225 147 72 147 10 147 60 147 20 343 18° Koeficienty BDFk metody lze získat prostřednictvím předpisů j" y2 ~^]Vn+l = T T-fn+l, ^k = y2~, Ck+1 = - , Okj-fJ ok j^J (k + l)ók podrobnosti viz např. [12] BDFk metody jsou implicitní a ve srovnání s implicitními AMk metodami mají značně větší chybové konstanty. Na druhé straně však metody zpětného derivování mají jednu ohromnou přednost, která plně ospravedlňuje jejich používání, a tou je neomezená oblast Ra absolutní stability. Pro metody BDF1 a BDF2 oblast Ra absolutní stability obsahuje 25 celou zápornou polorovinu komplexní roviny C, tj. R4 D {z G C [Rez < 0}. Takové metody se nazývají A-stabilní. Abychom mohli popsat oblast absolutní stability zbývajících BDFk metod, zavedeme si jeden nový pojem. Řekneme, že numerická metoda je A(a)-stabilní, a G (0, tt/2), jestliže její oblast absolutní stability R a obsahuje nekonečný klín Wa = {relíp G C I r > 0, \

1. (1.32) a 1, tj. jde o tuhý problém. Metoda DP54 má interval absolutní stability (—3,31; 0), takže podmínka stability (1.33) vyžaduje, aby —3,31 < lOOOr, tj. r < 0,00331. Na intervalu délky L je tak třeba til > L/0,00331 kroků délky menší než 0,00331, což je v souladu s tabulkou uvedenou v první části příkladu 1.1. Skutečně, na intervalu(10; 100) délky 90 je třeba alespoň 90/0,00331 = 27190 kroků délky 0,00331, ve skutečnosti program ode45 provedl přibližně stejný počet 30 071 — 2 953 = 27118 kroků proměnné délky. Protože přesná řešení y\ = —e~l a y2 = e~* jsou na intervalu (10; 100) téměř konstantní, rovná nule, je zřejmé, že délka kroku je omezena z důvodu stability a ne z důvodu přesnosti. □. Příklad 1.2 Uvažujme počáteční problém y' = y2-y3, t E (0,2/5), y(0) = 5. (1.34) 28 Diferenciální rovnice má dvě konstantní řešení y = 0 a y = 1: zvolíme-li počáteční podmínku y(0) < 0, pak y(t) —► 0 pro t —► oo, zatímco pro y(0) > 0 dostaneme y(t) —► 1 pro t —► oo. Zvolíme-li ô > 0 velmi malé, pak pravá strana í/2 — y3 diferenciální rovnice nabývá malých kladných hodnot, tj. funkce y(t) velmi pomalu roste a na poměrně dlouhém intervalu zůstávají její hodnoty blízké k 0. Konkrétně pro ô = 1CT4 je y(t) < 10~2 ještě pro t = 9 900, pak y(t) začíná prudce růst a pro t > 10 020 je už y(t) prakticky rovno 1. Spektrální poloměr g(fy) = \2y — 3y2\. Pro malé y « 0 je \2y — 3y2\ ~ 0, pro y blízké 1 je však \2y — 3y2\ ~ 1. Na intervalu (0,9 900) je výraz 9 900|2í/ — 3y2\ charakterizující tuhost poměrně malý, nejde zde proto o tuhý problém, takže délka kroku se řídí především přesností metody. V intervalu (9 900; 10 020) se řešení prudce mění, to mechanismus automatického řízení délky kroku zachytí a krok zkrátí. Důvodem zkrácení kroku je zde spíše prudká změna řešení než narůstající tuhost. Poté, co řešení nabude hodnotu rovnou přibližně 1, je délka kroku metody řízena stabilitou. 1 1 1 1 1 1 1 1 1 1 1 1 J i i i i i *-1-1-1-1-1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 104 Obr. 1.2 Příklad 4.3 řešený DP54 metodou, celý výpočet ■■>■ 0.98 1 1.02 1.04 1.06 1.08 1.1 1.12 x 104 Obr. 1.3. Příklad 4.3 řešený DP54 metodou, detail Úlohu (1.34) jsme řešili explicitní metodou DP54 (matlabovský program ode45) a implicitní TR metodou (matlabovský program ode23t). Oba programy jsme použili se stejnou přesností (er = 10~4, ea = 10~7). DP54 je explicitní Rungova-Kuttova metoda řádu 5 s intervalem absolutní stability (—3,31; 0). Délka kroku proto musí splňovat podmínku stability r\2y - 3y2\ < 3,31, což pro t > 10 020 znamená volit r = 3,31. TR metoda je A-stabilní metoda řádu 2, takže tato metoda délku kroku z důvodu stability nijak neomezuje. Průběh výpočtu znázorňuje pro každou z metod dvojice obrázků: horní zachycuje celý výpočet, dolní pak zvětšený výřez pro t G (9 800,11 200). V následující tabulce uvádíme 1.0001 1 0.9999 29 pro intervaly (0,£), kde £ je postupně rovno 9 900, 10 020 a 20 000, údaj pk/pf, kde pk je počet (úspěšně) provedených kroků a pf je celkový počet vyhodnocení pravé stany: i 9900 10 020 20 000 ode45 17/151 36/331 3 041/20 245 ode23t 85/169 184/382 192/396 Vidíme, že v intervalu (10 020, 20 000), kde přesné řešení je prakticky rovno 1, je TR-me-toda velmi efektivní, na zdolání intervalu (10 020,20 000) potřebovala jen 8 kroků. Zato metoda DP54 na tomto intervalu provedla 3 005 kroků délky r = 3,31, takže její použití rozhodně vhodné není. □ x 10 Obr. 1.4. Příklad 4.3 řešený TR metodou, celý výpočet 1.0001 0.9999 x 10 Obr. 1.5. Příklad 4.3 řešený TR metodou, detail Metody pro řešení tuhých problémů. Pro řešení tuhých problémů je třeba používat metody s neomezenou oblastí absolutní stability. V Matlabu jsou to metody NDF a BDF implementované v programu odel5s, TR metoda implementovaná v programu ode23t a dvě L-stabilní metody řádu 2: TR-BDF2 metoda (jde o kombinaci metod TR a BDF2) implementovaná v programu ode23tb a Rosenbrockova metoda implementovaná v programu ode23s, podrobnosti viz [15]. Programy odel5s, ode23t a ode23tb vyžadují řešení nelineárních soustav rovnic tvaru yn+i = lTÍ(tn+1,y, n+lj (1.35) 30 kde parametr 7 je charakteristická konstanta metody a ty je vektor nezávislý na yn+i. Například pro TR metodu (1.14) je 7 = I a ty = yn + irf (ín, yn). Rovnici (1.35) řešíme zjednodušenou Newtonovou metodou. Počáteční aproximaci yfXi získáme dostatečně přesnou extrapolací z hodnot yn,yn_i,.... Další aproximace y^+P získáme řešením soustav lineárních rovnic (I " r^J) (y^ " yiíi) = -nf(tn+1, ynsl) + V " yi+i, (1-36) kde I je jednotková matice a J je Jacobiova matice fy(£n+i_fc, yn+i-fc) Pro nějaké k > 0 (jedná se tedy o zjednodušenou Newtonovu metodu, pokud by J = fy(£n+i, y^-J, šlo by o klasickou Newtonovu metodu). Výpočet organizujeme takto: označíme G = I - r7J, ds = y(*+p - y21; gs = r-fí(ín+1, y^jj + ty - y^ a vypočteme nejdříve ds jako řešení soustavy lineárních rovnic Gds = gs a pak dopočítáme Yn+P = yi+i + ds. Je-li přírůstek ds dostatečně malý, klademe yn+i = y^i"^. Matice soustavy G obsahuje tři členy, které se mohou měnit: r při změně délky kroku, 7 při změně metody (třeba ve VSVO implementaci metod zpětného derivování) a J při přepočítání Jacobiovy matice. Pokud se žádný z těchto členů nezmění, zůstává matice G stejná. Toho je třeba využít: pouze při změně G provedeme výpočetně náročný LU rozklad matice soustavy G = LU, kde L je dolní trojúhelníková matice a U horní trojúhelníková matice. V následujících iteracích, kdy se matice soustavy G nemění, provádíme výpočtově nenáročná řešení dvou soustav lineárních rovnic s trojúhelníkovou maticí soustavy, tj. Lôs = gs a pak Uds = ôs. Program, který má pracovat efektivně, musí mít promyšlenou strategii, podle níž rozhodne, kdy změní J, což je výpočetně nej náročnější, a kdy jen r nebo 7, což znamená nový LU rozklad matice G. Používá se „konzervativní strategie": přepočítání Jacobiovy matice se provede až tehdy, když Newtonova metoda nekonverguje dostatečně rychle, délku kroku případně metodu změníme, až když očekávaný zisk takové akce převýší náklady spojené s LU rozkladem. Jacobiovu matici lze zadat přesně nebo ji lze spočítat přibližně numericky (Matlab užívá funkci odenumjac z privátní knihovny toolboxu matlab\funf un). Pokud je Jacobiova matice řídká a uživatel zadá pozice jejích nenulových prvků, výpočet Jacobiovy matice lze značně urychlit. Do dalších podrobností už zacházet nebudeme, zájemce odkazujeme na skvělou monografii [22] a na matlabovskou dokumentaci [15]. Významným zdrojem poučení je studium kódů jednotlivých matlabovských programů, příliš snadné čtení to však není. Program ode23s, založený na implementaci Rosenbrockovy metody, se od zbývajících programů odel5s, ode23t a ode23tb liší v tom, že se vyhýbá řešení nelineárních rovnic. V každém kroku Rosenbrockovy metody je třeba sestavit Jacobiovu matici fy(ín, yn) a vektor derivací ít(tn,yn) = {dfi(tn,yn)/dt}f=1, provést LU rozklad matice I - 7rfy(ín,yn) a užitím tohoto rozkladu vyřešit tři soustavy lineárních rovnic, podrobnosti viz [24]. 31 Statistika o činnosti matlabovských programů pro řešení ODR. Kvalitní programy uživateli vždy poskytují informaci o tom, jak úspěšně si při řešení konkrétního problému počínaly. V Matlabu se dodávají tato čísla: pk počet úspěšných kroků pn počet neúspěšných kroků pf počet vyhodnocených pravých stran f pj počet sestavených Jacobiových matic J pr počet LU rozkladů J = LU ps počet řešených soustav lineárních rovnic LUds = gs Pro ilustraci uvádíme Příklad 1.3 Robertsonův problém y'1 = -0,04y1 + l04y2y3, ž/i(0) = l, y'2 = 0,04 í/! - ÍO4^ - 3 • 107y22 , y2(0) = 0, y'3 = S-107y22, ž/s(0) = 0 popisuje koncentrace tří příměsí v chemické reakci, tj. 0 < y±,y2,y3 < 1, nezávisle proměnná t je čas, blíže viz [22]. Jacobiova matice této soustavy je /-0,04 104í/3 I04y2\ J= 0,04 -ÍO^-G-IO7^ -10S/2 . \ 0 6 • 107y2 0 / V čase t = 0, tj. pro yľ = 1, y2 = y3 = 0, má Jacobiova matice vlastní čísla {—0,04; 0; 0}. Z fyzikálních úvah plyne, že y1} y2 —► 0 a y3 —► 1 pro t —► oo. Vlastní čísla Jacobiovy matice pro yi = y2 = 0, y$ = 1, jsou { — 10 000,04; 0; 0}. Při řešení na intervalu (0,1010) je Robertsonův problém tuhý. O tom se lze ostatně snadno přesvědčit experimentálně: explicitní metody selhávají, metody pro řešení tuhých problémů zabírají. Numerickým výpočtem lze zjistit, že již pro t > 0,01 je spektrální poloměr ^(J) > 2 • 103, takže Robertsonův problém lze považovat za tuhý již na nepoměrně kratším intervalu délky řádově v jednotkách. Úlohu jsme řešili dvěma matlabovskými programy určenými pro tuhé problémy: programem ode23t (TR metoda) a programem odel5s (metody BDFk, k=l,2,... ,5). Délku kroku jsme řídili pomocí tolerancí er = 10~3, ea = 10~6, činnost programu odel5s jsme omezili tak, aby pracoval jen s BDF metodami řádů 1, 2 a 3. Do následující tabulky jsme zapsali „statistiku" výpočtu, tj. čísla pk, pn, pf, pj, pr, ps: pk pn Pf PJ pr ps ode23t 238 74 794 37 188 644 odel5s 245 15 504 11 67 458 Z tabulky plyne, že oba testované programy Robertsonův problém úspěšně vyřešily při zhruba stejných výpočetních nákladech. □ 32 2. Obyčejné diferenciální rovnice: okrajové úlohy V kapitole 2.1 ukážeme, jak řešit okrajový problém pro soustavu 0DR1 metodou střelby. V následujících kapitolách pak popíšeme tři nejznámější metody řešení okrajového problému pro ODR2: v kapitole 2.2 diferenční metodu, v kapitole 2.3 metodu konečných objemů a v kapitole 2.4 metodu konečných prvků. 2.1. Metoda střelby Okrajový problém pro soustavu ODR1 spočívá v určení funkce y (re), která v intervalu (a, b) splňuje diferenciální rovnici y'(x) = í(x,y(x)) (2.1) a v koncových bodech intervalu (a, b) vyhovuje okrajové podmínce r(y(a),y(6)) = o. (2.2) Stejně jako v kapitole 1.1 předpokládáme, že počet rovnic jakož i počet okrajových podmínek je roven d, tedy y(x) = (yi(x),y2(x),...,yd(x))T, y'(x) = (y[(x),y'2(x),...,y'd(x))T, í(x, y(x)) = (fľ(x, y(x)), f2(x, y(x)),fd(x, y(x)))T, r(y(a),y(6)) = (r1(y(a),y(6)),r2(y(a),y(6)),...,ríi(y(a),y(6)))T. Ve speciálním případě, když r(y(a),y(&)) = y (a) ~ V = o nebo-li y(a) = 77, (2.3) přechází problém (2.1)-(2.2) v počáteční problém (1.4). Pokud je však v podmínkách (2.2) obsažena alespoň jedna složka vektoru y (a) a současně také alespoň jedna složka vektoru y(b), jde o problém okrajový. Právě tímto případem se budeme v dalším zabývat. Podstatný rozdíl mezi počáteční a okrajovou úlohou spočívá v tom, že zatímco řešení úlohy s počátečními podmínkami existuje a je jediné pro dosti širokou třídu diferenciálních rovnic, u okrajové úlohy s velmi jednoduchou diferenciální rovnicí je možné, že řešení neexistuje nebo naopak, že řešení je nekonečně mnoho. Tuto skutečnost si ukažme na rovnici y' = y, jejíž obecné řešení je y = Cex. Odtud plyne, že pro okrajovou podmínku (a) y (ti) = y(l) existuje jediné řešení y = 0, (b) e y (ti) = y(l) existuje nekonečně mnoho řešení y = Cex, C libovolné, (c) e y (ti) = y(l) + 1 řešení neexistuje. 33 Metoda střelby je založena na numerickém řešení počátečního problému y'(x) = í(x,y(x)), y(a) = r), (2.4) za omezující podmínky r(77,y(&))=o. (2.5) Neznámá počáteční hodnota rj je určena implicitně rovnicí (2.5). Lichoběžníková metoda. Nechť a = x0 < X\ < ■ ■ ■ < x^ = b je dělení intervalu (a, b), hi = xi+1 — Xi je délka kroku a h = maxj h,i. Přibližné řešení yi? i = 0,1,..., N, dostaneme jako řešení soustavy d(N + 1) rovnic Yí+i =Yí + lhí[í(xí,yí) + f(a:i+i,yi+i)], r(yo,yiv) = 0. 0,l,...,iV (2.6) Soustava (2.6) je obecně nelineární. Její řešení se obvykle počítá pomocí nějaké varianty Newtonovy metody. Aby nastala konvergence, je třeba dodat dosti dobrou počáteční aproximaci y^ y(xi), i = 0,1,..., N. Lichoběžníková metoda je řádu 2, tj. je-li funkce f dostatečně hladká, pro chybu platí y(xi) ~ Yi = 0(h2), čímž se míní, že řádu 0(h2) je každá složka vektoru y{x,j) — yť. Ve speciálním případě, když í(x, y) je lineární v proměnné y a r(u, v) je lineární v obou proměnných u i v, bude soustava rovnic (2.6) lineární. Nechť tedy y' = A(x)y + q(x), Bay(a) + B6y(6) = o. Soustava (2.6) je pak tvaru [-1 - \híK(xí)\ yí + [I - \híK(xi+1)\ yi+1 Bay0 + BbyN = o , kde I je jednotková matice. Maticový zápis této soustavy je (2.7) \hi[q[xi) + q_{xi+1)} 1,2, N (Ro So Ri Si \ ( yo \ yi yjv-i V y^v / / vo \ Vl vjv-i v ° / (2.8) kde -\ - \hiK(x,i) , Si = I - \hiK(xi+1), \i = lhi[q(xi) + q(xi+1)] . 34 Simpsonova metoda. Matlab nabízí pro řešení okrajového problému (2.1)-(2.2) program bvp4c, který je založen na použití Simpsonovy formule y»+i = y» + ^i[f(^i>yi) + 4f (2^+1/2^+1/2) + f(a;i+i,yi+i)], (2.9) kde xi+1/2 = + yi+i/2 = + y»+i) + ^i[f(^i>yi) - f(zi+i>yi+i)] , Simpsonova formule (2.9) je řádu 4, tj. pro chybu platí y(xi) -yt = 0{hA). Další informace o Simpsonově formuli (2.9) lze načerpat v [23], [12] a také v [15]. V případě lineární úlohy (2.7) řešíme soustavu lineárních rovnic (2.8), matice Ri? Sj a vektory v j získáme z formule (2.9), v níž klademe í(x, y) = A(rr)y + q(rr). Odvovození vzorců pro Rj, S j a v j ponecháváme čtenáři jako cvičení. 2.2. Diferenční metoda Ve zbytku kapitoly 2 se budeme převážně zabývat lineární ODR2 -\p(x)u']' + q(x)u = f(x), x E (OJ). (2.10) Předpokládejme, že funkce p(x), p'(x), q(x) i f(x) jsou spojité, dále že p(x) > 0, q(x) > 0, a uvažujme nejdříve jednoduché okrajové podmínky u(0) = g0, (2.11a) u(í) = gl. (2.11b) Za uvedených předpokladů má úloha (2.10) - (2.11) jediné řešení. Úloha (2.10) - (2.11) může popisovat například problém tahu-tlaku prutu, tedy prutu namáhaného pouze tahem popřípadě tlakem. V tom případě je u(x) posunutí střednicové čáry prutu, p{x) = E(x)A(x), kde E(x) je Youngův modul pružnosti a A{x) je plocha průřezu prutu, q(x) je měrný odpor podloží, na němž prut spočívá, f(x) je intenzita zatížení a g0 a gi jsou předepsaná posunutí koncových bodů prutu. Jinou aplikací, popsanou stejnou rovnicí a stejnými okrajovými podmínkami, je například stacionární úloha vedení tepla v tyči. Pak u(x) je teplota, p(x) je koeficient tepelné vodivosti, f(x) — q{x)u{x) je intenzita tepelných zdrojů, g$ a gi jsou teploty koncových bodů tyče. Úlohu (2.10) - (2.11) lze přeformulovat do tvaru (2.1)-(2.2), stačí položit yi) = (u), f=f V2'PX r=(yim~g°)i (2-12) y2J \pu'J \qyi-fj \yiv)-9ij a následně řešit metodou střelby, viz kapitola 2.1. My si ale vysvětlíme jinou techniku, známou jako diferenční metoda (stručně FDM podle anglického finite difference method). FDM je klasická diskretizační metoda, zevrubně popsaná a analyzovaná v celé řadě publikací, viz např. [14], [27]. 35 Nechť O = x0 < X\ < ■ ■ ■ < Xn = l je dělení intervalu (0, l) a h,i = x,j_ — Xí_i je délka kroku. Body x,j_ nazýváme uzly a množinu {xi}f=Q uzlů nazýváme sítí. Pro jednoduchost se omezíme na ekvidistantní dělení, takže hi = h = Í/N a Xi = ih, i = 0,1,..., N. Splnění rovnice (2.10) budeme vyžadovat pouze ve vnitřních uzlech sítě, tj. -(pu')'\x=Xi + qíu(xl) = fl, i = l,2,...,N-l, (2.13) kde q,i = q(xi) a f i = f(xi). Clen —(pu')'\ nahradíme diferenčním podílem. Pomocí označení Xí-l/2 =Xi~\h, Xi+1/2 =Xi + \h, Pi-X/2 = p(^í-l/2) , Pí+1/2 = p(xi+1/2) lze přibližně položit ^ pu\x=xi+1/2 Pu \x=xi_1/2 _ pi+i/2u'(xi+1/2) -Pi-ipu'ixi-tp) \X-Xi fa fa u(xi+1) - u(xí) u(xí) - m(^_i) . Pt+1/2 h ~Pt-1/2 h h Užitím Taylorova rozvoje snadno ověříme, že tato aproximace je řádu 0(h2), tj. že platí / ,y| -Pi-l/2u{xi_l) + {pi_l/2+Pi+l/2)u{xi)-pi+l/2u{xi+l) 2 -i.Pu)\x=x% =-fa2-+ 0{h). (2.14) Po zanedbání chyby 0(h2) tak z (2.14) a (2.13) dostáváme pro přibližné řešení u,i « u(x,j) soustavu rovnic ~Pí-l/2Uí-l + (Pí-1/2 + Pí+1/2 + h2qi)ui - pi+1/2ui+1 _ — Ji i l^-l^J h2 pro i = 1,2,... ,N — 1. Z okrajových podmínek máme u0=g0, (2.16a) un = Qi ■ (2.16b) To nám umožní dosadit do první rovnice uq = Qq a člen —pi/29o/h2 převést na pravou stranu. Podobně v poslední rovnici položíme = ge a člen —pN-i/29e/h2 převedeme na pravou stranu. Matice takto upravené soustavy rovnic (2.15) je třídiagonální, pozitivně definitní, diagonálně dominantní (pro qri > 0 ryze). Tyto vlastnosti zaručují, že matice soustavy je regulární a soustava rovnic má jediné řešení. Jsou-li funkce p, q a / dostatečně hladké, pak pro chybu platí uix^-u^Oih2), (2.17) přesná formulace a příslušný důkaz viz [27]. 36 Okrajové podmínky s derivací. Okrajové podmínky (2.11) se nazývají Díríchletovy. V aplikacích se objevují ještě další typy okrajových podmínek. Podmínky p(0)u'(0) = a0u(0) - p0, (2.18a) -p(£)u'(e) = aeu(£) - fy (2.18b) se nazývají Newtonovy nebo také Robínovy. Pokud koeficient obsahující člen «o resp. on chybí, hovoříme o Neumannově okrajové podmínce. Aby byla zaručena jednoznačná existence řešení, předpokládejme a0 > 0, cti > 0. Jestliže uvažujeme současně obě podmínky (2.18a) i (2.18b), pak předpokládejme navíc buďto «o > 0 nebo on > 0 nebo q(x) > 0 alespoň na části intervalu (0,£). Interpretujeme-li Newtonovy okrajové podmínky v úloze tahu-tlaku prutu, je pu' normálová síla, a0, cti jsou tuhosti pružin v koncových bodech prutu a fy, fy jsou zadané síly působící na koncích prutu. V úloze stacionárního vedení tepla je pu' tepelný tok, a0, cti Jsou koeficienty přestupu tepla a fy, fy jsou zadané tepelné toky na okrajích tyče. Tepelné toky se často uvažují ve tvaru fy = QqUq, fy = a^, kde Uq resp. u\ je vnější teplota okolo levého resp. pravého konce tyče. V případě zadané okrajové podmínky (2.18a) eventuálně (2.18b) musíme sestavit rovnici umožňující výpočet neznámé u0 eventuálně u^. Ukažme si to třeba pro případ předepsané okrajové podmínky (2.18a). Pomůžeme si tak, že budeme požadovat, aby rovnice (2.10) platila také v levém krajním uzlu xq = 0, tj. aby rovnice (2.13) byla splněna rovněž pro i = 0. Clen — (pu')'\x_ vyjádříme nejdříve pomocí jednostranné diference a pak pomocí centrální diference a vztahu (2.18a), tj. pu'\ — pu'\ ( lx=x1/2 ľ \x=x0 - (PU ) \x=x0 =--p- + °(h) = u(xi) - u(x0) í \ a -Pi/2 h—+ aou{x0) - fy -+ 0{h). \h Zanedbáme-li chyby, dostaneme rovnici pro neznámou uq (pi/2 + ha0 + lh2q0)u0 - p1/2u1 11 -h*--hPo + 2f°- (2'19a) V případě okrajové podmínky (2.18b) postupujeme obdobně, tj. požadujeme splnění rovnice (2.13) také pro i = N, a po úpravách obdržíme rovnici pro neznámou u n -Pn-i/2UN-! + (pjv-1/2 + hae + \h2qN)uN 11 ,„1(lU h2 ti T Je-li předepsána okrajová podmínka (2.18a), zapíšeme jako první rovnici (2.19a), pak následují rovnice (2.15) pro i = 1, 2,..., N — 1, a je-li předepsána okrajová podmínka (2.18b), připojíme jako poslední rovnici (2.19b). Matice výsledné soustavy rovnic je opět 37 třídiagonální, pozitivně definitní a diagonálně dominantní. I když se při odvození rovnic (2.19) dopouštíme chyby řádu 0(h), pro chybu přibližného řešení platí zase vztah (2.17). Rovnice s konvekčním členem. V aplikacích často vzniká potřeba řešit poněkud obecnější rovnici - [p(x)u' - r(x)u}'+ q(x)u = f(x), xe(0,£). (2.20) Tato rovnice popisuje například transport chemické příměsi v kapalině. V tom případě je u(x) koncentrace příměsi v kapalině, r = gv, kde g je hustota kapaliny a v její rychlost, p(x) je koeficient difúze a f(x) — q(x)u(x) je intenzita objemového zdroje příměsi. Rovnici (2.20) lze použít také k popisu teplotního pole kapaliny. Pak u(x) je teplota, r = cgv, kde c je tepelná kapacita, g hustota a v rychlost kapaliny, p je tepelná vodivost a f(x) — q(x)u(x) je intenzita vnitřních tepelných zdrojů. S přihlédnutím k typické fyzikální interpretaci říkáme, že — (pu')' je difúzni clen, (ru)' je konvekční člen a / — qu je zdrojový člen. Newtonovy okrajové podmínky pro rovnici (2.20) uvažujeme ve tvaru (pu' - ru)\x=Q = a0u(0) - j3Q , (2.21a) -{pu'-ru)\x=e= aeu(i)-l3e. (2.21b) Předepisuje se tedy hodnota tzv. toku —(pu' — ru) veličiny u, člen —pu' je difúzni tok a člen ru je konvekční tok. Konvekční člen vyjádříme ve vnitřních uzlech pomocí centrální diference (ru)'\x=xi = rm/2M(^+1/2) " r^u^-^ + 0(h2) (2.22) a v koncových uzlech pomocí jednostranné diference ( V| 7-1/2^(2:1/2) -r-Qu(xQ) 2 (2.23) / VI rNu(xN) - rN-i/2u(xN-i/2) M L=XN =-p-+ °(h) ■ V (2.22) a (2.23) vyjádříme konvekční tok ru ve středech :Tj_i/2 úseček (xí-i,Xí) pomocí aritmetického průměru hodnot u v koncových bodech, stručně centrální aproximace, ri-i/2u(xi-i/2) = 2^-i/2K^-i) + u(xí)] + 0(h2), i = l,2,...,N. (2.24) Pomocí (2.21), a po zanedbání chybových členů v (2.22)-(2.24), obdržíme soustavu rovnic, jejíž tvar vyjádříme prostřednictvím dříve odvozených rovnic (2.15) a (2.19) takto: na levou stranu rovnice In^uo + u^/h, přidáme člen <{ \[ri+i/2(ui + ui+1) - rj_i/2(tíi-i + Ui)]/h , (2.25) -\rN-i/2(uN-i + uN)/h. 38 Matice takto vzniklé soustavy rovnic je třídiagonální a nesymetrická. Dá se ukázat, že když !%i-i/2| p. Pak je účelné vyjádřit konvekční tok ru ve středech :Tj_i/2 úseček {xí-i,Xí) takto: Íri-i/2u{xi-i) + 0{h) pro rj_i/2 > 0, (2.27) ri-i/2u(xi) + 0(h) pro rj_i/2 < 0. Aproximace (2.27) je z fyzikálního hlediska přirozená: informaci o řešení u v uzlu Xi čerpáme ze znalosti řešení proti „proudu", proti „větru": pro rj_i/2 > 0 „fouká zleva", proto použijeme hodnotu ■u(xi_1) v uzlu Xí_i ležícím nalevo od bodu rEj_i/2, pro rj_i/2 < 0 „fouká zprava" a proto použijeme hodnotu u(x,j) v uzlu x,j_ ležícím napravo od bodu :Tj_i/2. Jednostrannou aproximaci konvekčního toku podle (2.27) nazýváme upwind aproximací. Pomocí označení a+ = max(a, 0), aT = min(a, 0), kde a je libovolné číslo, lze (2.27) zapsat v kompaktním tvaru ri-x^u^i-xfr) = r+_1/2u(xi-i) + r7_1/2u(xi) + 0(h) . (2.28) Pomocí (2.21), a po zanedbání chybových členů v (2.22), (2.23) a (2.28), obdržíme výslednou soustavu rovnic, jejíž tvar vyjádříme prostřednictvím dříve odvozených rovnic (2.15) a (2.19) takto: na levou stranu rovnice (2.19a) j ( ir+2u0 + r-/2Ul]/h, (2-!5) \ přidáme člen < [rt+i/2ui + r7+i/2ui+i]/h ~ irti/2ui-i + Vi/2ííi]//i. (2-29) (2.19b) J 1^ -[r+_i/2uN_1+rN_i/2uN]/h. Matice takto vzniklé soustavy je regulární nezávisle na jemnosti dělení intervalu (0,£). Pro chybu však platí jen u(xí) -Ui = 0(h) . (2.30) Pro dosažení přesnosti řádu 0(h2) je třeba používat přesnější upwind aproximace konvekčního toku, viz např. [25], [6]. Příklad 2.1. Zabývejme se řešením modelové úlohy -£U" + U' = 0 pro a: G (0,1), u(0) = 0, u(l) = 1, 39 kde e > O je konstanta. Přesné řešení 1 - ex/£ u(x) ol/e je rostoucí funkce. Diskretizací konvekčního členu pomocí centrální diference (2.252), tj. pomocí druhého vztahu v (2.25), dostaneme rovnici Ui+1 - 2Ui + Ui+1 - Ui-x 0. h2 2h kterou lze při označení k := h/e zapsat ekvivalentně ve tvaru -(1 + ^KJUi-x + 2uí - (1 - \k)uí+1 = 0. Pro k = 2 dostaneme Ui = Ui-i, takže řešení je Ui = 0 pro i < N a u^ = 1. Pro c / 2 snadným výpočtem ověříme, že diferenční rovnici vyhovuje řešení k kde Ci a C2 jsou konstanty, které určíme z okrajových podmínek. Pro k > 2 přibližné řešení Ui osciluje okolo C±, což je v rozporu s chováním přesného řešení u, rostoucí posloupnost {ui}^L0 dostaneme jen pro \k\ < 2, což je podmínka (2.26) pro p = e, r = 1. Naši modelovou úlohu vyřešíme také upwind technikou. Konvekční člen nahradíme levostrannou diferencí podle (2.29) a dostaneme rovnici ui+1 - 2uí + Ui-i Ui Ui-l 0. h2 h Diferenční rovnici vyhovuje řešení Ui = C1 + C2(1 + k)í, které je rostoucí nezávisle na velikosti k. □ 2.3. Metoda konečných objemů (stručně FVM podle anglického finite volume method) se používá především pro řešení problémů proudění ve více prostorových proměnných, např. viz [25], [7]. Princip této metody však lze objasnit i pro jednodimenzionální úlohu popsanou diferenciální rovnicí (2.20) a okrajovými podmínkami (2.11) a (2.21). Ke každému uzlu Xi přiřadíme konečný objem E>i (stručně buňku) takto: pro vnitřní uzly Bi = (2^-1/2,^+1/2) , * = 1, 2,..., N - 1, a pro koncové uzly B0 = {0,x1/2), BN = {xN-i/2,ty- Zřejmě (0,£) = {jf=0Bi. Integrací rovnice (2.20) přes buňku Bi dostaneme bilanční rovnici / — \pu' — ru]' dx + / qudx = f dx. (2.31) J Bi J Bi J Bi 40 Předpokládejme nejdříve, že Bri přísluší vnitřnímu uzlu. Pak dostaneme _ rxi+i/2 rxi+i/2 -[pu1 - rw]z=^±í/2 + / qudx = f dx. (2.32) Difúzni tok aproximujeme pomocí centrální diference, Pí-i/2u'(xí_1/2) = pti/2u(x*)-u&-i) + 0^2)j i = 1,2,..., N-1, (2.33) a konvekční tok aproximujeme pomocí centrální aproximace (2.24) resp. upwind aproximace (2.28). Integrály v rovnici (2.32) spočteme obdélníkovou formulí, qudx = hqiu(xi) + O(tr) , / /dx = hfc + O(tr). Zanedbáme-li chyby, dostaneme aproximaci bilanční rovnice (2.32) ve tvaru (2.252) resp. (2.292). Je-li v koncovém bodě x = 0 resp. x = l předepsána Newtonova okrajová podmínka (2.21a) resp. (2.21b), je třeba uvážit bilanční rovnici (2.31) také pro buňku Bq resp. B^. Tak například pro buňku Bq máme — (pu' — ru) _ + / qudx = f dx. Difúzni tok v bodě xľ/2 aproximujeme centrální diferencí (2.33), konvekční tok v bodě Xi/2 aproximujeme pomocí (2.24) resp. (2.28), člen p(0)u'(0) — r(0)u(0) vyjádříme pomocí okrajové podmínky (2.21a) a integrály spočteme levostrannou obdélníkovou formulí, 1/2 ľxl/2 qu dx = 7;hq0u(x0) + 0(h2) , \ f dx = ^hf0 + 0(h2) . x0 Jx0 Zanedbáme-li chyby, dostaneme rovnici (2.25i) resp. (2.29i). Rovnici (2.25a) resp. (2.29a) odvodíme obdobně z bilanční rovnice (2.31) zapsané pro buňku B^. Metodou konečných objemů jsme tedy dostali stejné rovnice jako rovnice odvozené metodou diferenční. Přednosti metody konečných objemů ve srovnání s metodou diferenční vyniknou až při řešení parciálních diferenciálních rovnic ve více prostorových proměnných. 2.4. Metoda konečných prvků Nej univerzálnější metodou diskretizace okrajových úloh je metoda konečných prvků (stručně FEM podle anglického finíte element method, česky MKP). V úlohách mechaniky pevné fáze je to jednoznačně nej používanější metoda. I když přednosti MKP lze plně ocenit teprve u úloh ve dvou a třech prostorových proměnných, podstatu metody lze objasnit i na jednorozměrné úloze. Z řady publikací věnovaných MKP zmiňme např. [8], [28], [2]. Východiskem pro diskretizaci metodou konečných prvků je slabá formulace okrajové úlohy. Proto si teď ukážeme, jak z klasické formulace (2.10), (2.11), (2.18) formulaci slabou dostaneme. Nejdříve zavedeme následující 41 Označení. Symbolem C(0,£) budeme značit prostor všech funkcí, které jsou v intervalu (0,£) spojité, a symbolem Ck(0,£) pak prostor všech funkcí, které jsou v intervalu (0,£) spojité spolu se svými derivacemi až do řádu k včetně (C podle anglického continuous). Říkáme, že bod a je pro funkci f(x) bodem nespojítostí prvního druhu, existuje-li v a konečná limita zprava i zleva [označme tyto limity f (a + 0) resp. f (a — 0)] a je-li f(a + 0)rf(a-0). Funkce f(x) definovaná v intervalu (0, £) se nazývá po částech spojitá v intervalu (0, £), je-li v (0,£) spojitá s výjimkou konečného počtu bodů, v nichž má nespojitost prvního druhu. Prostor funkcí po částech spojitých v intervalu (0,£) označíme PC(0,£) (PC podle anglického piecewise continuous). Symbolem PCk(0,£) značíme prostor funkcí, které jsou v intervalu (0,£) spojité spolu se svými derivacemi až do řádu k — 1 včetně, a jejichž k-tá derivace je v intervalu (0,£) po částech spojitá. V dalším budeme používat zejména prostor C1(0, £) funkcí, které jsou v intervalu (0, £) spojité spolu se svou první derivací, a prostor PCľ(0, £) funkcí, které jsou v intervalu (0, £) spojité a mají v něm po částech spojitou první derivaci. Slabá formulace. Začneme tím, že zavedeme pojem testovací funkce: funkci v G Cľ(0,£) nazveme testovací, jestliže v = Ov tom krajním bodě intervalu (0, £), v němž je předepsána Dirichletova okrajová podmínka. Pro konkrétnost se omezíme na okrajové podmínky (2.11a) a (2.18b), takže testovací funkce v splňuje v(0) = 0. Násobme rovnici (2.10) testovací funkcí v a integrujme přes (0,£). Integrací per-partes členu J^[—(pu')']v dx a následným užitím okrajové podmínky (2.18b) a vztahu v(0) = 0 obdržíme r-l r-l fvdx = / [—(pu')' + qu] vdx = — pu'v^ + / [pu'v' + quv] dx = = [ct£u(£) — f3e]v(£) + / [pu'v' + quv] dx. Jo Odvodili jsme tedy, že řešení u úlohy (2.10), (2.11a) a (2.18b) musí splňovat kromě Dirichletovy okrajové podmínky -u(0) = Qq také rovnost [pu'v' + quv] dx + aeu(£)v(£) = fvdx + f3ev(£) (2.34) Jo pro každou funkci v G Cľ(0,£), v(0) = 0. Okrajová podmínka (2.18b) Newtonova typu, která se stala součástí integrální rovnice (2.34) a je tak automaticky splněna, se nazývá přirozenou okrajovou podmínkou. Dirichletovu okrajovou podmínku (2.11a), která součástí rovnice (2.34) není a jejíž explicitní splnění proto musíme vyžadovat, nazýváme podstatnou nebo také hlavní okrajovou podmínkou. Rovnice (2.34) je dobře definována i v případě, kdy funkce u a v jsou z prostoru X = PCľ(0,£). Testovací funkce pak volíme z prostoru V = {v G X\ v(0) = 0} testovacích funkcí a řešení u z množiny W = {v G X\ v(0) = Qq} přípustných řešení. Dále označíme a(u, v) = [pu'v' + quv] dx + aiu[£)v{£), L(v) = / fvdx + pev(£). Jo 42 Pak úlohu najít u EW splňující a(u,v) = L (v) \/v E V (2.36) nazýváme slabou formulací problému (2.10), (2.11a), (2.18b). Řešení úlohy (2.36) nazveme slabým řešením. Slabá formulace je obecnější než formulace klasická, neboť klade nižší nároky na hladkost dat: klasická formulace: p E C^O^r, / G C(0,£) , (2.37a) slabá formulace: p,q, f G PC(0, £) . (2.37b) Jestliže p > 0, q > 0, cti > 0 a když platí (2.37b), pak úloha (2.36) má jediné slabé řešení. Ukázali jsme si, že klasické řešení je vždy také řešení slabé, viz odvození rovnice (2.34). Opak obecně neplatí, tj. slabé řešení nemusí být řešení klasické. Jsou-li však funkce p, q a / dostatečně hladké, konkrétně platí-li podmínky (2.37a), pak lze dokázat, že slabé řešení u G C2(0,£) je řešení klasické. Slabá formulace má v úloze tahu-tlaku prutu význam principu virtuálních posunutí a samotné testovací funkce v E V mají význam virtuálních posunutí Su přípustných řešení u E W. Slabá formulace je tedy zcela přirozená, neboť konkrétně pro úlohu tahu-tlaku prutu popisuje základní fyzikální zákon mechaniky kontinua. Slabá formulace pro všechny kombinace okrajových podmínek. Uveďme si tvar V, W, a(u,v) a L (v) pro všechny možné kombinace okrajových podmínek. (DD) Okrajové podmínky (2.11a), (2.11b) V = {v E X | v(0) = v{£) = 0}, W = {v E X | v(0) = g0, v {i) = ge}, re ŕ (2.38-DD) a(u,v) = / [pu'v1 + quv] dx, L[v) = l fvdx. Jo Jo (DN) Okrajové podmínky (2.11a), (2.18b) V = {v E X | v(0) = 0}, W = {v E X | v(0) = g0}, re re (2.38-DN) a[u, v) = [pu'v' + quv] dx + ctiu[€)v{€), L (v) = fv dx + (3ev(£). Jo Jo (ND) Okrajové podmínky (2.18a), (2.11b) V = {v E X | v(£) = 0}, W = {v E X\v(£) =ge}, (2.38-ND) a(u, v) = j [pu'v' + quv] dx + a0u(0)v(0), L(v) = j fvdx + /3ov(0). (NN) Okrajové podmínky (2.18a), (2.18b) v = X, W = X, a{u, v) = j [pu'v' + quv] dx + a0u(0)v(0) + aeu(£)v(£), ^ ^ L(v) = í fvdx + p0v(0) + Pev(£). Jo 43 Ve všech čtyřech případech je zaručena jednoznačná existence slabého řešení úlohy (2.36), pokud p > 0, q > 0, a0 > 0, cti > 0 a když platí (2.37b). V případě úlohy (2.38-NN) je třeba navíc předpokládat «o > 0 nebo on > 0 nebo q > 0 alespoň na části (0,£). Diskretizace užitím lineárního prvku. Na intervalu (0,£) zvolíme dělení 0 = x0 < x\ < ■ ■ ■ < xn = £ a na každé úsečce {xí-i,Xí) délky hi = Xi — Xi-i hledáme přibližné řešení U(x) ve tvaru lineárního polynomu procházejícího body [xí-i, Ui-i] a [rr^iíj], takže U(x) = Ui_iWi_i(x) +UíWí(x), kde Wj_1(rr) Wi(x) hi h,, Funkce U (x) je tedy na celém intervalu (0, £) po částech lineární funkcí^ určenou předpisem n U(x) = ^2uíWí(x), (2.39) i=0 kde Wi(x) se jsou tzv. bázové funkce, lineární na každé úsečce (rrfc-ij^fc) a takové, že f 1 pro i = j, Wi(Xj) = < { 0 pro i ý j- 0 = x0 Xl Xi-1 Xi+1 XN-1 XN = l Obr. 5.1. Lineární Lagrangeovy bázové funkce Úsečku {xí-i,Xí), na které je definována lineární funkce určená svými hodnotami v uzlech rEj_i a xÍ7 nazýváme Lagrangeovým lineárním prvkem nebo také elementem. Délku největšího dílku dělení {xi}f=Q označíme jako h, tj. h = max^-K^v hi. Prostor všech po částech lineárních funkcí (nebo-li lineárních splajnů) označme X^. Zřejmě X^ C X. Nechť Vh = V H Xh a Wh = W fl Xh- Pak přibližné řešení U, tzv. MKP řešení, obdržíme z diskrétní slabé formulace najít U G Wh splňující a,h(U,v) = Lh(v) \/v G Vh- (2.40) Přitom index h v a,h(U, v) resp. Lh(v) značí, že integrál Jq[pU'v' + qllv] dx v a(U, v) resp. Jq fvdx v L (v) počítáme numericky kvadraturní formulí řádu alespoň jedna. 44 V dalším budeme pro konkrétnost uvažovat slabou formulaci (2.38-NN). Označíme-li Ql( i=0 n n \ / AT \ ^ AjWj, ©í^í ) - Lft í O^ J j=0 i=0 / V i=0 / ,i=o (2.42) 0' [KA - Fl kde © = (0O, ©i,..., Qn)t, K = {kíj}^j=0 pro % = ah(wj,Wi), A = (A0, Ai,..., AN) a F = (F0, F1}..., F^)T pro Fj = Lh(w,j). Protože © je libovolný vektor, musí platit KA = F. (2.43) Matice K bývá označována jako matice tuhosti a vektor F jako vektor zatížení. Toto pojmenování pochází z prvních aplikací MKP v pružnosti a stalo se univerzálním označením pro matici soustavy a pro vektor pravé strany v soustavě rovnic vzniklé diskretizací jakékoliv úlohy MKP. Soustavu rovnic (2.43) sestavíme pomocí tzv. elementárních matic tuhosti K* a elementárních vektorů zatížení YL příslušných elementům {xí-i,Xí), i = 1, 2,..., N. Pomocí obdélníkové formule vyjádříme Q\pU'v') =hlPl_1/2 A,t - A,_! 0, - [@i]TKilAi h,: h,: kde ©' ©,-i ©,: K Pí-l/2 a A* hi v -1 1 Dále pomocí lichoběžníkové formule obdržíme Q\qUv) = Ih^e^A^ + q^Aí) = [0fKí2A!, kde T^í2 _ u (Qí-i 0 2 M o qi A,_i A, 45 položíme K* = K*1 + K*2. Nakonec, opět pomocí lichoběžníkové formule, dostaneme Qi(fv) = ±hi(fi_1Qi_1 + fiQi) = [&i]TFi, kde P = i^ Z rovnice 0 =ah(U,v) - Lh(v) = ' n Ql(pU'v' + qUv) + a0U(x0)v(x0) + aíU(xN)v(xN) fi n Q\fv) + AM^o) + Pev(xN) Y ^[©^[KW - F] + 0o[«oAo - fa] + Gjvt^Ajv - Ä í=i a z rovnice (2.42) tak dostaneme rovnost 0T [KA - F] = ^[©^[KW - F] + GoKAo - A] + QN[^AN - fa], (2.44) í=i z níž plyne postup, jak pomocí elementárních matic K* = {klrs}2 s=l, elementárních vektorů F* = (Fl,Fi,)T a čísel a0, fa, a^, fa sestavit globální matici K a globální vektor F: stačí srovnat členy se stejnými indexy u parametrů 0 a A (pro určení prvků matice K) nebo jen u parametru 0 (pro určení prvků vektoru F) na levé a na pravé straně rovnice (2.44). Struktura matice soustavy K a vektoru pravé strany F je patrná z tabulky 2.1. k\x + ao k1 k1 ft21 ^22 ^11 k2 ft12 ^ + k2 h2 4- P ^22 ' "41 F| + i? h1^-1 4- kN ^22 ' "41 ft12 fn-i + fn kN &22 + ^ F2N + fa Tab 2.1: Matice soustavy K a vektor pravé strany F. Pro ekvidistantní dělení je výsledná soustava rovnic (2.43) stejná jako soustava rovnic (2.19a), (2.15) a (2.19b), kterou jsme odvodili diferenční metodou. Výpočet probíhá podle následujícího algoritmu: 1) Matici K a vektor F vynulujeme. 46 2) Postupně procházíme jednotlivé prvky {xí_i,Xí), i = 1, 2,..., N, na každém z nich vypočteme elementární matici K* = {klrs}^.s=1 a elementární vektor F* = {F^=1 a koeficienty klrs resp. F%r přičteme k odpovídajícím prvkům matice K resp. vektoru F v souladu s tabulkou 2.1. 3) Matici K a vektor F modifikujeme podle uvažovaných okrajových podmínek: a) Je-li v levém krajním bodě předepsána Dirichletova okrajová podmínka (2.11a), odstraníme první řádek matice K a první řádek vektoru F, pak od pravé strany odečteme první sloupec matice K násobený předepsanou hodnotou g§ a nakonec vynecháme také první sloupec matice K. b) Je-li v levém krajním bodě předepsána Newtonova okrajová podmínka (2.18a), přičteme k prvku v levém horním rohu matice K číslo «o a k prvnímu prvku vektoru F přičteme číslo fl0. c) Je-li v pravém krajním bodě předepsána Dirichletova okrajová podmínka (2.11b), odstraníme poslední řádek matice K a poslední řádek vektoru F, pak od pravé strany odečteme poslední sloupec matice K násobený předepsanou hodnotou ge a nakonec vynecháme také poslední sloupec matice K. d) Je-li v pravém krajním bodě předepsána Newtonova okrajová podmínka (2.18b), přičteme k prvku v pravém dolním rohu matice K číslo a^ak poslednímu prvku vektoru F přičteme číslo f}^. 4) Vyřešíme soustavu lineárních rovnic Ku = F. Podle zvolených okrajových podmínek tak získáme u = (ui,U2, ■ ■ ■ ,un-i)t v případě okrajových podmínek (2.11a), (2.11b), u = (ui,U2, ■ ■ ■ ,un)t v případě okrajových podmínek (2.11a), (2.18b), u = (u0,ui,... ,Mtv-i)t v případě okrajových podmínek (2.18a), (2.11b), u = (u0,ui,... ,Un)t v případě okrajových podmínek (2.18a), (2.18b). Pro chybu u — U a její derivaci platí za předpokladu u G C2(0,£) odhad u-U = 0(h2), u'-U' = 0(h). (2.45) 47 3. Parciální diferenciální rovnice Parciální diferenciální rovnice (stručně PDR) vyjadřuje vztah mezi funkcí několika proměnných a jejími parciálními derivacemi. Parciální rovnice a jejich soustavy jsou matematickým modelem mnoha technických úloh. Četné modely byly vytvořeny již v minulých staletích, jejich praktické řešení však umožnily teprve výkonné počítače. Prostředky klasické matematické analýzy se zkoumá existence, jednoznačnost, hladkost a další vlastnosti řešení v závislosti na koeficientech rovnice, okrajových a počátečních podmínkách a na oblasti, ve které má být rovnice splněna. Tuto oblast budeme standardně značit symbolem íl. Pro některé jednodušší úlohy lze těmito prostředky najít i přesné řešení, často ve tvaru nekonečné řady. Převážnou většinu těchto úloh však dovedeme řešit pouze přibližně, numericky. Označení eliptické, parabolické nebo hyperbolické získaly PDR na základě formální podobnosti s rovnicemi kuželoseček. Stacionární úlohy jsou na čase nezávislé, úlohy nestacionární na čase závisejí. K jednoznačnému určení řešení nestačí samotná diferenciální rovnice, je nutno zadat ještě okrajové podmínky a u nestacionárních úloh také počáteční podmínky. Ve třech následujících kapitolách uvedeme nej jednodušší rovnice druhého řádu. Omezíme se přitom na PDR ve dvou proměnných, tj. ve stacionárním případě jde o úlohu ve dvou prostorových proměnných x, y a v nestacionárním případě se uvažují úlohy s jedinou prostorovou proměnnou x, druhou proměnnou je čas t. Diskretizaci v prostorových proměnných provedeme pomocí diferenční metody, metody konečných objemů a metody konečných prvků. Metoda konečných prvků jednoznačně dominuje při řešení problémů mechaniky pevné fáze, zatímco pro proudění tekutin se více používají programy pracující na bázi metody konečných objemů. Několik pojmů. Uzávěr množiny M G IRd je sjednocením bodů množiny M a bodů ležících na její hranici dM. Uzávěr M značíme M, tj. M = M U dM. Oblastí rozumíme otevřenou souvislou množinu v IRd. Nechť íl je oblast. Prostor funkcí, které jsou v Cl spojité spolu se svými derivacemi až do řádu k, značíme Ck(Cl). Prostor C°(Cl) funkcí spojitých v Cl značíme stručně C (Cl). Řekneme, že funkce u je v íl po částech spojitá, jestliže Cl je sjednocením uzávěrů konečného počtu navzájem disjunktních podoblastí, tj. Cl = [jCli, Q i fl Qj = 0 pro i ^ j, a jestliže u je na každé z podoblastí ížj spojitá a spojitě prodloužitelná až do hranice, tj. existuje spojitá funkce u,i G C (Cli) s vlastností u,i = u,i v ížj. Prostor po částech spojitých funkcí značíme PC(ÍÍ). Symbolem PCk(íí) značíme prostor funkcí, které jsou v oblasti Cl spojité spolu se všemi svými derivacemi až do řádu k—l včetně a jejichž k-té derivace jsou po částech spojité. Tak třeba PCľ(íl) je prostor funkcí, které jsou v Cl spojité a jejichž první derivace jsou v íl po částech spojité. Definici funkce spojité a funkce po částech spojité na hranici lze prostřednictvím parametrického vyjádření hranice převést na definici funkce spojité a funkce po částech spojité na úsečce, viz kapitola 2.4. Pokud jde o značení, tak třeba PC {Vi) je prostor po částech spojitých funkcí na části C díl hranice oblasti íl. 48 3.1. Úloha eliptického typu 3.1.1. Formulace úlohy Buď Q omezená oblast v IR2. O hranící T = díl oblasti Q předpokládejme, že je sjednocením uzávěrů dvou navzájem disjunktních částí Ti a I^, tj. T = TilJl^, TiPlI^ = 0. Dále n = (ni,ri2)T nechť je jednotkový vektor vnější normály hranice a du du du — = n • y u = n-i--h rio— on dx dy je derivace ve směru vnější normály. Nechť p(x,y) > 0, q(x,y) > 0, f(x,y), g(x,y), a(x,y) > 0 a (3(x,y) jsou dané funkce. Naším úkolem je určit funkci u(x,y), která uvnitř Q vyhovuje diferenciální rovnici ô / du \ 0 / du \ na hranici Ti splňuje Dirichletovu okrajovou podmínku u = g(x,y) na Ti, (3.2) a na hranici T2 splňuje Newtonovu okrajovou podmínku du -p(x,y)— = a(x,y)u-P(x,y) na T2. (3.3) Je-li v (3.3) a = 0, dostaneme Neumannovu okrajovou podmínku V případě, že p je konstanta a q = 0, dělíme rovnici (3.1) číslem p a vznikne Poissonova rovnice -Au = f(x,y) ví), (3.4) kde Laplaceův operátor A aplikovaný na funkci -uje d2u d2u dx2 dy2 Rovnice Au = 0 se nazývá Laplaceva rovnice. Fyzikální význam. Úlohu (3.1)-(3.3) můžeme interpretovat jako stacionární vedení tepla v nekonečném hranolu o průřezu íl. Pak u(x,y) znamená teplotu, p(x,y) je tepelná vodivost, —pVu je tepelný tok, q(x,y) je měrný tepelný odpor, f(x,y) je intenzita vnitřních tepelných zdrojů, okrajová podmínka (3.2) předepisuje teplotu na povrchu a v okrajové podmínce (3.3) je —pdu/dn tepelný tok ve směru vnější normály, a(x,y) je koeficient přestupu tepla a (3(x,y) je předepsaný tepelný tok. Poissonova rovnice s homogenní Dirichletovou okrajovou podmínkou u = 0 vyjadřuje např. průhyb membrány upevněné na okraji a zatížené tlakem úměrným funkci /. Laplaceova rovnice s Neumannovou okrajovou podmínkou popisuje např. potenciální proudění: u je potenciál vektoru rychlosti, du/dx je rr-ová a du/dy y-ová složka vektoru 49 rychlosti, f3 je normálová rychlost s vlastností jTf3ás = 0 (reprezentující nestlačitelnost tekutiny). Potenciál u není určen jednoznačně: je-li u řešení, je také u + C řešení, kde C je libovolná konstanta. Rychlost [du/dx, du/dy] však už jednoznačně určena je. Existence a jednoznačnost řešení. Podmínky zaručující existenci a jednoznačnost klasického řešení u G C2(Q) problému (3.1)-(3.3) jsou komplikované a proto je zde uvádět nebudeme. V praktických aplikacích vystačíme s existencí tzv. slabého řešeni, které existuje za podmínek v aplikacích běžně splněných. V dalším budeme předpokládat, že Q je mnohoúhelník, p > 0, q > 0, / jsou po částech spojité v Q, g je spojitá na r1? a > 0, f3 jsou po částech spojité na r2, a pokud V = T2, pak buďto q > 0 na části Q nebo a > 0 na části r2. Za těchto předpokladů existuje jediné slabé řešení u G PC1 (Q) problému (3.1)-(3.3), viz např. [11]. Zesílením uvedených předpokladů lze docílit toho, že slabé řešení je také řešení klasické. Tyto zesílené předpoklady však obvykle odporují požadavkům praktických aplikací. 3.1.2. Diferenční metoda Diskretizace okrajové úlohy ve dvou dimenzích je analogická diskretizaci jednodimenzionální úlohy, viz kapitola 2.2. Princip metody. Abychom výklad nezatěžovali detaily nepodstatnými z hlediska numerické metody, začneme řešením Dirichletovy úlohy pro Poissonovu rovnici na čtverci, tj. řešíme rovnici (3.4) s okrajovou podmínkou (3.2) pro r\ = T, když Q = (0,£) x (0,£) je čtverec se stranou délky l. Na Q vytvoříme pravidelnou čtvercovou síť. Diferenční metoda se proto také často nazývá metoda sítí. Zvolme tedy N > 1 celé a definujme krok h = Í/N. Označme x,j_ = ih, i = 0,1,..., N, y j = jh, j = 0,1,..., N. Body [xh y j], i, j = 0,1,..., N, nazveme uzly sítě. Rovnice (3.4) má být splněna ve všech bodech [x, y] uvnitř Q, musí tedy být také splněna ve všech vnitřních uzlech, tj. d2u(xi,yj) d2u(xi,yj) dx2 dy2 Parciální derivace vyjádříme pomocí centrálních diferencí _ d2u(xi,yj) _ -ufa-uyj) + 2u(xi,yj) - u(xi+1,yj) dx2 h2 _ d2u(xi,yj) _ u(xj,yj—i) + 2u(xi,yj) -u(xí,yj+1) dy2 h2 f(xi,yj), i,3 = l,2,...,N-1. (3.5) 0(h2), (3.6a) 0(h2), (3.6b) dosadíme do rovnice (3.5) a chybové členy 0(h2) zanedbáme. Po vynásobení h2 dostaneme soustavu sítových rovnic - Mij-i + Auíj - uijj+1 - ui+1J = h2fij , i,j = l,2,...,N—l, (3.7) 50 ^ = y n Xi-i Xi Xi+1 XN = í Obr. 3.1. Síť kde Uij je aproximace u(xi,yj) a fy = f(xi,yj). Z okrajové podmínky (3.2) dostaneme u,j_j = Qij pro i = 0 nebo i = N nebo j = 0 nebo j = N, (3-8) přičemž g^- = g(x,nyj). Když z (3.8) dosadíme do (3.7) a na levé straně ponecháme pouze členy s neznámými u^, dostaneme soustavu (N — l)2 lineárních algebraických rovnic, kterou můžeme zapsat maticově ve tvaru Ku = F. Pro N = 4 má soustava rovnic (3.9) následující tvar: (3.9) 4 -1 0 -1 0 0 0 0 0 "\ / un \ / fc2/ii - ^9oi - 1-010 \ -1 4 -1 0 -1 0 0 0 0 h2fu - ^902 0 -1 4 0 0 -1 0 0 0 «13 h2fis - I" #03 - 1-014 - 1 0 0 4 -1 0 -1 0 0 «21 h2f2i - 1- 020 0 -1 0 -1 4 -1 0 -1 0 «22 = h2Í22 0 0 -1 0 -1 4 0 0 -1 «23 h2Í23 - 1- 024 0 0 0 -1 0 0 4 -1 0 «31 h2fsi - 1-041 - 1-030 0 0 0 0 -1 0 -1 4 -1 «32 h2Í32 - 1-042 0 0 0 0 0 -1 0 -1 4 y \ «33 / \ h2f33 - 1-043 " 1- 034 / Pravá strana rovnice odpovídající uzlu, který není sousedem hranice, obsahuje pouze člen h2fij, pro uzly nejbližší vrcholům čtverce přibudou dva členy s funkcí g a pro ostatní sousedy hranice jeden člen s funkcí g. 51 Soustavu (3.9) napsat v blokovém tvaru B -I O . . O -I B -I . . O O -I B . . O O O O . . -I O O O . . O B -I ° \ O O -I B / / Ui u3 uat-2 ( Fi F2 F3 FaT-2 \ Fat_i (3.10) kde I je jednotková matice řádu N —1, O je nulová matice řádu N — 1 a B je třídiagonální matice řádu N — 1 tvaru / 4 -1 0 B V 0 0 -1 4 -1 0 0 0 -1 4 0 o ,u 0 0 0 -1 o o o o 4 -1 0 o -1 4/ ,Af-l)T, * 1,2,...,N 1, F~i2 1, jsou části celkového vektoru u 1, jsou části celkového 1,2,. ..,N Vektory = (un,Ui2,.. neznámých, vektory F« = vektoru F pravých stran. Pro homogenní okrajovou podmínku g(x,y) - n ;" p - 12 okrajovou podmínku je v některých rovnicích příspěvek z hranice, přesně 0 je Fíj = h2fij, pro nehomogenní F h2h hz f- i íj + 9oj ) Fu = h2fn + g01 + g10 ,2. Fu — h2fn Fn-i,i — h f n-i.i + 9ni + 9 n-i.o 2 1 52 je počet dílků ve směru osy x a M > 1 je počet dílků ve směru osy y. Položíme h = a/N, k = b/M a definujeme x,j_ = ih, i = 0,1,..., N, a y^ = jk, j = 0,1,..., M. Síť je tvořena uzly [xí, yj] pro i = 0,1,..., N, j = 0,1,..., M. Při diskretizaci Poissonovy rovnice (3.4) postupuje analogicky jako dříve, jen místo (3.6b) použijeme d2u(xi,yj) _ -ufa^j-i) + 2u(xi,yj) - u(xhyj+1) 0(k2), (3.6b') dy2 k2 a tak místo rovnic (3.7) dostaneme pro i = 1, 2,..., N — 1 a j = 1, 2,..., M — 1 rovnice ~ui-i,j + 2iíjj — itj+ij — Uíj-i + 2iíjj — itjj+i /y • (3-7') /i2 A;2 Je-li řešení w dostatečně hladké, pro chybu platí m(^,%)-^ = 0(/i2 + A;2). (3.11') Obecnější rovnice. Při diskretizaci rovnice (3.1) aproximujeme členy —[pux]x a — [puy]y pomocí vzorce (2.14). Tak ve vnitřních uzlech dostaneme místo (3.7') rovnici -Pi-l/2,jUi-l,j + (Pi-l/2,j + Pi+l/2,j)Ujj - Pi+l/2,jUi+1j ^ h2 -Pi,j-l/2Ui,j-l + (Pi,j-1/2 + Pi,j+l/2)Uij — Pi,j+l/2Uij+l k2 -\- QijUij fij . Přitom index i ±1/2 znamená, že za argument x dosadíme x;t ± ^h, a podobně index j ±1/2 znamená, že za y dosadíme y j ±\k. Je-li řešení u dostatečně hladké, pro chybu opět platí (3.11'). Newtonova okrajová podmínka. Při diskretizaci postupujeme obdobně jako v jednodimenzionálním případě, viz odvození vztahů (2.19). Ukážeme si to na příkladu podmínky -P(a,y)^^- =ae(y)u(a,y)-pe(y), 0x)uNM 2k2 1 1 1 e 1 IQNMUnM = iÍnM + ~~faP 0, pokud [ri]i+i/2j < (3.19) u{xi+i,yj), { < 0. Clen (r2u)y(xi,yj) aproximujeme obdobně jako člen (riu)x(xi,yj), při rozhodování o hodnotě u(xi,yj+i/2) se řídíme analogií ke kritériu (3.17), tj. podmínkou \k\[r2]í,j+i/2\ NM Obdélník Bij = [(^-1/2,^+1/2) x {yj-1/2, y j+1/2)] n Q 0,1,. ..,N, j = 0,l,...,M, nazveme konečným objemem (stručně bunkou), viz Obr. 3.2. Integrací diferenciální rovnice (3.15) přes buňku dostaneme bilanční rovnici -(pux - riit)^ - (puy - r2u)y + qu] dxdy Bi f dxdy. (3.21) 56 Uvažujme nejdříve případ, kdy je vnitřní buňka, tj. když 0<2 0, pokud [r2]jvj+i/2 < u{xN,yj+1), { < 0. Čtvrtý integrál v (3.23) aproximujeme obdobně jako třetí integrál. Poslední dva integrály v (3.23) aproximujeme užitím jednostranné součinové obdélníkové formule qudxdy ^hkq^ju(xj^,yj), / / dxdy « ^hkfj^j . Po dosazení do (3.22) opět dojdeme ke stejné rovnici, jakou bychom dostali pomocí diferenční metody. Obecnější síť buněk. Mnohoúhelník Q vyjádříme jako sjednocení konečného počtu uzavřených trojúhelníků, z nichž každé dva jsou buďto disjunktní nebo mají společný vrchol nebo stranu. Vrcholy trojúhelníků nazveme uzly. Soubor všech trojúhelníků vytváří tzv. triangulací oblasti Q, v metodě konečných objemů označovanou jako primární síť, viz Obr. 3.3. Ke každému uzlu p přiřadíme buňku Pj. Sestavíme ji z přilehlých částí všech trojúhelníků s vrcholem p. Objasněme si to podrobněji. Jestliže je na celé hranici předepsána Dirichletova okrajová podmínka, stačí uvažovat jen buňky příslušné vnitřním uzlům P« ^ díl. Pro jednoduchost se omezme na tento případ. Nechť je trojúhelník s vrcholy p, Pj a Pk. Označme P^ střed strany pPj, Pik střed strany P{P\~ a P^ těžiště trojúhelníka T^. Pak do buňky Pj zahrneme čtyřúhelník 58 Obr. 3.6. Duální síť Cijk s vrcholy Pí, Píj, P^k a Pik, viz Obr. 3.4. Tento postup opakujeme pro všechny trojúhelníky T^, které obsahují uzel Pí, a sjednocením přilehlých částí C^k dostaneme buňku B>i, viz Obr. 3.5. Množina všech buněk se nazývá duální síť, viz Obr. 3.6. Diskretizace vychází opět z bilanční rovnice (3.21). Při integraci členů s derivacemi použijeme Gaussovu-Ostrogradského větu, viz např. [20], jn{%+^)dxdy=sdnPni+Qn2) ás, (3,24) kde P = — (pux — riu) a Q = — (puy — r2u), a dostaneme analog rovnice (3.22) dBi L du -p-—h (r • n)u on ds qudxdy = f dx dy . Bi J Bi (3.22') kde r • n = r\n\ + r2n2 a n = (ni,n2)T je jednotkový vektor vnější normály hranice dBi buňky Bi. P P ik ___—= _______s / s / y J Pijk / Cíjk Pj \ Pk Obr. 3.4. Cijk = B,n Tijk Obr. 3.5. Vnitřní buňka Bi 59 Klíčová je aproximace křivkového integrálu. Hranici dE>i vyjádříme jako sjednocení společných částí dB^ = dBi fl dBj buňky B i a sousedních buněk Bj, tj. dBi = [J. dBij. Protože jdB = jdB , stačí popsat aproximaci jen na dBij. To lze provést třeba takto: du -p-—h (r • nm on ds \dBtJ\ p(^)M(^,U^.(Pj')+(r-n)^J-) p.p, \rir3 I kde \dB,-Lj\ je délka ô-B^, |-Pí-P/| je délka úsečky PíPj, (r • n)íj = r (Pij) ■ = r^P^nf + r2(Píj)n2J a = (n^,n2J)T = (Pj - Pí)/\PíPj je jednotkový vektor ve směru vektoru PíPj. Zbývá aproximovat hodnotu u(Píj) ve středu úsečky PíPj. Jestliže je dominantní difúze, tj. když l\PlPJ\\(r-n)lJ\i(x,y), které jsou v Pí rovny jedné a v ostatních uzlech jsou rovny nule, tj. 61 y Wí(Pj O pro * = J, Obr. 3.7. Bázová funkce Wi(x,y) Každá funkce v G Xh může být pomocí svých hodnot v uzlech a pomocí bázových funkcí vyjádřena ve tvaru M v(x, y) = Y, v(xh Ví)wí{x, y). Dále definujme prostor testovacích funkcí Vh = {v G Xh | v(xi,yi) = 0 VPj G ľi} a množinu přípustných řešení Wft = G Xft | v(xi,yi) = g(xi,yi) MPi G Ti}. Trojúhelník, na němž je definována lineární funkce, jednoznačně určená svými hodnotami ve vrcholech, se nazývá Lagrangeův lineární trojúhelníkový prvek. Nyní už máme vše potřebné k dispozici: přibližné řešení U, tzv. MKP řešení, obdržíme z áiskrétní slabé formulace najít U G Wh splňující a,h(U,v) = Lh(v) \/v G Vh (3.29) kde ah(U,v) = ]T [QTe(p[Uxvx + UyVy]) + QTe(qUv)] + ]T Qs^aUv) Teei sees (3.30) Teei sees Přitom symbolem QTe(ip) jsme označili kvadraturní formuli pro výpočet jT Lpdxdy na trojúhelníku Te a symbolem QSe(f) formuli pro výpočet js ip ds na straně Se. Jako vhodné kvadraturní formule lze doporučit: 1) člen p[Uxvx + UyVy] integrujeme na trojúhelníku T formulí QT{ ^2 ®íwí ) ~ Lh ^ ®iWi \i=\ j=N+l i=l / \i=l N N N ^2 @i ^2 ah(wj,wi)Aj - y^Qj i=i j=i i=i m Lh(wi)- ^2 ah(wj,wi)9j j=N+l (3.35) e7 íka kde © = (0i, 02,..., 0at)t, k = {%}fi=i pro ktj = ah(wj, Wi), a = (Ai, A2,..., Ajv)t a F = (Fi, F2,..., FN)T pro = Lh(w,i) - YIjLn+i ah(wj, ™i)gj- Protože © je libovolný vektor, musí platit ka = F. (3.36) Matice K se často označuje jako matice tuhosti a vektor F jako vektor zatížení. Tato označení se přenesla z prvotních aplikací MKP pro řešení úloh pružnosti. Dá se ukázat, že matice tuhosti K je pozitivně definitní. Je také řídká a při vhodném očíslování uzlů je pásová. Vyřešením soustavy lineárních rovnic (3.36) získáme Aj = U(Pj), j = 1, 2,..., N. Za předpokladu dostatečné hladkosti slabého řešení u pro chybu metody platí u-U = 0(h2), Ux-Ux = 0(h), Uy-Uy = 0(h), (3.37) kde h je nej delší strana trojúhelníků triangulace 7. Algoritmus. Matici K a vektor F sestavíme pomocí elementárních matic KTe a elementárních vektorů FTe pro Te G 7 a elementárních matic K5e a elementárních vektorů FSe pro Se G §. Matici K budeme nazývat také globální maticí tuhosti a vektor F globálním vektorem zatížení. Elementární matice a elementární vektor na trojúhelníku. Nechť ELEM je tabulka typu NT x 3, která v řádku e obsahuje čísla vrcholů trojúhelníka Te. Uvažme jeden konkrétní trojúhelník Te triangulace 7 s vrcholy Pf = [xl,yf], P2 = [x2, y2] a P3e = y3]. Pro uzel Pre, r = 1,2,3, je r lokálním číslem uzlu na trojúhelníku Te. Globálním číslem uzlu Pre je číslo i = ELEM(e,r). Pre a P« jsou tedy jen různá označení téhož uzlu. Řešení U a testovací funkce v je na trojúhelníku Te tvaru U (x, y) = A\wl(x, y) + Ae2we2(x, y) + Ae3we3(x, y), Te (3.38) v(x, y) = Q\w\(x, y) + Q2we2(x, y) + 6|w|(a:, y), 63 kde A^ = U(Pre), Qer = v(Pre) a kde Wr(x, y) = aerx + bery + ďr je bázová funkce příslušná k uzlu Pre, r = 1,2,3. Z (3.38) a (3.39) dostaneme (3.39) a, Bea /Ai\ A^ WsJ kde Be 62 ^3 /©?\ W jsou vektory parametrů na trojúhelníku Te. Koeficienty ďr a ^ lze snadno spočítat z rovnic (P2e) + /(P3>(P3e)] kde JITeii/CPfje; + /(P2e)0^ + f{P!)Qi] = [e //(Ae)\ i IP /(ab: (3.46) (3.47) je elementární vektor na trojúhelníku Te. Elementární matice a elementární vektor na straně. Nechť SIDE je tabulka typu N s x 2, která v řádku e obsahuje čísla krajních bodů strany Se. Uvažme konkrétní stranu Se hranice r2 s koncovými body Pf = [rrf, ] a P| = , . Pro uzel Pre, r = 1,2, je r lokálním číslem uzlu na straně Se. Globálním číslem uzlu Pre je číslo i = SIDE(S,r). Pr5 a p jsou tedy jen různá označení téhož uzlu. Řešení U a testovací funkce v je na straně Se tvaru U(x,y) v(x,y) ^eiw\(x,y) + Ae2w*(x,y) eiwl(x,y) + Qe2we2(x,y) (3.4Í kde Aer = U(Pre), = v(P^) a kde w^(rr, y) je bázová funkce příslušná k uzlu Pre, r = 1, 2. Užitím formule (3.33) obdržíme kde QS°(aUv) =\\Se\[a{pz)U{P*)v{pz) + a{pz)U{pz)v{pz)} = \\Se\[a{P?)Q\&\ + a(PI)Qe2Ae2] = [05e]TK5ea5e a(Pf) K' o Op, 0 (3.49) (3.50) 65 je elementární matice na straně Se a 'al), = jsou vektory parametrů na straně Se. Podobně odvodíme qs^M =l\seMP>(PD + /?(p2bMAb)] = i|5e|[/3(pí)e? + /3(p2e)6|] = [&S*]TFS% (3.51) kde (3.52) je elementární vektor na straně Se. Sestavení soustavy rovnic. Zkombinujeme-li (3.35), (3.30), (3.44), (3.46) , (3.49) a (3.51) vidíme, že pro MKP řešení U a libovolnou testovací funkci v G Vh platí 0 = ah(U, v) - Lh(v) = &T [KA - F] = = [&Te]T [KTe ATe - FTe] + i@Sef [K5e A5e - F5e] . (3-53) TeeT 5ee§ Z této rovnosti lze odvodit pravidla pro sestavení globální matice K a globálního vektoru F z lokálních matic KTe, K5e a lokálních vektorů FTe, F5e. Postupujeme podle následujícího algoritmu: 1) Matici K řádu N a sloupcový vektor F délky N naplníme nulami. 2) Postupně pro každý trojúhelník Te G T sestavíme elementární matici KTe = {k^}f s=1, viz (3.45), a elementární vektor FTe = {FrTe}^=1, viz (3.47). Postupně pro r,s = 1, 2, 3 položíme i = ELEM(e, r), j = ELEM(e, s) a i < Naj < N, pokud <{ i < N a j > N, provedeme i < N, 3) Postupně pro každou stranu Se G S sestavíme elementární matici K5e = {k^} viz (3.50), a elementární vektor F5e = {F^e}2r=1, viz (3.52). Postupně pro r = 1, 2 položíme i = SIDE(e, r) a pokud i < N, provedeme ku <— ku + k'^, Fi <— Fi + Pr5e • 66 3.2. Úloha parabolického typu Formulace úlohy. Hledáme funkci u(x,ť) definovanou pro x G (0,£), t G (0,T), která vyhovuje diferenciální rovnici f)ii f) ŕ f)ii \ C{x)m ~ di Y^Yx) + q{x)U = f{X>th x E (0,£), í G (0, T), (3.54) Dirichletovým okrajovým podmínkám u(0,t) = g0(t), u(£,t) = ge(t), t e (0,T), (3.55) nebo Newtonovým okrajovým podmínkám p(o)^^ = "o«(o)-#,(*), ÍG(0,T), (3.56) a počáteční podmínce u(x,0) = ip(x), xe(0,£). (3.57) Možný je také případ, kdy je na jednom okraji předepsána Dirichletova podmínka a na druhém podmínka Newtonova. Úloha (3.54) - (3.57) popisuje nestacionárni vedeni tepla v tyči délky £, T je doba trvání děje. Proměnná x je prostorová, t má význam času, u(x, t) je teplota v bodě x a v čase t, funkce p, q, f, qq, ge, (3q, fy a konstanty «o, ae mají stejný význam jako v kapitole 2, c je součin tepelné kapacity a hustoty, ip je teplota tyče v čase t = 0. Tepelné toky se často uvažují ve tvaru flo(t) = qíoMq(í), fle(t) = aeuf(t), kde Uq resp. u\ je vnější teplota okolo levého resp. pravého konce tyče. Předpokládejme že funkce c, p, q, f, g0, ge, fl0, fy a ip jsou dostatečně hladké a že jsou splněny podmínky nezápornosti c > 0, p > 0, g > 0, «o > 0, > 0. Aby existovalo klasické řešení úlohy (3.55) - (3.57), je třeba pro Dirichletovy okrajové podmínky připojit ještě tzv. podmínky souhlasu V?(0) = 0o(o), = ge(0) ■ Tyto vztahy, vyjadřující soulad počáteční a Dirichletovy okrajové podmínky, nebývají v aplikacích obvykle splněny. Také funkce c, p, q, f, g0, ge, fy, fy bývají často jen po částech spojité. V tom případě má úloha (3.54) - (3.57) pouze tzv. slabé řešení (jehož přesnou definici zde uvádět nebudeme). Pro praktické účely je slabé řešení zcela vyhovující a numerické metody, které si uvedeme, budou poskytovat přibližné hodnoty takového slabého řešení. Úloha (3.54) - (3.57) je okrajovou úlohou druhého řádu vzhledem k proměnné x a počáteční úlohou prvního řádu vzhledem k proměnné t. Při její diskretizaci proto zkombinujeme postupy z prvních dvou kapitol. 67 Prostorová diskretizace se provádí metodami kapitoly 2. Pro jednoduchost budeme uvažovat úlohu s Dirichletovými okrajovými podmínkami, diferenční metodu a rovnoměrné dělení intervalu (0, l) s krokem h = £/N, tj. Xi = ih, i = 0,1,..., N. Budeme požadovat, aby rovnice (3.54) byla splněna ve vnitřních uzlech Xi, i = 1, 2,..., N — 1, tj. c(xí) du(xi, t) dt d f du dx \ dx (xí,t) + q(xi)u(xi,t) = f(xi,t). Derivaci podle x vyjádříme pomocí diference analogicky jako v (2.14), tj. d í du\ dx \ dx J (X{, í) h2 Pí-l/2u(xi_1,t) + + [Pí-i/2 +Pí+i/2]u(xi,t) - pi+1/2u(xi+uť) + 0(h2) Zanedbáme-li chybu, dostaneme soustavu rovnic CíÚiiť) h2 -Pi-l/2Ui-l{t) + (Pi-l/2 + Pi+1/2 + h2q,i)Ui(t)- (3.58) Pí+l/2Ui+l(t) = fi(t) i = 1,2,...,N- 1, ÍG(0,T) níž Ui(t) je aproximace u(xi,t), úi(t) dujjt) dt je aproximace dtu(xi,t) a fi(t) = f(xi,ť). Z okrajových podmínek (3.55) máme u0(t)=g0(t), uN(t)=ge(t), *G(0,T), (3.59) což dosadíme do (3.58) pro i = lai = N— 1. Z počáteční podmínky obdržíme Uí(0) = p(xi i = 0,l,...,N. (3.60) (3.58) je soustava obyčejných diferenciálních rovnic prvního řádu pro N — 1 hledaných funkcí Ui(t), i = 1, 2,..., N — 1, s počátečními podmínkami (3.60). Původní úlohu, tj. určení jedné funkce u(x,t), která na obdélníku (0,£) x (0,T) vyhovuje (3.54), (3.55) a (3.57), jsme nahradili počáteční úlohou pro soustavu N — 1 obyčejných diferenciálních rovnic s neznámými funkcemi Ui(t) definovanými na T 9o{t) Ui+1{t) Ui(t) Ui-!{t) ge(t) N — 1 úsečkách x = Xi, i = 1,2,..., N — 1, Xq Xi-i Xí xi+i XN Obr. 3.8. Metoda přímek t G (0,T). Tento postup se nazývá metoda přímek (anglicky method of Unes, stručně MOL). Postup, při němž se diskretizace provádí jen vzhledem k některým (nezávisle) proměnným, se nazývá semídískretízace. 68 Počáteční problém (3.58)-(3.60) jsme tedy získali semidiskretizací úlohy (3.54), (3.55) a (3.57) vzhledem k prostorové proměnné x. Úlohu (3.58) - (3.60) můžeme zapsat maticově ve tvaru Ců + Ku = F(ť), ÍG(0,T) u(0) = ^: (3.61) kde C je diagonální matice s kladnými prvky na diagonále, tzv. matice tepelné kapacity, K je třídiagonální pozitivně definitní matice známá jako matice tepelné vodivosti, F(t) je tzv. vektor tepelných zdrojů, u(t) = (u\(t), u2(t),..., UN-i(t))T je vektor neznámých funkcí a (p = ( 0 jsou kmity tlumené, pro c = 0 netlumené), p charakterizuje tuhost, q odpor okolního prostředí a / vnější síly. Dirichletovy okrajové podmínky předepisují deformaci, Newtonovy okrajové podmínky vnitřní síly: «o resp. on reprezentuje tuhost pružné podpory a fy resp. fy bodovou vnější sílu. Počáteční podmínky určují počáteční deformaci ip a její počáteční rychlost ip. Předpokládejme, že funkce p, c, p, q, f, g0, ge, fy, fy, ip a ip jsou dostatečně hladké, že jsou splněny podmínky nezápornosti p > 0, p > 0, g > 0, «o > 0, > 0 a že platí podmínky souhlasu ¥>(o) = 0o(o), ^(0) = ^(0), > > ditj(t) á2uAť) podmínku. Pro Ui(t) aproximující u{x,i,t) při označeni Ui(t) = ———, Ui(t) = ——— dostaneme soustavu rovnic PíUí(t) + Ciúi(t) + ^^-pj_i/2Mi-i(r) + [Pí-i/2 +PÍ+1/2 + h2qi]ui(t)- (3.71) - Pi+1/2Ui+1(tÝ) =fi(t), t = 1,2,..., N-1, t E (0, T) . Z okrajových podmínek (3.68) máme Mt) = 9o(t), uN(t) = ge(t), t G (0, T), (3.72) což dosadíme do (3.71) pro i = lai = N— 1. Z počátečních podmínek (3.70) dostaneme ut(0) = V(xt), ^(0) = ^), i = 0,l,...,JV-l. (3.73) (3.71) je soustava obyčejných diferenciálních rovnic druhého řádu pro N — 1 hledaných funkcí Ui(t), i = 1,2,... ,N — 1, s počátečními podmínkami (3.73). Počáteční úlohu (3.71)-(3.73) můžeme zapsat maticově ve tvaru Mú + Ců + Ku = F(ť), te(0,T), u(0) = (p, ů(0) = V, (3.74) kde M je diagonální matice s kladnými prvky na diagonále, tzv. matice hmotnosti, C je diagonální matice s nezápornými prvky na diagonále, tzv. matice útlumu, K je třídiagonální, pozitivně definitní-matice, tzv. matice tuhosti, F(t) je tzv. vektor (vnějšího) zatížení, cp = (if(x1),if(x2),...,(f(xN^1))T, if> = (^(xx),tfj(x2),...,il)(xN-1))T jsou vektory počátečních hodnot a u(t) = (ui(t), u2(t),..., UN-i(t))T je vektor neznámých funkcí. Časová diskretizace. Abychom mohli posoudit tuhost problému (3.74), zapíšeme ho jako počáteční problém prvního řádu, Rw + Sw = G(í), te(0,T), w(0) = k, (3.75) kde R-(ó m)--(u).^(k c).°W=U)).«=ft přičemž I je jednotková matice, O je nulová matice a o je nulový vektor. Prozkoumejme charakter počáteční úlohy (3.75) za zjednodušujících předpokladů. Pro p = p = 1, c = konst, q = f = 0, u(0, t) = u(í, ť) = 0, je M = I, K = — A a C = d, kde A je definována předpisem (3.63). Rovnice (3.75) se tak zjednoduší, dostaneme w = Pw, kde P=(° _c|V (3.76) Dá se ukázat, že vlastní čísla matice P 2N 7tí \ =-\c± ^\c2 - u2, kde ^ = ~sin^> « = l,2,...,iV-l. 72 Zřejmě maxj |Aj| —► oo pro N —► oo. Pro c = 0 jsou vlastní čísla ryze imaginární a pro c > 0 mají zápornou reálnou složku. Podle podle (1.32) je tedy problém (3.76) pro velké N tuhý. Protože reálné složky vlastních čísel jsou zdola ohraničené, Re(Aj) > —c, a velikost imaginárních složek je neomezená, říkáme, že problém (3.76) je oscílatorícky tuhy. Problém (3.75) se chová obdobně. Pro jeho řešení je proto vhodné použít metodu, jejíž oblast absolutní stability obsahuje celou zápornou komplexní polorovinu včetně imaginární osy. Teoretické úvahy i praktické zkušenosti potvrzují, že vhodnou metodou je lichoběžníková metoda. Z matlabovských programů lze doporučit programy ode23t, ode23tb a odel5s. Lichoběžníková metoda aplikovaná na úlohu (3.75) vede na předpis (R + irnS)w™+1 = (R - irnS)w™ + \rn{Gn + G™+1) , (3.77) viz (3.66). Pro efektivní výpočet složek un+1 a ůn+1 vektoru wn+1 je vhodné soustavu rovnic (3.77) zapsat po složkách, tj. u™+1 - \rniin+1 = un + \rniin , Můn+1 + ±rn[Kun+1 + Cůn+1] = Můn - ±rn[Kun + Cůn] + |rn(Fn+1 + Fn). Z prvé rovnice vyjádříme ůn+1, dosadíme do druhé a obdržíme rovnici pro výpočet un+1: Kun+1 = F, kde (3.78) K = 4tM + — C + K , F = f 4tM + —C — K j un + -Mu" + Fn + Fn+1. 2 ^i-----) \ 2 " Až vypočítáme un+1, dopočítáme ún+1 = — (un+1 - un) - ún . (3.79) Startujeme přitom z počátečních hodnot u° = (p , ú° = V • (3.80) Matice K soustavy (3.78) nezávisí na čase a je pozitivně definitní. Soustavu (3.78) je proto účelné řešit pomocí Choleského rozkladu, který stačí provést jen jednou. Metoda (3.78)-(3.80) je speciálním případem Newmarkovy metody hojně používané v inženýrské praxi, viz [2]. 3.4. Hyperbolická rovnice prvního řádu Formulace úlohy. Hledáme funkci u(x,t) definovanou pro x G (0,£), t G (0,T), která vyhovuje diferenciální rovnici ni i ni i — + a(x,t)— = f(x,t), xe(0,£),te(0,T), (3.81) 73 * G (O, T) (3.82) okrajovým podmínkám u(0,t) = g0(t), a(0,ŕ)>0 pokud u(£,t) = ge(t), a(£,t)<0, a počáteční podmínce u(x, 0) = ip{x) . (3.83) Nechť Q = (0,£) x (0, T) je obdélník, v němž hledáme řešení, a <9Q+ ta je část hranice dQ obdélníka Q, na níž je předepsána počáteční nebo okrajová podmínka, tj. dQ+ = {[x, t] G dQ 11 = 0 nebo x= 0 (pokud a(0, t) >0) nebo rr= £ (pokud a(£,t)<0)} . Pro každý bod Po = [^Oj^o] Q\ dQ+ určeme řešení x (t) počátečního problému dx(t) dt a(x(t),t) t < t o x (t0 X0. (3.84) du(x(t),t) Křivka x (t) se nazývá charakteristika příslušná bodu [rro,ŕo]- Předpokládejme, že charakteristika protíná dQ+ v jediném bodě P0* = [xq,Íq]- Tomuto bodu říkáme pata charakteristiky. Podle pravidla o derivování složené funkce du(x(t),ť) du(x(t),t) dx(t) du(x(t),t) du(x(t),t) du(x(t),t) ďt = di ďT+ dt = dt + ^> • Úlohu (3.81)-(3.83) proto můžeme na charakteristice x (t) zapsat ve tvaru f 0) , ip{£) = ge(0) (pokud a(£, 0) < 0). Pak úloha (3.85)-(3.84), a tedy rovněž úloha (3.81)-(3.83), má jediné klasické řešení. Rovnice (3.81) bývá označována jako advekční rovnice nebo také transportní rovnice. Úloha (3.81)-(3.83) popisuje třeba šíření příměsi prouděním: u(x,t) je koncentrace příměsi v proudícím médiu v bodu x a v čase t, a je rychlost proudění a / je zdrojový člen. Charakteristika x(t) je časoprostorová trajektorie, po které se částečka příměsi pohybuje. Pokud bychom připustili také difúzni šíření příměsi, dostali bychom o konvekčně-difúzni rovnici dQ+ 1 Obr. 3.9. Charakteristiky T ■t P2(x2M) A Pl(xi,h)r/ r P2*(x2' ^2) /Pľ(x*i,ťi) X 74 Ou O ô í ôu \ c(^)^ + ^(r(x,í)M)- — —J = f(x,t), x e (O,£), t g (O,T), v níž p je koeficient difúze a jejíž stacionární variantu jsme zavedli již dříve v kapitole 2.2, viz (2.20). Na rovnici advekce lze proto nahlížet jako na limitní případ konvekčně-difúzni rovnice, v níž zanedbáme účinky difúze. Jiná forma advekční rovnice tedy je c(x)— + —(r(x,t)u) = f(x,t), xe(0,£), íg(0,T). (3.81') Metoda přímek. Uvažujme rovnoměrné dělení intervalu (0, £), tj. Xi = ih, i = 0,1,..., N, kde h = £/N. Předpokládejme, že funkce a{x,t) v krajních bodech intervalu (0,£) nemění znaménko. Splnění rovnice (3.81) budeme požadovat ve vnitřních uzlech x1}x2,... ,xn_i a pro a(0, t) < 0 také v uzlu x0 a pro a(£, t) > 0 také v uzlu x^. V rovnici ^)+a(x,,ť)^)=/(x,,ř) aproximujeme člen ux(x,-nt) ve vnitřních uzlech centrální diferencí a v krajních uzlech jednostrannou diferencí. Pokud třeba a(0,t) > 0, a(£,t) > 0 pro všechna t g (0,T), dostaneme Ul(t) + ai(í)Mí) - g0(t)]/(2h) = hit), úl(t) + al(t)[ul+1(t)-ul_1(t)]/(2h) = f,(t), i = 2,3,...,N-l, (3.86) úN(t) + aN(t)[3uN(t) - 4ííjv-i(í) + uN_2(t)]/(2h) = fN(t) , kde a,i(t) = a(xi,t), fí(t) = f(x,nt) a Ui(t) je aproximace u(x,nt). Z (3.83) dostaneme počáteční podmínky Ul(0) = Vl, i = 1,2,...,N, (3.87) kde ifi = ip(xi). Počáteční problém (3.86)-(3.87) řešíme vhodnou numerickou metodou. K jejímu výběru nám poslouží zjednodušený model. Předpokládejme, že a = 1, f = 0, u(0,t) = 0 a uvažme periodické okrajové podmínky u(0,t) = u(£,t). Pomocí centrálních diferencí odvodíme ůiit) + ai(t)[u2(t) - uN(t)]/(2h) = hit), úi(t) + ai(t)[ui+1(t)-ui_1(t)]/(2h) = fi(t), i = 2,3,...,N-l, (3.86') úN(t) + aN(t)[Ul(t) - uN-1(t)]/(2h) = fN(t). Počáteční problém (3.86')-(3.87) je pak tvaru ů = Au, íg(0,T), u(0) = (p, (3.88) 75 kde u(t) (Ul(t),u2(t),...,uN(t))T, tp(t) /o 1 O ... O -1 O 1 ... o 0-10 ... o A 1 "2h 0 1 O O O o o -1 (Vl(*),V?2(í), -1 \ o 0 1 o ., fN(t)) a (3.89) Dá se ukázat, že vlastní čísla matice A Aí = /—srn—, 7 = V-1, í 1,2,...,iV, leží na imaginární ose mezi —N/í a N/í a maxj |Aj| —► oo pro A?" —► oo, tj. pro velké N je počáteční problém (3.86')-(3.87) tuhý. Problém (3.86)-(3.87) se chová podobně: vlastní čísla odpovídající matice A leží v záporné komplexní polorovině v blízkosti imaginární osy, největší vlastní číslo \max leží prakticky na imaginární ose, \max = IN/i, I = \f—\. Pro řešení počátečního problému (3.88)-(3.89) je proto třeba zvolit metodu, jejíž oblast absolutní stability obsahuje úsečku (—a, a) na imaginární ose pro nějaké a > 0. Vhodná je třeba A-stabilní TR metoda, pro kterou a = oo, takže délka časového kroku není z důvodu stability nijak omezována a řídí se jen požadavkem na přesnost. Použít lze ale také některé další metody. Třeba metodu BS32, pro kterou a = 1,72 nebo DP54, pro kterou a = 1, viz např. [12]. V případě a < oo bude délka časového kroku z důvodu stability omezena, toto omezení však nebude nijak dramatické, délka kroku bude řádu 0(h). Z matlabovských programů lze doporučit programy ode23t, ode23 a ode45. Stejné programy jsou vhodné i pro řešení úlohy (3.86)-(3.87). Metoda charakteristik je další vhodnou technikou pro řešení úlohy (3.81)-(3.83). Její podstatu vysvětlíme nejdříve pro zjednodušenou úlohu u(0,t) = g0(t), t E (0,T), u(x, 0) = 0 je konstanta. Diferenciální rovnici přepíšeme pomocí charakteristik na tvar du(x(t),t) dt 0. Zvolme rovnoměrné dělení intervalu (0,£) s krokem h = £/N, tj. Xi = ih, i a rovnoměrné dělení intervalu (0,T) s krokem r = T/Q, tj. tn = nr, n 0,1,. 0,1, (3.90) ..,iV, ..,Q. Charakteristika vycházející z bodu [:Tj,tn+1] je přímka x(t) = Xi + a(t — tn+í). V čase x(tr Xj ar. 76 Předpokládejme, že ar < h. Pak pro i = 1,2,... ,N bod xn G {xí_i,Xí), zejména tedy xn G (0,£), takže u(xn,tn) má smysl. Integrací rovnice (3.90) od tn do tn+1 obdržíme tn+l du(xi(t),t) dt = u(xí, tn+i) — u(x™, tn) = 0 , a odtud u(xí, £n+i) = u(x™, tn). Numerickou metodu dostaneme tak, že u(xn,tn) určíme přibližně interpolací, tj. počítáme u?+1 kde Pk{x) je interpolační polynom stupně k > 1 splňující Pk(Xi-l) = <_! , Pk(Xi) = < . Pro lineární polynom Px z podmínek (3.92) odvodíme Pi(rr) =< + (3.91) (3.92) ar /i — (x-Xi) a odtud Pi(x") = - - u?^) ln+l Dostali jsme tak upwind metodu x(t)/ J u n+1 ar (3.93) i i i i i i _i_i_ Xí—l xn X;L Xí+i Obr. 3.10. Upwind metoda Dá se ukázat, že když ar Název metody nám připomíná, že veličina u v uzlu Xi závisí jen hodnotách této veličiny x v uzlech ležících „proti proudu, proti větru", pro a > 0 tedy nalevo od Xi. V bodu xq = 0 užijeme okrajovou podmínku, takže urQl+1 = go(tn+i). (3.94) pak za předpokladu dostatečné hladkosti přesného řešení pro chybu platí u(xi,tn) - u? = 0(t) , (3.95) viz [14]. Říkáme, že metoda (3.93) je řádu jedna. Pro v = 1 dokonce un = u(xi,tn) je přesné! Číslo v = ar/h se nazývá Courantovo číslo a podmínka (3.94) se nazývá Courantova-Fríedríchsova-Lewyova podmínka, stručně CFL podmínka. Metodu řádu dva dostaneme tak, že v (3.91) použijeme interpolační polynom druhého stupně. Když k podmínkám (3.92) přidáme ještě podmínku P2(Xi+1) = <+1, (3.96) vypočteme P2{xn) a dosadíme do (3.91), dostaneme Laxovu-Wendroffovu metodu u. n+l _ n ar , ,„ „ , a2r2 (3.97) 77 Pro aproximaci v uzlu x^ můžeme použít interpolační polynom, který kromě podmínek (3.92) vyžaduje navíc splnění podmínky P2(xn_2) = mtv-2- Jestliže pro obecný uzel Xi přidáme k podmínkám (3.92) podmínku P2(rri_2) = <_2 , (3.98) vypočteme P2(xn) a dosadíme do (3.91), obdržíme Beamovu-Warmingovu metodu flT fl2T2 <+1 = < - ^(3< - 4<_x + <_2) - 2^2-« - 2<-i + <-2) • (3-99) Počítáme-li u™+1 podle (3.97) a podle (3.99), je-li splněna CFL podmínka (3.94) a je-li přesné řešení u dostatečně hladké, pro chybu platí u(xhtn) - u? = 0(>r). (3.100) Beamova-Warmingova metoda je upwind aproximací řádu 2, řešení v uzlu Xi závisí jen na „protiproudových" hodnotách v uzlech Věnujme se dále případu, kdy konstanta a < 0. Pak okrajová podmínka bude předepsána vpravo, tj. podmínku (3.892) nahradí podmínka u(£,t) = ge(t), te(0,T). (3.892) Charakteristika vycházející z bodu [xí, tn+i] bude směřovat doprava, takže za předpokladu platnosti CFL podmínky \a\r i/ = -LT-<1 3.101 h bude Xi(tn) G (xí,Xí+i) pro i = 0,1,..., N — 1. Numerické metody odvodíme opět ze vztahu (3.91), tentokrát však místo (3.92) požadujeme Pk(xt)=u?, Pk{xl+1)=unl+1. (3.92') Podobně jako dříve dostaneme upwind metodu <+1 = <-^W-<+i) (3.93') a Beamovu-Warmingovu metodu <+1 = < - ^(3< - 4<+1 + <+2) - ^«+2 - 2<+1 + O , (3.99') Laxova-Wendroffova metoda zůstane stejná, tj. postupujeme podle (3.97). V bodu x^ užijeme okrajovou podmínku, takže u^1 = g(tn+i). V obecném případě a = a(x,t) určíme xn numericky: xn ~ x(tn), kde = a(x(t),t), x{tn+l) = Xl. (3.102) dt 78 Pro upwind metodu použijeme EE metodu, takže x'1 = Xi — a™ t, kde a™ = a(x,n tn+1). Pro Laxovu-Wendroffovu metodu a Beamovu-Warmingovu použijeme EM2 metodu, tj. x'1 = Xi — a™t, kde a™ = \(ki + k2), k\ = a(xi,tn+1) , k2 = a(x,i — rki,tn). Vzorce (3.93), (3.93'), (3.97), (3.99) a (3.99') se změní jen v tom, že v nich místo a píšeme a™. Délka kroku rn = tn+1 — tn musí splňovat CFL podmínku h < 1. (3.101') Jestliže v rovnici (3.89i) uvažujeme nenulovou pravou stranu f(x, t), tj. řešíme-li místo rovnice (3.90) rovnici = Mt),t), (3.90-) přičteme k pravým stranám formulí aproximaci Q"(f) integrálu J^n+1 f(x(t),t) dt. Pro upwind metodu stačí použít jednostrannou obdélníkovou formuli QUf)=rf(xt,tn+1). Pro Laxovu-Wendroffovu metodu a Beamovu-Warmingovu metodu užijeme lichoběžníkovou formuli QUf) = 12r[f(xt,tn+1) + f(x1l,tn)}. Metoda konečných objemů. Pro diskretizaci rovnice (3.81') se výborně hodí metoda konečných objemů, na časoprostorových objemech E>i = (2^-1/2,^1+1/2) x (tn,tn+i) , více o tom viz [13]. 79 Literatura [1] R. Ashino, M.Nagase, R. Vaillancourt: A survey of the MATLAB ODE suite, Technical Report CRM-2651, Centre de recherches mathematiques, University of Ottawa, 2000. [2] K. J. Bathe: Finite Elements Procedures, Prentice-Hall, Upper Saddle River, NJ, 1996. [3] M. Brandner, J. Egermaier, H. Kopincová: Numerické metody pro řešení evolučních parciálních diferenciálních rovnic, učební text ZCU a VSB, Plzeň, 2012. [4] L. Čermák: Numerické metody II. Diferenciální rovnice, Akademické nakladatelství CERM, Brno, 2010. [5] L. Čermák, R. Hlavička: Numerické metody, CERM, učební text FSI VUT Brno, 2008. [6] D. R. Durran: Numerical Methods for Fluid Dynamics: with Application to Geophysics (2nd edition), Springer, New York, 2010. [7] J. H. Ferziger, M. Perič: Computational Methods for Fluid Dynamics (3rd edition), Springer, Berlin, 2002. [8] J. Fish, T. Belytschko: A First Course in Finite Elements, John Willey & Sons Ltd, Chichester, 2007. [9] M. T. Heath: Scientific Computing. An Introductory Survey, McGraw-Hill, New York, 2002. [10] R. Hlavička: Numerické metody při řešení diferenciálních rovnic: Průvodce softwarem a počítačová cvičení v prostředí MATLABu, učební text FSI VUT Brno. [11] P. Knabner, L. Angermann: Numerical Methods for Elliptic and Parabolic Partial Differential Equations, Springer, New York, 2003. [12] J. D. Lambert: Numerical Methods in Ordinary Differential Systems. The Initial Value Problem, J. Willey & Sons, Chichester, 1993. [13] R. J. LeVecque: Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, Cambridge, 2002. [14] R. J. LeVecque: Finite Difference Methods for Ordinary and Partial Differential Equations, Siam, Philadelphia, 2007. [15] The Math Works, Inc.: MATLAB® Mathematics, R2012b. http://www.mathworks.com/help/pdf_doc/matlab/math.pdf. [16] S.Mika, P. Přikryl, M. Brandner: Speciální numerické metody: Numerické metody řešení okrajových úloh pro diferenciální rovnice, Vydavatelský servis, Plzeň, 2006. [17] C.B.Moler: Numerical Computing with MATLAB, Siam, Philadelphia, 2004. Electronic edition: The MathWorks, Inc., Natick, MA, 2004. http://www.mathworks.com/moler. [18] S.Mika, P. Přikryl: Numerické metody řešení obyčejných diferenciálních rovnic, okrajové úlohy, učební text ZCU Plzeň, 1994. [19] A. Quarteroni, R. Sacco, F. Saleri: Numerical Mathematics, Springer, Berlin, 2000. 80 [20] K.Rektorys: Přehled užité matematiky 1,11, Prometheus, Praha, 1995. [21] L.F. Shampine: Some practical Runge-Kutta formulas, Math. Comp., Vol. 46, Num. 173 (1986), str. 135-150. [22] L.F. Shampine: Numerical Solution of ordinary differential equations, Chapman & Hall, New York, 1994. [23] L.F. Shampine, I. Gladwell, S. Thompson: Solving ODEs with MATLAB, Cambridge University Press, Cambridge, 2003. [24] L.F. Shampine, L. W. Reichelt: The MATLAB ODE suite, SIAM J. Sci. Comput., Vol. 18 (1997), No. 1, str. 1-22. [25] H. K. Versteeg, W. Malalasekera: An Introduction to Computational Fluid Dynamics: The finite volume method (2nd edition), Prentice Hall, Harlow, 2007. [26] E. Vitásek: Numerické metody, SNTL, Praha, 1987. [27] E. Vitásek: Základy teorie numerických metod pro řešení diferenciálních rovnic, Academia, Praha,1994. [28] O. C. Zienkiewicz, R.L.Taylor: The Finite Element Method, Volume I: The Basis, Butter -worth-Heinemann, Oxford, 2000. 81