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
i/(t) = f(t,y(t)) (1-1) a splňuje počáteční podmínku
y(a)=r). (1.2)
Je-li v nějakém okolí D bodu [a,r]] funkce f(t,y) spojitá a splňuje-li v tomto okolí Lip-schitzovu podmínku s konstantou L vzhledem k proměnné y, tj. platí-li
\f(t,u)-f(t,v)\<L\u-v\   V [t, u] ,[t,v] e D , (1.3)
pak bodem [a,r]] prochází jediné řešení y(t) rovnice (1.1). Jestliže funkce f(t,y) má v D omezenou parciální derivaci vzhledem k proměnné y, tj. když \dyf(t, y) \ < L, Lipschitzova podmínka (1.3) platí. Jiný standardní výsledek říká, že když je funkce f(t,y) spojitá na (a,b) x IR a splňuje tam Lipschitzovu podmínku, pak počáteční problém (1.1), (1.2) má jediné řešení definované v celém intervalu (a, b).
Počáteční problém pro soustavu ODR1 znamená určit funkce yi(t),..., y<i(f) splňující diferenciální rovnice
y'j{t) = fj{t,y1{t),...,yd{t)),   j = l,2,...,d, a počáteční podmínky
yj(a)=Vj,   j = 1,2,..., ti. Stručnější vektorový zápis je
y'(ŕ) = f(ŕ,y(ŕ)),  y(a) = r,, (1.4)
kde
y(í) = (yi(ŕ),y2(ŕ),...,y,(ŕ))T,     y\t) = (y[(t),y'2(t),... ,y'd(t))T,
f (t, y(t)) = (A (ŕ, y(í)), /2(í, y(í)),..., fd(t, y(ŕ)))T,      v = (Vl, V2,..., Vd)T.
Je-li v okolí D bodu [a,r]] funkce f (ŕ, y) spojitá a splňuje tam vzhledem k proměnné y Lipschitzovu podmínku s konstantou L , tj. platí-li
||f(í,u)-f(í,v)|| <L||u-v||   V[í,u], [ŕ, v] G D, (1.5)
5
pak bodem [a, rj] prochází jediné řešení počáteční úlohy (1.4). Má-li f v D spojité a omezené parciální derivace {dfi(t,y)/dyj}fj=1, pak Lipschitzova podmínka (1.5) platí. Je-li D = {a, b) x Rd, jediné řešení existuje v celém intervalu (a, b).
Rovnice vyššího řádu. Počáteční problém pro obyčejnou diferenciální rovnici řádu d,
y^(t) = F(t,y(t),y'(t),...,y(d-1\t)) (1.6) s počátečními podmínkami
y{°) =Vu   y'{a) =V2,---,y[d~1\a) = r)d, lze snadno převést na počáteční problém (1.4) pro d rovnic řádu prvního:
y[(t) = y2(t), y1(a)=r]1,
2/2 (ř) =Vs(t), ž/2(a) =??2,
y'd-i{t) = yd(t), Vd-i(a) = Vd-i ,
y'dit) = F(t, yi(t),y2(t),..., yd(t)),      yd(a) = Vd ,
kde yi(t) = y(t), y2(t) = y'(t), ..., ydit) = y^{ť).
V dalším budeme vždy předpokládat, že uvažovaná počáteční úloha má v intervalu (a, b) jediné řešení. Budeme také předpokládat, že funkce f (ŕ, y) má tolik spojitých derivací, kolik jich v dané situaci bude zapotřebí.
Jedna rovnice prvního řádu s jednou neznámou funkcí je v aplikacích méně významná než soustavy rovnic. Metody přibližného řešení se však snadněji odvodí pro jednu rovnici a lze je aplikovat bezprostředně i na soustavy. Také analýza numerických metod je pro jednu rovnici podstatně snadnější. Proto se v následujícím výkladu převážně omezíme jen na jednu rovnici. Z velkého množství metod uvedeme ty, které jsou pro své dobré vlastnosti široce používány. Mezi ně bezesporu patří metody implementované do Matlabu a právě na ně se v tomto textu zaměříme.
Numerickým řešením počáteční úlohy rozumíme výpočet přibližných hodnot hledaného řešení yit) v bodech tn dosti hustě vykrývajících interval (a,b). Nechť tedy
a = t0 < ri < • • • < tQ = b
je dělení intervalu (a, b). Body tn jsou uzly, vzdálenost rn = tn+i —tn dvou sousedních uzlů je délka kroku Jsou-li všechny kroky stejně dlouhé, tj. když rn = r = (b — a)/Q, hovoříme o rovnoměrném (ekvidistantním) dělení intervalu (a,b). V tom případě je tn = a + nr, n = 0,1, ..., Q. Hodnotu přesného řešení v uzlu tn budeme značit yitn) a hodnotu přibližného řešení yn. Jestliže se nám podaří najít přibližné řešení yn, n = 0,1, ..., Q, můžeme vypočítat přibližnou hodnotu řešení y(t) v libovolném bodě t E {a, b) interpolací.
Numerická metoda pro řešení počáteční úlohy (1.1), (1.2) je předpis pro postupný výpočet aproximací yi, y2, ..., yg, y o = 77 z počáteční podmínky. Metoda se nazývá
6
k-kroková, závisí-li předpis pro výpočet aproximace yn+i na předchozích aproximacích yn, yn-i, ■ ■ ■, Vn-k+i- Speciálně jednokroková metoda počítá přibližné řešení yn+i v uzlu tn+i jen pomocí znalosti přibližného řešení yn v uzlu tn, přibližná řešení yn-i, yn-2, ■ ■ ■ spočtená v předchozích uzlech ŕn_i, rn_2, • • • nepoužívá. Výpočet yn+i nazýváme krokem metody od tn do tn+i (stručně krokem). Při popisu kroku budeme u délky kroku rn = tn+i — tn vypouštět index, tj. píšeme rn = t.
1.2. Eulerovy metody
Explicitní Eulerova metoda. Nejjednodušší numerickou metodou pro řešení úlohy (1.1), (1-2) je explicitní Eulerova metoda (stručně EE metoda). EE metodu snadno odvodíme z Taylorovy formule
y(tn+1) = y{tn + t) = y{tn) + ry\tn) + \t2y"{Í-,n),   £„ G (tn, tn+1). (1.7)
Uvážíme-li, že y'(tn) = f(tn,y(tnj) a zanedbáme-li člen \r2y"{^n), dostaneme
y(tn+1) w y{tn) + t/(í„, y{tn)).
Výrazy yitn) a y(tn+i) nahradíme jejich přibližnými hodnotami yn a yn+i, znaménko přibližné rovnosti    nahradíme znaménkem rovnosti a obdržíme předpis EE metody
yn+1=yn + Tf(tn,yn). (1.8)
O explicitní metodě hovoříme proto, že pro určení yn+\ máme explicitní vzorec: dosazením známé hodnoty yn do pravé strany rovnice (1.8) obdržíme hledanou hodnotu yn+i. V anglicky psané literatuře se EE metoda označuje jako Euler method resp. explicit Euler method resp. forward Euler method.
Diskretizační chyby. Přesnost numerické metody měříme pomocí tzv. lokální diskreti-zační chyby (anglicky local truncation error)
lte„ = y{tn+1) - y(tn) - rf(tn, y(tn)).
Lokální diskretizační chyba je tedy chyba, které se dopustíme v jednom kroku metody za tzv. lokalizačního předpokladu, že yn = yitn) je přesné řešení počáteční úlohy (1.1), (1.2). Z (1.7) pro EE metodu plyne
lte„ = \t2y"{^n) a tedy  |lte„| < CV2,   kde C = \   max   \y"(t) \ ,
ín<í<ín + l
což lze stručně vyjádřit zápisem lte„ = Oir2). Následující poznámka připomíná význam Landauova symbolu 0(tp).
Poznámka. Nechť <p(t) je funkce definovaná v intervalu (0, r0) a jo je libovolné číslo. Řekneme, že funkce <p(t) je řádu 0(tp) a píšeme <p(t) = 0(tp), jestliže existuje kladné číslo C takové, že pro všechna 0 < r < r0 platí |<^(r)| < Ctp. □
Lokální diskretizační chyba při reálném výpočtu nevzniká, neboť obecně není splněn lokalizační předpoklad, tj. yn ^ y(tn). Lokální diskretizační chyba se uplatní jen při
7
analýze vlastností numerické metody, například konvergence yn —> yitn). Pro praktické účely, například pro řízení délky kroku, je nutné pracovat s tzv. lokální chybou (anglicky local error) definovanou předpisem
le„ = un(tn+1) - yn+1,
kde un(t) je tzv. lokální řešení počátečního problému
= f{t, un(tj) ,       un(tn) = yn .
Lokální chyba le„ je tedy chyba, které se skutečně dopustíme při reálném výpočtu v kroku od tn do tn+i. Dá se ukázat, že pro výpočet s dostatečně malými délkami kroků je rozdíl mezi oběma lokálními chybami prakticky zanedbatelný.
Obr. 1.1.  Diskretizační chyby
Hromaděním lokálních chyb vzniká globální áiskretizační chyba
e„ = y(t„) -y„. V případě rovnoměrného dělení lze dokázat, že
M = \y(tn)-yn\ <Ct,       n = 0,l,...,Q, (1.9)
kde C je konstanta nezávislá na r = (b — a)/Q. Tuto skutečnost stručně vyjádříme tvrzením, že globální áiskretizační chyba EE metoáy je řááu Oir). Říkáme také, že EE metoáa je řááu 1. Protože en —> 0 pro Q —> oo, numerické řešení získané EE metodou konverguje k řešení přesnému. Říkáme také, že rychlost (řáá) konvergence EE metoáy je rovna 1. Lokální chyby lte„, le„ a globální chyba era+1 jsou zakresleny v obrázku 1.1.
8
Tvrzení (1.9) lze snadno ověřit v případě, že f(t, y) = fit) nezávisí na y. Pak totiž
y(tk+1) = y(tk) + rf(tk) + ir2/(6), yk+1 = yk + rf(tk).
Odečtením obou rovnic dostaneme ek+\ = ek + \t2f'(^k), a protože 6q = 0, je
Označíme-li M = maxQ<g<č, |/'(£)|, pak
\en\ < \r2nM < \ [tQ\Mt = \{b- a)Mr,    neboť tQ = b - a.
Zaokrouhlovací chyby. Dopustíme-li se v kroku od tn do tn+i zaokrouhlovací chyby, jejíž velikost |Ayra+1| = \ýn+i — yn+i\ nepřesáhne e, pak lze dokázat, že po Q krocích délky r velikost zaokrouhlovací chyby nepřesáhne Ker-1, kde K je konstanta nezávislá na e a t. Pro celkovou chybu EE metody pak platí
max \y(tn) - yn\ <Ct + eKr~ľ,
0<n<Q
kde C, K jsou konstanty nezávislé na r a e. Protože konstanta e je malá, vliv zaokrouhlování se projeví až pro extrémně velký počet kroků Q (tj. pro velmi malé r). Tato situace při řešení běžných úloh nenastává a tudíž vliv zaokrouhlovacích chyb bývá nepodstatný.
Stabilita. Řekneme, že počáteční problém (1.1), (1.2) je stabilní vzhledem k počáteční podmínce, jestliže malá změna počáteční hodnoty 77 vyvolá malou změnu řešení. Elementárním příkladem takového problému je testovací úloha
y' = Xy,   y(0) = l, (1.10)
kde A = a + í/3 je komplexní číslo se zápornou reálnou složkou, tj. Re(A) = a < 0. Skutečně, jestliže
uit) = Xuit),      u(0) = 1 ==>   u(t) = ext
v'(t) = Xv(t),      v(0) = l + S v(t) = (1 + S)eXt,
pak
\u(t)-v(t)\ < \ó\ ■ |e(Q+l/3)í| = |5|eQÍ • | cosf3t +1 sin/3í| = |5|eQÍ < \b\ .
Pro řešení y(t) = eXt testovací úlohy (1.10) rovněž platí
\y(t)\ = \ext\ = eat \ cos f3t +i sin f3t\ -> 0   pro   t -> 00 .
Je proto přirozené požadovat, aby na rovnoměrném dělení tn = nr, n = 0, 1, .. ., numerické řešení yn testovací úlohy (1.10) splňovalo analogickou relaci, tzv. podmínku stability
\yn\ ~^ 0   pro   n —> 00. (1-11)
9
Aplikujeme-li EE metodu na testovací rovnici (1.10), dostaneme
yn+1 =y„ + r\yn = (1 + r\)yn = (1 + rA)2y„_i = • • • = (1 + rX)n+1y0 .
Podmínka stability (1.11) bude splněna, právě když |1 + t\\ < 1, neboli když r A leží v tzv. oblasti absolutní stability Ra-
t A E Ra = {z E C| \z + 1| < 1}.
Oblast absolutní stability EE metody je tedy vnitřek jednotkového kruhu |z + l| < 1 komplexní roviny C se středem v bodě [—1,0]. Průnik oblasti absolutní stability se zápornou částí reálné osy je interval absolutní stability Ia- Pro EE metodu Ia = (—2, 0). Pro reálné A < 0 podmínka stability (1.11) vyžaduje volit krok r < 2/|A|.
Tvar a velikost oblasti absolutní stability metody je spolu s řádem metody základní charakteristikou kvality numerické metody. EE metoda z tohoto pohledu příliš kvalitní není: je pouze řádu 1 a oblast její absolutní stability je malá. EE metoda se používá jen výjimečně.
Implicitní Eulerova metoda. Vyjdeme opět z Taylorova rozvoje
y(tn) = y(tn+i) ~ Ty'{tn+1) + ±T2y"(£n),   £„ G (tn, yn+1). (1.12)
Vypuštěním členu \y"{^n) a užitím rovnosti y'(tn+i) = f(tn+i,y(tn+i)) obdržíme implicitní Eulerovu metodu (stručně IE metodu) jako předpis
Vn+i =Vn + rf(tn+1, yn+1). (1.13)
V anglicky psané literatuře se IE metoda označuje jako implicit Euler method resp. back-ward Euler method. O implicitní metodě mluvíme proto, že neznámá yn+i je rovnicí (1.13) určena implicitně. Určit yn+i znamená řešit obecně nelineární rovnici. To je ve srovnání s EE metodou problém navíc. Aby mělo použití IE metody nějaký smysl, musí mít IE metoda oproti EE metodě také nějakou přednost. Pokusme se ji najít.
Nejdříve prozkoumáme přesnost IE metody. Lokální diskretizační chyba IE metody je definována předpisem
lte„ = y{tn+1) - y(tn) - rf(tn+1, y{tn+1)),
kde y(t) je řešení úlohy (1.1), (1.2). Z (1.12) plyne
lten = -|rV'(O = 0(r2).
IE metoda je tedy řádu 1 stejně jako EE metoda. Při podrobnějším zkoumání lze zjistit,
že
\y"{tn)T2 + 0(r3) pro EE metodu, -\y"(tn)T2 + 0(r3)   pro IE metodu.
Hlavní členy lokálních diskretizačních chyb (v případě Eulerových metod to jsou členy obsahující r2) jsou co do absolutní hodnoty stejné, liší se jen znaménkem. Můžeme proto oprávněně soudit, že obě metody jsou stejně přesné.
10
Podívejme se také na stabilitu IE metody. Pro testovací úlohu (1.10) je
1 / !
yn+1 = yn + r\yn+1,   odtud   yn+1 = y—-^ž/n = • • • = \^—^\J y°
a tedy podmínka stability (1.11) platí, právě když |1 — rA| > 1, neboli když r A leží v oblasti absolutní stability Ra'
t A g Ra = {z g C| \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+i získáme jako přibližné řešení rovnice g (z) = 0, kde
g{z) = z-yn-Tf(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: yn+1 =yn + \r[f{tn, yn) + f(tn+1, yn+1)}. (1.14)
Pro lokální diskretizační chybu
lte„ = 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 = -±T3y"'(tn) + 0(T4).
Lichoběžníková metoda (stručně TR metoda podle anglického trapezoidal rule) je tedy implicitní metoda řádu 2. Stabilitu TR metody zjistíme řešením testovací úlohy (1.10):
2 + r\in+1
y o ■
yn+1 =yn + \r\[yn + yn+i],   odtud yn+1 = 2 ^ TJ^yn = ■■■
Není těžké ověřit, že podmínka stability (1.11) platí, právě když
rA g Ra = {z g 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 + t(&i*4 + b2k2 H-----h bsks), (1.15)
kde koeficienty kl, i = 1,2,..., s, jsou určeny předpisem
h = f(tn,yn),
k2 = f(t„ + tc2, y„ + ra2ik1),
h = f(tn + tc3, yn + r(a3iA;i + 032^2)), (1-16)
ks = f(tn + rcs, yn + r(asi/ci + as2k2 H-----h aS)S_i*;s_i)),
a kde bt, cl5 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 k% pomocí k\, k2 a tak dále, až nakonec spočteme ks pomocí k\, k2,..., /cs_i. Vypočtené koeficienty kt, 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+i potřebujeme znát jen yn, předchozí hodnoty yn-i, yn-2, • • • v kroku od tn do tn+i nepoužijeme.
Koeficient kt je směrnicí lokálního řešení procházejícího bodem [t*,y*], kde
[tliVi] = [tn,yn], t*=tn+TCi, y* = y„+r(ali/ci+al2/c2H-----ha^-i^-i), i = 2, ?>,..., 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ľi=i b% = 1, viz (1.17).
Abychom dostali konkrétní metodu, musíme určit stupeň s a dále konstanty cl5 bt, aiy Konstanty RK metod je zvykem zapisovat do tabulky známé jako Butcherova tabulka:
C2	a2\		
c3		a-32	
cs	asi	as2 ...	
	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
lte„ = y(tn+1) - y(tn) - ^ bjki(tn),
i=i
12
kde
ki(t„) = f(tn,y(tn)),
l-l
h{tn) = f{tn + tc15 yn + r ^ aiókó{tn)),   i = 2, 3,..., s.
J=l
Lokální diskretizační chyba je chyba, které se dopustíme v jednom kroku za lokalizačního předpokladu yn = yitn). 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:
s
řád 1:    ^2 h = 1,
i=i
s s
řád 2:    ^6i = l,   J^&lCl = ±, (1.17)
i=l i=2
s    i — l
řád 3:    ^ ^ = 1,    ^ btc% = \,    ^ ^c- = \ ,    ^ ^ fcja.
i=l i=2 i=2 i=2 j=2
V C3 6
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 = 1, 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(tp). 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(ts), 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,...
í>(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 yl+1 =yl + rf(U, y i).
Metody řádu 2. Pro s = p = 2 má explicitní RK metoda Butcherovu tabulku
Podmínky (1.17) pro metodu řádu 2 stanoví
C2	Q21
	h b2
h + b2 = 1,    b2c2 = \ , 13
a protože ve shodě s (1.18) předpokládáme 021 = c2, dostáváme tabulku
kde ab = |. Parametry a, b jsou tedy svázány jednou podmínkou, takže zvolíme-li a 7^ 0, je b = l/(2a).
a	a
	1-b b
Pro a =	= \ je b = 1
Vn+l = Vn +	
k2 = f{tn + \t, yn + \rh), h = f{tn, yn), (1.19)
známou pod názvem modifikovaná Eulerova metoda. Budeme ji značit EM1 jako první modifikace Eulerovy metody. V anglicky psané literatuře je metoda (1.19) známa jako midpoint Euler formula.
Pro a = ljefo = |a dostáváme metodu
yn+1 = yn + t;t(/ci + k2),    kde   kx = f{tn,yn), k2 = f(tn + r,yn + r/ci). (1.20)
Budeme ji značit EM2 jako druhou modifikaci Eulerovy metody. Metoda (1.20) se také často uvádí pod názvem Heunova metoda. Pro a = |jefo = |a dostáváme metodu
yn+1 = yn + |r(A;i + 3/c2),    kde   kx = f{tn,yn), k2 = f(tn + |r, y„ + §r/ci),
známou jako Ralstonova metoda řádu 2 (stručně R2 metoda). Metody řádu 3. Pro s = p = 3 dostáváme Butcherovu tabulku
a 4 podmínky pro metodu řádu 3:
b1 + b2+b3 = l,   b2c2 + &3c3 = \ , &2c2 + &3c| = § ,    b3a32c2 = i .
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 metodu řádu 3 (stručně R3 metodu):
I = y„ + \t (2/ci + 3/c2 + 4fc3),
0  f /ci = /(í„, y„),   /c2 = /(ín + \t, yn + ér/ci)>
c2	c2	
c3	c3 - 032	Q32
	bi	&2 &3
2 14 9     3 9
ks = f(tn+'ÍT,yn+'±Tk2).
Ralstonova metoda řádu 3 je základem Runge-Kutta-Bogacki-Shampine 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 metoda
yn+1 = yn + \t (h + 2k2 + 2/c3 + k4), k\ = f(t„, y„), k2 = f(tn + \t, yn + |r/ci),
h = f(tn + \t, y„ + \rk2),   k4 = f(tn + r,yn + rk3).
1	1		
2	2		
1	0	1	
2		2	
1	0	0	1
	1	1	1 1
	6	3	3 6
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-Prince, 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 est„ lokální chyby le„ nabývala pořád přibližně stejné hodnoty e. Krok od yn do yn+\ je úspěšný, když
|est„| < e. (1.21)
Určení odhadu est„ lokální chyby se věnuje následující odstavec.
Je-li podmínka (1.21) splněna, krok je úspěšný a pokračujeme výpočtem yn+2- Pokud podmínka (1.21) splněna není, krok je neúspěšný a výpočet yn+\ opakujeme. Novou délku kroku r* určíme v případě úspěchu i neúspěchu stejným postupem, který si teď vysvětlíme.
Předpokládejme, že yn+i počítáme metodou řádu p, takže
le„ = Ctp+1 = est„    ==>    C = estn/rp+1.
Novou délku kroku r* zvolíme tak, aby velikost |le*| lokální chyby le* = C(t*)p+1 byla přibližně rovna zadané toleranci e, tj.
\K\ = \C\(r*)p+1 = \estn/rp+1\(r*Y+1 = e    =>    r* = r (e/\estn\)1/{p+1) .
Nová délka kroku r* se ještě redukuje pomocí parametru 6 < 1, takže
t* = 0t(e/\estn\)1/{p+1) . (1.22)
V matlabovských programech ode23 a ode45 se bere 6 = 0,8. Současně se ještě uplatňují následující zásady.
1) Tolerance e se uvažuje ve tvaru
e = max{er max{|y„|, |yn+i|},eQ},
kde er je relativní tolerance a ea je tolerance absolutní. Matlabem přednastavené hodnoty jsou er = 10~3 a ea = 10~6.
2) Označme rmm resp. tmax minimální resp. maximální povolenou délku kroku. Jestliže t* < tmm, výpočet končí konstatováním, že danou diferenciální rovnici program neumí s požadovanou přesností vyřešit. Přitom rmm = 16em(r„), kde emitn) je tzv. relativní přesnost aritmetiky pohyblivé řádové čárky, emitn) = 2,2 • 10_16|r„|, viz funkce eps v Matlabu. Jestliže r* > Tmax, položí se r = Tmax, kde je Tmax = 0,l(b—a).
3) V případě neúspěšného kroku navrženou délku r* redukujeme:
r* = max(r*, <?mmr) ,
kde qmm = 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* = mm(T*,qmaxT), 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 = Oe]/^ [max(M,eQ/£r)]/\f(a,V)\, (1.23) přičemž pro r < rmm změníme r na rmm a pro r > rmax změníme r na Tmax.
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'n*+1 přesnější metodou a Vn+i méně přesnou metodou. Použitelný odhad lokální chyby je
esln = Vn+1 ~ Vn+1-
Obě metody používají tutéž množinu koeficientů {k0\sJ=l. V t°m 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ž
y*+i = yn + T Yf]=i tfkj   je metoda řádu p, y*+i = y-n + r Y^j=i b*j*k-j   je metoda řádu p + 1, est„ = r Yľj=i Ejkj   je odhad lokální chyby,
takže Ej = b** -6* i = 1, 2,..., s.
Pokud ve výpočtu pokračujeme přesnější metodou, tj. když yn+i = yn*+i = y^+i+ est„, ří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+\ = 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.
c2	«21		
cs	asi	as2 ■	
	bl	h*	■ K
	bl*	b*2* ■	. b*s*
	Eí	E2 ■	■ Es
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 4	0	3 4		
1	2	1	4	
	9	3	9	
	7	1	1	1
	24	4	3	8
	2 9	1 3	4 9	0
	5	1	1	1
	72	12	9	8
