Numerická matematika 1 Parabolické rovnice Budeme se zabývat rovnicí du d2u , N äi = Där* (1) tato rovnice určuje chování funkce u(t, x), která závisí na dvou proměnných. První proměnná t mívá význam času, druhá x bývá prostorová souřadnice. Předchozí rovnice lze snadno zobecnit pro více prostorových proměnných např. du nf92u d2u\ ~dt " + ch/2) ( ' Příkladem parabolické rovnice je rovnice vedení tepla d /, dT\ dT cTx [kdx-j = PC9ľ (3) tato rovnice popisuje časový vývoj teploty T(t, x) uvnitř nekonečné stěny, k je tepelná vodivost materiálu, c je tepelná vodivost a p je hustota. Jiným příkladem je rovnice difúze pro hustotu částic n(t,x) dn 1 dn ~ ~k?~ďt = ( ' kde k je konstanta popisující rychlost difúze. Parabolická rovnice nám vlastně říká, že tam kde je druhá derivace záporná, hledaná funkce v čase klesá a naopak. Dochází tak k vyhlazování průběhu funkce - teploty se vyrovnávají, hustoty částic v různých bodech se srovnávají. Parabolické rovnice většinou řešíme na omezeném intervalu. Budeme předpokládat, že se jedná o interval x G (0, Ľ). Např. při řešení rovnice vedení tepla, bude x — 0 reprezentovat vnitřní povrch stěny a x — L vnější povrch stěny. Abychom mohli rovnici řešit je třeba znát hodnoty hledané funkce v počátečním čase t — 0 tzv. počáteční podmínku. Například je třeba zadat počáteční hodnoty ve stěně. Počáteční podmínka má obecně tvar u(0, x) — p(x) (5) kde p{x) je známá funkce. Zadání počáteční podmínky však pro výpočet nedostačuje. Je třeba také vědět co se odehrává na okrajích studované oblasti. Je tedy třeba zadat chování funkce pro x — 0 a x — L, tomuto říkáme okrajové podmínky. Např. je třeba zadat teplotu na vnitřním a vnějším povrchu stěny. Okrajová podmínka může mít dva základní tvary. První možností je přímo zadání hodnot na hranici oblasti tj. u(t,0) = ff0(t) nebo u(t,L)=gi(t) (6) takovéto podmínce se říká Dirichletova podmínka. Druhou možností je zadání prostorové derivace funkce tj. du du — (t,0) = <72(t) nebo — (í, L) = g3(t) (7) takovéto podmínce se říká Neumannova podmínka. Například u rovnice vedení tepla této podmínce odpovídá zadání tepelného toku, např. nulová derivace odpovídá dokonale izolovanému povrchu stěny. Naším cílem tedy bude nalezení funkce u(t,x) pro x G (0, L) a í G (0,oo), která splňuje rovnici (1), počáteční podmínku (5) a okrajové podmínky (6) či (7). Numerická matematika 2 Příklad Řešme diferenciální rovnici du jj^u dt dx2 na intervalu (0, L) s okrajovými podmínkami u(í,0) = 0 u(t,L) = l a s počáteční podmínkou u(0, x) — x + sin \ j (8) (9) (10) Snadno se přesvědčíme, že počáteční a okrajové podmínky jsou v souladu tj. podle obojího je w(0,0) — 0 a w(0,L) = 1. Dosazením se lze přesvědčit, že řešením dané rovnice s těmito podmínkami je u(t, x) — x + sin y — J e -c-2 Pro T —> oo se toto řešení plíží k ustálenému stavu it(oo, x) — X Metoda sítí pro parabolické rovnice (11) (12) Ukážeme si základní metodu na numerické řešení PDR, tzv. metodu sítí přesněji zvanou metodu konečných diferencí. V této metodě se místo spojité funkce u(t,x) hledají pouze odhady řešení v konečném počtu bodů. Tyto body tvoří v oblasti řešení síť, odtud je název metody. Existuje celá řada druhů sítě, my se omezíme na nejjednodušší případ pravoúhlé rovnoměrné sítě. Body této sítě mají souřadnice [í",Xj]. Časové okamžiky tn jsou rovnoměrně rozmístěny s časovým krokem t a prostorové body jsou rovnoměrně rozmístěny s krokem h tj. Xj — j h J = 0,1, n = 0, 1, h = L/M (13) (14) tnO- ť i-o tu6ôôôôôôôôg>g>g>6 o o o o o o Xq Xj X|_[ X, X.4.1 Xm x Obrázek 1: Konstrukce rovnoměrné ortogonální sítě. Místo spojité funkce u(t, x) tedy budeme hledat odhady tohoto řešení w™, tak aby v ideálním případě platilo uf — u(tn,Xi). Tento vztah nebude platit přesně, protože Numerická matematika 3 numerická metoda bude opět zatížena chybou. Místo diferenciálních rovnic pro u(t, x) je třeba získat tzv. diferenční rovnici pro uf. Tuto diferenční rovnici získáme tak, že parciální derivace v diferenciální rovnici nahradíme některým přibližným vzorcem. Explicitní metoda V nejjednodušší tzv. explicitní metodě použijeme následující přibližné vzorce tíľ+1 - ur- du ~dt = t d2u ^ <+1-2<+<_! dx2 h2 Dosazením do (1) dostaneme diferenční rovnici pro uf i?+i-2<+<_1 u?+1 - u? D- h2 Z této rovnice můžeme vyjádřit u n+l n+l (15) (16) (17) (18) Pokud známe řešení v časový okamžik tn, můžeme z této rovnice vypočítat řešení v následující časový okamžik fn+1. Výpočet začneme v čase t — to — 0, kde máme zadáno řešení počáteční podmínkou. Poté postupujeme časem a ze vztahu (18) vypočítáme řešení v čase t1, t2, .... Vztah (18) však nelze použít v krajních hodnotách, protože například k výpočtu Mq+1 potřebujeme znát hodnotu která ale leží mimo studovanou oblast. Pro výpočet krajních hodnot Mq+1 a u™+1 je tedy třeba použít okrajových podmínek. L n+l ♦-» X|-l *i xi4l explicitní u n+l Xj_] X, Xj4] inplicitni Xi-l X| xi+] Crank-Nicholsonovo Obrázek 2: Schemata pro řešení parabolické rovnice. Explicitní metoda je celkem jednoduchá, má ale jedno podstatné omezení. Pokud je časový krok příliš velký metoda není tzv. stabilní To se projeví po několika krocích metody, kdy se objeví na průběhu funkce zákmity, které se postupně zesilují, až řešení diverguje. Takové chování je samozřejmě nežádoucí. Aby toto chování nenastalo tj. aby metoda byla stabilní, je třeba, aby byla splněna následující podmínka Takováto podmínka, která zajišťuje stabilitu metody, se nazývá podmínka stabilita. Metody se z tohoto hlediska dělí na tři druhy - metody stabilní, nestabilní Numerická matematika 4 a podmíněně stabilní. Metoda se nazývá podmíněně stabilní pokud potřebuje k zajištění stability nějakou podmínku, jako výše uvedená metoda. Podmínka (19) omezuje velikost časového kroku. Tato podmínka je bohužel dost přísná. Pokud je například tato podmínka přesně splněna a my chceme zmenšit prostorový krok na polovinu je třeba zmenšit časový krok čtyřikrát, v důsledku toho vzroste časová náročnost osmkrát. Řešení viz program pdr-parabol-expl. f 90. Implicitní metoda Implicitní metoda vylepšuje explicitní metodu tím, že odstraňuje podmínku stability. Je tedy stabilní, nikoli jen podmíněně stabilní. Tato metoda používá jinou aproximaci druhé prostorové derivace než je (16) a to dx2 -2u n+1 ,n+l h2 (20) Použijeme tedy stejný derivační vzorec, ale vztažený k jinému časovému okamžiku. Dostaneme tak tíľ+1 - uv = D .n+1 - 2u n+1 ,n+l h2 (21) Z této rovnice nemůžeme přímo vypočítat řešení v novém čase tn+1 protože k výpočtu u™+1 je třeba již znát řešení v sousedních bodech w™^1 a . Rovnice (21) představuje tedy vztah mezi třemi hledanými hodnotami v novém čase tn+1 ve třech sousedních bodech. z předchozího vztahu dostaneme Označme si a — ^ u n+1 u n a«+1 ~3 - -v^j+1 -ait"+11 + (l + 2aK+1 2^+1 n+1 (22) (23) Pro různá i máme různé rovnice, dohromady tvoří tyto rovnice soustavu lineárních rovnic pro neznáme V maticovém zápisu 1 - 2a -a 1 -a -2a -a -a - 2a 1 + 2a —a —a 1 + 2a u0 n+1 -i 0 n+1 1 n+1 2 u0 n+1 J-1 n+1 J J . -ľ (24) Řešením této soustavy dostaneme řešení v novém časovém okamžiku tn+1. Crankovo-Nicholsonovo schéma Předchozí implicitní metoda je vždy stabilní. Z hlediska stability lze v této metodě volit libovolně velký krok. Pokud však zvolíme větší krok diskretizační chyba se zvětší. Chyba implicitní i explicitní metody je 0(r,h2). Chyba je tedy úměrná pouze první mocnině časového kroku, to je způsobeno tím, že obě metody jsou nesymetrické v čase. Tuto asymetrii odstraňuje Crankova-Nicholsonova metoda. Tato metoda místo prostorové derivace v čase tn (jako explicitní metoda) nebo v čase tn+1 (jako implicitní metoda) používá průměr z obou možností tj. Numerická matematika 5 neboli <+1 - = f [K+í - K+1 + + K+i - + (26) Z této rovnice opět nelze přímo vypočítat řešení v novém čase, rovnici tedy nejprve upravíme -\u^\ + (1 + a)u^ = u] + |(^+1 - 2M™ + u]_,) (27) Máme tedy opět soustavu lineárních rovnic pro řešení v čase tn+1. V maticovém zápisu je tato soustava "1+a -§ -f 1+« -f -f 1 + a Crankova-Nicholsonova metoda je nepodmíněně stabilní stejně jako implicitní metoda. Její chyba je ale menší je řádu 0(t2, ti2). Tato metoda je téměř stejně časově náročná jako implicitní metoda, je ale přesnější, a proto bychom jí měli dávat přednost. Řešení viz program pdr-parabol-impl. f 90. „,"+! Ul „,"+! — u2 «ô + !K-2«ô) (28)