[§*1
Vn+l ~ Vn
řádu 3. Pomocná metoda
y*n+i = yn + r[^ik1
3fc2
= Vn+l
řádu 2 používá kromě koeficientů ki, k2 a /c3 navíc ještě koeficient k± = 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 k± z kroku předchozího, takže nově se počítají jen koeficienty k2, k% a k±. Výjimkou je případ neúspěšného kroku, kdy se hodnota yn+i 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.
Vn+l = Vn+1-
Hodnoty přibližného řešení pro t E {tn,tn+i) spočteme dostatečně přesně pomocí kubického Hermitova polynomu H^(t) určeného podmínkami
H3(tn) =
yn,        Há(tn) = ki, = Vn+l,   H'z (tn+i) = ki .
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	9 017 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 57600	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
	57600		16 695	1920	339 200	525	40
Metoda DP54 je typu FSAL, neboť k7 = f(tn+1, yn+i)- Proto se v každém úspěšném kroku
17
metody počítají jen koeficienty &2,..., k?, koeficient k\ byl už vypočítán jako koeficient k'j v předchozím kroku. Zdůrazněme, že DP54 metoda se používá jako metoda s lokální extrapolací, tj. yn+1 = y*n*+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 = iki, &2, k^,     k5, ke, kj),
B
183	37	145 -i
64	12	128
0	0	0
1500	1000	1000
371	159	371
125	125	375
32	12	64
9 477	729	25 515
3 392	106	6 784
11	11	55
7	3	28
3 2	-4	5 2 J
q = (q, q2, g3, q4)T, přičemž q = (t — tn)/r.
Metoda DP54 je rovněž vynikající: je řádu 5, v každém úspěšném kroku se pravá strana počítá jen 6-krát, a to stačí jak pro řízení délky kroku tak pro výpočet řešení v intervalu
{tni tn+l)-
Stabilita. Rešíme-li testovací rovnici (1.10) RK metodou na rovnoměrném dělení s krokem r, dostaneme
yn+1 = Ps(r\)yn,
kde Ps je polynom stupně s určený pomocí konstant bl, al3 definujících RK metodu. Podmínka stability (1.11) tedy platí, právě když |Ps(rA)| < 1, neboli když rA leží v oblasti absolutní stability R a.
t\eRa = {z G C I \Pa(z)\ < 1}.
Explicitní R K metody mají omezenou oblast absolutní stability, neboť pro \z\ —> oo také |Ps(z)| —> oo. Dá se ukázat, že pro p = s < 4 je Ps(z) = 5ľi=ozV*'- Pr°t° každá explicitní s-stupňová RK metoda řádu p = s < 4 (stručně RKs metoda) má stejnou oblast absolutní stability. Oblast absolutní stability RKs metod pro s = 1, 2, 3,4 a metody DP54 je zobrazena např. v [12]. Intervaly absolutní stability těchto metod jsou Ia = (a, 0), kde
metoda	RK1	RK2	RK3	RK4	DP54
a	-2	-2	-2,51	-2,79	-3,31
Zdůrazněme, že z pohledu stability je BS32 metoda ekvivalentní s RK3 metodou, obě metody tedy mají stejnou oblast a stejný interval absolutní stability.
18
1.4. Lineární mnohokrokové metody
V této kapitole se budeme zabývat metodami, které počítají přibližné řešení yn+\ v uzlu tn+i pomocí dříve spočtených aproximací yn, yn-i, yn-2, • • • a odpovídajících hodnot f(tn, Vn), f{tn-i, Vn-i), f(tn-2, Vn-2), ■ ■ ■ pravé strany diferenciální rovnice. Tyto hodnoty jsou znovu použity tak, abychom získali yn+i s vysokou přesností pomocí jen několika málo nových vyhodnocení funkce fit, y). Nejznámějšími metodami tohoto typu jsou Adamsovy metody a metody zpětného derivování. Obě skupiny metod patří do obecné třídy metod známých jako lineární mnohokrokové metoáy, stručně LMM.
1.4.1. Obecná lineární mnohokroková metoda
LMM je předpis
a0yn+i + axyn H-----Y akyn+1_k = r[/30/(ŕ„+i, yn+i) + Pif n H-----h (3kfn+i-k\, (1-24)
ze kterého počítáme yn+i. Přitom a3 a (3j jsou číselné koeficienty, které formuli jednoznačně určují, a f j je zkrácený zápis pro f(tj,yj). V dalším budeme předpokládat, že platí normalizační poámínka a0 = 1. Jestliže alespoň jeden z koeficientů ak nebo [3k je různý od nuly, metoáa je k-kroková.
Pro [3q 0 je nová hodnota yn+\ určena implicitně, hovoříme proto o implicitní metoáě, pro [3q = 0 máme metoáu explicitní. Abychom v implicitní metodě určili yn+i, musíme vyřešit obecně nelineární rovnici.
LMM lze použít, jen když jsou zadány startovací hoánoty yo, yi, ..., yk-\. yo = f] určíme z počáteční podmínky, zbývající startovací hodnoty je však třeba získat jinou vhodnou metodou, yr metodou nejvýše r-krokovou.
Lokální diskretizační chyba je chyba, která vznikne, když do formule (1.24) dosadíme místo přibližného řešení yn+i_j přesné řešení y(tn+i-j), tedy
k k
lte„ = ^^(ín+i-j) -t^ft/(Ul-J,!/(UH))'
3=0 J=0
Jestliže
lte„ = Cp+1Tľ+1y^+1\tn) + 0(t^+2),
řekneme, že metoda je řádu p. Clen Cp+iTp+1y(p+1\tn) se nazývá hlavní člen lokální áiskretizační chyby, konstanta Cp+i je tzv. chybová konstanta. LMM je tím přesnější, čím je vyššího řádu. Z několika metod téhož řádu je pak nejpřesnější ta metoda, pro kterou je velikost chybové konstanty |Cp+i| nej menší.
D-stabilita. Řekneme, že LMM je stabilní ve smyslu Dahlquista (stručně D-stabilní), jestliže všechny kořeny prvního charakteristického polynomu
Q(Z) = Šk + "i?""1 + • • • + + dk
leží uvnitř jednotkového kruhu \z\ < 1 komplexní roviny C a pokud některý kořen leží hranici \z\ = 1, pak je jednoduchý.
19
Význam D-stability lze vysvětlit na rovnici y' = 0. Její řešení pomocí LMM vede na předpis Y^=o ajyn+i-j = 0. Zvolíme-li startovací hodnoty y3 = er°, j = 0,1,..., k — 1, kde gir) = 0 a e je libovolné číslo, pak yn = ern pro každé n. Skutečně,
k k
a3ern+1-J = ern+1~k ^ ajrk-j = ern+1~kg(r) = 0.
3=0 3=0
Pro \r\ > 1 a e ^ 0 je lim^^ \yn\ —> oo, což je nepřijatelné: pro e = 0 je yn = 0 přesné řešení rovnice y' = 0, avšak pro e ^ 0, |e| 1, tj. již pro nepatrnou poruchu startovacích hodnot y3, 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(tp), pak globální diskretizační chyba je rovněž řádu 0(tp).
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(ts).
Absolutní stabilita. Rešíme-li testovací úlohu (1.10) LMM na rovnoměrném dělení s krokem r, dostaneme
k
^2(aj - rX^y^-j = 0 . (1.25)
3=0
Řešení hledejme ve tvaru yn = rn. Po dosazení do diferenční rovnice (1.25) obdržíme
k k
^2(aj - T\f3j)rn+1-J = rn+1~k ^(^ - tA/^)^-3' = rn+1-kir(r, t A) = 0 ,
3=0 3=0 k
kde ^,z) = ^2(aj-zPj)e-j
3=0
je tzv. polynom stability LMM. Jestliže 7r(r, t A) = 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 Ra 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ž
rA g Ra = {z g C I tt(£, z) = 0 => |£| < 1}.
1.4.2. Adamsovy metody
Integrací diferenciální rovnice (1.1) od tn do tn+i dostaneme
y(tn+1) - y(tn) = [n+1 /(í,y(í))dí.
Jtn
20
Funkci f (t, y (t)) aproximujeme pomocí interpolačního polynomu Pk-i(t) stupně k — 1, tj.
ľtn + l
y{tn+1) !=t y{tn) +/      Pfc_i(ŕ)dŕ,    kde   Pfc_1(ŕíl+1_;;-) = f(tnJrl_j, y^+i^-)).
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
ľtn+l
yn+i=yn+ Pfc_i(ŕ)dŕ,   kde   Pfc_1(ŕíl+1_;;-) = f{tn+1_v yn+1-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
yn+1 =yn + r^2 Pk,j •
Stručně ji budeme označovat jako ABk metodu. ABk metoda je explicitní, /c-kroková, řádu k, D-stabilní. Pro konstantní délku kroku, tj. když
ín+l-j = tn+1 - JT,    j = 1, 2,..., k,
jsou koeficienty ABk metod pro k = 1,2,..., 6, spolu s chybovými konstantami Cl+1 a dolními mezemi a*h intervalů absolutní stability (aj*, 0), uvedeny v následující tabulce:
k	Pl,i	Pl,2	Pl,3	ftk,4	0k,5		^fc+l	
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+1 =yn + r^2 7j* VJ/n,
J=0
21
kde V-7 je operátor zpětné diference,
V°/„ = /„,   V/„ = /„ - • • • V7™ = V1"1/™ - V1"1/™-!,   * = 2, 3,... .
a koeficienty 7* jsou definovány rekurentním předpisem
70 = 1, 7; = i-É^ ; = i,2,...
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
t	tn+l			tn+2-k
Pk-l{t)	fn+1	f n		fn+2-k
Adams-Moultonovu metodu lze zapsat ve tvaru fc-i
yn+1 =Vn + T1^2 (3kj fn+l-j ■ 3=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 ak intervalů absolutní stability (ak, 0), jsou uvedeny v následující tabulce:
k	Pk,0	Ä,l	Ä,2	Ä,3	Pk,4	Ä,5	Ck+i	
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 + T^2 li V* fn+1,
3=0
22
kde koeficienty 7^ jsou definovány předpisem 7o = l,    7J=7*-7*_1,    i = 1,2,.... Chybová konstanta Ck+i = jk- 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
yn+1 = ip(yn+i),   kde   ip(z) = yn + r(3kfi f(tn+1, z) + r
Použít můžeme metodu prosté iterace: zvolíme počáteční aproximaci y^+i a postupně počítáme y^+i = (p(yns+i^), s = 1, 2,.... Pro dostatečně malé r metoda konverguje. Jestliže počáteční aproximaci y^+i určíme pomocí AB metody, provedeme jen několik málo iterací
Ä = v{yns+i),  s = 1,2,...,s,
a nakonec položíme
Vn+l = Vn+1 ,      fn+1 = f(tn+l, Un+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 f(tn+i, yns+i)- Zvolíme-li jako prediktor metodu ABk a jako korektor metodu AMk, dostaneme metodu, kterou značíme ABk-AMk-P(EC)'sE. 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 est„ lokální chyby. Jestliže = y^+i spočteme prediktorem ABk a y'n*+i = yi+i korektorem AMk, pak tzv. Milnův odhad lokální chyby dává
est„ =     Ck+1r    {y*:+í - (1.27)
odvození viz [12]. Nakonec položíme
Vn+l = y*n*+l + estn- (1-28)
Korekce y'n*+1 pomocí odhadu lokální chyby est„ 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+i spočteme metodou AM(k+l) a položíme
est„ = yn+1 - y*n*+l. (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/„ -
fn+l = f{tn+l, Vn+l)
AM2:      y*n*+1=yn + lT(f:+1 + fn)
est„ = -|(y*+i - y*n+1),   yn+i = yn*+i + est™
fn+l = f{tn+l, Vn+l)-
V případě AB2-AM3-PECE metody nahradíme řádek L řádkem
C:      AM3:      yn+1 = yn + ^r(5/*+1 + 8/„ - estn = yn+1 - y™+1.
Startovací hodnotu yi 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 (variable step variable order). Základní myšlenka je jednoduchá. Předpokládejme, že jsme vypočetli yn+i metodou ABk-AM(k+l)-PECE s krokem délky t. Určíme odhad estjj 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ů t^_1, t£ a t£+1 stanovíme podobně jako pro jednokrokovou metodu,
-r/ = 6>t(e/|es<|)1/(m)   pro í = k - 1, k, k + 1,
a největší z čísel t^_15 t£, t£+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* = rjľ, je-li největší Tfc+u zvětšíme k o jedničku a pokračujeme s krokem délky r* = rjľ+1 a je-li největší t^_15 k o jedničku snížíme a pokračujeme s krokem délky r* = t£_v
Pro krok délky r* 7^ r je třeba vypočítat hodnoty f*+í_j pro - = tn+1 — jr*.
To lze snadno provést interpolací pomocí fn+i, 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 /c-krokovou metodu zpětného derivování použijeme zkrácený zápis BDFk metoda. Dostaneme ji tak, že v rovnici
y'(fn+i) = f(tn+1,y(tn+1))
nahradíme derivaci y'(tn+i) pomocí derivace P^(rra+i) interpolačního polynomu Pk(t) stupně k procházejícího body [tn+1,y(tn+1)], [tn, y(tn)], [tn+i-k, y(tn+1_k)]. Když pak nahradíme y(tn+i_j) přibližnými hodnotami yn+i_j, j = 0,1, ..., k, dostaneme metodu
Pkifn+l) = f(tn+1,yn+1). Přehledné vyjádření aproximujícího polynomu Pkit) uvádí následující tabulka
t	tn+l			tn+i-k
Pkit)	Vn+1	yn		yn+i-k
BDFk metodu zapíšeme ve tvaru
k
&k,jyn+l-j = TPkfifn+l ■
3=0
Metoda BDFk je implicitní, /c-kroková, řádu k, pro k < 6 je D-stabilní. Pro konstantní délku kroku jsou koeficienty akj, [3kfl a chybové konstanty Ck+i BDFk metod uvedeny v následující tabulce:
k	«fc,0		«fc,2	«fc,3	«fc,4	«fc,5	«fc,6	Pk,0	Ck+1	
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	36	16	3			12	12	73°
		25	25	25	25			25	125	
5	1	300 137	300 137	200 137	75 137	12 137		60 137	10 137	52°
6	1	360	450	400	225	72	10	60	20	18°
		147	147	147	147	147	147	147	343	
Koeficienty BDFk metody lze získat prostřednictvím předpisů
k    ^ ^ k    ^ ^
7)]-VVl = t — fn+1,      4 = V-,      Ck+1 =
Ókj^J ók j^J 0 + 1)4
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 R& absolutní stability obsahuje
25
celou zápornou polorovinu komplexní roviny C, tj. D {z E 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 E (0,7r/2), jestliže její oblast absolutní stability Ra obsahuje nekonečný klín
Wa = {relv E C I r > 0, \(p - tt| < a}.
BDFk metody jsou (pro k < 6) A(a)-stabilní, příslušné úhly (pro větší názornost ve stupních) jsou uvedeny jako poslední sloupec výše uvedené tabulky.
BDFk metody (pro k < 6) jsou dokonce L(a)-stabilní. L(a) stabilní metodu přitom definujeme jako A(a)-stabilní metodu, pro kterou z podmínek 7r(£, z) = 0 a Re(z) —> —oo plyne £ —> 0. Pro a = f dostáváme L-stabilní metodu. BDF1 a BDF2 jsou tedy L-stabilní.
Uhel 18° metody BDF6 je příliš malý a proto se tato metoda obvykle nepoužívá. V matlabovském programu odel5s jsou implementovány metody zpětného derivování řádů 1 až 5. Program odel5s je typu VSVO, tj. volí optimální délku kroku i řád metody. Základní metodou programu odel5s je metoda NDF (podle numerical differentiation formula). NDFk metody jsou modifikace BDFk metod, mají o něco menší chybové konstanty (o 26% pro k = 1, 2, 3 a o 15% pro k = 4 ) a poněkud menší úhly A(a) stability (o 7% pro k = 3 a o 10% pro k=4) než odpovídající BDFk metody. Podrobnosti týkající se NDFk metod lze najít v [24].
Nelineární rovnice pro výpočet yn+\ se řeší pomocí několika málo kroků Newtonovy metody.
1.5. Tuhé problémy
Tuhý počáteční problém (v anglicky psané literatuře stiff problém) se vyznačuje několika charakteristikami, z nichž dvě si postupně uvedeme. Praktická a snadno ověřitelná je
Charakteristika 1. Počáteční problém je tuhý, když počet kroků potřebných k jeho vyřešení metodou s omezenou oblastí absolutní stability je podstatně větší než počet kroků, který k jeho vyřešení potřebuje metoda s neomezenou oblastní absolutní stability.
Platnost této charakteristiky ukážeme na příkladu počátečního problému
Příklad 1.1
)©•     Gffi=(-0- ™
Řešení je yi = — e_í, y2 = e~l. Úlohu (1.30) jsme řešili matlabovským programem ode45 (DP54 metoda řádu 5 s omezenou oblastí absolutní stability) a matlabovským programem odel5s (BDF metody VSVO řádů 1-5 s neomezenými oblastmi absolutní stability). Oba programy jsme použili se stejným požadavkem na přesnost: er = 10~3 a ea = 10~6. Efektivnost obou metod lze přibližně porovnat podle počtu úspěšně provedených kroků pk a počtu pf vyhodnocení pravé strany. V následující tabulce jsou uvedeny hodnoty pk/pf pro několik délek í intervalu integrace.
26
í	10"2	ío-1	10°	101	102
ode45	10/61	22/151	269/1747	2 953/18 919	30 071/192 475
odel5s	10/24	10/24	12/28	42/88	71/146
Pro malé í = 10~2 se na intervalu (0; 10~2) řešení poměrně rychle mění. Délku kroku zde určuje požadavek na přesnost a protože obě metody jsou téhož řádu, potřebují přibližně stejný počet kroků. Jak však í vzrůstá, délku kroku stále více začíná ovlivňovat podmínka stability. □
Abychom tomuto efektu lépe porozuměli, potřebujeme několik dalších poznatků.
Stabilní problém. Počáteční problémy' = f (ŕ, y), y (a) = r\, je stabilní, když malá změna dat f, a, r\ způsobí malou změnu řešení y. Zabývejme se speciálně stabilitou vzhledem k počáteční podmínce, konkrétně jak změna počáteční hodnoty r\ ovlivní řešení y. Pro jednoduchost uvažujme počáteční problém
y' = Ay,       y(a) = r,, (1.31)
kde A je číselná matice řádu d. Nechť u je řešení problému (1.31) a v je řešení téže diferenciální rovnice, avšak s porušenou počáteční podmínkou, tj.
u' = Au,       u(a) = r), v' = Av,       v(a) = r) + ô.
Pro w = u — v dostaneme problém
w' = Aw,       w(a) = ô,
který popisuje šíření počáteční poruchy ô. Předpokládejme, že matice A má navzájem různá vlastní čísla {A-,}^=1. Pak
d
w(t) = J>e^-<S,
kde ~Vj jsou vlastní vektory příslušné vlastním číslům \0 a Cj jsou konstanty, které určíme z počáteční podmínky:
Vc = ô,   kde   V = (vi, v2,..., vd),   c = (cľ, c2,..., cd)T.
Můžeme tedy vyslovit tyto závěry:
1) Jestliže Re(\j) < 0 pro všechna j, pak ||w(r)|| —> 0 exponenciálně pro t —> oo, tj. porucha se velmi rychle zmenšuje.
2) Jestliže Re(Xj) < 0 pro všechna j, pak ||w(r)|| < ||V|| • ||V_1|| • \\S\\, tj. porucha bude omezená, přitom ||w(r)|| —> 0 pro ô —> o.
27
3) Jestliže nějaké vlastní číslo \j má kladnou reálnou složku a počáteční porucha ô je taková, že Cj 7^ 0, pak ||w(r)|| —> 00 exponenciálně pro t —> 00, tj. porucha se velmi rychle zvětšuje.
Problém (1.31) je tedy stabilní (vzhledem k počáteční podmínce), jestliže všechna vlastní čísla matice A mají nekladnou reálnou složku. □
Spektrální poloměr. Označme symbolem g(A) spektrální poloměr matice A definovaný jako velikost nej většího vlastního čísla A, tj.
g(A) =   max kde \j, j = 1, 2,..., d, jsou vlastní čísla A.
j=l,2,...,d
Nyní již můžeme zformulovat další charakteristika tuhého systému.
Charakteristika 2. Stabilní problém je tuhý, jestliže součin spektrálního poloměru Jaco-biovy matice fy(ŕ, y) = {dfi(t,y)/dy~j}fj=1 a délky intervalu integrace b — a je velký, tj. když
max£(fy(ŕ,y(ŕ)))(&-a)>>l. (1.32)
a<t<o
Vraťme se nyní k úloze (1.31). Pro f = Ay je fy(r, y) = A. Problém (1.31) je tedy tuhý, pokud všechna vlastní čísla matice A mají nekladnou reálnou složku a g(A)(b — a) je velké číslo.
Jsou-li všechna vlastní čísla matice A jednoduchá a mají zápornou reálnou složku, pak pro řešení y (ŕ) úlohy (1.31) platí y (ŕ) —> o pro t —> 00. Rešíme-li takovou úlohu numericky, pak snadno dokážeme, že podmínka stability yn —> o pro tn = nr —> 00 platí, právě když
{Air, A2r,..., Xdr} G RA, (1.33)
kde Ra je oblast absolutní stability uvažované numerické metody. Příklad 1.1 — pokračování. Pravá strana diferenciální rovnice je
0 1
Ay,       kde   A = ,
1   -1000 -1001
Vlastní čísla A jsou Ai = —1, A2 = —1000, takže g(A) = 1000. Jde o stabilní problém (vlastní čísla A jsou záporná) a pro větší í je g(A)(b — a) = 1000í 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 < 1000t, tj. r < 0,00331. Na intervalu délky L je tak třeba nL > 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í yi = —e_í 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,      te (0,2/5),      y(0) = 6. (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 yit) —> 0 pro t —> oo, zatímco pro y(0) > 0 dostaneme y(r) —> 1 pro t —> oo. Zvolíme-li 5 > 0 velmi malé, pak pravá strana y2 — 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 ó = 10~4 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|2y — 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
0.8 0.6 0.4 0.2 0
x 10
Obr. 1.2 Příklad 4.3 řešený DP54 metodou, celý výpočet
Ú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,11200). V následující tabulce uvádíme
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:
í	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í. □
1.0001
x 10
Obr. 1.4. Příklad 4.3 řešený TR metodou, celý výpočet
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
y„+i = 7Tf (tn+1, yn+1) + if),
;i.35)
30
kde parametr 7 je charakteristická konstanta metody a rj) je vektor nezávislý na yn+i-Například pro TR metodu (1.14) je
1= \ a     = yn + \ri(tn, y„).
Rovnici (1.35) řešíme zjednodušenou Newtonovou metodou. Počáteční aproximaci y^+1
získáme dostatečně přesnou extrapolací z hodnot yra,yra_i,.... Další aproximace y^^ získáme řešením soustav lineárních rovnic
(i - r7J) (yJJi15 - yiti) = r7f(Wi,yiti) + </> - yiti, (1.36)
kde I je jednotková matice a J je Jacobiova matice {y(tn+i_k,yn+i-k) W° nějaké k > 0 (jedná se tedy o zjednodušenou Newtonovu metodu, pokud by J = fy (tn+i, y^+i), šlo by o klasickou Newtonovu metodu). Výpočet organizujeme takto: označíme
G = I - r7J,       ds = yi!^ - yit,       gs = r7f (tn+1, y%) + ^ - y%
a vypočteme nejdříve ds jako řešení soustavy lineárních rovnic Gds = gs a pak dopočítáme yi+i"^ = yi+i + ds. Je-li přírůstek ds dostatečně malý, klademe yn+i = Yn+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\funfun). 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 íy(tn, yn) a vektor derivací ft(r„,yra) = {dfi(tn,yn)/dt}f=1, provést LU rozklad matice I — jTÍy(tn,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
^ = -0,04yi + 104y2y3, Z/i(0) = l,
y'2 = 0,04 yi - 104y2y3 - 3 • 107 y\ ,      y2(0) = 0 , y3 = 3-io7y22, y3(0) = 0
popisuje koncentrace tří příměsí v chemické reakci, tj. 0 < yi,y2,y3 < 1, nezávisle proměnná t je čas, blíže viz [22]. Jacobiova matice této soustavy je
/-0,04 104y3 104y2\
J =      0,04  -104y3 - 6 • 107y2   -104y2 . V   0 6-107y2 0 /
V čase t = 0, tj. pro yi = 1, y2 = y3 = 0, má Jacobiova matice vlastní čísla {—0,04; 0; 0}. Z fyzikálních úvah plyne, že yl5 y2 —* 0 a y3 —> 1 pro t —> oo. Vlastní čísla Jacobiovy matice pro yi = y2 = 0, y3 = 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 g(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 nej zná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{x), 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),..., yd{x))T, i{x,y{x)) = (fi(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)), .. ., rd(y(a), y(fo)))T. Ve speciálním případě, když
r(y(a), y (6)) = y (a) - r] = o   nebo-li   y (a) = r], (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(0) = y(l)       existuje jediné řešení y = 0,
(b) e y {tí) = y{l)        existuje nekonečně mnoho řešení y = Cex, C libovolné,
(c) e y {tí) = y{l) + 1 řešení neexistuje.
33
Metoda střelby je založena na numerickém řešení počátečního problému
ý{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 r\ je určena implicitně rovnicí (2.5).
Lichoběžníková metoda. Nechť a = xq < x\ < ■ ■ ■ < x^ = b je dělení intervalu (a, b), ht = xl+i — Xj je délka kroku a h = max, hl. Přibližné řešení yl5 i = 0,1,..., N, dostaneme jako řešení soustavy d(N + 1) rovnic
yl+i = yl + \h%[í(x%,y%) r(yo, Yn) = 0 .
+ f(xl+i,yl+i)
i = 0,1,... ,iV — 1,
(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 yf^1 a; y(xj), 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 Oih2) je každá složka vektoru y(xj) — yl.
Ve speciálním případě, když f(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
(2.7)
y' = A(x)y+ q(x), BQy(a) + B6y(fo) = o. Soustava (2.6) je pak tvaru
[-1 - \fnA{xi)\ yl + [I - ±hiA(xi+1)] yl+i = ^[qfc) + q{xi+1)], i = 1,2,... ,N - 1,
Bay0 + B6yAr = o , kde I je jednotková matice. Maticový zápis této soustavy je
/ R-o So
Ri Si
B6 J
í yo \
yi
yjv-i
í vo \
Vl
VAT-1
v ° /
(2.8)
kde
R,
I 2/ijA(xj
S, = I
\hiA{xi+1) ,        = i/ij[q(xj) + q(xl+i)] . 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
yl+i = yl + \hi[i{xuyi) + 4f(xl+i/2,yl+i/2) + f (xi+1, yi+1)],
(2.9)
kde   xl+1/2 = Xi + \hu    yl+i/2 = |(yl +y"i+i) +        (x,, y i) - í(xi+1,yi+1)] , Simpsonova formule (2.9) je řádu 4, tj. pro chybu platí y(xi)-yi = 0(h*).
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 Rl5 Sj a vektory Vj získáme z formule (2.9), v níž klademe f(x,y) = A(x)y + q(x). Odvovození vzorců pro Rl5 Sj a Vj 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),       xe{0,£). (2.10)
Předpokládejme, že funkce p(x), p'(x), q(x) i fix) jsou spojité, dále že pix) > 0, q(x) > 0, a uvažujme nejdříve jednoduché okrajové podmínky
u{0)=g0, (2.11a) u(£)=ge. (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 uix) posunutí střednicové čáry prutu, pix) = Eix)Aix), kde Eix) je Youngův modul pružnosti a Aix) je plocha průřezu prutu, qix) je měrný odpor podloží, na němž prut spočívá, fix) je intenzita zatížení a g0 a    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, pix) je koeficient tepelné vodivosti, fix) — qix)uix) je intenzita tepelných zdrojů, go 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
y=(yi) = (U),       í=( m/Pf),       r=(yi?l~9°), (2-12)
W   Wy \<m-fj \yÁZ)-9i)
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ť 0 = xo<a;i<---<a;7v = ^je dělení intervalu (O, £) a ht = xt — je délka kroku. Body xl nazýváme uzly a množinu {x^^Lq uzlů nazýváme sítí. Pro jednoduchost se omezíme na ekvidistantní dělení, takže ht = h = £/n a xt = íh, i = 0,1,..., n.
Splnění rovnice (2.10) budeme vyžadovat pouze ve vnitřních uzlech sítě, tj.
-tim')'\x=Xi + qM*i) = /i>   * = 1, 2,..., JV - 1, (2.13)
kde ql = q{xi) a fl = /(xj). Clen —{pu')'\ nahradíme diferenčním podílem. Pomocí označení
Xi-l/2 =Xi-\h,  xl+1/2 = xl + \h, pl-l/2 = P{xl-l/2) , Pi+1/2 = P(^i+l/2)
lze přibližně položit
^u)\x=xi- h h
uixl+i) — u{xj) uixj) — uixl_i)
^ Pí+1/2    h    ~PÍ-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{Xj-l) + (Pi-l/2 + Pi+l/2)u{Xj) ~ pl+l/2u{xl+1)   | 2
I x=xí fa2
Po zanedbání chyby Oih2) tak z (2.14) a (2.13) dostáváme pro přibližné řešení ul m u(xl) soustavu rovnic
+ (Pi-1/2 + Pi+1/2 + ^2g»K ~ K+l/2^+1 _ f , ,
h2 -Ji, (Z.Lb)
pro i = 1, 2,..., n — 1. Z okrajových podmínek máme
u0 = #o, (2.16a) uN = gí. (2.16b)
To nám umožní dosadit do první rovnice u0 = g0 a člen — Pi/2go/h2 převést na pravou stranu. Podobně v poslední rovnici položíme = j( a člen — PN-i/2ge/h2 převedeme na pravou stranu. Matice takto upravené soustavy rovnic (2.15) je třídiagonální, pozitivně definitní, diagonálně dominantní (pro qt > 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í
u(xí) -u, = 0(h2), (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í Dirichletovy. V aplikacích se objevují ještě další typy okrajových podmínek. Podmínky
p(0)u'(0) = a0u(0) - /30, (2.18a) -p{t)u\t) = aeu(£) - pe (2.18b)
se nazývají Newtonovy nebo také Robinovy. Pokud koeficient obsahující člen oíq resp. chybí, hovoříme o Neumannově okrajové podmínce. Aby byla zaručena jednoznačná existence řešení, předpokládejme «o > 0,     > 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     > 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, jsou tuhosti pružin v koncových bodech prutu a f30, (3g jsou zadané síly působící na koncích prutu. V úloze stacionárního vedení tepla je pu' tepelný tok, «o, OLi jsou koeficienty přestupu tepla a (3q, [3^ jsou zadané tepelné toky na okrajích tyče. Tepelné toky se často uvažují ve tvaru [3q = oíqUq, [3^ = 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ě uN. 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_x vyjádříme nejdříve pomocí jednostranné diference a pak pomocí centrální diference a vztahu (2.18a), tj.
pu' _     — pu' _
- (p"')'U =--1,2ih       ° + Oih) =
-Pl/2—-'—-*-L + 010U(X0) - Po
--\-0(h).
\h
Zanedbáme-li chyby, dostaneme rovnici pro neznámou u0 (pi/2 + ha0 + \h2q0)u0 - py2U\ 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 í = N, a po úpravách obdržíme rovnici pro neznámou u^
-ř»Ar_i/2iíAr-i + (pn-1/2 + ha? + \h2qN)uN     11 ,ninl.
= tPí + 7.jn ■ (2.19b)
h2 h' 1 2*
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),       x G (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 fix) — 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 uix) 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 člen, (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) - (30 , (2.21a)
— (pu — ru) |     = oi(u(Í) — Pí . (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^ = r-+1/2"(x-+1/2)~r-1/2"(x-1/2) + 0(h-) (2.22)
a v koncových uzlech pomocí jednostranné diference
,   yi r1/2u(x1/2) -r0u(x0)
2 (2.23)
/   x/i rNu(xN) - rN_1/2u(xN_1/2)
M l=XN =-p-+ °W •
V (2.22) a (2.23) vyjádříme konvekční tok ru ve středech xt_i/2 úseček pomocí aritmetického průměru hodnot u v koncových bodech, stručně centrální aproximace,
Ti-yziiixi-yz) = \ri-1/2[u{xi-.1) + u(xí)] + 0(h2),       t = 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
^1/2(^0 + ^1)/^,
přidáme člen   ^     \[rl+i/2(u% + ui+1) - r^/^u^ + u%)]/h , (2.25) -1^-1/2(^-1 + uN)/h.
38
Matice takto vzniklé soustavy rovnic je třídiagonální a nesymetrická. Dá se ukázat, že když
\h\rt_l/2\ <Pi-1/2,       2 = 1,2,..., N, (2.26)
pak je matice soustavy regulární. Pro dostatečně jemné dělení intervalu (0, £) bude podmínka (2.26) jistě splněna. Pro chybu opět platí (2.17).
Dominantní konvekce. Podmínka (2.26) může být značně omezující v případě, kdy konvekční koeficient výrazně převažuje nad koeficientem difúzním, tedy pro \r\ ^> p. Pak je účelné vyjádřit konvekční tok ru ve středech xl_i/2 úseček takto:
{7^-1/2^(^-1) + O (ti)   pro rvi/2 > 0, (2.27) rl-i/2u(xl) + O (ti)      pro ri_1/2 < 0.
Aproximace (2.27) je z fyzikálního hlediska přirozená: informaci o řešení u v uzlu xt čerpáme ze znalosti řešení proti „proudu", proti „větru": pro rl_i/2 > 0 „fouká zleva", proto použijeme hodnotu u(x^i) v uzlu Xj_i ležícím nalevo od bodu Xj_i/2, pro rj_i/2 < 0 „fouká zprava" a proto použijeme hodnotu u(xl) v uzlu xt ležícím napravo od bodu Xj_i/2-Jednostrannou aproximaci konvekčního toku podle (2.27) nazýváme upwind aproximací. Pomocí označení
a+ = max(a, 0),       a~ = min (a, 0),    kde a je libovolné číslo,
lze (2.27) zapsat v kompaktním tvaru
^-1/2^(^-1/2) = r+_ 1/2u(xl_1) + r~_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 r    [rt/2u0 + r-/2Ul]/h,
(2.15) l přidáme člen < [rt+1/2ui + r^i/2^+i]/^ ~ [rt-1/2^-1 + r1-i/21IJ/řl' (2.29) (2-19b) J 1^ -[r+_i/2UN_1 + r-_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 O (ti1) 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
-eu" + u' = 0   pro   x G (0,1),   u(0) = 0, u(l) = 1,
39
kde e > O je konstanta. Přesné řešení
1 - ex/£
u\x) =
1
,1/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 - 2-Uj + iíj_i    iíj+i - iíj_i
= 0,
h2 2h kterou lze při označení n := h/e zapsat ekvivalentně ve tvaru
-(1 + \k)u%-i + 2u% - (1 - \k)uí+1 = 0.
Pro k = 2 dostaneme -Uj = iíj_i, takže řešení je -Uj = 0 pro í < N a = 1. Pro řt / 2 snadným výpočtem ověříme, že diferenční rovnici vyhovuje řešení
li, = Ci + C2
kde C\ a C2 jsou konstanty, které určíme z okrajových podmínek. Pro n > 2 přibližné řešení -Uj osciluje okolo C±, což je v rozporu s chováním přesného řešení u, rostoucí posloupnost {tiij^o dostaneme jen pro \n\ < 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+i - 2-Uj + iíj_i
Ui-l
= 0.
h2 h Diferenční rovnici vyhovuje řešení
Ui = C1 + C2{1 + k)\
které je rostoucí nezávisle na velikosti n. □
2.3. Metoda konečných objemů
(stručně FVM podle anglického finite volume methoď) 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 xt přiřadíme konečný objem Bl (stručně buňku) takto: pro vnitřní uzly Bl = {xt_i/2, Xi+1/2) , i = 1, 2,..., N — 1, a pro koncové uzly B0 = (0, £1/2), BN = (xn-i/2, ty■ Zřejmě (0,£) = {J^B^ Integrací rovnice (2.20) přes buňku Bl dostaneme bilanční rovnici
J B,t J Bi J Bi
40
Předpokládejme nejdříve, že Bl přísluší vnitřnímu uzlu. Pak dostaneme
X=X- ľxi + l/2 ľxi + l/2
— [pu' — ru]x=x^_\,íi2Jr j        qudx = f dx. (2.32)
Difúzni tok aproximujeme pomocí centrální diference,
Pi-i/2u'(xi_1/2)=Pi_1/2U(Xi) ~^Xl~l} +0(h2),   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í,
i + l/2 ľxi + l/2
qu dx = hqiu{xi) + 0(h3) ,        /        / dx = hjl + 0{k?).
-1/2 ^xi-l/2
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 = í 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 B0 resp. BN. Tak například pro buňku Bq máme
rxi/2 rxi/2
\X=Xi ' '
{pu — ru) I     1/2 + /      qudx = f dx .
IX—xn f I
=x0
I XQ -JX0
Difúzni tok v bodě xi/% aproximujeme centrální diferencí (2.33), konvekční tok v bodě x\j2 aproximujeme pomocí (2.24) resp. (2.28), člen p(0)u'(0) — r(0)ií(0) vyjádříme pomocí okrajové podmínky (2.21a) a integrály spočteme levostrannou obdélníkovou formulí,
xl/2 fXl/2
qu dx = Tjhqou^xo) + 0{h2),        /      f dx = hhfo + 0{h2) .
xo j x0
Zanedbáme-li chyby, dostaneme rovnici (2.25i) resp. (2.29i). Rovnici (2.253) resp. (2.293) odvodíme obdobně z bilanční rovnice (2.31) zapsané pro buňku BN.
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 finite 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 fix) bodem nespojitosti 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 C1(0,£) nazveme testovací, jestliže v = 0 v 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')']vdx a následným užitím okrajové podmínky (2.18b) a vztahu v(0) = 0 obdržíme
fvdx = /  [—(pu)' + qu] v dx = — pu'v\^    + /   [puv' + quv] dx = o Jo x~ Jo
= [a.£u(£) — (3i]v(£) + /   [puv' + 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) = go také rovnost
puv' + quv] dx + oí£u(£)v(£) = / fvdx + (3^v(£) (2.34) o 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 = PC1(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) = go} přípustných řešení. Dále označíme
aiu,v) = I   [puv' + quv] dx + aeu(£)v(£),
L(v) = / fvdx +(3ev(£). Jo
42
Pak úlohu
najít u EW splňující a(u, v) = L(v)   \/v G 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 G C1 {Q,£),r,fe C(0,i) , (2.37a) slabá formulace:      p, q, f e PC(0, £) . (2.37b)
Jestliže p > 0, q > 0, a>e > 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í óu přípustných řešení u G 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 G X | v(0) = v(í) = 0},    W = {v G X | v(0) = g0, v(í) = ge},
fe ŕ (2.38-DD)
a(u,v) = / [puv + quv] dx,    L(v) = / fvdx. Jo Jo
(DN) Okrajové podmínky (2.11a), (2.18b)
V = {v G X\v(0) = 0},   W = {v eX\v(0) = g0},
f* ŕ (2.38-DN)
a(u,v) = / [puv' + quv] dx + a>eu(£)v(£),    L(v) = /  fvdx + (3^v[£).
JO Jo (ND) Okrajové podmínky (2.18a), (2.11b)
V = {v G X\v(£) = 0},    W = {v G X\v(£) = ge}, (2.38-ND)
pt pí a(u, v) =     [puv' + quv] dx + a>ou(0)v(0),    L (v) =      fvdx + (3qv(Q) . Jo Jo
(NN) Okrajové podmínky (2.18a), (2.18b)
V = X,   W = X,
a(u,v) = j [puv' + quv] dx + aou(0)v(0) + aeu(£)v(£), (2 38 NN)
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, ae > 0 a když platí (2.37b). V případě úlohy (2.38-NN) je třeba navíc předpokládat «o > 0 nebo     > 0 nebo q > 0 alespoň na části (0,£).
Diskretizace užitím lineárního prvku. Na intervalu (0, £) zvolíme dělení 0 = xq < x\ < ■ ■ ■ < xn = i a na každé úsečce (xl_i,xl) délky ht = xt — hledáme přibližné řešení Uix) ve tvaru lineárního polynomu procházejícího body [xj_i,iíj_i] a [xl,ul], takže
U(x) = ul_iwl_i(x) + ulwl(x),    kde wl_i(x) =
hi hi
Funkce Uix) 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 wl(x) se jsou tzv. bázové funkce, lineární na každé úsečce {xk-i,Xk) a takové, že
1 pro i = j, 0 pro i 7^ j.
	y			
1 -				
	W0(x)	s i^V wáx)	wn(x) j	' x —1->■
0 = x0
x\
Xi-1
Xi+i xN_i      xN = £
Obr. 5.1. Lineární Lagrangeovy bázové funkce
Úsečku na které je definována lineární funkce určená svými hodnotami
v uzlech a xl5 nazýváme Lagrangeovým lineárním prvkem nebo také elementem.
Délku největšího dílku dělení {xj}^0 označíme jako h, tj. h = maxi<l<N hl. Prostor všech po částech lineárních funkcí (nebo-li lineárních splajnů) označme Xh- Zřejmě Xh C X. Nechť Vh = V n Xh a Wh = W n 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í ah{U,v) = Lh{v)   \/v G Vh-
(2.40)
Přitom index h v ah(U, v) resp. Lh{y) značí, že integrál JQ [pU'v' + qUv] dx v a(ř7, -u) 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(tp) přibližně spočtenou hodnotu J** ^ ipdx, je
N
ah(U, v) = Y2 Ql(pU'v' + qUv) + a0U(x0)v(x0) + a^U{x^)v{x^),
i=i
N
Lhiv) =      Ql(fv) + /30v(x0) + Pev(xN).
i=i
Nechť v(x) = ®íwí{x) £ Vh je libovolná testovací funkce (tj. Gj =
(2.41)
v{Xi) je
libovolné číslo)   a U(x) = Ylf=o^jwj(x) Je MKP řešení (tj. A j = U (x j) = u j). Pak
Jj=0 ^3 3
z (2.40) pro úlohu (2.38-NN) plyne
0 =ah(U,v) - Lh(v) = ah
'N N \ / N \
v3=0 i=0 / \ i=0 /
N
i=0 l]=0
N
Y a^w^WijAj - Lh{w%
(2.42)
GT [K A - F],
kde G = (60, 9i,..., GN)T, K = {kl0}fJ=o pro kl3 = ah(Wj, Wi), A = (A0, AU...,AN) a F = (Fq, Fi, ..., F/v)T pro Fj = Lh(wi). 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 K1 a elementárních vektorů zatížení Fl příslušných elementům Xj), i = 1, 2,. .., N. Pomocí obdélníkové formule vyjádříme
a, - Aj_x e, - eť_i
7{pU'v>) =/^-1/2~* r~L = [©l]TKllAl,
hi
kde
1 -1 -1 1
a   A1 =
A,_i
h,   V-l     1/    "   ~~     V A,
Dále pomocí lichoběžníkové formule obdržíme
Q\qUv) = ^(ft.xei-iAi-i + qíQíAí) = [Gl]TKl2Al,
kde
w2 = \hjq%-x °
0     qi '
45
fi-1
a položíme Kl = Kl1 + Kl2. Nakonec, opět pomocí lichoběžníkové formule, dostaneme
Ql(fv) = lhl(fl.1ei.1 + m) = [@l]TFl,   kde   F = i^^ Z rovnice
0=aft(ř7,W)-Lft(«) =
" N
Ql(pU'v' + gřJ-u) + ctoř7(xo)u(a;o) + a>eU(xN)v(xN)
,i=i
' AT
i=1 V
J^[ei]T[KlAl - Fl] + e0[a0A0 - Po] + eN[aeAN - (3,
i=i
a z rovnice (2.42) tak dostaneme rovnost
AT
GT [KA - F] = ^[ei]T[KlAl _ Fi] + e0[a0A0 - A,] + Ojv^Ajv - fa], (2.44)
i=i
z níž plyne postup, jak pomocí elementárních matic Kl = {klrs\2 s=í1 elementárních vektorů F1 = (F^F^)7" a čísel a0, (30, c^, (3^ sestavit globální matici K a globální vektor F: stačí srovnat členy se stejnými indexy u parametrů G a A (pro určení prvků matice K) nebo jen u parametru G (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.
kh + a0	k1 ^12					Ä + F}
^21	^22 + ^11	k2 ^12				F,1 + F2
	k2 K21	k2   4- P "22   ' "11				F2 + Ff
						
				k1^-1 4- kN "22       ' "11	I.at K12	FN-1 + fat
				^21	^22 + ^	F2N + Pe
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 Xj), i = 1, 2, ..., N, na každém z nich vypočteme elementární matici K1 = {klrs}ls=l a elementární vektor Fl = {Frl}2=1 a koeficienty /c*s resp. F* 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 go 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 a0 a k prvnímu prvku vektoru F přičteme číslo [3q.
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 gi> 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 [3^.
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,..., "Uat_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 = (uq, ui, ..., -uat_i)t v případě okrajových podmínek (2.18a), (2.11b), u = (uq, 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 E 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 fl. 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 nejjednodušší 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 Rd 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 M.d.
Nechť fl je oblast. Prostor funkcí, které jsou v fl spojité spolu se svými derivacemi až do řádu k, značíme Ck(fl). Prostor C°(fl) funkcí spojitých v fl značíme stručně C(fl).
Řekneme, že funkce u je v fl po částech spojitá, jestliže fl je sjednocením uzávěrů konečného počtu navzájem disjunktních podoblastí, tj. fl = {Jflt, flt Pl íl, = 0 pro i 7^ j, a jestliže u je na každé z podoblastí flt spojitá a spojitě prodloužitelná až do hranice, tj. existuje spojitá funkce u% G C(fžj) s vlastností u% = ií,v fl%. Prostor po částech spojitých funkcí značíme PC (Cl). Symbolem PCk(fí) značíme prostor funkcí, které jsou v oblasti fl 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 PC1(fl) je prostor funkcí, které jsou v fl spojité a jejichž první derivace jsou v fl 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 PCiTp) je prostor po částech spojitých funkcí na části Tg C dfl hranice oblasti fl.
48
3.1. Úloha eliptického typu
3.1.1. Formulace úlohy
Buď íž omezená oblast v IR2. O hranici T = <9íž oblasti íž předpokládejme, že je sjednocením uzávěrů dvou navzájem disjunktních částí Ti a r2, tj. T = T1UT2, I^nl^ = 0. Dále n = (ni, n2)T nechť je jednotkový vektor vnější normály hranice a
du du du
— = n • \u = ni— + n2 — 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 f3(x,y) jsou dané funkce. Naším úkolem je určit funkci u(x,y), která uvnitř íl vyhovuje diferenciální rovnici
d /        du \     d /        du \
na hranici T± 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 — /3(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)   vfi, (3.4)
kde Laplaceův operátor A aplikovaný na funkci u je
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, —pS/u 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 x-ová a du/dy y-ová složka vektoru
49
rychlosti, (3 je normálová rychlost s vlastností JT(3á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(Ú) 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šení, které existuje za podmínek v aplikacích běžně splněných.
V dalším budeme předpokládat, že íl je mnohoúhelník, p > 0, q > 0, / jsou po částech spojité v íl, g je spojitá na Ti, a > 0, [3 jsou po částech spojité na r2, a pokud Y = Y2, pak buďto q > 0 na části íl nebo a > 0 na části T2. Za těchto předpokladů existuje jediné slabé řešení u G PC1 (Cl) 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 Ti = T, když íž = (0, £) x (0, £) je čtverec se stranou délky £.
Na íl 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 xt = íh, i = 0,1,..., N, y3 = jh, j = 0,1,..., N. Body [xt, y3], í, j = 0,1,..., N, nazveme uzly sítě. Rovnice (3.4) má být splněna ve všech bodech [x, y] uvnitř íl, musí tedy být také splněna ve všech vnitřních uzlech, tj.
d2u(xl,yJ)     d2u(xl,yJ)       .... , --Q^ŕ---QyT^ = f^yi)^   hJ = l,2,...,N-l. (3.5)
Parciální derivace vyjádříme pomocí centrálních diferencí
d2u{xl, yj) _ -ií(xj_i, yó) + 2u(xl, yj) - u{xl+1, %•)
dx2 h? d2u{xl,yj) _ -u(xj, yj-j) + 2u(xl, y^ - ujx,, yJ+1) dy2 h2
0(h2), (3.6a) 0(h2), (3.6b)
dosadíme do rovnice (3.5) a chybové členy Oih2) zanedbáme. Po vynásobení h2 dostaneme soustavu síťových rovnic
-Ui-1:j - -uM_i + Auíj - ui:j+1 - ui+1:j = h2fi:j,    í,j = l,2,...,N—l, (3.7)
50
ľ = VN
Xi-Í
Obr. 3.1. Síť
kde Uij je aproximace ií(xj,yj) a f^j — Z okrajové podmínky (3.2) dostaneme
gij   pro i = 0 nebo i = N nebo j = 0 nebo j = N, (3-8)
přičemž ^ = g(xl,yJ). Když z (3.8) dosadíme do (3.7) a na levé straně ponecháme pouze členy s neznámými ul01 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					( h2^ -	^goi -	^9io \
-1	4	-1	0	-1	0	0	0	0				h2fi2 ~	^902	
0	-1	4	0	0	-1	0	0	0		Ul3		h2fi3 ~	^ 903 -	V 914
- 1	0	0	4	-1	0	-1	0	0		U21		h2f2i ~	{-	920
0	-1	0	-1	4	-1	0	-1	0		U22	=	h2f22		
0	0	-1	0	-1	4	0	0	-1		U23		h2f23 -	{-	924
0	0	0	-1	0	0	4	-1	0		U31		h2.f3i -	^941 -	V 930
0	0	0	0	-1	0	-1	4	-1		U32		h2f32 -	^942	
0	0	0	0	0	-1	0	-1			\ U33 J		\ h2h3 -	^ 943 -	y 934 J
\
Pravá strana rovnice odpovídající uzlu, který není sousedem hranice, obsahuje pouze člen h2 fij, 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	O	°\		(U1 )		
-I	B	-I .	. O	O	O				F2
o	-I	B .	. O	O	o				F3
o	O	O .	. -I	B	-I		u7v-2		
V o o o
o
B/
\ Uat_i /
(3.10)
kde I je jednotková matice řádu N — l, Oje nulová matice řádu N — 1 a B je třídiagonální matice řádu N — l tvaru
B
/ 4
-1
0
0
V o
-1
4 -1
0
o
o
-1
4
0
o
o o o
-1 o
o 0\
o o o o
4
-1
-1
4/
Vektory Uj = (un, -ul2,..., uhN_i)T, i = 1, 2,..., N — 1, jsou části celkového vektoru u
neznámých, vektory Fj = (Fll5 Fl2,..., FlAr-i) , i = 1, 2,..., N — 1, jsou části celkového vektoru F pravých stran.
Pro homogenní okrajovou podmínku g{x,y) = 0 je Fl3 = h2fl0, pro nehomogenní okrajovou podmínku je v některých rovnicích příspěvek z hranice, přesně
Fij = h2jl3 -
■9o j
Fa = h2 fn + <7jo ) Fn = h2fn + #oi + #io , Fn-i,i = ^2/tv-i,i + g ni + 9n-i,o
2<i<N - 2, 2<j<N-2, FN_1:J = h2fN_1:J +gNj, 2<j<N-2, Fl>Ar_i = /i2/líAr_i + glN , 2 < i < N - 2 , Ft ,N-1 — hz.fi ,n-i + go,N-i + giN i
Fn-i,n-i =     ÍN-l,N-l + gN,N-l + 9N-1,
Matice K soustavy (3.9) má řadu vynikajících vlastností. Některé z nich jsou patrné okamžitě: K je symetrická, diagonálně dominantní, pásová a pětidiagonální. Matice K je pro velké N řídká, neboť počet jejích nenulových prvků je malý: z celkového počtu (N — l)4 je jich nenulových jen 5(iV — l)2 — A(N — 1). Dá se dokázat, že K je pozitivně definitní a tedy regulární.
Jsou-li funkce f a g dostatečně hladké, pak pro chybu metody platí
u(xi,yj) - Uíj = 0(h2) .
(3.11)
Obdélníková oblast. Je-li íž = (0, a) x (0, b) obdélník a na něm chceme zvolit pravidelnou síť, musíme ve směru osy y vybrat obecně jiný krok než ve směru osy x. Nechť tedy N > 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 xt = íh, i = 0,1,..., N, a y3 = 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(xu yj) _ -u(xu yj-i) + 2u{xl, y3) - u(xu yj+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
~uí-i,j + 2iiij — líi+ij    — uí,j-i + 2iiij — líij+i
h2 k2 Je-li řešení u dostatečně hladké, pro chybu platí
= k ■ (3-7')
u
(Xl,yj) - ul0 = Oih2 + k2). (3.1ľ)
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)uij — Pi+l/2,jui+lj ^ ^„ ^
h2
-Pi,j-l/2Ui,j-l + (pv-l/2 +P i,j+l/2)uij ~ Pi,j+l/2ui,j+l
k2
-\- QíjUij j\j .
Přitom index i ± 1/2 znamená, že za argument x dosadíme xl ± \h, a podobně index j ± 1/2 znamená, že za y dosadíme y j ±
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,l/)^T^       =ď(y)u(a,y)-ny),       0<y<b, (3.12e)
Oy x=a
na východní (east) straně obdélníka íl. Splnění rovnice (3.5) požadujeme i v uzlech [xn, yj] na straně x = a. Při diskretizaci —\pux]x{x^,yj), 1 < j < M — 1, postupujeme zcela analogicky jako v jedné dimenzi, tj.
-lpuM?K, V,) = - + 0(h) = (3.13e)
-p-h--+0{h) ■
53
Pomocí (3.13e) a užitím standardní aproximace členu — \puy]y(xN, y3) dostaneme ve vnitřních uzlech [xn, Vj] východní strany x = a, tj. pro j = 1,2,. .., M — 1, rovnici
-Pn-i/2,jUn-i,j + {PN-i/2,j + haf)uNj
'/. »7/ - 11 Ar .• , i mil. »r     , 1
(3.14e)
/i2
-PN,j-l/2^N,j-l + (PAT ,j —1/2 "i" PN,j+l/2)uNj PN,j+l/2uN,j+l
2k2
\qN3uNj - \ fNj + -f3] .
Předpokládejme, že Newtonova okrajová podmínka je předepsána také na severní (north) straně straně íl, tj. že platí
du(x, y)
= an(x)u(x,b) - (3n(x),       0<x<a. (3.12n)
y=b
dy
Podobně jako při odvození (3.13e) dostaneme
n nn   , UiM ~ ^i,M-l
d2u(xl,yM)     a* UlM ~ & + Pi'M-i/2-7-
--=-p-"-+ (3-13n)
Odtud a užitím standardní aproximace členu — [pux]x(xl,yM) dostaneme pro vnitřní uzly horní strany y = b, tj. pro i = 1,2,. .., N — 1, rovnici
—Pi-l/2,mui-l,m + {Pi-1/2,m + Pi+1/2,m)uiM — Pi+l/2,mui+l,m t 1 , x
-2h2-+ (3-Mn)
-Pi,M-l/2Ui,M~l + {Pi,M-í/2 + ka>?)uiM
k2
\q%MU%M = \ fiM + -^fii ■
Pokud je Newtonova okrajová podmínka předepsána současně na východní i severní straně, dostaneme v severovýchodním rohu [xn, Vm] (pomocí (3.13e) pro j = M a (3.13n) pro í = N) rovnici
—Pn-i/2,mUn-i,m + ÍPn-i/2,m + haeM)uNM
(3.14ne)
2h2
—PN,M-1/2UN,M-1 + {PN,M-1/2 + kď^)uNM
2k2
Při diskretizaci Newtonovy okrajové podmínky na ostatních stranách a v rozích postupujeme podobně. Je-li řešení u dostatečně hladké, pro chybu opět platí (3.1ľ).
54
Rovnice s konvekčním členem je tvaru
d í        du \     d í        du \
~dx V (x>y^Q^~ri(x>y^u) ~ ]fy{p(X,y^~dy~ -rÁxiV)u) +q(x,y)u = f (x,y) (3.15)
a odpovídající Newtonova okrajová podmínka je du
—p(x,y) ——h [ti(x, y)ni + r2(x, y)ri2\u = a(x, y)u — (3(x, y)      na ľ2. (3.16) on
Konvekční členy (riu)x a (r2u)y aproximujeme podobně jako v jednorozměrné úloze, viz. kapitola 2.2. Ukažme si to pro člen (riu)x. Omezíme se přitom jen na vnitřní uzel [xl1 y0]. Pomocí centrální diference dostaneme
{riu^x^yj) & ([ri]í+i/2)jií(xí+i/2,%-) - [r^^^j^x^^, yj))/h.
Dalším krokem je aproximace členů u(xl+i/2, y3) a u(xl_i/2, y-j). Protože aproximace obou těchto členů je založena na stejných pravidlech, věnujme se podrobně jen aproximaci členu u(xl+i/2,y-j). Ta závisí na dvourozměrné analogii podmínky (2.26): pokud platí
|/i|[ri]l+i/2j| <pi+i/2,j, (3.17) použijeme centrální aproximaci
u(xi+1/2,yj) « \[u{xi,yj) + u{xi+1,yj)}, (3.18) v opačném případě použijeme jednostrannou upwind aproximaci
{u(xí, yj), |  > 0,
pokud  [ri]l+i/2íJ   < (3.19) u{xi+i,yj), { < 0.
Clen ir2u)y{xl1 yj) aproximujeme obdobně jako člen (riu)x(xl, y3), při rozhodování o hodnotě u(xl,yJ+i/2) se řídíme analogií ke kritériu (3.17), tj. podmínkou
|/c|[r2]líJ+i/2| <pi:j+1/2. (3.20)
Aproximujeme-li všechny konvekční členy pomocí centrální aproximace, pak za předpokladu dostatečné hladkosti řešení u pro chybu platí (3.11'). Jestliže však konvekční členy aproximujeme užitím jednostranné upwind aproximace, řád chyby se o jednotku sníží, pro chybu platí jen
u(xí, yj) - Uíj = 0(h + k) .
Pro dosažení přesnosti řádu Oih2 + k2) je třeba používat přesnější upwind aproximaci konvekčního toku, viz např. [25], [6], [11].
Dirichletova okrajová podmínka. Standardní postup je tento: je-li v uzlu [xl,y]]
předepsána hodnota řešení u(xl,yJ) = glv rovnici pro tento uzel vůbec nesestavujeme a v rovnicích obsahujících ul3 položíme ul3 = giy
55
Podmínky ul3 = gl3 lze vynutit i jinak. Nejdříve sestavíme soustavu rovnic jako kdyby na části Y\ hranice byla předepsána Newtonova okrajová podmínka (3.3) s a = [3 = 0. Následně modifikujeme rovnice příslušné uzlům ležícím na Ti. Předpokládejme, že uzlu [xí, y-j] s předepsanou hodnotou u(xl, y3) = gl3 přísluší r-tá rovnice soustavy Ku = F. Pak provedeme krr := K>krr,Fr := K,krrglv kde k je velké číslo, např. k = 1020. To způsobí, ze mimodoagonální koeficienty v r-té rovnici budou oproti velkému diagonálnímu koeficientu prakticky zanedbatelné, takže r-tá rovnice nabude přibližně tvaru K,krrur = Kkrrgl0 nebo-li ur = gl3, což jsme potřebovali zajistit.
3.1.3. Metoda konečných objemů
Metodu vysvětlíme pro konvekčně-difúzní úlohu (3.15), (3.2), (3.16).
Pravidelná síť. Uvažujme nejdříve případ, kdy íž = (0, a) x (0, b) je obdélník. Na íž zvolme uzly [xu     = [ih,jk], i = 0,1,..., N, j = 0,1,..., M, kde h = a/N, k = b/M.
b = yM
Vj+l/2
V]
26-1/2
yo
, h/2 , i*—■-n
	i	i i i i	i	i i i	Bnm
------		i i i i -----«----- i i i i	—:------	i i i i i	------1
------		i		i i i i i	
					
		+ 1 1 1			
	i	1 1 1	i	i i	
		w i i i i		T 1 1	
	i	i i i i	!	1 1 1	
k/2
XQ
Xi-l
-1/2
^1+1/2 Xl+1
Obr. 3.2. Vnitřní buňka B^, hraniční buňka Bnj a rohová buňka B nm
Obdélník
Bij = [(^-1/2,^+1/2) x {vj-i/2, yj+1/2)] HO,       i = 0,1,..., N,   j = 0,1,..., M,
nazveme konečným objemem (stručně bunkou), viz Obr. 3.2. Integrací diferenciální rovnice (3.15) přes buňku Bl3 dostaneme bilanční rovnici
(pUy
r2u
I Bi,
qu] dxdy = / fdxdy.
J Ba
(3.21)
56
Uvažujme nejdříve případ, kdy Bl3 je vnitřní buňka, tj. když 0<i<iVaO<j< M. Pak      = {xt_i/2, a^+1/2) x {vj-1/2, yj+1/2), viz Obr.3.2. Integrací per-partes obdržíme
Vj+1/2 rvj+i/2
[pux - riu] (xj+i/2, y) dy + /        [pux - nu] {xi_1/2,y) dy-
Vj-1/2 -hj-1/2
+1/2 r^i+1/2
[puy - r2u] (x, yj+1/2) dx + [puy - r2u] (x, yj-1/2) dx+ (3.22)
Xi-1/2 Jxi-1/2
qudxdy = i    f dx dy .
J Bij
Ukážeme si, jak provést aproximaci prvního členu v (3.22), zbývající jednoduché integrály se aproximují podobně. Pomocí obdélníkové formule dostaneme
ľVj + 1/2
- /        [pux - nu] (xi+1/2,y) dy t=t k[-pi+1/2dux(xi+1/2,yj) + [^+1/2,^(^+1/2, %')] > Jyj-1/2
derivaci ux(xl+i/2,yJ) aproximujeme pomocí centrální diference
ux{xi+1/2,yj) « [u(xi+1,yj) -u(xi,yj)]/h
a u(xl+i/2,y3) aproximujeme stejně jako v diferenční metodě, viz (3.17)-(3.19). Dvojné integrály v (3.22) aproximujeme pomocí součinové obdélníkové formule:
qudxdy    hkqlJu(xl, y j) ,        /   f dxdy & hkfij .
Bij
Po dosazení do (3.22) dospějeme ke stejné rovnici, jakou bychom dostali pomocí diferenční metody.
Věnujme se nyní případu, kdy uzel [xl,yJ] leží na hranici íž, tj. když i = 0 nebo í = N nebo j = 0 nebo j = M. Jestliže uzel [xl,yJ] leží na části T± hranice a je v něm tedy předepsána hodnota u(xl,yJ) = g(xl,yJ), pak bilanční rovnici (3.21) pro takový uzel vůbec nesestavujeme. Uvažujme tedy případ, kdy uzel [xl,yJ] je vnitřním bodem hranice T2. Pro konkrétnost budeme stejně jako v kapitole 3.1.2 uvažovat vnitřní uzel pravé strany obdélníka íž, tj. uzel [xw,yj], 0 < j < M, pro který Bnj = {xn-i/2,xn) x (y?-i/2iyj+1/2)i viz Obr. 3.2. Bilanční rovnici (3.21) upravíme integrací per-partes,
Vj+1/2 rvj+i/2
[pux-r1u](xN,y)dy + /        [pux - rľu] {xN_1/2, y) dy-Vj-i/2 Jyj-1/2
■xN pxn
[puy - r2u] (x,yj+1/2) dx + /        [puy - r2u] (x, y j-1/2) dx+ (3.23) XN-1/2 JxN-1/2
qudxdy = /    / dx dy .
J Bij
První integrál v (3.23) upravíme užitím Newtonovy okrajové podmínky (3.16),
yj+1/2 rVj+i/2
[pux - r±u] (xN, y) dy = /        [ď(y)u(xN, y) - f3r{y)] dy , yj-1/2 Jyj-1/2
57
a pak použijeme obdélníkovou formuli, takže
ľVj + l/2
I [pux - nu] (xN, y) dy « /cfo^i^, -yj-i/2
Druhý integrál v (3.23) aproximujeme stejně jako obdobný integrál ve vnitřní buňce. Třetí integrál v (3.23) aproximujeme užitím jednostranné obdélníkové formule
derivaci Uy(xN,yJ+i/2) aproximujeme centrální diferencí
uy{xN,yj+1/2) ~ [u(xN,yj+1) - u(xN, yó)]/h,
aproximace členu u(xN,yj+i/2) závisí na poměru velikostí konvekčního členu [t2]n,j+i/2 a difúzního členu Pn,j+i/2, yiz (3.20): jestliže
|^|[r2]7V,j+l/2| < PN,j+l/2 i
použijeme centrální aproximaci
u(xN,yj+1/2) ~ \[u{xN,yj) + u(xN,yj+1)] , v opačném případě použijeme jednostrannou upwind aproximaci,
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 íž 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. triangulaci oblasti Q, v metodě konečných objemů označovanou jako primární síť, viz Obr. 3.3. Ke každému uzlu Pt přiřadíme buňku Bl. Sestavíme ji z přilehlých částí Cljk všech trojúhelníků s vrcholem Pl. 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 Pt (£ dfl. Pro jednoduchost se omezme na tento případ. Nechť Tljk je trojúhelník s vrcholy Pl5 P3 a Pk. Označme Pl3 střed strany PiP3, Pik střed strany PlP\z a Pljk těžiště trojúhelníka Tljk- Pak do buňky Bl zahrneme čtyřúhelník
\h\puy - r2u](xN,y:i+1/2)
Č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
58
Obr. 3.6. Duální síť
Cijk s vrcholy Pl, Pl3, Pl3k a P^, viz Obr. 3.4. Tento postup opakujeme pro všechny trojúhelníky Tljk, které obsahují uzel Pl, a sjednocením přilehlých částí Cljk dostaneme buňku Bi, 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],
d P    dQ \ f — + — ) dxdy = /   (Pni + Qn2) ds,
kde P = —{pux — riu) a Q = —{puy — r2u), a dostaneme analog rovnice (3.22) du
(3.24)
dBi
-P
dn
r • n )u
ds + /   qudxdy = \   f dxdy,
'Bi J Bi
(3.22')
kde r • n = r\n\ + r2n2 a n = (ni, n2)T je jednotkový vektor vnější normály hranice dBl buňky -Bj.
Pí
P	ik ——^
	. Pij k /
&ij k	
P \ ri3 \	
Pk
Obr. 3.4. Cijk = BinTijk
Obr. 3.5. Vnitřní buňka Bi
59
Klíčová je aproximace křivkového integrálu. Hranici dBl vyjádříme jako sjednocení společných částí dB%] = dB% n dBj buňky B% a sousedních buněk B j, tj. <9Pj = [J. dB%]. Protože JaB = ^2j JdB--i s^ačí popsat aproximaci jen na dBiy To lze provést třeba takto:
dBi
du ^ dn
r • n )u
ds
\dBy
p[- 13 > \-p-p\
kde |<9Pv,-| je délka dBl3, |pP,| je délka úsečky P%P3,
r • n)ťj- = r(Pťj) ■ nv = + r2(P>1/   a   n, = (n1^,nf)T = (P,- - P)/|PP
je jednotkový vektor ve směru vektoru PjPj. Zbývá aproximovat hodnotu u(Píj) ve středu úsečky P%P3. Jestliže je dominantní difúze, tj. když
\\PlPJ\\(Y-^)lJ\<V{PlJ),
zvolíme aritmetický průměr
^p^iMPO + ^p,-)),
v opačném případě užijeme upwind aproximaci,
í u(Pi),        .   .   .     .      í >0,
«(Pj
pokud (r • n)
< 0.
Dvojné integrály aproximujeme jako součin obsahu \Bl3 \ buňky Bl3 a hodnoty integrandu v bodu Pj, tj.
/ quáxáy^\Bl\q{Pl)u{Pl), Í f dxdy « ^(/(P) .
Více podrobností o metodě konečných objemů, včetně přesnější varianty upwind aproximace konvekčního členu, lze najít např. v [11].
3.1.4. Metoda konečných prvků
Metodu vysvětlíme pro úlohu (3.1)-(3.3).
Slabé řešení. Stejně jako v kapitole 2.4 úlohu převedeme na tvar vhodný pro nasazení MKP, tj. odvodíme slabou formulaci naší úlohy. K tomu účelu násobíme rovnici (3.1) testovací funkcí v G C1 (Ú) s vlastností v = 0 na Ti a integrujeme přes oblast ŕž, tj. provedeme
/ [(pUx)x + (Puy)y]v dxdy + / quvdxdy= / f v dxdy.
(3.25)
60
První člen na levé straně upravíme pomocí Gaussovy-Ostrogradského formule (3.24) (v níž zvolíme P = puxv, Q = puyv) na tvar
(pux)x + (puy)y]v dxdy = - /   p[uxni + uyn2]v ds + / p[uxvx + uyvy] dxdy . o Jao Jo
Křivkový integrál dále upravíme: využijeme toho, že v = 0 na T± a že na T2 platí okrajová podmínka (3.3). Tak dostaneme
f    du f
p[uxni + uyn2]v ds = — /  p —— v ds = /  [au — (3\v ds. ao Jr2   dn JT2
Dosadíme-li z posledních dvou rovností do (3.25) vidíme, že řešení u úlohy (3.1)-(3.3) splňuje rovnici
p{uxvx + uyVy) + quv\dxdy + j  auvds= / f v dxdy + /  (3v ds (3.26)
pro každou funkci v G C1(Q) s vlastností v = 0 na IV Rovnice (3.26) má zřejmě smysl i v případě, když funkce u a v jsou jen z prostoru X = PC1(Ů). Testovací funkce tedy volíme z prostoru V = {v G X \ v = 0 na Ťi} testovacích funkcí a řešení u z množiny W = {v G X | v = g na Ťi} přípustných řešení. Označíme-li
a(u, v) =     [p(uxvx + uyVy) + quv] dxdy + /   auv ds , Jo Jr2
(3.27)
L (v) = / f v dxdy + / (3vds, Jo Jt2
pak úlohu
najít u EW splňující a(u, v) = L(v)   \/v G V (3.28)
nazveme slabou formulací úlohy (3.1)-(3.3) a funkci u nazveme slabým řešením.
Diskretizace. Omezíme se na případ, že Cl je mnohoúhelník. Cl vyjádříme jako sjednocení konečného počtu uzavřených trojúhelníků Te, z nichž každé dva jsou buďto disjunktní nebo mají společný vrchol nebo společnou stranu, viz Obr. 3.3. Množinu 7 všech trojúhelníků nazveme triangulací oblasti Cl. Trojúhelníky budeme značit Tl5 T2,..., TNt. Vrcholy trojúhelníků budeme nazývat uzly, značíme je Pi, P2,..., Pm- Předpokládejme, že společné body částí Ťi a T2 hranice V jsou uzly tringulace 7. Množinu stran (trojúhelníků T E 7), jejichž sjednocení je T2, označíme §. Strany budeme značit S±, S2,..., Sns.
Funkci, která je v Cl spojitá a která je na každém trojúhelníku T E 7 lineární, nazveme spojitou po částech lineární funkcí. Každá taková funkce v(x,y) je jednoznačně určena svými hodnotami v(xt, y i) ve vrcholech Pt = [xt, yt] triangulace 7. Prostor všech spojitých po částech lineárních funkcí označíme X^. Speciálním případem funkcí z Xu jsou bázové funkce w^x, y), které jsou v Pl rovny jedné a v ostatních uzlech jsou rovny nule, tj.
61
0
pro
i = J, i + j-
Obr. 3.7. Bázová funkce Wi(x,y)
Každá funkce v E 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, V) = Y1 v(xn Ví)wí{x, y).
1=1
Dále definujme prostor testovacích funkcí
Vh = {v E Xh | v{xu yi) = 0 VP, E fi}
a množinu přípustných řešení
Wh = {v E Xh | v(xí, yi) = g(xi, y i) VPj E 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 diskrétní slabé formulace
najít U E Wh splňující ah(U, v) = Lh{v)   \/v E Vh ,
(3.29)
kde
ah(U,v) =       [QTe(p[Uxvx + UyVy\) + QT*{qUv)] + ^ Qs'{olUv) ,
TeeT
Lh(v) = Y QTeW + Y QSe(Pv) ■
5eeS
(3.30)
TeeT
SeeS
Přitom symbolem QTe(ip) jsme označili kvadraturní formuli pro výpočet JT ipdxdy na trojúhelníku Te a symbolem QSe(<p) 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(tp) = \T\<p(P0), (3.31)
kde \T\ je plocha trojúhelníka T, Po = ^(Pi + Pí + P'i) Je těžiště T a Pi, P2, P3 jsou vrcholy T;
2) členy qUv a f v integrujeme na trojúhelníku T formulí
QT(ip) = §|T| [^(Px) + <p(P2) + <p(P3)] ; (3.32)
3) členy aUv a [3v integrujeme na straně S lichoběžníkovou formulí
Qs(<p) = ±\S\[<p(P1)+<p(P2)], (3.33) kde 15*1 je délka strany S a Pl5 P2 jsou koncové body S.
62
Formule (3.31), (3.32) a (3.33) jsou řádu 1, tj. jsou přesné, když ip je polynom stupně 1.
Předpokládejme, že uzly jsou očíslovány tak, že Pi, P2,..., Pn leží buďto uvnitř oblasti íl nebo uvnitř hranice T 2 a že Pn+i, Pn+2, ■ ■ ■, Pm leží na hranici Y\. To není žádné omezení, neboť uzly lze vždy přečíslovat tak, aby tento předpoklad byl splněn. Pak
N M N
U(x,y) = ^ Ajwj(x,y) +  ^2 9j^j(x,y),       v(x, y) = ^ Qlwl(x, y) , (3.34)
j=l j=N+l i=l
kde Aj = U(Pj), (/-, = g(Pj) a0t = v (Pí). Dosazením do (3.29) obdržíme
(N M N \ / N \
AjWj + ^2 9jwv ^2 ®iWi) ~Lk [ ®iWi)= 3=1 J=N+1 i=l / \ i=l )
N N N
i=l       J—1 i=l
QT [K A - F],
M
Lh(wi) -  ^2 ah(wj,Wi)gi
3=N+1
(3.35)
kde G = (6i, 62,..., QN)T, K = {h,}fJ=1 pro % = ah(Wj, Wi), A = (Au A2,..., AN)T
a F = (Fi, F2,..., FN)T pro F, = Lh(wi) - Y!^=n+i ah(wv w%)gr Protože O 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ě deŕmitní. 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 nejdelší strana trojúhelníků triangulace T.
Algoritmus. Matici K a vektor F sestavíme pomocí elementárních matic KTe a elementárních vektorů FTe pro Te e 7 a elementárních matic K.Se a elementárních vektorů FSe pro Se £ §. 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 P£ = [x\, yf], P2e = [x%, y2] a -F3 = [x\, yf]-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\w\{x, y) + Ae2we2(x, y) + AB3wB3(x, y),
te
v(x, y)    = ®\w{{x, y) + Qe2we2(x, y) + ®%w%{x, y),
te
63
(3.38)
kde Aer = U(Pre), Oer = v(Pre) a kde
wer(x,y) = aerx + bery + cer (3.39) je bázová funkce příslušná k uzlu Pre, r = 1, 2, 3. Z (3.38) a (3.39) dostaneme
u„
BeA
AI
w
e
Te
= BeeJ
/ef\
W
kde   Be =
b\ b\
jsou vektory parametrů na trojúhelníku Te. Koeficienty aer a ^ lze snadno spočítat z rovnic
nebo-li maticově
( x\   y{   1 \ Ví 1
V
pro
/ a?
6f 6|
r = s , r/s,
a3 \
61
^ z/I 1 J
Označíme-li
( x\   y{ 1\
x\   ye2 1
7
r,s = 1,2,3,
/l  0 0\
0 1 o 0  0 1
J
De
pak   Be = [De
\ 4  Vl  1 /
-i
[1:2,1:3] '
tj. matice Be je rovna prvním dvěma řádkům matice inverzní k matici De. Užitím kvadraturní formule (3.31) tak obdržíme
QTe(p[Uxvx + Uyvy]) = [eTe]T[Be]T|Te|p(P0e)BeATe = [0T«]TKT-1ATe, kde elementární matice
K^1 = |Te|KPo)[Be]TBe, jo(P0e) je hodnota funkce p v těžišti P0e = |(-Pf + P2e + ^í) trojúhelníka Te a
|Te| = i|det(De)|
je plocha trojúhelníka Te.
Užitím kvadraturní formule (3.32) dostaneme
QT*(qUv) 4|Te|[g(Pf)tf(PfMP^ =
l\Te\[q(P^eiAl + q(PI)eiAl+q(Pi)eiAl] = [@Te]TKe'2 ATe,
(3.40)
(3.41)
(3.42)
64
kde elementární matice
o
TfTeč _ l\rp
3\1e
V o
o o
o
q(PÍ)J
(3.43)
Celkem tak
QTe(p[UXVx   + UyVy])-
kde
-Te,2
;(9£/t>) = [eTe]TKTeATe
KTe = KTe>1
(3.44)
(3.45)
je elementární matice na trojúhelníku Te.
Užitím kvadraturní formule (3.32) obdržíme
(fv) =i\Te\[f(PZ)v(P?) + / (P2>(P2e) + f(PI)v(Pi)] =
i\Te\[f(PZ)Qt + f(PŠ)Qe2 + f(PŠ)Ql] = [&Te]TFTe
kde
■pTe  _ liji
/(Ae)
(3.46)
(3.47)
je elementární vektor na trojúhelníku Te.
Elementární matice a elementární vektor na straně. Nechť SI DE je tabulka typu Ns 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 P-f = [x\,y^\ a P2e = [x2,y2]. Pro uzel Pf, 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). P^ a Pl 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)     =Ae1wl(x,y) + AZwl(x,y),
Se
v(x,y)     = Qe1we1(x,y) + Qe2we2(x,y) ,
Se
(3.48)
kde Aer = U(P^), Qer = v(P^) a kde wer{x, y) je bázová funkce příslušná k uzlu Pre, r = 1,2. Užitím formule (3.33) obdržíme
kde
f*(aUv) =±\Se\[a(Pl)U(P?)v(Pl) + a(P2e)U(P2e)v(P2e)] =
l\Se\[a(Pl)QlAl + a(P2e)ee2Ae2] = [QSe]TKSe ASe,
K5e = hsP.
a(P
(3.49)
(3.50)
65
je elementární matice na straně Se a
jsou vektory parametrů na straně Se. Podobně odvodíme
Qs*(f3v) =i|5e|[/3(P1>(P1e) + /3(P2>(P2e)] =
!|Se|[/?(Pie)e? + /3(P2e)6^] = [e5e]TF^, (3.51)
kde
je elementární vektor na straně
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 E Vu platí
0 = ah(U, v) - Lh{v) = QT [KA - F] = =      [®Te]T [KTe ATe - FTe] +      [&Se]T [K5e A5e - F
-,&]. (3-53)
TeeT 5eeS
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, K.Se a lokálních vektorů FTe, FSe. 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 E T sestavíme elementární matici KTe = {kJl}fs=1, viz (3.45), a elementární vektor FTe = {Fje}l=í, viz (3.47).
Postupně pro r, s = 1, 2, 3 položíme i = ELEM(e, r), j = ELEM(e, s) a
{i < N a j < N, i < N a j > N, provedeme i < N,
3) Postupně pro každou stranu Se E § sestavíme elementární matici K.Se = {kfs}2.s=l, viz (3.50), a elementární vektor FSe = {Fr5e}y=1, viz (3.52).
Postupně pro r = 1, 2 položíme í = SIDE(e, r) a
pokud i < N,   provedeme ku^—ku + k^, Fl^Fl + F^e.
66
3.2. Úloha parabolického typu
Formulace úlohy. Hledáme funkci u(x,t) definovanou pro x G {0,£), t G (0,T), která vyhovuje diferenciální rovnici
du     d í      du \
c^m~^[p^^)+q^u = f^,t^ xe(0,£), te(0,T), (3.54) Dirichletovým okrajovým podmínkám
u{0, t) = g0(t), u{£, t) = ge(t), t G (0, T), (3.55) nebo Newtonovým okrajovým podmínkám
P(°)dUg^ = a0u(0) -/30(t),
te(0,T), (3.56)
Q = n.n(í t) - /?,ít)
dx
a počáteční podmínce
u{x,0) = (p(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ární vedení 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, g0, g^, f30, (3e a konstanty a0, 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 (3o(t) = aou^t), [3^{t) = a^u^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, go, gg, (3q, (3e a ip jsou dostatečně hladké a že jsou splněny podmínky nezápornosti c > 0, p > 0, q > 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
<p(0) = g0(0),   <p(£) = 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, g^, f30, (3e 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, £) s krokem h = £/N, tj. x% = ih, i = 0,1,..., N. Budeme požadovat, aby rovnice (3.54) byla splněna ve vnitřních uzlech xt, i = 1, 2,..., N — 1, tj.
c(xí)
du(xl, t)
dt
d ( du
dx \^ dx
(xí, t) + q(xl)u(xl, t) = f(xi, t).
Derivaci podle x vyjádříme pomocí diference analogicky jako v (2.14), tj.
" d	í du\	i    ^ 1		_l/2lí(Xj_i,r) +	
dx	dx J	(Xi,t) =	—Pi-		
		+ [Pi-	l/2 +	Pi+1/2]u(xi,t) - pi+1/2u(xi+1,t)	+ 0(h2
Zanedbáme-li chybu, dostaneme soustavu rovnic
CiŮi(t)
1
Ji?
-pj_i/2iíj_i(r) + {pi-i/2 +p1+1/2 + h2qi)ui(t)-
(3.58)
pl+1/2ul+1(t) =fi{t),
1,2,..., N-1, te{0,T),
v níž Ui(t) je aproximace u(xl, t), ů^t) = —-^-^ je aproximace dtu(xl, t) a /i(ŕ) = f(xl, t). Z okrajových podmínek (3.55) máme
M*) = 9o(t),       uN(t) = ge(t),       t G (0, T), (3.59)
což dosadíme do (3.58) pro í = Iělí = N — 1. Z počáteční podmínky obdržíme
Ui{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 u^t) definovanými na N — 1 úsečkách x = x%, i = 1, 2,..., N — 1, 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á semidiskretizace.
9o(t)
Ul+Í(t)
Ui(t)
Ui-!(t)
9e(t)
XQ 2^-1      Xl      Xl+l XN
Obr. 3.8. Metoda přímek
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(ŕ),   re(0,T),      u(0) = cp,
(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(r) je tzv. vektor tepelných zdrojů, u(r) = (ui(t), «2 (ŕ), ..., uN_i(t))T je vektor neznámých funkcí a cp = (ip(xi), (pfaz),..., ipixN_i))T je vektor počátečních teplot.
Newtonova okrajová podmínka. Protože Newtonovy okrajové podmínky se v úloze vedení tepla používají nejčastěji, uvedeme si rovnice reprezentující diskterizaci Newtonovy okrajové podmínky vzhledem k proměnné x v levém resp. pravém koncovém bodu intervalu (0, £). Bude-li tedy Newtonova okrajová podmínka předepsána v bodě x = 0, zařadíme před rovnice (3.58) jako první rovnici
^c0ú0(t) + -\ [(p1/2 + ha0 + \h2q0)u0{t) - p1/2ui(t)] = ^-(30 + ^/o(*) ,
(3.58a)
a bude-li Newtonova okrajová podmínka předepsána v bodě x = £, přidáme za poslední z rovnic (3.58) rovnici
^cNůN(t) + -\ [-pAr_i/2Wjv-i(t) + (pn + hae + \h2qN)uN{ť)] = ^/3e + ^fN(t) . (3.58b)
Časová diskretizace. Prozkoumejme nejdříve charakter počáteční úlohy (3.61) za zjednodušujících předpokladů. Pro c = p = l,q = f = 0, u(0, t) = u(£, t) = 0, lze úlohu (3.61) zapsat ve tvaru
ů = Au,   re(0,T),       u(0) = <p,
(3.62)
kde
1
"h?
( 2
-1 0
V
o o
-1
2 -1
0 0
0
-1 2
0 0
-1 0
0 0 0
2 -1
0\
0
o
-1
2
(3.63)
Vlastní čísla matice A 4N2
7TI
£2
sin
2N '
i = \,2,...,N-\,
jsou reálná záporná a max, |Aj| —> oo pro N —> oo. Podle (1.32) je tedy problém (3.62) pro velké N tuhý. Totéž platí také pro úlohu (3.61), tj. jde o tuhý problém. Počáteční problém (3.61) je tedy vhodné řešit metodou, jejíž oblast absolutní stability obsahuje celou
69
zápornou poloosu komplexní roviny. Metody s touto vlastností se nazývají Ao-stabilní. Patří mezi ně všechny metody doporučené v kapitole 1.5 pro řešení tuhých problémů, zejména tedy implicitní Eulerova metoda, lichoběžníková metoda nebo metody zpětného derivování. V Matlabu lze použít některý z programů ode23s, ode23t, ode23tb a odel5s. Při řešení úloh vedení tepla se často používá metoda časové diskretizace známá jako
Theta metoda, kterou v následujícím textu stručně popíšeme. Nechť tedy 0 = ro < ti < • • • < íq = T je dělení intervalu (0, T), rn = tn+i — tn je délka kroku a u™ je přibližné řešení v čase tn, tj. uf ~ Ui(tn) ~ u(xl,tn). Nechť 6 je pevně zvolené číslo z intervalu (0,1). Označme tn+e = tn + 6rn. Rovnici (3.61) zapíšeme pro t = tn+e a dostaneme
Cún+d + Kun+d = Fn+e , (3.64)
kde Fn+e = F(tn+g), un+d = u(tn+g), ún+d = ú(tn+d). Předpokládejme, že u(í) je na intervalu {tn,tn+i) lineární funkce, tj.
u(í) = un+t-^(un+1 -un).
un+l _ un
Pak iin+d =-a un+d = (1 — 9)un -\-6un+1. Po dosazením do (3.64) a malé úpravě
dostaneme 6-metodu
(C + rn#kk+1 = (C - r„(l - 6)K)un + rnFn+e . (3.65)
Snadno ověříme, že pro | < 6 < 1 je ^-metoda A-stabilní a tedy také ^o-stabilní, a že pro 0 < 6 < | je oblast absolutní stability ^-metody omezená a interval absolutní stability je interval {2/(26 — 1),0). Pokud jde o přesnost, ^-metoda je řádu 1 pro 6 7^ \ a řádu 2 pro 6 = \. Všimněte si, že z ^-metody dostaneme pro 6 = 0 EE metodu, pro 6 = 1 IE metodu a pro 6 = 1 TR metodu. Metoda (3.64) pro 6 = i je známa jako Crankova-Nicolsonova metoda. Často se zapisuje v analogickém tvaru
(C + irnkk+1 = (C - ir„kk + \Tn{Fn + F™+1). (3.66)
Přesnost i stabilita metody (3.65) pro 6 = | a metody (3.66) je stejná.
Matice C + Tn6J£ soustavy (3.65) nezávisí na čase a je pozitivně definitní. Soustavu (3.65) je proto účelné řešit pomocí Choleského rozkladu, který stačí provést jen jednou.
Nelineární úloha vedení tepla vznikne tehdy, když některé z funkcí c, p, q, f, a, [3 závisejí na teplotě u. Odpovídající počáteční úlohu
C(u)ú = F(í,u) -K(u)u,   te(0,T),      u(0) = (p, (3.61')
lze v Matlabu vyřešit programem ode23t nebo odel5s. Řešení jednodimenzionálního parabolického problému, a to i nelineárního, umožňuje matlabovský program pdepe.
70
3.3. Úloha hyperbolického typu
Formulace úlohy. Hledáme funkci u(x,t) definovanou pro x G (0, £), t E (0, T), která vyhovuje diferenciální rovnici
p(x& + c(x)iz ~ i- (p(x)it) +v(x)u = /(*>ŕ)'   x G (°. *)> ŕ G (°.T)' (3-67)
dt2 dt    dx \     dx r
Dirichletovým okrajovým podmínkám
u(0,t) = g0(t),       u(£,t) = ge(t),       t G (0,T), (3.68
nebo Newtonovým okrajovým podmínkám
te{0,T), (3.69)
P(°)dUg^ =a0u(0,t)- /30(t),
-p{e)^^ = atu{i,ť)-t3t(ť), a počátečním podmínkám
u(x, 0) = <p(x), °^ = ip(x),       x G (0, £). (3.70)
Možný je také případ, kdy na jednom okraji je předepsána Dirichletova podmínka a na druhém Newtonova. Tato úloha vyjadřuje příčné kmitaní tenké struny (nebo podélné kmitaní prutu) délky £. Proměnná x je prostorová, t má význam času, u(x, t) je deformace (tj. příčná výchylka pro strunu nebo podélné posunutí v případě prutu) v bodě x a v čase t, p je hustota, c udává tlumení (pro c > 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. reprezentuje tuhost pružné podpory a [3q resp. [3^ 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, g^ f30, (3g, ip a ip jsou dostatečně hladké, že jsou splněny podmínky nezápornosti p > 0, p > 0, q > 0, «o > 0, > 0 a že platí podmínky souhlasu
V(0)=g0(0),   V(0)=^(0), <p(t)=9e(0), = 9e(t),
vyjadřující soulad počátečních a Dirichletových okrajových podmínek. Pak má úloha (3.67) - (3.70) jediné klasické řešení.
Podmínky souhlasu v aplikacích někdy nebývají splněny. Také funkce p, c, p, q, f, go, ge, /3o, Pí 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 ovšem 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.67) - (3.70) je okrajovou úlohou druhého řádu vzhledem k proměnné x a počáteční úlohou druhého řádu vzhledem k proměnné t. Při její diskretizaci proto opět použijeme postupy z prvních dvou kapitol.
71
Prostorová diskretizace se provede stejně jako u parabolického problému metodou
přímek, viz kapitola 3.2. Opět předpokládáme rovnoměrné dělení a Dirichletovu okrajovou
t-. , ,   ,     \ > dií j (ŕ) A2uAt)
podmínku. Pro Ui(t) aproximující u{xl,t) při označeni ůi(t) = ———, ůi(t) = ———
dostaneme soustavu rovnic
PiUiit) + Ciůiit) + ^(-^-1/2^-1 (ŕ) + [Pi-i/2 +P1+1/2 + h2qi]ui(t)- (3.71)
- pi+1/2ui+1(tfj =fi(t),       i = l,2,...,N-l, te(0,T). Z okrajových podmínek (3.68) máme
M*) = 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
Ui(0) = <p(xi),      ůi(0) = ip(xi),      i = 0,1, ...,N-1. (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(r),   te{0,T),      u(0) = (p,   ů(0) = tp, (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 (ŕ) je tzv. vektor (vnějšího) zatížení, cp = (ip(xi), (p(x2),•••, tp(xN_i))T, xp = (ifj(xi), /4,{x2)i ■ ■ ■, ^(%-i))T jsou vektory počátečních hodnot a u(r) = (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(r),   te{0,T),      w(0) = k, (3.75)
kde
R^(Ó  £)-=(u)'S=(£ c)-G(t»=(Fw)-«=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(£, t) = 0, je M = I, K = -A a C = cl, kde A je definována předpisem (3.63). Rovnice (3.75) se tak zjednoduší, dostaneme
w = Pw,       kde   P = (a   -cl) ' (3'76) Dá se ukázat, že vlastní čísla matice P
2N   . 7TÍ
Xl = -\c± \j\c2 - íú\,   kde   ^ = ~sin^y'   * = 1, 2,..., iV - 1.
72
Zřejmě max, |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 oscilatoricky tuhý. 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 + ÍTnS)w™+1 = (R - ±rnS)wn + \rn{Gn + Gn+1), (3.77)
viz (3.66). Pro efektivní výpočet složek ura+1 a iin+1 vektoru wra+1 je vhodné soustavu rovnic (3.77) zapsat po složkách, tj.
Mún+1 + ir„[Kun+1 + Cún+1] = Mtin - ±r„[Kun + Cún] + \rn(Fn+1 + Fn). Z prvé rovnice vyjádříme ůra+1, dosadíme do druhé a obdržíme rovnici pro výpočet ura+1: Kun+1 = F,      kde (3.78)
K = 4m + — C + K ,      F = f 4m + —C - K | un + —Mú" + Fn + Fn+1.
Až vypočítáme ura+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
^ + a(x,t)^ = f(x,t),   x E (OJ), ŕG(0,T), (3.81)
73
t E (O, T),
(3.82)
okrajovým podmínkám
u(0,ť) = g0(ť) , a(0,í)>0, pokud
u(£,t) = ge(t), a(£,t)<0,
a počáteční podmínce
u(x,0) = <p(x). (3.83)
Nechť Q = (0, £) x (0, T) je obdélník, v němž hledáme řešení, a dQ+ 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 | ŕ = 0 nebo x=0 (pokud a(0,í)>0) nebo x = £ (pokud a(£,ŕ)<0)}.
Pro každý bod Po = [xo, to] E 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)
Křivka x (t) se nazývá charakteristika příslušná bodu [xo,ŕo]- Předpokládejme, že charakteristika protíná dQ+ v jediném bodě P0* = [x^ŕo]- Tomuto bodu říkáme pata charakteristiky. Podle pravidla o derivování složené funkce
du(x(t),t)     du(x(t), t) dx(t) t du(x(t),t)     du(x(t),t) ^ ^du(x(t),t)
dt dx        dt dt dt
Úlohu (3.81)-(3.83) proto můžeme na charakteristice x(t) zapsat ve tvaru
ip(x*0)   pro ť0 = 0 ,
dx
du(x(t), t) ďí
= f(x(t),t), te(ť0,t0),   u(x(ť0),ť0) = { g0(ť0)  pro x*0 = 0, (3.85)
9í(t*o)   pro x% = £.
Předpokládejme, že funkce a, f, go, ge a ip jsou spojité, že funkce a(0,t) i a(£,t) nemění znaménko a že okrajová podmínka je kompatibilní s počáteční podmínkou, tj. platí
V?(0) = g0(0) (pokud a(0, 0) > 0),    <p(£) = 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čni 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 konvekčně-difúzni rovnici
<9<2H
•í		
P2(x2,t2) A		
	Pi(xiM)/	
r P£(x'2, Í2)		
	/p{(*\A)	X
dQ+
Obr. 3.9. Charakteristiky
74
du     d d í      du \
c{x)— + — {r{x,t)u) -— \P{x)fa) =f(x,t),   x e {OJ), te{0,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úzní rovnice, v níž zanedbáme účinky difúze. Jiná forma advekční rovnice tedy je
Hlt H
<x)^r + ^~(r(x,t)u) = f(x,t),   x e (OJ), t e (0,T). (3.8ľ) ot ox
Metoda přímek. Uvažujme rovnoměrné dělení intervalu (0, £), tj. xt = íh, í = 0,1,..., N, kde h = i/N. Předpokládejme, že funkce aix, t) v krajních bodech intervalu (0, í) nemění znaménko. Splnění rovnice (3.81) budeme požadovat ve vnitřních uzlech x\, x2, • • • , Xjv_i a pro a(0, t) < 0 také v uzlu xq a pro a(í, t) > 0 také v uzlu x^. V rovnici
du(xij) du(xij) —3--+a(xiJ)—3-— = f(xi,t)
aproximujeme člen ux(xl,t) 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
ůi(í) + ai(í)Mí) - 9o(t)]/(2h) = fjt),
ůi(t) + ai(t)[ui+1(t)-ui-1(t)]/(2h) = fi(t), i = 2,3,...,N-l, (3.86) ůN(t) + aN(t)[3uN(t) - 4«Ar_i(í) + uN_2(t)]/(2h) = fN(t),
kde (Zj(r) = a(xij), fl(t) = f(xl,t) a u^t) je aproximace u(xl,t). Z (3.83) dostaneme počáteční podmínky
u.
■i(0) = <pi,   i = 1,2,..., N, (3.87)
kde (fi = ipixi). 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
ůi(í) + ai(í)Mí) - uN(t)]/(2h) = /x(í) ,
ul{t) + al{t)[ul+1{t)-ul_1{t)]/{2h) = f^t),   i = 2,3,...,N-1, (3.86') ůN(t) + aN(t)[Ul(t) - uN.1(t)]/(2h) = fN(ť). Počáteční problém (3.86')-(3.87) je pak tvaru
ů = Au,   te{0,T),      u{0) = (p, (3.88)
75
kde u (ŕ)
(u1(t),u2(t),...,uN(t))T, tp{ť)
/   O     1     O   ... O -1     O     1   ... O 0-1     O   ... O
A =---
1
27i
V
0
1
o o
o o
o
-1
(^l(í),^2(í),
-1\
o
0
1
., (pN(t))T a
(3.89)
Dá se ukázat, že vlastní čísla matice A
N 2irí A, = I — sin-,
í      N '
I =
-1,
i = 1,2,..., N,
leží na imaginární ose mezi —N/i a N/l a max, |Aj| —> oo pro N —> 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árni osy, nej větší vlastní číslo AmQX leží prakticky na imaginárni ose, AmQX = IN/í, I = y/—l.
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 Oih). 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
du      du , mN
u(0,t) = g0(t), te(0,T), u{x, 0) = ip{x),   x G (0, í),
(3.89)
kde a > 0 je konstanta. Diferenciální rovnici přepíšeme pomocí charakteristik na tvar du(x(t), t)
dt
= 0.
(3.90)
Zvolme rovnoměrné dělení intervalu (0, £) s krokem h = Í/N, tj. xt = íh, i = 0,1,. .., N, a rovnoměrné dělení intervalu (0, T) s krokem r = T/Q, tj. tn = m, n = 0,1,..., Q. Charakteristika vycházející z bodu [xl,tn+i] je přímka x(t) = xt + ait — tn+i). V čase t = tn je
x,
x(tTl
Xj
ar.
76
Předpokládejme, že ar < h. Pak pro i = 1, 2,..., N bod x™ G (xj_i,Xj), zejména tedy x™ G (0, £), takže u(x™,tn) má smysl. Integrací rovnice (3.90) od tn do tn+i obdržíme
ín + l
+1 dii(xl(ŕ),ŕ) ~~ďí
dŕ = ií(xj, rn+i) — u(x™, tn) = 0 , a odtud -u(xj, rn+i) = u(x™, tn).
Numerickou metodu dostaneme tak, že u(xll,tn) určíme přibližně interpolací, tj. počítáme
<+1 = Pfe«), kde Pfc(x) je interpolační polynom stupně k > 1 splňující
Pk(Xi-l) = K-i,    pk{xi) = < • Pro lineární polynom P\ z podmínek (3.92) odvodíme
(3.91)
(3.92)
Pí (x) = ur:
ar
ln+l
x(t)/	j-
——■	
-j^(x-Xí)   a odtud   Px(x") =<-^-«-<_i)
Dostali jsme tak upwind metodu ar
iíľ+1 = u?
h
« " <_i) •
(3.93)
i i i i i i _i_i_
Obr. 3.10. Upwind metoda
Název metody nám připomíná, že veličina u v uzlu Xj závisí jen hodnotách této veličiny ;   x ^   v uzlech ležících „proti proudu, proti větru", pro xl+i       a > 0 tedy nalevo od Xj. V bodu xo = 0 užijeme okrajovou podmínku, takže Uq+1 = go(tn+i).
Dá se ukázat, že když ar
(3.94)
pak za předpokladu dostatečné hladkosti přesného řešení pro chybu platí
u{xi,tn) -< = O(r), (3.95)
viz [14]. Říkáme, že metoda (3.93) je řádu jedna. Pro v = 1 dokonce = u(xt,tn) je přesné! Číslo v = ar/h se nazývá Courantovo číslo a podmínka (3.94) se nazývá Courantova-Friedrichsova-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(xl+1) =<+1, (3.96)
vypočteme P2(x") a dosadíme do (3.91), dostaneme Laxovu-Wendroffovu metodu
<+1 = <-^«+i-<-i)
2h2
«+1-2< + <_!).
(3.97)
77
Pro aproximaci v uzlu můžeme použít interpolační polynom, který kromě podmínek (3.92) vyžaduje navíc splnění podmínky ^2(^-2) = ^Ar-2-
Jestliže pro obecný uzel Xj přidáme k podmínkám (3.92) podmínku
^-2) =<_2, (3-98)
vypočteme P2(xf) a dosadíme do (3.91), obdržíme Beamovu-Warmingovu metodu
2 2
<+1 = < " 2^(3< - 4<_x + <_2) " 2^"« " 2<-i + <-2) • (3-99)
Počítáme-li <+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{xi,tn) - u? = 0{t2) . (3.100)
Beamova-Warmingova metoda je upwind aproximací řádu 2, řešení v uzlu xl 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),      t E (0, T). (3.89á)
Charakteristika vycházející z bodu [xl5 tn+i] bude směřovat doprava, takže za předpokladu platnosti CFL podmínky
v=^-<1 (3.101)
bude Xi(tn) E {xl1xl+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(Xi) = < ,    Pk{xl+Í) = <+1. (3.92')
Podobně jako dříve dostaneme upwind metodu
<+1 = <-^«-<+i) (3-93')
a Beamovu-Warmingovu metodu
11 22 <+1 = < " ^(3< - 4<+1 + <+2) - ^«+2 - 2<+1 + <), (3.99')
Laxova-Wendrofřova metoda zůstane stejná, tj. postupujeme podle (3.97). V bodu xN užijeme okrajovou podmínku, takže u^1 = g(tn+i).
V obecném případě a = a(x,t) určíme x™ numericky: x™ a; x(tn), kde
^l = a(x(t),t),       x(tn+í)=xl. (3.102)
78
Pro upwind metodu použijeme EE metodu, takže
x™ = xl — a™T ,    kde   a™ = aixl1 tn+i). Pro Laxovu-Wendroffovu metodu a Beamovu-Warmingovu použijeme EM2 metodu, tj.
x™ = xl — a™T,    kde   a™ = \{ki + k2),    k\ = a(xl5 tn+i) , k2 = a(xl — 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+i — tn musí splňovat CFL podmínku
max, \a\ \rn
h
< 1. (3.101')
Jestliže v rovnici (3.89i) uvažujeme nenulovou pravou stranu /(x, ŕ), tj. řešíme-li místo rovnice (3.90) rovnici
^J)=/Wí).t). (3.90-)
přičteme k pravým stranám formulí aproximaci Q?(/) integrálu ft^+1 f(x(t), t) dt. Pro upwind metodu stačí použít jednostrannou obdélníkovou formuli
Q7(f)=rf(xi,tn+1).
Pro Laxovu-Wendroffovu metodu a Beamovu-Warmingovu metodu užijeme lichoběžníkovou formuli
Qnf) = ^[f(xi,tn+1)+f(x^,tn)].
Metoda konečných objemů. Pro diskretizaci rovnice (3.81') se výborně hodí metoda konečných objemů, na časoprostorových objemech Bl = {xl_i/2, xl+i/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 MathWorks, 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 I,II, 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