4. Obyčejné diferenciální rovnice - počáteční úlohy
V této kapitole se budeme zabývat problematikou numerického řešení úloh s počátečními podmínkami pro obyčejné diferenciální rovnice.
4.1. Formulace, základní pojmy
Počáteční problém pro ODR. Obyčejná diferenciální rovnice (ODR) se nazývá rovnice n-tého řádu (ODRn), jestliže nejvyšší derivace neznámé funkce jedné proměnné je n—tého řádu. Obecný tvar rovnice prvního řádu je g(x,y(x),y'(x)) — 0. Budeme předpokládat, že z rovnice je možné explicitně vyjádřit y*{x), tj. že rovnici lze psát ve tvaru
y\x) = f{x,y{x)).
(4.1)
y{a) = rj.
(4.2)
Podmínka (4.2) se nazývá počáteční a, znamená, že integrální křivka y = y(x) prochází bodem [a, n].
Lipschitzova podmínka. Je-li v nějakém okolí D bodu [a, r/] funkce f(x, y) spojitá a splňuje-li v tomto okolí Lipschitzovu podmínku s konstantou L, tj. platí-li
\f{x,u)-f{x,v)\ < L\u-v\   V[x,u], [x,v] e D,
(4.3)
pak bodem [a, 7y] prochází jediné řešení y(x) rovnice (4.1). Jestliže má funkce f(x, y) v okolí D ohraničenou parciální derivaci \dyf(x,y)\ < L, Lipschitzova podmínka (4.3) platí. To je přímý důsledek věty o střední hodnotě:
f(x,u)-f(x,v) = dyf{x,Q(u-v),
kde [x, £] je nějaký vnitřní bod úsečky s krajními body [x, u] a [x, v].
Jiný standardní výsledek říká, že když je funkce f(x, y) spojitá na (a, 6) x R a splňuje tam Lipschitzovu podmínku, pak počáteční problém (4.1),(4.2) má jediné řešení definované v celém intervalu (a, b).
Rovnice (4.1) spolu s počáteční podmínkou (4.2) se nazývá počáteční úloha nebo také Cauchyova úloha.
Soustava obyčejných diferenciálních rovnic 1. řádu je soustava n rovnic pro n neznámých funkcí yi(x),...,yn(x) tvaru
Vj = fÁX> 2/1: ■ ■ ■ : Vn),     j = 1, 2, . . . , TI.
Počáteční podmínky jsou tvaru Vj(a) = Vi,   J = 1,2, ...,n.
Připomeňme, že obecné řešení rovnice (4.1) obsahuje integrační konstantu, která může nabýt libovolné hodnoty. Proto k jednoznačnému určení y(x) nestačí sama rovnice. Potřebujeme znát hodnotu řešení pro jisté x — a, tj.
Počáteční úlohu zapíšeme stručněji vektorově: y' = f(x,y), y (a) = i),
kde
y = (l/i, V2, ■ ■ ■, VnY, y' = (y[, yí, ■ ■ ■,y'JT, f(x,y) = (fi(x,y), f2(x,y),.
(4.4)
Existence a jednoznačnost řešení. Je-li v okolí D bodu [a,Tj] funkce t(x,y) spojitá a splňuje tam Lipschitzovu podmínku s konstantou L, tj. platí-li
\\f(x,u)-f(x,v)\\<L\\u-V\\
(4.5)
pro všechny body [x,u] a [x,\] z okolí D, prochází bodem [a,rj\ jediné řešení počáteční úlohy (4.4). Je-li D = (a, b) xř, jediné řešení existuje v celém intervalu (a, b). Z věty o střední hodnotě plyne
i(x, u) — f(x, v) = 3(x, *)(u — v) ,
kde J(x, *) = {dfi(x, Zi)/dyj}fj=1 je Jacobiova matice funkce í(x, y), jejíž složky v i-tém řádku jsou vyhodnoceny v nějakém vnitřním bodě [x,^] úsečky s koncovými body [x,u] a [x, v]. Jsou-li všechny parciální derivace dijdy v D ohraničené, tj. platí-li
dfi(x,y)
9y3
< L   V[i,y]eDa pro i, j = 1,2,
je Lipschitzova podmínka (4.5) zřejmě splněna.
Rovnice n-tého řádu. Každou rovnici n-tého řádu tvaru
j/<n> = F{x,y,y',...Jn-V)
s počátečními podmínkami
y{a) = ni, y'{o) =
y(n'l\o.)=nn
(4.6)
yn = y^\ dále
lze snadno převést na soustavu (4.4). Stačí položit y\ = y,yi = y'-. f i = Vj+i, j = 1,2,- ■ ■ ,n - 1 a /„ = F(x,yuy2, ■ ..,yn).
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 í(x,y) splňuje Lipschitzovu podmínku a 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 někdy omezíme jen na jednu rovnici. Z velkého množství metod uvedeme jen některé a to takové, které jsou pro své dobré vlastnosti široce používány.
49
4.2. Úvod do numerických metod řešení
Numerickým řešením počáteční úlohy rozumíme výpočet přibližných hodnot hledaného řešení y(.r.) v bodech a-, dosti hustě vykiývajících interval (a,b). Nechť tedy a = x0 < xi < ■ ■ ■ < x*? = b je dostatečně jemné děleni intervalu (a, b). Body Xi nazýváme uzly, vzdálenost hi = xi+í — xi dvou sousedních uzlů nazýváme délkou kroku (stručně krokem), a množinu {x0, x%,..., Xjs/} všech uzlů nazýváme sítí. Jsou-li všechny kroky stejně dlouhé, tj. h{ — h := (b — a)jN, hovoříme o rovnoměrném (nebo také ekvidistantním) dělení intervalu (a, b). V tom případě je Xi = a + ih, i = 0,1,..., N. Hodnotu přesného řešení v uzlu x, budeme značit y(xi) a hodnotu přibližného řešení yt. Jestliže se nám podaří najít přibližné řešení yh i = 0,1,..., AT, můžeme přibližně vypočítat hodnotu řešení y(x) v libovolném bodě x e {a, b) interpolací.
Numerická metoda pro řešení počáteční úlohy (4.1),(4.2) je předpis pro postupný výpočet aproximací yi, y2,..., yN, y0 = ij vždy. Metoda se nazývá k-kroková, závisí-li předpis pro výpočet aproximace yi+\ na předchozích aproximacích j/j, yi-\,... Speciálně jednokroková metoda počítá přibližné řešení v uzlu jen pomocí znalosti přibližného řešení yi v uzlu přibližná řešení yt-i, yi-?, ■ ■ - spočtená v předchozích uzlech Xi-i, Xi-2, ■ ■ • nepoužívá.
Explicitní Eulerova metoda. Nejjednodušší numerickou metodou pro řešení úlohy (4.1),(4.2) je explicitní Eulerova metoda, stručně EE metoda. EE metodu snadno odvodíme z Taylorovy formule
y{xi+i) = y(xi+h) = y{xi)+hy'{xi)+\h2y"{íi),   & € (xhxi+1).
(4.7)
Uvážíme-li, že y'(xi) = f(xi, y{xi)), a zanedbáme-li člen j/iV'(£>), dostaneme formuli
y(xi)+hf(xi,y(xi)).
(4.8)
Tento výsledek je však nepoužitelný, protože neznáme hodnotu y(xi) a navíc se v tomto vztahu vyskytuje znaménko přibližné rovnosti. Proto zaměníme ve formuli (4.8) výrazy y(Xi) a y{xi+i) jejich přibližnými hodnotami y,, yi+i a znaménko přibližné rovnosti as nahradíme znaménkem rovnosti. EE metoda je tedy předpis
Vi+i = yi + hf(xi,yi)
(4.9)
pro i = 0,1,..., N - 1. ?/o = y(xo) = V je přesné.
Diskretizační chyby. Do bodu [xi, y{\, kde y: = y0 + hf(x0, ya), se dostaneme posunem po tečně sestrojené k přesnému řešení y(x) ve startovacím bodě [a;0,i/o]- Posouváme se doprava tak dlouho, až xx-x0 = h, kde h je zvolený krok. Nový bod [xi, y\\ už na integrální křivce y{x) neleží, podle (4.7) je y(x,) - yi = f/i2j/"(£o)- Výraz \h2y"(£0) označíme jako ře0 a nazveme lokální chybou (le jako local error) v uzlu xi.
Dále spočteme y2 = yi + hf(x±, y{). Vidíme, že do bodu [x2, y2] se dostaneme posunem po přímce y — y1 = f(xi,yL)(x-xi). Tato přímka je tečnou k integrální křivce procházející bodem [xi,y-i], tj. k funkci Ui(x), kde
u', = f{x,ui), u1{x1)=yi.
50
i
Funkci ui(x) nazveme lokálním řešením úlohy (4.1),(4.2) příslušným uzlu x\. Do bodu [^2,2/2] se tedy dostaneme po tečně sestrojené k lokálnímu řešení u\(x) v bodě Nový bod [x2,?/2] na integrální křivce u\{x) zase neleží,
ui(x2) = yi + hf(x1,y1) + -h2u"(n-í),
kde r}i € (xi, x2),
takže u-ifa) = y2 + lex, kde lei = \h2u'[(r)i) je lokální chyba v uzlu x2. Celková chyba y{x2) — V2 P° dvou krocích se označuje e2 a nazývá se globální diskretizační chyba (globál truncation error). Globální diskretizační chyba je důsledkem hromadění účinků vzniklých chyb lokálních. Říkáme, že globální diskretizační chyba vzniká kumulací lokálních chyb, Zřejmě e\ = leo avšak leo + le% ~ e2 je pouze aproximace e2.
Xi Xi+l
Obr. 4.1. Diskretizační chyby
Stejně postupujeme dál. Předpokládejme, že jsme už spočetli přibližné řešení yt v uzlu Xi. Sestrojíme lokální řešení Ui(x) definované předpisem
u'i = f(x,ui),   Ui(xi) = yi,
a bod [xi+i, j/í+i] dostaneme posunem po tečně sestrojené k lokálnímu řešení Ui(x) v bodě [xi,yi\. Přitom vznikne lokální chyba lei = Uí(zí+i) - 2/i+i, pro kterou platí
Ui(xi+1) = Ui(xi)+hf(xi,Ui{Xi))+lei, kde leA = |/i2ií"(íJí) a ^ g (xí,xí+1) . (4.10)
Celková chyba ei+1 = y(xi+1)-yt+i vznikne kumulací lokálních chyb lej pro j = 0,1,... ,i.
Možná jste si všimli zdánlivé nesrovnalosti: u lokální chyby Ze; jsme nepoužili slůvko diskretizační. To byl úmysl, spojení lokální diskretizační chyba (local truncation error) je rezervováno pro výraz Itei = |fcV'(&) definovaný rovností (4.7), tj.
y(xi+i) = y(xi)+hf(xi, y(xi))+ltei, kde Ita = |/i2y"(&) a & e (xh xi+l). (4.11)
51
Obě chyby že* a Itei jsou téhož řádu řádu 0(h2), stejné však nejsou. Lokální chyba let je chyba, které se dopustíme v jednom kroku metody. Lokální diskretizační chybu lte{ lze charakterizovat podobně: je to chyba, které se dopustíme v jednom kroku metody, ovšem za tzv. lokalizačního předpokladu, že j/j = y(xi) je přesné. Lokální chyba je tedy chyba, která při výpočtu danou metodou v každém kroku skutečně vzniká. Lokální diskretizační chyba vyjadřuje chybu, které se dopustíme, když do formule dosadíme přesné řešení. V dalším si však ukážeme, že pro malá h je rozdíl mezi oběma chybami prakticky zanedbatelný. Obě lokální chyby le; a lte,, spolu s globální chybou ei+1, jsou zakresleny v obrázku 4.1.
Rád, konvergence. Řekneme, že metoda je řádu p, jestliže pro její lokální diskretizační chybu lte, platí
Hei = 0{hP+v). (4.12) EE metoda je tedy řádu 1. Pro globální diskretizační chybu e; EE metody lze dokázat
e; = y{xi) - Ví = 0{h). (4.13) Rovnici (4.13) je třeba chápat tak, že existují kladné konstanty h0 a K takové, že
max IějI < K h   pro všechna 0 < h < /ín .
Zejména tedy
max \y{x%)-Vi\
0   pro h    0 .
Metody s touto vlastností se nazývají konvergentní metody. Protože globální diskretizační chyba je řádu O(h) říkáme, že konvergence EE metody je řádu 1.
Důkaz konvergence EE metody lze najít například v [34], [35]. Tvrzení (4.13) však lze snadno ověřit v případě, že f(x,y) = f(x) nezávisí na y. Pak totiž
y(xi+l) = y(Xj) + hf(xj) + i/i2/'fó), Vj+i = Ví + hf(xj). Odečtením obou rovnic dostaneme e,- .H i = e j eí = ^2[/'(6) + /'(6) + --- + /'te-i)].
\h2 ľ {Ši), a protože e0 = 0, je
Označíme-li M = max^j, pak |e;| < \h2iM < \[hN]Mh = \{b - a)Mh, neboť
Nh = b- a.
x) = (x2 + c)2, počáteční podmínka
Příklad 4.1. Obecné řešení rovnice y' = Ax^/y je y 2/(1) = 4 určí partikulární řešení y(x) = (x2 + l)2. Řešme tuto počáteční úlohu EE metodou na intervalu (1,3) s různým krokem h. Nejprve zvolme N = 10, tj. h = 0, 2, x0 = 1, X\ = 1,2, x2 = 1,4,..., xw = 3. Z počáteční podmínky máme y0 = 4, dále podle (4.9)
52
počítáme
J/i=2/o + hf{x0, y0) = yo + hixo^yíi = 4 + 0,2-4-l-V/4 = 5,6
?/2 = 2/i + = 2/i + MxiVížľ= 5,6 + 0, 2 ■ 4 ■ 1,2 ■ V/M = 7,8718
2/jv = ž/jv-i + Ä/(a:jv-i,2/a'-i) = 81,826 ,
zatímco přesná hodnota y(3) = 100, tedy absolutní chyba pro x = 3 je 18,2 a relativní chyba je 18,2%.
Pro N = 20 je h = 0,1. Z počáteční podmínky máme opět 2/0 = 4, dále 2/i=2/o + ft/Oro, J/o) = 2/o + Ä4x0v/ž?ô = 4 + 0,1 • 4 • l\/í = 4,8 2/2 = 2/i + fc/(ari,lh) = 2/i + AteiVSh = 4, 8 + 0,1 • 4 • 1,1 • sfcŽ = 5,763992
atd. až konečně yN = 90,40. Absolutní chyba i relativní chyba se zmenšily přibližně na polovinu, což odpovídá tomu, že EE metoda je řádu 1.
Vztah mezi lokální chybou a lokální diskretizační chybou. Jak už jsme uvedli, chyby let a Itei jsou pro malá h prakticky k nerozlíšení. Dokažme si to. Z Taylorova rozvoje dostáváme
y(xi+i) = y{xi) + hy'{xi) + \h2y"{xi)+0(hz), takže lokální diskretizační chybu říe; lze vyjádřit ve tvaru
(4.14)
Člen \h2y"{xi) se nazývá hlavní člen lokální diskretizační chyby lte,. Podobně lze lokální chybu lei vyjádřit ve tvaru
lei = kh2u'l{xi)+0{h?).
Odečtením obou vyjádření dostaneme
lte,-Jcí = lAV(^)-«f(*i)]+0(*S)-Derivováním rovnice (4.1) obdržíme
y"{x) = fx(x,y)+fv{x,y)ý = fx(x, y)+fy(x,y)f(x,y),
takže
y"{xi) - u"{xí) =fx{xi:y(xi)) + fy{xi,y{xi))f(xi,y{xi))--[fx{xhUi(xi)) + fv{Xi,Ui(Xi))f(Xi,Ui{Xi))] ,
a odtud užitím věty o střední hodnotě po snadné úpravě dostaneme
\lf{xl)-v!l(xt)\<M\y(zi)-Mzi)\=M\v{x,)-vt\,
53
kde konstanta M závisí na prvních a druhých parciálních derivacích funkce f(x, y). Avšak y{xi) -Vi = O(h) podle (4.13), takže y"(xi) - uj'fo) = O(h) a
Ite, = lei + 0(h3
(4.15)
Protože jak Itei tak ier jsou řádu 0(h2), lze pro malá h rozdíl mezi oběma chybami prakticky zanedbat.
Zaokrouhlovací chyby. Doposud jsme vůbec neuvažovali zaokrouhlován chyby. Ukážeme si, že jejich vliv nelze zanedbat. Pro jednoduchost předpokládejme, že
2/0 = ??,   2/í+i =Vi+hf(xi,yi)+£i,   i = 0,1,...,N-l, a nechť \s,\ <e. Pro r i := yt - yt proto platí
r,+i = U + h[f(xi, yi) - f(xi,y{)] -e{, takže užitím Lipschitzovy podmínky (4.3) dostáváme
< H + hLIyi-yil+či < (l + hL)\n\+e. Protože r0 = 0, je |ri| < e, |r2| < (1 + hL)s + £,...,
\n\ <e[l + (1 + hL) + (1 + hLf + ■■■ + {! + hL)1-1] =
^ + hVý-1<£é^-i<Kph_,   VApK = u^)L
hL hL ~
Pro celkovou chybu EE metody tedy platí
max \y(xi)-yi\ < Koh + sK^hT1,
0<i<N
<K eh-\   kde K
1]/L.
(4.16)
kde h je dostatečně malý krok a KD,KR jsou konstanty na /i nezávislé. Protože konstanta e je malá, vliv zaokrouhlování se projeví až pro značné velký počet kroků N (tj. pro velmi malé h), v tom případě může dojít dokonce k převážení účinků zaokrouhlovacích chyb a vypočtené hodnoty y i tak mohou být i přes ohromné „výpočetní náklady" zcela bezcenné.
Stabilita. Abychom minimalizovali náklady, snažíme se při reálném výpočtu měnit délku kroku tak, aby přibližné řešení bylo dostatečně přesné a současně aby celkový objem výpočtů byl co nejmenší. Řízení délky kroku je podrobně popsáno v dalších odstavcích. V zásadě však jde o to, jak volit kroky co nejdelší. Vzniká proto přirozená otázka: jak velký krok je ještě možné volit, aby přibližné řešení bylo přijatelné? Z řady dobrých důvodů je účelné pro odpověď na tuto otázku použít testovací úlohu
v' = xy, y(o)
(417)
kde A < 0 je daná konstanta. Co je na této úloze tak zajímavého? Přesné řešení y(x) = eXx je klesající funkce, eXx —> 0 pro x —» +00. Počáteční hodnotu y(0) = 1 lze považovat za jednotkovou poruchu a nezávisle proměnnou x za čas. Úlohu (4.17) lze proto interpretovat
54
jako proces stabilizace: s rostoucím časem dochází k útlumu počáteční poruchy neustálým přibližováním se k nulovému stacionárnímu stavu. Takové chování vykazuje celá řada dějů reálného světa, úloha (4.17) je proto rozumným modelem. Co bychom tedy měli očekávat od numerické metody aplikované na řešení této modelové úlohy? Přijatelné numerické řešení se musí blížit ke stacionárnímu stavu, tedy j/j —> 0 pro i —¥ 00. Nemusí být nutně yi > 0, jde nám totiž přednostně o zmenšování velikosti poruchy, požadujeme však
\Vi+i\ < \Vi\,
což zaručuje, aby utlumování bylo nepřetržité. EE metoda dá pro testovací úlohu (4.18)
(4.18)
Vi+i =Vi+hXyi = (l + h\)yi = ■•• = (l+hX)i+1yQ . Aby platilo (4.18), musí být \l + h\\ < 1.
Předpokládejme nyní, že v úloze (4.17) je A komplexní číslo a označme h := hX. Množina všech h, pro která je při řešení testovací úlohy (4.17) zvolenou numerickou metodou splněna podmínka (4.18), se nazývá oblast absolutní stability numerické metody. Odvodili jsme, že pro EE metodu je oblast absolutní stability jednotkový kruh \z + 1| < 1 komplexní roviny se středem v bodě [-1,0]. Průnik oblasti absolutní stability s reálnou osou se nazývá interval absolutní stability. Pro EE metodu je interval absolutní stability interval (—2,0). Pro reálné A < 0 je tedy třeba vybrat h < 2/|A|. Řekneme, že metoda je absolutně stabilní pro zvolené komplexní číslo A a krok h, jestliže h = hX leží v oblasti absolutní stability této metody. EE metoda je tedy absolutně stabilní pro A a h taková, pro která \h\ + 1| < 1.
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ě. Zabývali jsme se jí především proto, že je velmi jednoduchá a proto na ní lze snadno ukázat, s jakými problémy se při numerickém řešení počátečních úloh pro ODŘI setkáváme a jaké pojmy se při popisu numerických metod běžně používají.
Implicitní Eulerova metoda. Vysvětlení jednoho pojmu jsme ale zatím přeskočili: proč mluvíme o explicitní metodě, co tím máme na mysli? To nejlépe pochopíme, když zavedeme implicitní Eulerovu metodu (stručně IE metodu) jako předpis
Vi+i = yi + hf{xi+L,yi+1).
(4.19)
Metoda je implicitní proto, že neznámá yí+i je argumentem funkce f{x, y), takže z rovnice (4.19) nelze odvodit explicitní vzoreček j/j+i = „něco známého". je proto třeba určit jako řešení z rovnice
yi-hf{xi+1,z) =0.
(4.20)
55
Je-li f(x, y) lineární funkcí proměnné y, tj. f(x, y) = a(x)y + b(x), je řešení takové rovnice snadné,
2/i+i
y{ + hb(xi+i) 1 - ha(xi+i)
Obecně je však třeba řešit rovnici (4.20) přibližně pomocí vhodné iterační metody. 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 Zíěj IE metody je definována předpisem
yfa+i) = yfa)+hffa+í,yfa+l))+ltei,
kde y{x) je přesné řešení rovnice (4.1). Z Taylorova rozvoje y(xi) = yfa+i) — hy'(xi+i) + \h2y"fa+x) + 0(h3), a protože y"(xi+l) = y"fa) + 0(h), dostaneme
ltei = -\h2y"fa) + 0[h3).
(4.21)
IE metoda je tedy řádu 1. Srovnáním (4.21) a (4.14) vidíme, že hlavní člen lokální diskretizační chyby IE metody je co do velikosti stejný jako hlavní člen lokální diskretizační chyby EE metody, má však opačné znaménko. Obě metody jsou tedy stejně přesné. Podívejme se dále na stabilitu IE metody. Pro testovací úlohu (4.17) je
í/i+i:
yi+h\yi+i,   odtud yj+1
hX"
1
1 - h\
yo
a tedy podmínka (4.18) platí pro |1 - h\\ > 1. Oblast absolutní stability IE metody je tedy obrovská, je to celý vnějšek \z — 1| > 1 jednotkové kružnice komplexní roviny se středem v bodě [1,0]. Pro Re(A) < 0 lze proto volit krok h > 0 libovolně velký a podmínka stability (4.18) bude pořád splněna. Metoda, jejíž oblast absolutní stability obsahuje celou zápornou polorovinu Re(,z) < 0, se nazývá A-stabilní metoda. IE metoda je tedy A-stabilní.
Jak řešit nelineární rovnici? 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 (4.20). To neumíme jinak než užitím nějaké iterační metody. Počáteční aproximace se určí snadno, například pomocí EE metody položíme
z0 = yi + hf(xi,yi).
Jak ale počítat další aproximace? Zkusme metodu prosté iterace
Zs+i = <pfa),   kde ip(z) =y{ + hffa+1,z). Konvergence zs -> yi+l nastane, když
\tp'(z)\ — h\fv(xi+i, z)\ < a < 1,   kde a je nějaká konstanta .
Tato podmínka představuje omezení délky kroku h. Konkrétně pro testovací rovnici (4.17) dostaneme \hX\ < 1, což je podmínka ještě přísnější, než jakou vyžaduje podmínka absolutní stability EE metody. Metoda prosté iterace je tedy naprosto nevhodná.
Protože zo je dosti dobrá aproximace dá se očekávat rychlá konvergence Newtonovy metody. Praktická zkušenost potvrzuje, že tomu tak skutečně je. Počítáme proto
Z3+l
g(z3) = z3-yj- hf(xi+hzs) g'{zs)      s      1 - hfvfa+uzs)
s = 0,1,..., Si,
a nakonec položíme y,+i = Zs(+i. Počet iterací Si se volí, zhruba řešeno, tak, aby odhad chyby \z — Zs;+ij Newtonovy iterace byl menší než qe, kde e je požadovaná velikost lokální chyby le; IE metody a q < 1, např. q = 0,5. Přibližné splnění podmínky |le;| < e zajistíme výběrem vhodné délky h kroku.
Jestliže derivaci fy nepočítáme v každé iteraci, dostáváme modifikaci Newtonovy metody, kterou budeme nazývat zjednodušenou Newtonovou metodou. Úspory, které zjednodušená Newtonova metoda přináší, vyniknou zejména při řešení soustav ODR. Více si o tom dozvíme v odstavci 4.5.2.
4.3. Jednokrokové metody
se vyznačují tím, že přibližné řešení y;+i v uzlu = Xi + h se počítá ze vztahu, v němž kromě neznámé yi+1 vystupuje dále už jen předchozí uzel x,, přibližné řešení y, v uzlu %i, krok h a samozřejmě také pravá strana f(s,y) diferenciální rovnice. Jednokrokovou metodu lze tedy formálně zapsat ve tvaru
yi+i = y> + h$(xi, y<, yi+1, h; f).
(4.22)
kde $ je funkce čtyř proměnných x,, y;, ym a h, závislá na funkci f(x, y). Pokud funkce $ na yi+i nezávisí, jde o metodu explicitní, pokud $ na yi+i závisí, jde o metodu implicitní. Pro explicitní Eulerovu metodu (4.9) je * = f(£i,yi) a pro implicitní Eulerovu metodu (4.19) je * = f^j + h,yi+l).
Lokální diskretizační chyba obecné jednokrokové metody (4.22) je definována rovnicí
y (xi+-i) = y (i,-)+/í* (xí , y fa), y (xť+1), h; f)+lte,-,
(4.23)
kde y(x) je řešení úlohy (4.4). Ze srovnání (4.22) a (4.23) plyne, že lokální diskretizační chyba explicitní jednokrokové metody je chyba, které se dopustíme v jednom kroku za předpokladu, že y, = y fa) je přesné. Vskutku, označíme-li
9i+i = yfa)+h$fa,yfa),yfa+i),h;f),  pak ltei = y(x»+i)-ýi+i■
Předpoklad yť = yfa) bývá označován jako lokalizační předpoklad. Největší celé číslo p, pro něž platí
ltež = 0(/r+1),
(4.24)
se nazývá řád metody. Rovností (4.24) přitom míníme, že každá složka vektoru ltej je řádu 0(hp+1). Víme již, že EE metoda a IE metoda je řádu 1. Dokázali jsme to sice zatím
56
57
jen pro případ jedné rovnice, tvrzení však platí také pro soustavy. Tak například pro IE metodu dostaneme
lte :=y(x + h) — y(x) - hf(x + h, y(x + h}) = hy'{x + h) + 0(h2) - hy'(x + h)= 0(h2).
4.3.1. Metody Taylorova rozvoje
Explicitní Eulerovu metodu jsme odvodili z Taylorova rozvoje
y(x+h) = y(x)+hý(x)+±h2y"(x)+lh3y"'(x)+- ■ ■+^y^(x)+0(h"+1) (4-25)
tak, že jsme položili p = 1 a zanedbali jsme chybu 0(h2). Protože y'(x) = {(x,y(x)), dostali jsme EE metodu yi+1 = yt + hí{Xi,yi). Podobně můžeme pro p = 2 dostat metodu řádu 2. Stačí jen vyjádřit druhou derivaci y"{x). Platí
y"(x) =±f(x,y(x)) = ttxjtxJl+MxjWly'W = ío;{x,y(x)) +íy{x,y{x)){{x,y(x)), kde fy(x,y) = {dfj{x,y)/dyt}^e=1 je Jacobiova matice funkce f(:r,y). Metoda
yi+i = yi+hí{xi,yi) + ^h2[í:r(xi,yi)+íy(xi,yi)í{xhyi)]
(4.26)
je proto řádu 2. Za zvýšenou přesnost však platíme značnou daň: v každém kroku metody (4.26) musíme oproti EE metodě počítat navíc n parciálních derivací dfj(xi,yi)/dx a n2 parciálních derivací dfj(x,yi)/dyi. Metoda (4.26) se nazývá metoda Taylorova rozvoje řádu 2. Podobně lze zavést metodu Taylorova rozvoje řádu p předpisem
y>+i
: yi+hí(xi; yj+m^ixi,yi)+- • .+1 ftPffr-Dfo, y;;
(4.27)
kde í^(xi,yi) dostaneme tak, že do výrazu &f{x,y(a;))/dxJ dosadíme xt místo x a y; místo y{x). Výrazy pro vyšší derivace funkce f(x,y(x)) jsou obecně velmi složité. Pro některé zajímavé speciální případy však tyto výrazy mohou být docela jednoduché, takže pak lze metodu Taylorova rozvoje vysokého řádu použít velmi dobře. Jsou-li metody Taylorova rozvoje kvalitně implementovány, mohou být velmi efektivní i pro obecné problémy.
4.3.2. Rungovy-Kuttovy metody
V dalším půjdeme jinou cestou a odvodíme jednokrokové metody vyšších řádů, v nichž vystupují jen hodnoty funkce f. Takové metody se nazývají Rungovy-Kuttovy metody. Obecný tvar s-stupňové Rungovy-Kuttovy metody je
Yi+i = yi + h^2 bjkj
(4.28)
58
kde koeficienty k, jsou určeny předpisem
8
kj = i(Xj +hcj,y<+flj/kf),   j = 1,2,.,.jí
(4.29)
Abychom dostali konkrétní metodu, musíme určit stupeň s a dále reálné konstanty bj, Cj a a,jt pro j,£ = 1,2,..., s. Tak například pro s = 1, b\ = 1, c\ = 0, on = 0 dostaneme EE metodu a pro s = 1, &i = 1, C\ = 1, an = 1 dostaneme IE metodu. Označíme-li
Xj = Xi+hcj, Y j = yj+hy^ajtki,   j = 1,2,
vidíme, že kj = f(Xj, Yj), takže
y>+i =yi + hJ2W{Xj,Yj
(4.30)
(4.31)
3=1
určujeme pomocí s hodnot pravé strany i(x,y) ve vhodně vybraných bodech [X,-,Yj]. Rungovy-Kuttovy metody včleníme do obecného schématu (4.22) jednokrokových metod tak, že položíme <ř = $ľj=i bjí(Xj,Yj).
Obecná Rungova-Kuttova metoda s koeficienty (4.29) je implicitní: abychom určili vektory kj, musíme řešit soustavu ns obecně nelineárních rovnic, neznámé jsou složky všech vektorů ki, k2,..., ks. Dá se ukázat, že pro malé h má soustava (4.29) vždy jediné řešení. Situace se zjednoduší v případě semi-implicitních Rungových-Kuttových metod, kdy a jí = 0 pro l > j, takže
kj = f(Xi+hcj,yi + h^ajŕkŕ),   j = 1,2, ...,s.
(4.32)
V tomto případě lze nejdříve vyřešit soustavu n obecně nelineárních rovnic pro neznámé složky vektoru ki, pak další soustavu pro složky vektoru k2 a tak dále, až nakonec vyřešíme soustavu pro neznámé složky vektoru ks. Tedy místo abychom řešili jednu soustavu o ns neznámých, řešíme postupně S soustav pro n neznámých. Nejznámější Rungovy-Kuttovy metody jsou však metody explicitní, pro které cijt = 0 pro i > j, takže
ki = i(xi,yi),   kj = f(xi+hcj,yi+h^raj{ke),   j = 2,3,.
(4.33)
Pak se totiž žádné soustavy rovnic řešit nemusí, stačí jen postupně počítat ki,k2,... ,ks. Poznamenejme, že pod pojmem Rungova-Kuttova metoda se v literatuře běžně rozumí právě metoda explicitní, zatímco skutečnost, že metoda je semi-implicitní nebo plně implicitní se zdůrazňuje označením semi-implicitní Rungova-Kuttova metoda nebo implicitní Rungova-Kuttova metoda.
Koeficienty Rungových-Kuttových metod je zvykem zapisovat do tabulky známé jako Butcherova tabttlka:
59
Cl		012 •	■ als
C2	021	022 ■	■ a2s
Cs	ftsl	as2 .	
	f>l	h ■	. bs
U explicitních metod se do tabulky zapisují koeficienty jen pro l < k. Protože všechny prakticky používané metody splňují podmínku
cj - E ail ■■ l=\
3 = 1,2,...,s,
(4.34)
budeme i my předpokládat, že podmínka (4.34) platí.
4.3.2.1. Podmínky řádu
Věnujme se nyní otázce, jak určit koeficienty bj, Cj a tak, aby metoda byla zvoleného řádu. Naše úvahy se zjednoduší, když problém (4.4) přepíšeme do autonomního tvaru
y(t)\'_ (í(x(t),y(t)) x(t)    - { 1
y(o)
x(0)
(4.35)
Aplikujeme-li na tento systém Rungovu-Kuttovu metodu, dostaneme
(
yi+i
xi+l
3 = 1
kde body [Yj,Xj] splňují
1
3 = 1,2,
Vzhledem k podmínce (4.34) je Xj =Xi + hcó, takže body [X,, Y,] splňují (4.30). To ale znamená, že odvození Rungovy-Kuttovy metody pro autonomní systém (4.35) je ekvivalentní odvození Rungovy-Kuttovy metody pro obecný problém (4.4).
Předpokládejme tedy, že pravá strana f(i,y) je funkcí jen proměnné y, řešíme tedy rovnici y' = f (y). Rungova-Kuttova metoda je řádu p, jestliže pro její lokální diskterizační chybu platí
lte ;=y{x + h) - y(x) - h^bjk, = 0{hp+1).
3 = 1
(4.36)
kde
k, =f (y(x) +       ajikt),   j = 1,2,..., a.
Podmínky nutné pro dosažení řádu p dostaneme tak, že ve výrazu lte rozvineme funkci y(x + h) okolo x, kj okolo y(x) a anulujeme koeficienty u mocnin fý pro j < p.
60
Odvození podmínek řádu 3. Začneme tím, že v rozvoji (4.25) funkce y{x + h) vyjádříme derivace řešení y(x) pomocí pravé strany f(y). Derivováním rovnice y'(a;) = f{y(x)) po složkách dostaneme
m = ±fr(y(x)) = £ ~(y(x))y'a(x) = E |r(y(*))/«»(y(*)) ■
QX a=l °y" a=l °Va
Abychom zjednodušili zápis, r-tou složku vektoru označíme indexem r vpravo nahoře a parciální derivaci podle složky a označíme čárkou a indexem a vpravo dole. Tedy fr je r-tá složka vektoru f, ya je složka a vektoru y, f*a je derivace r-té složky vektoru f podle složky ya vektoru y. Vyšší parciální derivace zapíšeme podobně, například druhou parciální derivaci r-té složky vektoru f podle proměnných ya a ys označíme frag. Užitím tohoto označení dostaneme
n
(y"r = E/y-
Q = l
Elementární diferenciály. V dalším se nám bude hodit, když definujeme vektor Dj = f a vektor ~D\, jehož r-tá, složka (D^)r = YZ=i Pak totiž V'=D} a Ý' = Dľ Další
derivace je
(y"T = E
f-' dyf> dx 1     J>a^ dyP dx
.3=1    * 3=1 *
■ E [ľ^rľ+W]
a,a--i
Položíme-li (DIY = E^=i fofr a (D|)r = £^=1 &f$ľ, P»k Ý" = Df + D|. Vektory D|, T)\, T5\, Ti\ jsou tzv. elementární diferenciály.
V tomto postupu můžeme pokračovat a vyjádřit obecně y^' pomocí elementárních diferenciálů T)\, £ = 1,2,..., \j. S růstem j počet Aj elementárních diferenciálů D3e rychle roste. Abychom výklad příliš nekomplikovali, spokojíme se s odvozením postačujících podmínek pro metody řádu p < 3 a k tomu nám již odvozená vyjádření pro prvou, druhou a třetí derivaci řešení y(x) stačí.
Pokračujme proto rozvojem složek
kr](y) = My(x) + hJ2^í) koeficientu kj okolo bodu y(x). Z Taylorova rozvoje dostaneme
(4.37)
Q=l l=\
8,7=1 1=1
V tomto výrazu se opět vyskytují složky koeficientů k^. Za ně znovu dosadíme z Taylorova rozvoje. Protože člen kf je ve výrazu (4.37) násoben h, stačí použít zkrácený rozvoj
n s
k? = /Q + '*E /,í E a^kSm + '
{=1 7/1=1
61
Dosadíme-li toto vyjádření do (4.37) zjistíme, že se v něm vyskytují členy ksm, kf a kj násobené h2. Proto je nahradíme ještě kratším rozvojem
fc£=/'+0(A),   k} = r+0{h),   ksm = f+0(h).
Jestliže provedeme všechna naznačená dosazení, dostaneme z (4.37)
*j =r+/> E &É (/"+h E É a^\ŕ+°^+°^ j+
a=l       1=1        l 6=1      m=l J
^ E   É fl*^+°w] É a^f++°^3) ■
0,7=1        £=1 £=1
Pomocí (4.34) poslední výraz upravíme na tvar
a=l a,S=\ 1=1
a odtud pomocí elementárních diferenciálů
B,i=i
kj = D}+/icíD2+ä!!
^D? + E°^D;
+ 0(ft3).
Protože
y(x+/i)-y(x) = ftDJ + I^DHI^ÍDj+D^+O^4) dosazením do (4.36) dostaneme
lte =y(x + h) - y (x) - h E &jkj = feDi
.7=1
i-E6^
ft2D2
I-E6^
(4.38)
I-IE^
3=1
i - E haHct
3,1=1
+ 0(fe4),
nebo stručněji lte = M?Dj + ha2?D? + /ŕpfD? + T|D|] + 0(ň.4), kde T/ jsou tzv. koeficienty lokálni diskretizační chyby.
Podmínky pro metody řádu 1 až 3. Aby Rungova-Kuttova metoda byla řádu p < 3, musí být tzv. podmínky řádu T/ = 0 pro j = 1,2,... ,p a i = 1.2,..., Aj, kde
Ai = l,	2?	= 1	" Ei=l 6i .
A2 = l,		_ i 2	- Eí=i 6tCi,
A3 = 2,	T?	_ 1 2	
02
Obecně, má-li funkce f (y) spojité derivace všech řádů, můžeme lokální diskretizační chybu Rungovy-Kuttovy metody vyjádřit ve tvaru
lte = E^Et/d2-
j=i i=i
(4.39)
Podmínky pro metodu řádu 4. Podmínky řádu Tj = 0, £ — 1, 2,..., Áj, j = 1, 2,... ,p,
nutné a postačující pro metody řádu p, lze najít ve většině učebnic numerické matematiky nejméně pro p < 4. Uveďme si proto koeficienty Tf lokální diskretizační chyby pro p = 4 také my:
r24
^3 — 2 [12    Eij=l ftfOy^] i TA
' Eí,j=i biCiOijCj, - E°j,*=i biaaajkCk ■
V [13] je popsán algoritmus pro generování koeficientů T/ a tyto koeficienty jsou tam explicitně uvedeny pro p < 7.
Čtenáře možná napadlo, proč jsme podmínky řádu odvozovali s poměrně velkým úsilím pro soustavu y' = f (x, y), když pro jednu rovnici y' = f(x, y) by to jistě bylo jednodušší, Důvod je ten, že podmínky řádu pro soustavu se od podmínek řádu pro jednu rovnici mohou lišit. Lze ukázat, že pro p < 4 tato situace nenastane, takže všechny Rungovy-Kuttovy metody, které jsou řádu nejvýše 4 pro jednu rovnici, jsou téhož řádu i pro soustavu rovnic. Pro p > 5 se však podmínky řádu pro soustavu a pro jednu rovnici už liší a existují metody, které jsou řádu p pro jednu rovnici a řádu menšího než p pro soustavu. Když byla tato skutečnost odhalena, vyvolala rozčarování a obavy, neboť v té době bylo známo několik metod vyšších řádů, odvozených pro jednu rovnici, ale používaných i pro systémy. Naštěstí se ukázalo, že řád těchto metod na počtu rovnic nezávisel.
Konvergence. Pro rychlost konvergence Rungovy-Kuttovy metody platí následující
Věta (o rychlosti konvergence Rungovy-Kuttovy metody). Rungova-Kuttova metoda řádu p > 1 má globálni diskretizační chybu er = y(xi) — y, řádu 0(hp) .
Důkaz věty je uveden např. v [38], [17]. Předpokladem platnosti věty je dostatečná hladkost pravé strany f, konkrétně je třeba, aby funkce i(x,y) měla spojité derivace až do řádu p včetně. Pokud f má spojité derivace jen do řádu s < p, pak lze pro globální chybu dokázat pouze řád 0(hs).
4.3.2.2. Explicitní metody
Označme p(s) maximální dosažitelný řád s-stupňové explicitní Rungovy-Kuttovy metody abychom zdůraznili, že závisí na s. O funkci p(s) je známo, že p(s) —¥ co pro s —> oo. Kromě toho platí
p(s) = s   pro s = 1, 2,3, 4. P(5) = 4, p(6) = 5, p(7) = 6,
p(8) = 6, p(9) = 7 p(s) < s - 2
pro   s = 10,11,.
63
Explicitní p-stupňové Rungovy-Kuttovy metody řádu p tedy existují jen pro 1 < p < 4, tj. pouze pro tyto řády existují metody, které počítají pravou stranu jen p-krát. Uveďme si několik nejznámějších metod tohoto typu.
Metoda řádu 1. Pro s = p = 1 existuje jediná explicitní metoda a tou je nám již známá explicitní Eulerova metoda yi+i = yť + hí(xi,yi).
Metody řádu 2. Pro s = p = 2 má explicitní metoda Butcherovu tabulku Podmínky pro metodu řádu 2 stanoví
61 + 62 = 1,   62c2 = a protože ve shodě s (4.34) předpokládáme 021 = 02, dostáváme tabulku
kde ob = |. Parametry a, b jsou tedy svázány jednou podmínkou, takže zvolíme-li a / 0, je b = l/(2a).
0	
Cl	0-21
	h b2
1
2 '
0	
a	a
	1-b b
Yi+i = Yi+hki,   kde   k2 = ffo+jft, y;+|/iki), k! = f(xh yt
(4.40)
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 (4.40) známa jako „midpoint Euler formula".
Pro a = ljer> = |a dostáváme metodu
yt+i = yi+5Ä[ki+k2],   kde  kx = f(xi,yj), k2 = f(xi+/i,yi+/ik1)
(4.41)
uváděnou též pod názvem modifikovaná Eulerova metoda. Budeme ji značit EM2 jako druhou modifikaci Eulerovy metody. Metoda (4.41) se také často uvádí pod názvem Heunova metoda.
Dosazením do (4.38) dostaneme pro metodu řádu 2 lokální diskretizační chybu
lte = Ä3[7?D?+2fDä]+0(fc*) = h3 [(\ - \ba2) D? + iD»]+0(A4).
Odtud předěsím vidíme, že metodu řádu 3 nemůžeme dostat pro žádnou volbu a, b, neboť člen |D| není čím anulovat. Pomocí vhodné volby koeficientu o můžeme ovlivnit g - 56a2 = i - \a. Pro EM1 metodu je Tx3 = i a pro EM2 metodu ^. Pokud bychom tedy mohli vliv druhého sčítance ^D3, zanedbat, mohli bychom tvrdit, že EM1 metoda je dvakrát přesnější než EM2 metoda. Pro a = § lze dokonce docílit anulování koeficientu T3. Pak b = l/(2a)
jen koeficient T3
|, takže dostáváme metodu
yť+i = yi+j/i[k!+3k2],   kdekj = f(zj,yť), k2 = ffo+f/i, yi+f/iki),
(4.42)
která je optimální v tom smyslu, že pro ni jsou absolutní hodnoty koeficientů Tf a Tf nejmenší možné. Taková metoda se nazývá metoda s vyladěnou chybou. I když metoda (4.42) příliš rozšířená není, nejlepší používané metody vyšších řádů mívají chybu vyladěnu, tj. jde-li o metodu řádu p, pak jsou absolutní hodnoty \T?+l\ koeficientů T?+1 hlavního
64
členu 5Zj=T T?+1DJ+1 lokální diskretizační chyby lte malé. Je však třeba uvést, že kvalitu metody ovlivňuje ještě řada dalších faktorů, které jsme doposud nezkoumali a o nichž se zmíníme později.
Metody řádu 3. Pro s = p-
0		
r-1	c2	
Ci	c3 - 132	1.12
	h	í>2 h
■ 3 dostáváme Butcherovu tabulku
a 4 podmínky pro metodu řádu 3 bi + b2 + h = 1, b2c2 -i- &3c3 = \ b2c\ + b3cl = i ,   63a32C2 = i .
Ukazuje se, že když zvolíme dva parametry 0 < c2 < c3, jsou tím všechny koeficienty
metody jednoznačně určeny. Volba c2 0
, c3
vede na populární Ralstonovu metodu:
0
y«+i = y< + \h [2ki + 3k2 + 4k3] . k,_ = f(xť,yť),   k2 = f(xť + \h,yi + |Äki), k3 =f(xi + $h,yi + lhk2).
Ralstonova metoda je řádu 3 a má vyladěnou chybu: dva koeficienty T24 a T3 hlavního členu lokální diskretizační chyby jsou nulové, koeficient T44 = ^ pro jakoukoliv volbu parametrů c\, c2 a koeficient T4 = 2gg je malý.
Metody řádu 4. Pro s = p = 4 je nejznámější klasická Rungova-Kuttova metoda
Yi+l = y« + \h [ki + 2k2 + 2k3 + k4] , ki = t{xhyi), k2 = t(xi + íh,yi + |/iki)
k3 = f{Xi + \h, yi + \hk2), k4 = f {x i + h, yť + ftk3).
0			
1	1		
2	2		
1	0	1	
2		2	
1	0	0	1
	1	1	1 1
	6	3	3 6
Klasická Rungova-Kuttova metoda 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. Jednu velice kvalitní dvojici metod řádu 4 a 5 uvedeme v odstavci 4.3.2.4.
4.3.2.3. Řízení délky kroku
Doposud jsme předpokládali, že délka h kroku metody je konstantní. Avšak ani z odvození metod ani z jejich tvaru nevyplývá nutnost takové volby. V každém kroku můžeme délku h měnit, otázkou zůstává jak. Ideální přání uživatele je donutit program, aby mu spočítal řešení s požadovanou přesností, tj. aby pro zadanou toleranci e > 0 platilo ||ej|| = ||y; - y(zi)|| < £■ To se však neumí, současné programy pro řešení ODR dovedou zajistit jen to, aby pro dostatečně malé e byla globální chyba ||e,|| dosti malá. Docilují toho tím, že délku kroku h vybírají tak, aby velikost ||lei|| lokální chyby nabývala pořád přibližně stejné hodnoty e popřípadě he. Velikost ||ei|| globální chyby lze jen odhadnout (řídit ji zatím neumíme), v odstavci 4.3.2.6 si ukážeme, jak se to dělá. Předesíláme však, že monitorování velikosti globální chyby je nákladné a že většina programů tuto službu
65
nenabízí.
Strategie EPS a EPUS. Připomeňme si, že lokální chybou tej rozumíme chybu, které se dopustíme v jednom kroku metody yí+] = y» + hi&fa, yť, yi+i, ht; i), tj.
leť = u(xi+hi)-yi+i,   kde   u' = f(z,u), u(ar;) = yť
je lokální řešení. Velikost lokální chyby se obvykle ovládá jedním ze dvou následujícících způsobů. Podle kritéria chyby v kroku, zkráceně EPS podle anglického error per step, se krok z Xi do xi+i = Xi + K akceptuje, jen když
IM <e. (4-43)
Podle kritéria chyby v jednotkovém kroku, zkráceně EPUS podle anglického error per unit step, se krok akceptuje, jen když
||le,-1| <hi£. (4-44)
Není-li podmínka (4.43) resp. (4.44) splněna, krok je odmítnut a je třeba ho opakovat pro menší délku h* kroku takovou, aby tato podmínka chyby splněna byla. Délku h\ kroku však nevolíme ani zbytečně malou, naopak, vybereme ji co možná nej větší.
Když program provede krok z i; do spočte odhad est; vzniklé lokální chyby le». Jak se to provede, to si ukážeme v odstavci 4.3.2.4. Nyní je pro nás podstatné, že lokální chyba Rungovy-Kuttovy metody řádu p je
íe^v-^y.^+otO, (4-45)
kde funkce ip{x,y) je tzv. hlavní funkce chyby (z (4.39) plyne ^ = Y,e£i 27+1Df+1), takže
esti^let = h^(Xi,yi)+0{h^2). Pokud bychom z bodu zť postupovali krokem délky h*, dostali bychom lokální chybu
le« = u(xi+h*)-y*+1 = (ft*rV(^,yi)+0((^r2) ■ Srovnáním obou výrazů dospíváme k závěru, že lokální chyba
le* « (h'/hiY+'esti (4.46)
je aproximována výrazem, který- umíme spočítat. Jestliže je norma odhadu estj lokální chyby le, větší než tolerance e, tedy platí-li ||estj|| > e, použijeme přibližnou rovnost (4.46) pro určení optimální délky kroku, tj. maximální délky, pro kterou je zvolené kritérium chyby ještě splněno. Jestliže k řízení délky kroku použijeme kritérium EPS, vybereme h\ Lak, aby
e=||laJ||»(fti7Äír1ll«rti Zvolíme tedy
h\ = /í,(£/||estd|)1/(p+I),
(4.47)
6G
[
a výpočet yi+1 opakujeme s krokem délky h*. Nová délka h* kroku je zřejmě menší než původní délka Ať, jde o to, jestli to stačí: při odvozování h* jsme se totiž dopustili řady nepřesností, sice malých, řádu 0(h?+2), ale dopustili. Po výpočtu yi+1 s krokem délky h\ se proto může stát, že test chyby znovu neprojde, i když třeba jen těsně. Opakované krácení kroku je „drahé" a proto je na místě být opatrný a za h* zvolit jen jistou část 0 optimální délky. V [34] se uvádí, že reprezentativní hodnota používaná v řadě programů je 0 = 0,9. Proto zvolíme
h\ = ©^(e/Hestdl)1^1'.
Zatím jsme se zabývali jen tím, co máme dělat, když je krok zamítnut. Jestliže však je krok z xt do úspěšný, je třeba rozhodnout, jak vybrat délku h,i+1 dalšího kroku. Pro lokální chybu platí
lej+1 = u(a;í+i + /zi+i)-yi+2 = h^^{xi+1,yi+1) + 0{h^).
Tuto chybu můžeme dobře aproximovat pomocí odhadu estj lokální chyby ve stávajícím kroku. Vskutku, protože xi+r = xt + hj a yi+1 = yť -f O(foj), užitím Taylorova rozvoje dostaneme ip(xi+1, yi+1) = ij)(xi, y;) + 0(/i;), takže výraz pro lokální chybu
ten.! = AJíX*i,yí)+0(AEÍ?)
je formálně stejný jako výraz pro lokální chybu lej. Proto volíme délku
hi+1 = eMe/l|estd|)1/(p+,)
kroku z xi+i do xi+2 stejně jako délku h* zamítnutého kroku.
Jestliže k řízení délky kroku použijeme kritérium EPUS, nahradíme (4.47) pomocí
fe*e=||le*||«(^Ar+1||estd|. Odtud vyjádříme /i* a dostaneme předpis
Ä* = 0Äŕ(£VI|estdl)1/p.
Tentýž předpis použijeme také pro určení délky kroku hi+i. Následuje hrubý Algoritmus řízení délky kroku
Krok 1. Spočteme yí+I = y, + h${Xi, yit yi+1, h\ f).
0/i(£/||estd|)1/{,,+1)  pro EPS, XEPS, e^e/i/Hesídl)1/''     pro EPUS, XEPUS. Je-li ||estj|| > £, následuje krok 3, v opačném případě pokračujeme krokem 4.
Krok 3. Položíme h := M^(h,h*) a pokračujeme krokem 1.
Krok 4. Položíme h := M4(/i, h*), zvýšíme i := i + 1 a pokračujeme krokem 1.
67
Krok 2. Určíme estj a položíme h*
XEPS je modifikace kritéria EPS a podobně XEPUS je modifikace kritéria BPUS, obě modifikace jsou vysvětleny v následujícím odstavci 4.3.2.4. Zbývá vysvětlit význam příkazů h := M3(h, h*) v kroku 3 a h := M4(/í.ft*) v kroku 4. m3(A,/i*) a M4(h, h") jsou modifikační funkce, pro J\/li(h,h*) = M,4(h,h*) = h* dostáváme algoritmus, který plně odpovídá předchozímu výkladu. V programech pro řešení ODR však bývá navržená délka h* kroku ještě modifikována postupy, které jsme si formálně označili pomocí funkcí M3 a M4. Abychom byli konkrétní uveďme si, o jaké modifikace jde v implementacích metod BS32 resp. DP54 v programech ODE23 resp. ODE45 v MATLABU, viz [19]. Popis metod BS32 a DP54 uvedeme v odstavci 4.3.2.4. Následují poznámky specifikující zmiňované modifikace m3 a M4.
1. Označme hmin resp. hmax minimální resp. maximální povolenou délku kroku. Pak je třeba algoritmus modifikovat tak, aby bylo zajištěno splnění podmínky
hmin í h < /tmaa; ■
Pokud tedy v kroku 3 algoritmu nastane h* < hmm, je třeba výpočet ukončit konstatováním, že danou diferenciální rovnici program neumí s požadovanou přesností vyřešit. Přitom hmin = 16em, kde sm je tzv. strojová přesnost. Jestliže v kroku 4 algoritmu nastane h* > hmax, položí se h = hmax. Použito je hmax = 0,1(6 - a).
2. Požadujeme, aby se nová délka M3(h, h") v kroku 3 resp. Mn(h, h*) v kroku 4 příliš nelišila od stávající délky h. Konkrétně, při prvém zkracování v kroku 3
M3(h, h*) = max(ť, qmin h).
kde
f 0,5 v l 0,1 v
0, 5   v ODE23 (metoda řádu 3), ODE45 (metoda řádu 5),
při opakovaném zkracování téhož kroku pak Maf/i, h*) = 0,5 h. Při prodlužování v kroku 4
Mi(h, h") = mm(h", qmax h),
přičemž qmax = 5.
3. Bezprostředně po zkrácení kroku nesmí následovat jeho prodloužení. □
Stranou našeho zájmu zatím zůstala otázka, jak smysluplně volit toleranci e. V [34] se doporučuje zadávat ji ve tvaru e = max-fe^A^Y), «rA}, kde Sr je relativní tolerance, iVj(y) reprezentuje velikost řešení, například iV»(y) = maxdlyiH, ||yi+i||), £a je absolutní tolerance, jako standardní hodnota se doporučuje sR = 10-3 a    = 10~6.
Na začátku výpočtu je třeba vybrat délku počátečního kroku. Spokojíme se s hrubým odhadem založeným na následující myšlence: metoda ýi = rf rí y(a + h) dává lokální chybu le = y(a + h) — r) w hi (a, ti) řádu 0(h) a my předpokládáme, že při výpočtu yi
68
metodou řádu p dostaneme lokální chybu le « (le)p+l . Odtud, například podle kritéria EPS, obdržíme ||le|| w (/?,[|f(a, -ř?)!!)^1 sa £ a proto klademe
h = e[max(ea||f|||,ex)]1/(^1)/l|f(o.l)ll.
přičemž pro h < hmin změníme h na hmin a pro h > hmax změníme h na hmax.
Při programovaní algoritmu pro řízení délky kroku je třeba být opatrný. Příkazy programu je nutné uspořádat tak, aby nemohlo dojít k dělení nulou (v kroku 2 algoritmu a také při určování délky počátečního kroku).Podrobnější informace týkající se řízení délky kroku lze najít například v [34], [19], [36].
4.3.2.4. Odhad lokální chyby
je založen na použití dvou vhodně zvolených metod, z nichž jedna je řádu p a druhá řádu p + 1. V dalším uvedeme dva nejčastěji používané postupy odhadu lokální chyby a sice metodu polovičního kroku a metodu vnořených odhadů chyby.
Metoda polovičního kroku. Vysvětleme si nejdříve hlavní myšlenku. Z výchozí hodnoty yi spočteme v jednom kroku délky h přibližné řešení yj+1 a kromě toho spočteme ze stejné výchozí hodnoty y j pomocí dvou kroků poloviční délky h/2 přibližné řešení ý,+i. Odhad lokální chyby pak dostaneme vhodnou kombinací hodnot yi+1 a yi+1. Následuje podrobnější výklad.
Odvození. V kroku od xt do xi+x = Xi + h spočteme yi+1 metodou řádu p, takže pro lokální chybu platí
u(^+i)~yi+i = ftP+1'/'(^,yi)+0(^+2) = ^-
Postupujeme-li z x, do xi+i/2 := Xi + h/2 krokem poloviční délky, dostaneme přibližné řešení, které označíme yi+1/2 • Při tom vznikne chyba
/ h \ p+1 1 —
u(zi+1/2)-yj+1/2 = ^-J    V(z,,y.)+0(^+2) = ^lei+0(h^)
V dalším půlkroku získáme z výchozí hodnoty ýi+1/2 přibližné řešení yi+i s chybou 'hV+1
v(si+i)-yi+i
2)
VK^+i/2,yi+1/2)+0(/jp+2),
kde v(x) je lokální řešení procházející bodem [Xi+1/2, yi+1/2]. Protože ýj+1/2 = Yi + 0(h), z Taylorova rozvoje dostaneme ip(xi+ij2,yi+1/2) = ^/K^yi) + 0(h), takže
v(zi+1)-yí+i = (2J    ^i,Yi)+0{h^) = _lei+0(tf+2). Odečteme-li od rovnice u(si+1) — yi+1 = le; předchozí rovnici, po malé úpravě obdržíme
vi+i-y;+i
2p+i J
le,-[u(zi+1)-v(zl+1)]+0(/r+2),
69
Dále upravujme rozdíl u(xj+i) — v(x;+i). Označíme-li w = u(a;i+i)-v(a;i+i)-[u(xj+1/2)-v(a;i+i/2)],
pak
u(xm)-v(xi+1) = u(xi+i/2)-v(a:i+i/2)+w. Protože v(xi+1/2) = ýi+1/2, je
u(xi+i)-v(xi+i) = u(xí+1/2)-ýi+1/2+w = ^ríei+w+0(^+2). Naším dalším cílem je dokázat, že u; = 0{hP+2). Protože
w= /I,+'[u'(í)-v'(í)]dí= r+1[f(í,u(í))-f(í,v(ť))]dí,
užitím věty o střední hodnotě a Lipschitzovy podmínky dostaneme ||w|| < il/imax||u(í)-v(í)||,   kde  I = {xi+1/2,xi+1).
Abychom postoupili dále, použijeme známý výsledek, který tvrdí, že pro řešení u(x) a v(x) diferenciální rovnice y' = f (a;, y) platí
||u(x)-v(x)|| < ||u(a)-v(a)||e
í.(i-u)
pro x > a,
důkaz viz např. [34]. Zvolíme-li v poslední nerovnosti a := xi+1/2, i := í E / a dosadíme do předposlední nerovnosti, dostaneme
H| < ±LheLh/2\\u(xi+l/2)-v(xi+1/2)\\ = \LheLh/2\\u(xl+1/2)-y(x,+1/2)\\,
takže w = 0(/ŕ+2) a
u(xj+i)-v(xí+1)
Vrátíme-li sc zpět a dosadíme do výrazu pro ýj+i — y~j+1, obdržíme
yí+i-ýí+1 = (1~^i) ^-t^íí.+o^2) = ^íi.+o(^+2) ■
Odtud dostáváme odhad chyby, který umíme spočítat, neboť
u{xi+i)-yi+1 = kj
2"
2P-1'
(ýi+1-yi+1)+0(/r+2),
Je přirozené, abychom pro další výpočet použili přesnější řešení yi+1. V tom případě nás bude zajímat odhad chyby lej := u(xi+1) - ýť+1. Ten dostaneme snadno, neboť
u(zi+i)-ýť+i = K^»+i)-v(a;i+i)]+[v(xi+i)-yi+1],
70
a protože výraz v každé ze složených závorek na pravé straně je roven lei/2p+1 + 0(hp+2), u(xi+1)-yi+1 = ler = Iíii+0(^+2) = _Í_(yi+1-y.+1)+0(^+2).
Pro chybu metody s krokem h/2 tedy použijeme odhad 1
est,-
2" - 1 Když definujeme
(yi+i-yi+i) ■
9i+i = yi+i + estj - yi+1-
2f- 1
(yi+i-yi+i)
vidíme, že u(x,-+i)—ýj+i = 0(hp+2), tj. metoda výpočtu ýi+1 je vyššího řádu p+1. Způsob, jakým jsme z přibliných řešení yi+1 a ýj+1 řádu p zkonstruovali přibližné řešení ýj+l řádu p + 1 není nic jiného než Richardsonova extrapolace. Vyšší přesnosti lze dosáhnout také v Xi+i/2- Stačí si uvědomit, že
u(xi+1/2)-ýť+i/2 = j^+oíw*2) = \ ^(j^-y^+O^2),
a položit
yi+i/2 = yi+i/2 +j expresnější řešení ýi+i použijeme jako výchozí v dalším kroku, tj. položíme yi+l = ý,+1. Definujeme-li ještě y»+i/2 = ýi+1/2, dostaneme výsledný
Algoritmus metody polovičního kroku
ýť+i/2 = yi + \h<&{xu yt, yi+i/2, \h\ f),
yi+i = ýi+1/2 + |/i*(xj+i/2, ýi+i/2, yi+i, \h; f),
yi+i = yi + h*(x;,y;,yi+1,/i;f), (4.48)
esti = ^—^(yí+i -yi+1),
ýi+1/2 = ýi+1/2 + \ est;,   yi+1 = yi+i + est* •
Lokální extrapolace. Metoda (4.48) je ukázkou univerzálního postupu, který bývá označován jako lokální extrapolace: k základní metodě řádu p, v našem případě jde o vzorec pro výpočet yi+i, se sestaví předpis pro odhad její lokální chyby
estj = íe, + 0(//+2) = u(xi + /i)~yi+i+0(/í)3+2)
a definuje se nové řešení yi+i = yť+1 + estj, které je aproximací řádu p + 1, neboť
u(xi + h)-yi+1=0(^+2).
71
V algoritmu (4.48) první dva řádky představují výpočet základní metodou, další dva odhad lokální chyby základní metody a v posledním řádku je popsána lokální extrapolace. Výpočet bez lokální extrapolace bychom dostali, kdybychom poslední řádek v (4.48) nahradili příkazy
y.+i/2 = y<+i/21 y«+i = yt+i •
Odhad estj lokální chyby základní metody použijeme v algoritmu pro řízení délky kroku.
V případě lokální extrapolace tak vlastně řídíme délku kroku extrapolované metody řádu p + 1 pomocí odhadu lokální chyby základní metody řádu p, postupujeme tedy podle modifikovaných kritérií, která značíme XEPS místo EPS a XEPUS místo EPUS. Původní označení EPS a EPUS se ponechává pro případ, kdy se lokální extrapolace nepoužije. Řízení délky kroku extrapolované metody podle kritérií XEPS nebo XEPUS se běžně používá, funguje velmi dobře a je také teoreticky zdůvodněno.
Zamysleme se ještě nad výpočtovou náročností metody (4.48). Omezíme se na explicitní Rungovy-Kuttovy metod}7, p značí řád a s je počet koeficientů kj, tj. v každém kroku metody se funkce f počítá s-krát. Pak výpočet yi+t/2, y.+i a y«+i vyžaduje 3s — 1 vyhodnocení funkce f. Protože počítáme ve dvou krocích délky h/2 hodnoty yj+i/2 a yi+1, připadá na jeden krok (3s — l)/2 vyhodnocení funkce f. To představuje oproti základní metodě, kdy se počítá jen yi+1/2 a yi+i, zvýšení nákladů o 100 (s - l)/(2s)%, získáváme však odhad lokální chyby a navíc metodu vyššího řádu p + 1. Uveďme si dva konkrétní příklady.
a) Uvažme jako základní Rungovu-Kuttovu metodu řádu p = s = 2, například metodu EM1 nebo EM2. Pak metoda (4.48) řádu 3 potřebuje 5 vyhodnocení funkce f, tedy najeden krok připadá 2,5 vyhodnocení. Avšak pro Rungovu-Kuttovu metodu řádu 3 je s > 3.
b) Je-li základní metodou klasická Rungova-Kuttova metoda řádu 4, tj. p = s = 4, pak metoda (4.48) řádu 5 potřebuje 11 vyhodnocení funkce f, tedy na jeden krok připadá 5,5 vyhodnocení. Avšak pro Rungovu-Kuttovu metodu řádu 5 je s > 6.
K největším nedostatkům metody (4.48) patří ta skutečnost, že délku kroku lze měnit vždy až po provedení dvou půlkroků. To je neefektivní, řízení provádíme vlasně jen na 50%, velmi drahé je z toho důvodu také odmítnutí kroku. Proto se v současných programech dává přednost metodě, která bývá označována jako
Metoda vnořených odhadů 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 y, spočteme y£1 přesnější metodou a y*+1 méně přesnou metodou. Pak lokální chyba méně přesné metody splňuje
le* = u(si+1)-y*+1 = (y^-y^HW^-y^) = (y'^-y'+^+O^2) a protože le* = 0(hp+1), je esti = y"! - y*+1
72
použitelný odhad lokální chyby méně přesné metody. O vnořeném odhadu chyby mluvíme tehdy, když obě metody sdílejí společnou množinu koeficientů {k>}J=1. V tom případě je totiž získání odhadu laciné. Pokud ve výpočtu pokračujeme přesnější metodou, tj. pro y.+i = y**i = y*+1 + esti, ří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 yi+1 = y*+1. V tom připadě 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.
Dvojice Rungových-Kuttových metod se popisují pomocí rozšířené Butcherovy tabulky,
Cl	On	«12 •	au
C'l	021	022 ■	
Cs	así	«s2 ■	Q>ss
	k	k ■	■ k
	k'	k ■	■ k-
	E,	i?2	■ Es
ve které jako obvykle kj j = 1,2,...,s, přičemž
■ f (xí + hej, yi + h Yle=i ayťk«)
yí+i = y. + h XľJ=i bjkj   je metoda řádu p, y*+i =Yi + h Ej=i b"kj  je metoda řádu p + 1 a estj = h Y?j=t Ejkj  je odhad lokálni chyby,
takže E,
1,2,..
Jestliže se používá lokální extrapolace, tj. yi+i = y"x, připojíme k názvu metody značku (p + l,p) a pokud se lokální extrapolace nepoužívá, tj. yj+1 = y*+1, připojíme značku (p, p + 1). Existuje celá řada osvědčených a používaných dvojic metod. Dvě takové si uvedeme.
Bogacki-Shampine (3,2) metoda, stručně BS32 metoda. Rozšířená Butcherova tabulka BS32 metody je
Přesnější z obou metod páru je Ralstonova metoda
0				
1	1			
2	2			
■i 4	0	3 1		
1	2 9	1 3	1 9	0
	7	1	1	1
	24	4	3	8
	2 9	1 3	4 9	0
	5	1	1	1
	72	12	:')	a
L|k! + |k2 +
:yi+i
řádu 3 s vyladěnou chybou. Pomocná metoda y*+1 = yi+h [£k! + ik2 + |k3 + |k4]
řádu 2 používá kromě koeficientů ki, k2 a k3 navíc koeficient kj = f(£;+i, yj+i)- V každém kroku se tedy počítají jen 3 nové hodnoty funkce f, neboť koeficient ki ve stávajícím kroku je roven koeficientu k4 z kroku
předchozího, takže nově se počítají jen koeficienty k2, k3 a k4. Výjimkou je případ, kdy se hodnota yi+1 neakceptuje a krok se krátí. Tyto případy však nebývají časté. Zařazení koeficientu k4 do metody řádu 2 nás tedy nic nestojí, umožní však zlepšit vlastnosti této metody a tím i celého páru. Tento postup bývá označován FSAL podle anglického First Same As Last.
Hodnoty přibližného řešení pro x G (xí,Xj+i) spočteme dostatečně přesně pomocí
73
kubického Hermitova polynomu H3(x) určeného podmínkami
n3(xi)=yi,       Hi(au) = ki, H3(ii+i) = y«+i, H!,(zi+i) = k4.
Metoda BS32 je tedy skvělá: je řádu 3, v každém úspěšném kroku se pravá strana f počítá jen 3-krát a to stačí jak na řízení délky kroku tak na výpočet řešení mezi uzly Xi a xi+i. Dormand-Prince (5,4) metoda, stručně DP54 metoda je definována rozšířenou But-cherovou tabulkou
0							
1 5	1 5						
3 10	3 40	9 40					
4 5	44 45	56 15	32 9				
8 9	19372 6 561	25 360 2 187	64448 6561	212 729			
1	9 017 3 168	355 33	46 732 5 247	49 176	5103 18 656		
1	35 384	0	500 1113	125 192	2187 6 784	11 87	
	5179 57 600	0	7 571 16 695	393 640	92 097 339 200	187 2100	1 10
	35 384	0	500 1113	125 192	2187 6 784	11 84	0
	71 57 600	0	71 16 695	71 1920	17253 339 200	22 525	i 10
Metoda DP54 je typu FSAL, neboť k7 = f{xi+1, yi+i). Proto se v každém úspěšném kroku metody počítají jen koeficienty k2,..., k7, koeficient ki.byl už spočítán jako koeficient k7 v předchozím kroku. Metoda DP54 má vyladěnou chybu.
Hodnoty přibližného řešení pro x e {xí,xí+t) spočteme dostatečně přesně pomocí interpolačního polynomu H4(s) = y; + /iKBs, kde K = (k1,k2,k3,k4,ks,k6,k7).
	" 1	1 93	37	145 -i
		64	12	128
	0	0	0	0
	0	1500	1000	1000
		371	159	371
	0	125	125	375
		32	12	64
	0	9 177 3 392	729 106	25 515 6 784
	0	11 7	11 3	56 28
	. 0	3 2	-4	5 2 J
s2,	s3,	s4)T a	s — (x —	Xi)/h.
s — [s
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
[xi,Xi+l).
74
4.3.2.5. Implicitní a semi-implicitní metody
V úvodu kapitoly 4 jsme ukázali, proč má smysl používat implicitní Bulerovu metodu přesto, že je ve srovnání s explicitní Eulerovou metodou výpočetně velmi náročná. Jediným důvodem je to, že EE metoda má malou oblast absolutní stability zatímco IE metoda je A-stabilní. Ukážeme si, že u Rungových-Kuttových metod je to podobné: explicitní metody mají oblast absolutní stability omezenou, takže A-stabilní metody je třeba hledat jen mezi implicitními nebo semi-implicitními metodami. A-stabilní Rungovy-Kuttovy metody se používají prakticky výhradně k řešení tzv. tuhých problémů. Co to tuhé problémy jsou, s jakými obtížemi se při jejich řešení setkáváme a jaké metody se pro jejich řešení hodí, to vše jsou témata, se kterými se seznámíme podrobněji v samostatném odstavci 4.5. Zde si jen uvedeme několik známých a přitom poměrně jednoduchých A-stabilních Rungových-Kuttových metod.
Začneme tím, že si připomeneme, a případně nově zavedeme, potřebné pojmy. Uvažujme tedy testovací rovnici
y' = Xy,   kde A je komplexní číslo se zápornou reálnou složkou, tj. Re(A) < 0.
Aplikujeme-li na tuto rovnici obecnou Rungovu-Kuttovu metodu, dostaneme diferenční rovnici
Vi+i = R(h)yi,   kde h = h\. Funkce R(h) se nazývá funkce stability metody. |j/j+i| |Ä(A)|<1,
zřejmě platí, právě když
(4.49)
tj. metoda je absolutně stabilní pro ty hodnoty h, pro které platí (4.49). Oblast absolutní stability metody je oblast S$a komplexní roviny taková, že pro h £ &a platí (4.49), Jestliže M a obsahuje celou zápornou polorovinu, metoda je A-stabilní. Na řešení velmi tuhých problémů ani A-stabilní metody nestačí a používají se metody ještě silnějšího kalibru, tzv. L-stabilní metody. Jsou to A-stabilní metody s vlastností |iž(/í)| —> 0 pro Re(/i) -> —oo. Průnik oblasti absolutní stability s reálnou osou je interval absolutní stability.
Snadno se ověří, že pro explicitní s-stupňovou Rungovu-Kuttovu metodu je R(h) polynom stupně s a proto oblast M a nemůže být neomezená. Odtud plyne závěr: explicitní Rungova-Kuttova metoda není A-stabilní. Pro s = p < 4 funkce stability R(z) nezávisí na volbě konkrétní metody. Skutečně, lokální řešení testovací rovnice splňuje u(xi+h) = ehxyi a protože u(xj + h) - yi+í = [ehx — R(h\)]yi = 0(hp+1), R(z) je Taylorův polynom funkce ez stupně p. Oblasti absolutní stability se zjišťují numericky a jejich obrázky lze najít v monografiích věnovaných numerickému řešení ODR, např. v [17], [34]. Intervaly absolutní stability lze prezentovat jednodušeji. Pro s = p = 1,2, 3,4 dostáváme postupně intervaly (-2; 0), (-2; 0), (-2,51; 0), (-2,78; 0) a pro DP54 metodu interval (-3,30; 0).
Pro implicitní s-stupňové Rungovy-Kuttovy metody je R(h) racionální lomená funkce, čitatel stejně jako jmenovatel je polynom stupně nejvýše s. Tak například pro IE metodu
yi+i = Yi+htíi,   kde   kx = f(aľj+ft,yť+/iki),   tj. pro  yi+i = yi+hí(zi+i,y,
75
dostaneme R{z)
1 - z
takže pro z = a+ib, a < 0, je |Ä(«)|2 =
l + \a\Y + b2:
tj. \R(z)\ < 1 a \R{z)\ —► 0 pro a —► -oo. IE metoda je tedy L-stabilní.
Lichoběžníková metoda. Velmi známou metodou, která nechybí v žádné učebnici věnované numerickému řešení ODR, je lichoběžníková metoda, stručně TR metoda podle aglického trapezoidal rule,
yi+i = y i + 3 h[í (xí , yť) + f (xi+1, yť.
(4.50)
Za svůj název vděčí způsobu odvození: rovnici u'(ai) = f(x,u(x)) integrujeme v mezích od Xi do xi+x a        i(x,u(x)) dx spočteme přibližně lichoběžníkovou metodou, tj.
U(xi+i) - u(Xi) = \h[í (Xi, u(x;)) + f (Xi+U u(.Ti+, ))] + 0(/j3) .
TR metoda je jednokroková implicitní metoda řádu 2. kterou lze zapsat také jako semi-implicitní Rungovu-Kuttovu metodu:
y«+i =yť+|/i(ki+k2),   kde  ki = f(xr,yť), k2 = f(x«+/i,yž+|/i(ki+k2)), Funkce stability TR metody je 2 + z
--,   takže pro z = a + ib, a < 0, je |ii(z)|2
R(z)
(2-H)2 + &2 (2+|a|)2 + 62!
tj. |iZ(z)| < 1, \R{z)\ -+ 1 pro a -» -oo. TR metoda je tedy A-stabilní, L-stabilní však není. TR metoda je základem programu ODE23T MATLABu.
TRX2 metoda. Jestliže provedeme lichoběžníkovou metodou dva kroky stejné délky h/2, tj.
yi+i/2 = y i + \h[i(xi,yi) + f(sj+1/2,ym/2)],
y;+i = y.'+i/2 + \h[f(xi+i/2,yi+i/2) + f(aľi+i,yi+i)],
dostaneme semi-implicitní Rungovu-Kuttovu metodu, známou jako TRX2 metoda, kterou popíšeme pomocí rozšířené Butcherovy tabulky a vzorců
ki = f(xi,yť),
k2 = í(xí + \h, y{ + i/i(ki + k2)),
k3 = f {xí + h,yi + |A(ki + 2k2 + k3)),
yi+i =yi + iA(k1 + 2k2 + k3),
esti = ift,(-k1 + 2k2-k3).
76
0	0	0	0
1	1	1	o
2	4	4	
l	1	1	1
	4	2	4
	1	1	1
	4	2	4
	1	2	1
	6	3	6
	1	1	1
	12	6	12
TRX2 metoda je A-stabilní metoda řádu 2, používá se jako (2, 3) metoda, tj. bez lokální extrapolace, přesnější metoda řádu 3 se použije jen k získání odhadu est, lokální chyby, k výpočtu přibližného řešení ji použít nelze, neboť není A-stabilní. TRX2 metodu je třeba kvalitně implementovat, viz [14],
TR-BDF2 metoda. TRX2 metoda není L-stabilní. Existuje však podobná metoda řádu 2, která L-stabilní je. Je to metoda označovaná jako TR-BDF2: tato metoda používá nejdříve lichoběžníkovou metodu, odtud TR v jejím názvu, a pak metodu zpětného derivování řádu 2, zkráceně BDF2 metodu, viz odstavec 4.4.3. Výpočet je popsán vzorci
TR:      yi+7 = ys + \^h[f {xh yť) + f (si+7, yi+7)],
2—t              1 I-7 BDF2:---yi+1 - -^ryť+7 +        y« = Af {xi+1, yi+1),
1-7"
7(1-7)'
7
0	0	0	0
7	d	d	0
1	w	w	d
	w	w	d
	1(1 ' t")	\{3w + 1)	¥
	|(1 -Aw)	1 3	-\d
kde 7 = 2 - s/2 + hrf. TR-BDF2 metodu lze zapsat jako semi-implicitní
Rungovu-Kuttovu metodu pomocí rozšířené Butcherovy tabulky a vzorců
ki = i{xi,yi),
k2 = f (xí + 7A, yt + hd{ki + k2)), k3 = f(xi + h,yi + h(wki + wk2 + dk3)), yi+i = Yi + h(wk! + wk2 + dka), est, = |A[(1 - 4tf)ki H- k, - 2dk»|,
kde d = 7/2 a w = V2/A. Metoda řádu 3 se opět používá jen pro odhad est; lokální chyby. Podrobnosti o TR-BDF2 metodě lze najít v [14]. Metoda je implementována v MATLABu jako program ODE23TB.
Je známo, že existují implicitní a semi-implicitní Rungovy-Kuttovy metody, které jsou A-stabilní nebo i L-stabilní a jsou přitom libovolně vysokého řádu, viz například [17]. To je zajímavé z teoretického hlediska, v praxi se však takové metody vyšších řádů příliš nepoužívají, neboť algoritmy na nich založené jsou výpočtově velmi náročné.
4.3.2.6. Odhad globální chyby
je založen na obdobné myšlence jako odhad lokální chyby metodou polovičního kroku, viz odstavec 4.3.2.4. Stejné je to, že yť+1 spočteme v jednom kroku délky h z výchozí hodnoty y*. Stejné je také to, že yi+1 spočteme pomocí dvou kroků poloviční délky h/2, jediný rozdíl je v tom, že tentokrát vyjdeme z výchozí hodnoty y;. Položíme tedy yo = yo = v a počítáme
yi+i = y i + h^{xi,Yi,yi+1, A; f),
yV+1/2 = ýi + |A*(xt, ýi, ýi+i/i, \h; f),
y«+i = y>+i/2 + |A$(xi+i/2,yť+1/2,yi+1, \h;í) .
77
Dále využijeme toho, že pro metodu řádu p platí yi+1=y(z1+0 + ciŕP + O(ŕP+1),
yi+i=y(^+i) + ci(f)?+0(í/p+1),
kde H = max h je nejdelší krok a c; je konstanta nezávislá na H (jde o hluboký teoretický výsledek, jehož důkaz lze najít např. v [34]). Odečtením dvou posledních výrazů obdržíme
yi+i-yi+i = (y-i)c r|Y+0(tf'+1)
a odtud pro globální chybu e;+1 := y(xi+i) — yj+1 přesnější metody dostaneme
ě,+i = y(^+1)-yí+i = -oí (jj+o^1) = ^(yM-fi^+oiH"^),
takže jako odhad globální chyby ei+1 lze použít výraz 1
Ai+1 :=
2" - 1
(y>+i -yi+i)
Považujeme-li za výsledky přibližné hodnoty yi+i/2, y.+i vidíme, že odhad globální chyby zvyšuje výpočtové náklady přibližně o 50%. Postup, který jsme použili ke stanovení odhadu Aj+1 globální chyby, se nazývá globální extrapolace. Na závěr připojíme dvě poznámky:
1) Přibližné řešení y,+i lze zpřesnit: položíme-li
ýV+i = yi+i + Ai+i,
dostaneme y(xi+1) - yi+i = 0{HP+1).
2) Délku kroku řídíme pomocí odhadu est; lokální chyby le, = u(xi+i) - yi+1 méně přesné metody.
4.4. Lineární mnohokrokové metody
V této kapitole se budeme zabývat metodami, které počítají přibližné řešení yi+i v uzlu Xí+i pomocí dříve spočtených hodnot y;,yi-i,y;_2> • ■ • a odpovídajících hodnot f(xi;yi), f (.Ti_i,yi_i),f(aj_2,yt^2), ■ ■ ■ pravé strany diferenciální rovnice. Tyto hodnoty jsou znovu použity tak, abychom získali yi+1 s vysokou přesností pomocí jen několika málo nových vyhodnocení funkce f(x, 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é metody, stručně LMM.
78
4.4.1. Obecná lineární mnohokroková metoda
Pro konstantní délku h kroku je lineární mnohokroková metoda předpis aoy;+i+aiyjH-----Vakyi+1_k = h[j30{(xi+i, yi+i)+/3ifH-----h^fi+i-fc],
(4.51)
ze kterého počítáme yi+1. Přitom o,- a Bj jsou reálné koeficienty, které formuli jednoznačně určují, a ij je zkrácený zápis pro í(xj,yj). V dalším budeme předpokládat, že a0 / 0. Jestliže alespoň jeden z koeficientů ak nebo fik je různý od nuly, metoda je k-kroková.
Pro /3o Ý 0 Je nová hodnota y;+1 určena implicitně, hovoříme proto o implicitní metodě, pro /30 = 0 máme metodu explicitní. Abychom v implicitní metodě určili yi+i, musíme vyřešit rovnici
8 1 ^
yi+i = V(yi+i),   kde   ip(z) = h—f(xi+1. z)H--Y] [hfijU+i-j - ^jVi+i-j] ■
J=l
Pro dostatečně malé h má tato rovnice jediné řešení. Skutečně, užitím Lipschitzovy podmínky (4.4) dostaneme ||¥>(u) - <r(v)|| < /i(/3o/ao)£||u - v||, takže y(z) je pro h(Po/cto)L < 1 kontrakce a jednoznačná existence yi+i je důsledkem věty o pevném bodu.
LMM lze použít, jen když jsou zadány startovací hodnoty y0, yi,..., yk-i- ya 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.
Všimněte si, že výsledky získané formulí se nezmění, když formuli normalizujeme tak, že všechny koeficienty vydělíme nenulovým číslem. Používají se různé typy normalizací. Často se volí a0 = 1. My použijeme jinou normalizaci zahrnující všechny koeficienty aj. K tomu účelu se nám bude hodit (první) charakteristický polynom q{&) definovaný předpisem
o{9) = a0ek + oi!©*"1 + -.. + ak.
V dalším budeme za normalizovanou považovat takovou formuli, pro kterou j/(l) = 1.
Lokální diskretizační chyba. Chybu, které se dopustíme v jednom kroku, měříme pomocí lokální diskterizační chyby lte; definované rovnicí
3=0
nebo ekvivalentně
k
E«jy(^+i-j)
j=0
k
h^2^y'{xi+i-j)+ltei.
3=0
Lokální diskretizační chyba je tedy chyba, která vznikne, když do normalizované formule (4.51) dosadíme místo přibližného řešení yt+i-j přesné řešení y(xi+í-j).
79
Chybu ltej rozvineme do mocnin h tak, že provedeme Taylorův rozvoj všech funkcí okolo bodu xi+l. Pro řešení samotné dostaneme
y(xť+i-j) =y(xi+i) - jhy'(xi+1)
Hh)'< 2!
p+1
= y(xi+1) + £ ^yW(ím) + 0(0
a podobně pro derivaci
y'(zi+i-j) = y'(xi+i) - jhy"{xi+1) + ■ ■ ■ = ]T -
-i/>)ä
-yW(xi+1)+0(Ap+2).
Když tyto rozvoje dosadíme do vzorce pro lokální diskretizační chybu, obdržíme
ltei = ^Cs^yW(xi+1) + 0(^+2);
kde
Co = £ ai   a   C.= (-l)'5]
j=0 j=o
J «j , 3° Pi L s! (s-1)!
pro s = 1, 2,... ,p + 1.
Řád metody. Stejně jako u jednokrokových metod řekneme, že LMM je řádu p, jestliže ltej = 0(/ip+1), tj. pokud C0 = Ci = ■ • • = Cp = 0 a Cp+i ^ 0. Všimněte si, že výraz pro lte* zůstane v platnosti, když uzel xi+í nahradíme bodem x* = + 0(h), například uzlem Xj. Pro metodu řádu p se člen Cp+i/tp+1y^+1)(a;ť) nazývá hlavní člen lokální diskretizační chyby a konstanta Cp+i je tzv. chybová konstanta. Ta závisí na normalizaci LMM, takže pokud chceme srovnávat přesnost dvou formulí, musíme si být jisti, že jsou stejně normalizovány.
Podmínky C, = 0 pro s < p jsou ekvivalentní podmínkám pro dosažení řádu p u Rungových - Kuttových metod. Zatímco u Rungových - Kuttových metod představovaly podmínky řádu nelineární rovnice pro koeficienty metody, u LMM dostáváme v důsledku linearity jen rovnice lineární.
Pro LMM lze definovat také lokální chybu lej jako lokální diskretizační chybu lokálního řešení U;(z) splňujícího = i(x,Ui), u^Xi) = y,. Význam takto definovaného lokálního řešení je však omezen tím, že nelze předpokládat Uj(:ri+1_3) = y»+i-j, j = 2,3,..., k. Na rozdíl od jednokrokových metod tedy lokální chyba lej v případě LMM nemá bezprostřední názornou interpretaci a proto ji nebudeme používat.
Metoda konzistententní s ODR. Dá se ukázat, že každá konvergentní LMM musí být řádu p > 1, viz např. [38]. LMM řádu alespoň 1 se nazývá metoda konzistentní (s diferenciální rovnicí). V dalším budeme uvažovat jen konzistentní LMM. Všimněte si, že podmínku C0 = 0 lze ekvivalentně zapsat jako g(l) = 0. Abychom podobně vyjádřili také druhou podmínku konzistence, tj. podmínku Cj, = 0, definujeme druhý charakteristický polynom
(7(0) = 8o0k + Piek
Z podmínky
* Pí;
Ci = - ^2[jUj+Pj] = 0   dostaneme   <j(1) = ^ft- = - y] ja.
3=0
3=0
3=0
a dále
= X^-JK-e*-''-1   a odtud   ^(l)=*5^aJ-53ioj = fce(l)+íT(l).
3=0 j-o j=l
Konzistentní LMM je tedy ekvivalentně charakterizována dvěma podmínkami na charakteristické polynomy: ^(1) = 0 a g'(l) = <r(l).
D-stabilita. Jednokroková metoda řádu p > 1 je vždy konvergentní, pro LMM metodu to však neplatí. Abychom si objasnili v čem je problém, zabývejme se řešením triviální úlohy
y'(x) = 0  pro x e (0,1) s počáteční podmínkou y(Q) = 0.
Pokud startovací hodnoty zvolíme přesně, tj. yr = 0 pro 0 < r < k, pak yi = 0 pro všechna i. Prozkoumejme co se stane, když startovací hodnoty nepatrně pozměníme, tj. řešíme-li úlohu
k
Y^otjZi+i-j = 0
3=0
pro 0 < r < k pro  i > k — 1.
kde Ar jsou malé poruchy. Zi je řešení diferenční rovnice s konstantními koeficienty. Taková rovnice se řeší podobně jako diferenciální rovnice s konstantními koeficienty: řešení navrhneme ve tvaru zt = 7O, kde 7 a C jsou zatím neurčené konstanty, a dosadíme ho do diferenční rovnice. Tak dostaneme
5>7ci+1-J' = 7Ci+1-fc£^
:7Ci+1~*í?(C) = o,
3=0
3=0
takže je-li ( kořen charakteristického polynomu g{0), je zť řešení diferenční rovnice.
Vybereme-li startovací hodnoty Ar = 7Č7, r = 0,1,..., k — 1, pak zt = 7O je řešení pozměněné úlohy. Nechť TV je počet stejných dílků, na které je interval (0,1) rozdělen. zN je tedy aproximace přesné hodnoty y{l) = 0, která by měla být tím lepší, čím je N větší. Avšak \yN - zN\ = \zN\ = -y\C\N, takže pro |C| > 1 je lim^^ \yN - zN\ = 00. To je ale zcela nepřijatelné: z malých počátečních poruch vznikla koncová porucha, jejíž velikost s rostoucím N neomezeně roste.
Je-li C, = reIlp, I = -J—í, komplexní kořen velikosti r > 1, předchozí analýza zůstává v platnosti, pokud bychom se však chtěli vyhnout počítám s komplexními čísly, stačí uvažovat z, = 71-' cos(i(p). Vyloučit je třeba také případ, kdy C = eIv je násobný kořen velikosti 1, pak totiž z{ = 7i cos(itp) je řešením diferenční rovnice a opět limjv^oo \zN\ = 00.
81
Předchozí analýza nás opravňuje zavést nový pojem: řekneme, že LMM je stabilní ve smyslu Dahlquista nebo stručně D-stabilní, jestliže pro všechny kořeny G charakteristického polynomu g{0) platí |C,| < 1, a pokud |C| = 1, pak je kořen Q jednoduchý.
Konvergence. Nyní již můžeme vyslovit nutnou a postačující podmínku zaručující konvergenci LMM: LMM je konvergentní, tj. max0<j<w ||y(xj) - y»|| -+ 0 pro h -> 0, právě když je D-stabilní řádu alespoň 1. Důkaz viz např. [38]. Toto tvrzení lze stručně vyjádřit schématem
konzistence + D-stabilita konvergence Pro rychlost konvergence LMM platí následující
Věta (o rychlosti konvergence LMM). Uvažujme D-stabilní LMM řádu p > 1. Jestliže startovací hodnoty zadáme s chybou řádu 0{hp), pak metoda je řádu p, tj. pro globální diskretizační chybu e, = y(ij) — y i platíc = 0{hp).
Precizní formulaci a důkaz této věty lze najít např. v [38], [17]. Předpokladem její platnosti je dostatečná hladkost pravé strany f, konkrétně je třeba, aby funkce f(x,y) měla spojité derivace až do řádu p včetně. Pokud f má spojité derivace jen do řádu s < p, pak lze pro globální chybu dokázat pouze řád 0(h").
Absolutní stabilita. D-stabilita souvisí s konvergencí a tedy s délkou kroku h -¥ 0, v anglicky psané literatuře se proto místo D-stabilita často používá ekvivalentní termín „zero-stability". Při praktickém výpočtu však počítáme s konkrétní nenulovou délkou kroku a proto je třeba požadovat od numerické metody více než jen D-stabilitu. Vhodným nástrojem pro formulování takového požadavku je testovací úloha
y' = \y   pro x e (0, oo), y(0) je libovolné, A je komplexní číslo, Re(A) < 0 .
V tom případě totiž y(x) -> 0 pro x -> +oo a proto je žádoucí, aby při pevně zvolené délce h kroku y, -¥ 0 pro xt = ih -> oo. Testovací úlohu řešíme LMM a dostaneme lineární diferenční rovnici
Předpokládáme y i — 7c1, dosadíme do diferenční rovnice a obdržíme
Y,(*3-hW3hC+1-j = 7ť+1-*l>*-/lAftK*"i = 7C+1~Mc,£) =0.
3=0 j=o kde h = hX a
Jfc
■K(zth) = '^2l{o-j-h0j)zk~'i - e{z)-ha{z) i=o
je charakteristický polynom diferenční rovnice, označovaný také jako polynom stability LMM. Je-li (s kořen polynomu stability ír(z, h), pak Q je zřejmě řešení diferenční rovnice.
82
1
Je-li kořen (s násobný, pak další řešení diferenční rovnice jsou inQs pro n = 1,2,..., m, -1, kde m,s je násobnost kořene £s.
Vidíme tedy, že yt 0 pro i ->• 00, právě když |£s| < 1. Označme 31 a množinu těch bodů h komplexní roviny, pro které všechny kořeny f4 polynomu stability in (z, h) leží uvnitř jednotkového kruhu \z\ < 1. Pak h € 3ČA => Vi -> 0 pro i -> 00. Oblast MA nazýváme oblastí absolutní stability LMM. Jsou-li všechny kořeny polynomu stability ir(z,h) jednoduché, pak zřejmě < takže definice absolutní stability LMM je ve shodě s definicí absolutní stability zavedené pro jednokrokové metody, viz (4.18) a (4.49).
Velikost oblasti !%A absolutní stability je významnou charakteristikou LMM. Jak uvidíme v odstavci 4.5, metody s velkou oblastí absolutní stability se používají při řešení tuhých problémů. Pro řešení tuhých problémů se výtečně hodí metody zpětného derivování, seznámíme se s nimi v odstavci 4.4.3. Metody prediktor-korektor založené na Adamsových metodách, viz následující odstavec 4.4.2, mají poměrně malou oblast absolutní stability a proto se pro řešení tuhých problémů nehodí. To ale vůbec neznamená, že to jsou metody špatné, naopak, jsou to velmi dobré a hojně používané metody pro řešení problémů, které tuhé nejsou.
Přibližné řešení mezi uzly Xi a xi+i získáme interpolací z hodnot yj+i,yi,... ,yi+i-*.
4.4.2. Adamsovy metody
Integrací diferenciální rovnice (4.4) od X; do xi+1 = X; + ht dostaneme
y{x
i+i)~y(xi) = í
Jxí
f(x,y(x))dx.
Funkci f(x, y(x)) aproximujeme pomocí interpolačního polynomu Pv~i(x) stupně v - 1 a dostaneme metodu
y»+i
n+ P„-i(í)dí.
J Xi
(4.52)
4.4.2.1. Adamsovy-Bashforthovy metody
Jestliže interpolační polynom určíme podmínkami
P»-i(zť+i-j) = U+i-j ,   j = 1> 2,..., v,
dostaneme ^-krokovou Adarnsovu-Bashforthovu metodu, stručně ABľ metodu. Když vyjádříme interpolační polynom v Lagrangeově tvaru
3=1
obdržíme ABv metodu ve tvaru
v
yi+i = yi+in E Pij fi+i-j
k = 1
X — Xj+i-k Xi+l-j   ' ^i+l-fc
j = \,2,...,u,
(4.53)
83
kde
1   fx*+i  JL f
■ ^i+i-*
fc = i
■dť,   j = 1,2,
.,v.
Koeficienty /3*j závisejí na vzájemných vzdálenostech uzlů £j+i_j a proto se musejí počítat v každém kroku. Jsou však známy postupy, jak tyto koeficienty počítat efektivně, navíc výpočtové náklady lze snížit tím, že délka kroku se změní, jen když se vyplatí spočítat novou sadu koeficientů. Druhou, tradiční možností je pracovat s krokem konstantní délky h,
Xi+i-j = xi+1-jh,   j = 1,2,...,v,
pak koeficienty f3*j na h nezávisejí,
k-
ľ1 "
s-=í n
•'O L_,
k = 1
k-j
ds,   j = 1,2,
To ale neznamená, že bychom měli celý výpočet provádět s pevnou délkou kroku. Jestliže se rozhodneme pokračovat z x, do Xi + h* s novou délkou kroku h*, pak přibližné hodnoty tf+i-k Prayé strany f (x, y{x)) v uzlech iť - h*, Xi - 2h",..., xt — {v — l)h* stačí spočítat interpolací z původních hodnot fj+i-rs, tj. určíme
ÍT+1-* = P*-i(xi-(k-l)h*) = l3{xi-{k-l)h*),   k = 2,3,...,v.
i=i
Jeden krok Adamsovy-Bashforthovy metody vyžaduje jen jediné vyčíslení funkce f. Skutečně, předpokládejme, že jsme v uzlu xt, právě jsme spočetli yj a máme k dispozici y»-i,yi-2,... v předchozích uzlech xí-i,Xí-2, ... a dále také hodnoty fi_i, fj-2, — Abychom mohli provést další krok, musíme spočítat í{xi,yi), pokud je to zapotřebí tak také koeficienty    •, a nakonec užitím metody (4.53) vypočteme yi+1.
Snadno ověříme, že nej jednodušší Adamsova-Bashforthova metoda, kterou dostaneme pro v = 1, tedy AB1 metoda, není nic jiného než nám dobře známá jednokroková explicitní Eulerova metoda. První vícekrokovou metodou je AB2 metoda. V tom případě
J i-Xi+hi
t Xi.
-dť :
1    hi _   1 ľ
1 + 2XT' ^-hjxi
Xi+hi £
Xi
dt =
i h
"2Vi
kde hi = xi+i tvaru
Yi+1 = yi + h
Xi—l Xi
- Xi-i. Pro konstantní délku kroku nabývá AB2 metoda
3f *f
Formule pro konstantní délku kroku. V následující tabulce uvádíme koeficienty ABv metody pro v = 1,2,..., 6 a pro konstantní délku kroku:
v		ä.2			ä,5	
1	1					
2	3 2	1 2				
3	23	16	5			
	12	12	12			
4	55	59	37	9		
	24	24	24	24		
5	1901 720	2 774 720	2 616 720	1274 720	251 720	
6	4277	7 923	9982	7298	2 877	475
	1440	1440	1440	1440	1440	1440
Tvar (4.53) AB metody je důsledkem použitého Lagrangeova tvaru interpolačního polynomu. Pokud používáme krok konstantní délky h, je účelné použít Newtonův tvar interpolačního polynomu
pĽ_1(x) = fi+^ivfi+---+
(x     Xj)(x     Xj-i)...(x — Xj+2-v) yjv-\f
hr-\v-\)\ i!
kde operátor zpětné diference je
V% = U,   Vfi = íi-fi_1,   a obecně   VJ'+1fj = V(V'fi) ■.
Adamsova-Bashforthova metoda, která vznikne integrací Newtonova tvaru interpolačního polynomu, je tvaru
yi+^yi+^Tl-iV^fj, i=i
kde 7j = 1 a pro j > 1
^+lí (t-xi)...(t-xi+1.j)
(4.54)
J Xi
h>j\
dť
s+l)...(s+j-l)ds.
Všimněte si, že 73* nezávisí na v. Toho lze využít: označíme-li yf+1 řešení získané ^-krokovou metodou a y^1 řešení získané metodou (ľ + l)-krokovou, pak
y£} = y\+1+hyyii.
Dá se ukázat, že koeficienty 7J lze počítat pomocí rekurze
7o* = l,   t; + ^T*-i + ^7;-2+--- + t^I7o = 1,   i = 1,2,.... Z ní snadno spočteme několik prvních koeficientů
3	0	1	2	3	4	5	6
	1	1	5		251	95	19 087
u	1	2	12	s	720	288	6D480
84
85
T
Adamsova-Bashforthova formule (4.53) resp. (4.54) je normalizovaná. Skutečně, charakteristický polynom g(Q) = 6" - 0"-1, takže g'(l) = 1.
Lokální diskretizační chyba. Přesnost metody měříme pomocí lokální diskretizační chyby lte* definované rovností
v
y(ii+i)-y(xi) = fcs^Z^. f(sť+1_í,y(a;m_í))-(-ltef.
3=1
Jestliže do polynomu P„-i(x) dosadíme í(xi+i-j,y(xW-j)) místo íi+i-j vidíme, že
Piz-ito+i-j) = f(xi+1_j,y(a;í+1_j)) = y'(zi+i_j),   j = 1,2,...,v.
Ze známé věty o Lagrangeově tvaru zbytku interpolačního polynomu plyne, že pro každé í a každou složku m existuje číslo £m(r) e (xí+i-v,Xí) takové, že
y[m(t) = P^m(t)+^y%+1)(Ut))f[(t~xi+i-k) ■
To zapíšeme stručněji ve vektorové formě
y'(í) = P^(í) + ~ y<"+1>(*) n(t-o;i+1_fc),
V' k=l
hvězdička přitom vyjadřuje, že každá složka vektoru y'"+1' je vyhodnocena pro jiný argument. Pak
y(xi+1) = y(x4)+ ľ'+1 Pu-x(t) át+ ľ'1 i y<"+1>(*) Il(í-^m-
I dt.
Protože člen niUií* ~~ xi+i-k) nemění v intervalu (xí,xí+i) znaménko, užitím první věty o střední hodnotě integrálu dostaneme
lte* = P''ň(*-*m-»)
dť.
Jestliže H je nejdelší z kroků vystupujících ve formuli, pak pro lokální diskretizační chybu platí odhad
llltenmiy^ll//"*1,
takže ABi/ metoda je řádu v a protože v > 1, je ABi> metoda konzistentní. Je-li délka kroku konstanta h. výraz pro lokální diskretizační chybu nabude tvaru
lte* =y("+1)(*)/l"+1i /1n(ť+fc)dí = 7;y("+1)(*)^+1.
8C
Protože y("+1)(*) = y("+1)(xj)+0(/i), je 7; yí"+1'(a;i)/i'/+1 hlavní člen lokální diskretizační chyby a 7* je chybová konstanta. Lokální chybu lze snadno odhadnout. Přidáme-li ještě uzel £;_„ a sestrojíme polynom P„(x), dostaneme chybu
y(z*+i
ř(a;ť)+ /
J Xi
Pv{t)át+0(H"+2).
Odečteme-li tuto rovnici od rovnice y(xi+i) = y{xi) +       P„_i(r) dť + lte*, dostaneme
lte* = f '+1pv(t)dt- ľ'^ Pv^{t)át+0{Hv
J Xi J Xi
2)-
Tyto výrazy se spočtou snadno, když délka kroku je konstanta h a interpolační polynom je vyjádřen pomocí zpětných diferencí. V tom případě
lte* = ft7;Vf(^,y(a:1)) + 0(0.
Proto je přirozené odhadnout lokální diskretizační chybu výrazem A7* Vfj.
Další vlastnosti ABi/ metod. Interpolační polynom vyjádřený pomocí zpětných diferencí je zřejmě neobyčejně vhodný jak pro sledování lokální chyby, a následně k řízení délky kroku, tak ke změně řádu v metody. Tyto přednosti se neomezují jen na konstantní délku kroku. Na obecné síti lze použít formálně stejný Newtonův interpolační polynom, rozdíl spočívá jen v tom, že zpětné diference se nahradí zpětnými poměrnými diferencemi, viz např. [35]. Integrační koeficienty, ekvivalenty 7*, závisejí nyní na vzájemné poloze uzlů Xi,Xi-i,..., a proto se musejí v každém kroku počítat znovu. Výpočet na nerovnoměrné síti je proto v detailech jiný, podstata však zůstává stejná. Příslušné algoritmy jsou důkladně propracovány a úspěšně implementovány.
Adamsova-Bashforthova metoda je D-stabilní: charakteristický polynom ABp metody o(0) = 6" — 6"_1 má (v - l)-násobný kořen 0 a jednoduchý kořen 1. ABu metoda je tedy konvergentní.
Adamsova-Bashfortova metoda je explicitní: do pravé strany formule (4.53) dosadíme známé hodnoty a přibližné řešení yi+1 je určeno. Jednu věc jsme však zatím zamlčeli, totiž jak určit startovací hodnoty yo, y1;..., y„_i. Z počáteční podmínky vždy máme yo = r}. Zbývající startovací hodnoty můžeme určit třeba tak, že je spočteme dostatečně přesnou Rungovovou-Kuttovou metodou, teoretické výsledky říkají, že by měla být řádu alespoň 0(h"). To je však spíše klasický přístup, který měl oprávnění v dobách, kdy programy pracovaly s pevnou délkou kroku. Jestliže však délku kroku řídíme tak, aby lokální chyba nepřesáhla zadanou toleranci, můžeme dostatečně přesné startovací hodnoty získat i pomocí méně přesné metody, pokud délky kroků volíme dostatečně malé. Tomuto přístupu se v moderních programech dává přednost a používají se tzv. samostartující algoritmy: startovací hodnoty pro v-krokovou Adamsovu-Bashforthovu metodu se spočtou opět pomocí Adamsových-Bashforthových metod, yi pomocí AB1, y2 pomocí AB2 až nakonec y„_i pomocí AB(t/ — 1).
Adamsovy-Bashforthovy metody se samostatně téměř nepoužívají. To však neznamená, že by se nepoužívaly vůbec, nej významnější uplatnění nacházejí, spolu s Adam-sovými-Moultonovými metodami, v algoritmu známém jako prediktor-korektor. Kvalitní
87
implementace Adamsových metod technikou prediktor-korektor přitom patří v současnosti k základnímu programovému vybavení pro řešení ODR. V následujícím textu si proto nejdříve zavedeme třídu Adamsových-Moultonových metod a pak si vysvětlíme, jak pracuje algoritmus prediktor-korektor.
4.4.2.2. Adamsovy-Moultonovy metody
jsou metody implicitní. Odvodí se podobně jako Adamsovy-Bashforhovy metody, rozdíl je jen v tom, že interpolační polynom P„_i(a;) stupně v - 1 je nyní určen podmínkami
P*-l(xí+i-j) = fi+i-j,   3 = 0,1, • •., v-1 ■
Dosazením Lagrangeova tvaru polynomu P„_i do (4.52) dostaneme Adamsovu-Moulto-novu metodu, stručně AMv metodu
yi+i = yi+ftj/9„,of(ii+i,y.+i)+'»iXl^'fi+1-í ■ ^4'55^
Člen hiP„r0f(xi+i,yi+1) jsme úmyslně vypsali zvlášť mimo sumu, neboť právě tento člen „je zodpovědný" za to, že metoda je implicitní: počítaná hodnota yj+i je nejen vlevo, ale také vpravo jako druhý argument funkce f.
Pro v = 1 dostáváme implicitní Eulerovu metodu a pro v = 2 lichoběžníkovou metodu. Obě tyto metody jsou jednokrokové, vícekrokovou metodu dostaneme až pro v > 3, AM3 je už dvoukroková. Obecně, AMv metoda je max(l, v — l)-kroková. Formule pro konstantní délku kroku. V následující tabulce uvádíme, koeficienty AMv metody pro v = 1,2,..., 6 a pro konstantní délku kroku:
v	Ayo	ÄM	&,2	&,3		A*
1	1					
2	1 2	1 2				
3	5 12	8 12	1 12			
4	9 24	19 21	5 24	1 24		
5	251 720	616 720	264 720	106 720	19 720	
6	475 1440	1427 1440	798 1 440	482 1 440	173 1440	27 1440
Jestliže interpolační polynom P„_i vyjádříme v Newtonově tvaru
D I \ f , X ~ TI t , {X - Xi+i){x - Xj)...(x - Xi+3-v) Vv-lr , P„_1(x) = fi+l + —^- Vfí+1+- ■ ■ +-jpj^ _ ľy-■- V 1,4-1
a dosadíme do (4.52), dostaneme AMv formuli ve tvaru vhodném pro odhad chyby a změnu řádu:
yi+i = yi + h ]T li-\V3 lfi+i.
(4.56)
kde 7o = 1 a pro j > 1
li
-i:
(t-xi+1)...(t-xi+2_i) dt_l f1
hij\
i ŕ
= 7Í / («-l)*-(»+J-2)d«.
Koeficienty 7,- lze opět počítat pomocí rekurze
7o = l, 7i + -7j-i + -7;_2+-"+— 7o = 0, j = l,2, několik prvních koeficientů uvádí následující tabulka
j	0	1	2	3	4	5	6
1}	1	1 2	1 12	1 24	19 720	3 160	863 60 480
Z definice koeficientů 7, a 7? odvodíme, že pro jejich vzájemný vztah platí
7í = 7J-7?-i   Pro  i>l. (4.57)
Lokální diskretizační chyba lte, Adamsovy-Moultonovy metody se definuje předpisem
v-l
y{xi+i)-y(xi) = k^&jf(zm-í>y(z»ti-j))+ltej
3=0
a podobně jako dříve pro Adamsovu-Bashforthovu metodu odvodíme
iteí=iy(^)(*) r+in(í-zi+1_fc)dt
Jx> k=Q
Jestliže H je nejdelší z kroků vystupujících ve formuli, pak pro lokální diskretizační chybu platí odhad
||ltei||<||y<"+1>||ir+1,
takže AMv metoda je řádu v. Je-li délka kroku konstanta h,
ltof = y^1^*)^1- /■1ff(ť+fc-l)dí = 7!,y(''+1)
1/1 JO l_n
(*)h"+1,
takže 7», y("+1) (xi)hv+1 je hlavní člen lokální diskretizační chyby a 7„ je chybová konstanta. Pomocí zpětných diferencí lze, podobně jako jsme to provedli pro Adamsovu-Bashforthovu metodu, odvodit
lte, = H, V"f(Ij+1,y(ii+1)) + 0(fc"+2),
takže lokální diskretizační chybu můžeme aproximovat výrazem est; = hyv V"fj+1,
Snadno se ověří, že AMv metoda je normalizovaná, D-stabilní, konzistentní a tedy konvergentní.
4.4.2.3. Metody prediktor-korektor
Adamsova-Bashforthova a Adamsova-Moultonova formule je pro zvolené v téhož řádu v. Pro v = 1 je 71 = -TÍ a tedy implicitní Eulerova metoda IE = AM1 je přibližně stejně přesná jako explicitní Eulerova metoda EE = AB1. Pro v > 1 je však |7„| značně menší než a proto je Adamsova-Moultonova metoda výrazně přesnější než Adamsova-Bashforthova metoda téhož řádu v > 1.
Abychom určili yi+i, musíme řešit rovnici (4.55), tj. rovnici
v-1
Vi+i = v(yť+i). kde      =yi+^,of(aľi+i,z)+/»J3^'fi+1--'-
Použít můžeme metodu prosté iterace: zvolíme počáteční aproximaci y,-°\ a postupně
počítáme y£\ = v»(y.
= 1,2,.... Pro dostatečně malé h metoda konverguje, jak jsme ukázali v odstavci 4.4.1. Jestliže počáteční aproximaci yí°\ určíme pomocí Adamsovy-Bashforthovy metody, provedeme jen několik málo iterací
y^ = v{y^]), s = i,2,..,s,
a nakonec položíme
y»+i = y<+i, f,-+i = f(xi+u y«+i),
dostaneme metodu prediktor-korektor, kterou schématicky označujeme P(EC)SE. Přitom P značí předpověď počáteční aproximace explicitní Adamsovou-Bashforthovou metodou. C korekci implicitní Adamsovou-Moultonovou metodou a E vyhodnocení pravé strany f(ij+1,yj|)1). Někdy se klade fm = f^^y,^"11), takové schéma se označuje jako P(EC)S. Nejčastěji se používá schéma PECE, kdy se korekce provede jen jednou a pravá strana se počítá dvakrát. Konkrétně pro dvojici AB2-AM2 a PECE schéma organizujeme výpočet následovně:
P      AB2     y*+1 = yť + \h{3Si - S-i)
e ff+i = f(Ki+i.yí+i)
C     AM 2    yi+i =y, + |ft(f*+1 + fi)
E fí+1 =f(a:i+1,yj+i)
Efektivní tvar metody prediktor-korektor. Zabývejme se dále obecným případem, kdy prediktor je ABu metoda
(4.58)
.7 = 1
90
a korektor je AM;/ metoda
y.+i = yi + /iX^73-iVJ'-1fi+1.
3=1
Diference v Adamsově-Moultonově metodě obsahují členy fť+1 = f(xi+i,yi+1), které jsou v PECE schématu nahrazeny členy f*+1 = f(xi+1,y?+1), takže
V
yi+i = yi + /i^7hv)'"1t41.
3=1
Diference VJ_1fř+1 se snadno spočtou pomocí diferencí Ví'~"1fi, které máme uloženy z předchozího kroku:
V%i    = f^i-f.-Vf.^-Vfi,
Tyto diference se použijí jen pro výpočet y<+1 a k odhadu lokální chyby, takže se nemusejí
schovávat. Až spočteme yi+1 a fi+1, vytvoříme diference
V°f;+i     =  fí+i!
V!fi+1     = fi+1-fj = V0fi+i-V0fž,
V'+1íi+l  = V^fi+1 - W,.
Nové diference VJfí+i mohou být uloženy na místo starých diferencí V^fj, čímž se ušetří paměťové místo.
Pokud použijeme dvojici ABv-AUv, dostaneme pro určení yi+1 sympatický vztah Yi+i = yľ+i + Ä7Í-1 Vfá*+1. (4.59) Abychom to dokázali, odečteme formuli prediktoru od formule korektoru a obdržíme
Prvé členy z obou sum sloučíme do tvaru A7oV°f<;i-/l7o*V°fi = H*V1f^1,
přičemž jsme využili toho, že 70 = 7J. Když k tomu přidáme druhé členy z každé sumy, dostaneme
/^v^Vi + /V/iV1^ - fcrív1* = fc(70* + 71)v1^ - ätív1* =
a7í(v1fÍ1-v1fi) = /l7l*v2fj;
\
91
Zde jsme použili rovnost (4.57), tj. 7? = 7J_i+7j- Opakováním tohoto postupu dostaneme nakonec formuli (4.59).
Adamsova-Moultonova metoda řádu v nepoužívá hodnotu fi+i^„, kterou Adamsova-Bashforthova metoda řádu v používá. To nás přivádí k myšlence použít Adamsovu-Moultonovu metodu řádu ^ + 1 jako korektor pro Adamsovy-Bashforthovu metodu řádu v, Výraz pro yj+i dostaneme, když k pravé straně rovnosti (4.59) přičteme člen hyvVi^+v Užitím vztahu 7* = 7*_j + 7^ dostaneme
yi+i = y?+i+h-yl V"f;*+1.
(4.60)
Ze srovnání vzorců (4.60) a (4.59) je patrné, že použití korektoru vyššího řádu nás nic nestojí a proto je výpočet yť+1 podle vzorce (4.60) dosti populární.
Lokální diskretizační chyba. V dalším prozkoumáme lokální diskretizační chybu metody prediktor-korektor pro dvojici ABp*-AMp a schéma PECE. Připomeňme si, že
p*
lte? = y(xi+1)-y{xi)-h^2l3*.jf{xi+1^,y(xi+1.j))
3=1
je lokální diskretizační chyba prediktoru ABp* a
pí
ltej = y(xi+i)-y{xi)-hl5Pfli{xi+i,y{xi+i))-h^PPjf(xi+i-jly(xi+1^j))
je lokální diskretizační chyba korektoru ABp. Lokální diskretizační chyba lte; predikto-ru-korektoru je proto určena rovností
lte; = y(zí+i) - y{xi) - h0PiOi{xi+i,y(xi+i) - ltef) -h^Pvé f fo+W. yfa+i-j)) =
j'=i
= lte? + A0P)O[f(íCj+i, y(zť+i)) - f(xi+1, y(xi+i) - ltef)].
Odtud užitím věty o střední hodnotě dostaneme df{xi+u*)
lte; = lte^ + hppfi-
dy
-ltef,
(4.61)
přičemž hvězdička vyjadřuje, že prvky Jacobiovy matice jsou vyhodnoceny pro různé hodnoty argumentů.
Protože prediktor je řádu p* a korektor řádu p, je ltef = 0(/ŕ*+1) a lte,c = 0(hp+l), takže lte; = 0(/i,+1), kde q = min(p,p* + 1), tj. formule prediktor-korektor je řádu q. To je důvod, proč se nejčastěji volí prediktor i korektor stejného řádu p* = p: pak je totiž řád q páru prediktor-korektor roven řádu p korektoru a navíc, hlavní člen lokální diskretizační chyby páru je stejný jako hlavní člen lokální diskretizační chyby korektoru. To je skvělé, stačí jediná korekce a máme přesnost korektoru! Je-li řád prediktoru o jedničku nižší než řád korektoru,tj. p* = p - 1, je řád páru stále roven řádu korektoru, hlavní člen lokální diskretizační chyby páru je však již jiný než hlavní člen lokální diskretizační chyby korektoru.
Všechny úvahy jsme prováděli pro PECE schéma, které se v praxi zdaleka nejvíc používá. Není těžké zjistit, že každá korekce redukuje vliv prediktoru o jeden řád. Tak například pro prediktor řádu p* = p — 1 a dvě korekce, tj. výpočet podle schématu P(EC)2E, je hlavní člen lokální diskretizační chyby páru opět stejný jako hlavní člen lokální diskretizační chyby korektoru.
Odhad lokální diskretizační chyby. Vyjádření lokální diskretizační chyby lte, ve tvaru (4.61) má význam teoretický, pro praktický výpočet však potřebujeme odhad est; lokální diskretizační chyby, který se dá snadno spočítat. Ukažme si proto, jak se takový odhad dá zkonstruovat.
Pro odvození odhadu est; je účelné přijmout lokalizační předpoklad, totiž že dříve spočtené hodnoty yť,..., yi+i_„ jsou přesné, tj. y, = y{xt),yi+i_„ = y(xi+i_„). Tento předpoklad je (pro v > 1) zcela nerealistický, umožní nám však uhodnout, jak odhad est; může vypadat. Dodatečně se pak ukáže, že jde o odhad korektní, plně použitelný pro praktické výpočty.
Nechť ýi+1 je řešení získané formulí řádu i/+la yť+1 je řešení získané formulí řádu v, pak lokální diskretizační chybu y(z;+i) — y;+i můžeme přibližně odhadnout pomocí
yi+i-yi+i = (y(2i+i)-yi+i)-(y(Ti+i)-yVt-i) = y(^+i)-y;+i+0(fe"+2),
neboť lokální diskretizační chyba y(x;+i) -y;+i je 0(/i1/+2). Přesnější řešení y;+i můžeme spočítat pomocí Adamsova-Moultonova korektoru řádu v + 1 podle formule (4.60), tj.
Jestliže y;+1 spočteme pomocí ABi/-AMí/ páru metodou PECE, dostaneme pomocí (4.59) a (4.57)
yj+i-yi+i = ^V^rAl.V"^, = fcfcVfk ■ Odhad estj lokální diskretizační chyby proto volíme jako
est; = h-y„Vf*+
(4.62)
V [34] je dokázáno, že tento odhad je správný i bez lokalizačního předpokladu.
Jestliže jako přibližné řešení v vezmeme přesnější hodnotu y;+i = y;+i + est; říkáme, že jsme použili metodu prediktor-korektor s lokální extrapolací. Tedy výpočet podle (4.58), (4.59) je bez lokální extrapolace a výpočet podle (4.58), (4.60) je s lokální extrapolací.
Řízení délky kroku a řádu metody. Kvalitní programy založené na metodách prediktor-korektor mění jak délku kroku tak řád metody. Tak například program ODE113 v MATLABu používá ABf-AMi/ prediktor-korektor typu PECE s lokální extrapolací pro i/ = l,2,...,12.
Změna řádu se provádí současně se změnou délky kroku. Algoritmy tohoto typu se se označují jako VSVO algoritmy (variable step variable order). Základní myšlenka je jednoduchá. Předpokládejme, že jsme spočetli yi+1 metodou řádu v s lokální extrapolací a krokem délky h. Určíme odhad E„ := hy^Vf^ lokální diskretizační chyby. Jestliže
92
93
I
||E„|| > e, jde o neúspěch a výpočet yi+1 je třeba opakovat, v opačném případě pokračujeme výpočtem přibližného řešení yj+2- V každém případě je však třeba stanovit novou délku kroku a nový řád. Pokud jde o řád, připustíme jen řády v — 1, u nebo v + 1. Dá se ukázat, že chybu každé z těchto metod lze odhadnout pomocí E„_i, EĽ nebo E„+i. Novou délku kroku stanovíme podobně jako pro jednokrokovou metodu. Podle kriteria XEPS dostaneme pro každou ze tří kandidujících metod
h\ = eA(£/||Et||)1/(*+1)   pro k = v - 1, v a v + 1,
a největší z čísel A*, A*+1 určí jak novou délku kroku tak nový řád. Tedy, je-li
největší h*v, řád zůstává a jako novou délku kroku vezmeme h* = /t*, je-li největší , zvětšíme řád o jedničku a pokračujeme s krokem délky h* = a je-li největší řád o jedničku snížíme a pokračujeme s krokem délky h* = h*v_x.
Znovu upozorňujeme, že jsme zde vyložili jen základní ideu algoritmu VSVO, detaily je třeba hledat v dokumentaci úspěšných implementací, odkazy na ně lze najít např. v [17], [34].
Start metody není žádný problém: začneme metodou řádu 1, počáteční délku kroku určíme např. tak, jak je to popsáno v odstavci 4.3.2.3, a algoritmus VSVO se už sám rychle vyladí na správné hodnoty jak délky kroku tak řádu metody.
Absolutní stabilita. Doposud jsme se nezmínili o absolutní stabilitě Adamsových formulí. A-stabilní jsou jen první dvě Adamsovy-Moultonovy metody, totiž AM1, což je implicitní Eulerova metoda, a AM2, což je metoda lichoběžníková. Adamsovy-Moultonovy metody mají výrazně větší oblast absolutní stability než Adamsovy-Bashforthovy metody téhož řádu. S rostoucím řádem se oblasti absolutní stability zmenšují. Pro ilustraci uvádíme v následující tabulce dolní meze a intervalů absolutní stability (a, 0) pro AB^ a AMi/, v = 1,2,3,4:
v	1	2       3 4
ABi/	-2	i          6 3 '■        U 10
KVív	—oo	-oo   —6 —3
Pro metody prediktor-korektor typu P(EC)SjE oblast absolutní stability roste, když zvětšujeme počet S iterací: pro S = 0 je stabilita nejmenší, odpovídá stabilitě ABv metody, a pro S —► oo se stabilita blíží svému maximu, tj. stabilitě AMv metody. Tak například interval absolutní stability pro AB4-AM4 PECE metodu je (-1,25; 0), tedy větší než pro AB4 metodu, avšak menší než pro AM4 metodu. Interval absolutní stability AB4-AM4 PECE metody je více než čtyřikrát větší než interval absolutní stability AB4 metody, proto si (za jistých okolností) můžeme dovolit až čtyřikrát delší krok, a vzhledem k tomu. že PECE metoda v každém kroku vyhodnocuje pravou stranu dvakrát, můžeme oproti AB4 metodě docílit až 50 % úspory v počtu vyhodnocení pravých stran.
Oblasti absolutní stability pro metody AB^-ÁMi/ PECE metody pro v < 4 jsou uvedeny např. v [17]. Obecně lze konstatovat, že oblasti absolutní stability ABľ-AMi/ PECE metod nejsou příliš veliké a s rostoucím řádem v se dále zmenšují. Proto se programy založené ABi/-AMi/ PECE metodách, jako třeba program ODE113 v MATLABu, nehodí
pro řešení úloh se silným tlumením (stiff problémů), viz odstavec 4.5. Existuje však celá řada jiných úloh, které nejsou stiff, a pro ně jsou Adamsovy metody vhodné.
4.4.3. Metody zpětného derivování
Při řešení tuhých problémů je třeba pracovat s metodami, které se vyznačují velkou oblastí absolutní stability. Metody zpětného derivování (stručně BDF podle anglického backward differentiation formulas) tuto vlastnost mají. Pro označení ^-krokové metody zpětného derivování použijeme zkrácený zápis BDFi/ metoda. BDFf metodu odvodíme tak, že v rovnici
y'(xi+l) = ff>j+i,y(xj+1))
nahradíme derivaci ý(xi+1) pomocí derivace Pl(xi+1) interpolačního polynomu P„(x) stupně v procházejícího body [xi+í, ym], [x„ y,:], ..., [xi+l.v, y,-+1_„]. Tak dostaneme metodu
K(zí+i) = f{xw,yi+1). Vyjádříme-li interpolační polynom v Lagrangeově tvaru
P,(s) *,(*),   kde   ^.(x)=  J]     *-"+'-*    , j = i>2.....v,
j=0 k = 1   xi+l~J Xj+i_fc
dostaneme BDFí/ metodu jako předpis
ov,oy1+i+ov,iy<+- • •+cwyi+i„Ľ = hif(xi+uyi+l),   kde   avó = /iiij(a:i+1).
Metoda BDF^ je implicitní jA-kroková metoda. Speciálně BDFÍ je jednokroková implicitní Eulerovu metoda, jak se snadno ověří.
Formule pro konstantní délku kroku. V následující tabulce uvádíme koeficienty BĽFu metody pro ľ=l,2,... ,6 a pro konstantní délku kroku:
v	°V,0			tti/,3		«!/,5 Q'i/,6
1	1	-1				
2	3 2	-2	1 2			
3	11 6	-3	3 2	1 3		
4	25 12	-4	3	4 3	1 4	
5	137 60	-5	5	10 3	5 4	1 5
6	147 60	-6	15 2	20 3	15 4	6 1 5 6
Interpolační polynom P„(x) můžeme zapsat také v Newtonově tvaru n h" v\
94
95
Pro konstantní délku kroku
_d_ (x - Xi+i)(x - Xj)...{x - Xj+2-j)
á~x h> j\
1_
h j'
takže BDFť metoda zapsána pomocí zpětných diferencí je tvaru
(4.63)
j=l
Tento zápis je sice formálně velmi atraktivní, pro výpočty se však příliš nehodí. yi+i se totiž počítá iterací, tj. vyjde se z počáteční aproximace y|+\ a postupně se počítají aproximace y^\. To ale znamená, že se v každé iteraci musejí přepočítat zpětné diference VJy,-+i Pr0 všechna j = 1,2,..., v. Ukážeme si, jak se tomu můžeme vyhnout.
Efektivní tvar BDFi> metody. Začneme počáteční aproximací yj+1 = y|°\. Určíme ji extrapolací z předchozích hodnot yi+i-j v uzlech xi+i-j, j = 1,2,..., v + 1, do uzlu Newtonův tvar interpolačního polynomu je
-  , . (a - Xj)(x - Xj-l)...(x - sm-„)
Q*(z) = yi-t—-— Vyi+- • •+--V yj,
takže předpovězená hodnota
/í Ä(2h)... (ľ/l) _„ V^r,i
y*+1 = Q,(*m) = yt+^Vyi+... v  ^   Vyi = £vyi ■
j=0
To je velmi vhodná počáteční aproximace, neboť zpětné diference máme schovány z předchozího kroku. V dalším se nám bude hodit vyjádření rozdílu počítané hodnoty yi+i a její počáteční aproximace y*+1:
yi+i - y*+i = (yi+i - yO - vy. - vVi-----v»yi =
= (Vyi+i-Vyi)-V2yi-----V"yi =
= (V2yi+i - V2y.)-----V"y; = ■ ■ • = V+1yi+i ■
Je účelné stanovit také počáteční aproximaci f?+1 pro hodnotu pravé strany f(xi+1,yi+i) jako Q|,(aľi+i). Po jednoduché úpravě odvodíme
1 " 3 1
fT+i = QLtet+i) = r E 5ívjyi - kde nové konstanty 5i= Yl i ■
" j=l h=l
Všimněte si, že 6i = 1 a že 5j - &j-i = - pro j > 1. Užitím tohoto vztahu obdržíme
f(x,+1, yi+i) =i £ ±V'yM = iVyi+1 + \ J2{5j ~ ^-i)VVí+i =
i=2
96
= jt^yw - Jí>-i(V-W - ^'"Vi) =
j'=] j=2 1
Kombinací tohoto a předchozích výsledků dostaneme
f(zi+i,yi+i)-f*« = \svvyi+1~svvvyi = i^v+Vi+i = Ji-frm-yW ■
Výsledná algebraická rovnice pro určení yť+1 proto může být zapsána ve tvaru
I(x,+l,yl+l)-f*+, -I<5„(yi+1_y'+1) = 0. (4.64)
BDF metody se používají prakticky výhradně pro řešení tuhých problémů a jak uvidíme v odstavci 4.5.2, yi+1 se počítá řešením rovnice
g(z) := f(xi+1, z) -f*+1 - i í„(z-y?+1) = 0
pomocí zjednodušené Newtonovy metody. Pro takový výpočet se rovnice g(z) = 0 skvěle hodí: z0 = y*+1 je dobrá počáteční aproximace, y*+1 a f*+1 se určí snadno pomocí zpětných diferencí uchovaných z předchozího kroku, v rovnici se nevyskytují zpětné diference s neznámou z.
Lokální diskretizační chyba ltej je definována vztahem
lte« = £a.,y(zi+i-j)-hif(z;+i,yi+i) = /íť(PÍ,(:rj+i)-y'(:rť+i)),
kde P„(x) interpoluje y(x): P^Xj+i-j) = y(x<+1_,), j — 0,1,...,v. Podle (2.4) pro derivaci Lagrangeova interpolačního polynomu P„(x) v uzlu xi+i platí
y'(xi+l) = ^M+T^^'Wn^-^i-t) ■ Jestliže H je maximální délka kroku, pak zřejmě
\\K(^i)-y,M\\<j^\\y^1)\\s'.
Dosadíme-li odtud do výrazu pro lte; vidíme, že BDFv je metoda řádu v. Speciálně pro konstantní délku h kroku dostaneme
Ite^-^y'^H*)/^1.
97
Lokální diskterizační chybu můžeme také určit pomocí vztahu (4.63), do kterého místo yi+i-j dosadíme y(xi+w), tj.
V -Vjy{xi+1) = hf{xi+uy(xi+i))+ltei.
3=1 J
Jestliže horní mez v sumě zvětšíme o jedničku, dostaneme metodu řádu v + 1, takže
£ -VJy(%!) = hf(xi+1,y(xi+i)) + 0(h"+2). Odečtením obou rovnic pak získáme
ite, = --L-v-'+Vte+O+o^2) ■
v +1
Proto můžeme jako odhad est; lokální diskretizační chyby ltej použít výraz
est, = —t-V+1yi+1- (4-65) v +1
Další vlastnosti BDFf metod. Formule zpětného derivování jsou normalizované. Skutečně, druhý charakteristický polynom er(6) = Q" a protože BDFv metoda je řádu v>l, je g'(l) = (j(l) = 1. Formule zpětného derivování jsou D-stabilní jen pro v < 6, pro v > 7 jsou BDFv metody nepoužitelné. V dalším se proto omezíme na BDFľ metody pro v < 6. Ty jsou vždy konvergentní.
Všimněte si, že chybová konstanta -l/(v+l) je pro v > 1 podstatné větší než u Adam-sových metod stejného řádu, tj. BDFi/ metody jsou značně méně přesné než Adamsovy metody. Metody zpětného derivování však mají jednu ohromnou přednost, která plně ospravedlňuje jejich používání, a tou je skutečnost, že mají neomezenou oblast &a absolutní stability. Pro metody BDF1 a BDF2 oblast 3&A absolutní stability obsahuje celou zápornou komplexní polorovinu, St a 2 {ž|Rez < 0}, tj. tyto metody jsou A-stabilní. Abychom mohli popsat oblast absolutní stability zbývajících BDFť metod, zavedeme si jeden nový pojem. Řekneme, že numerická metoda je A{a)-stabilní, a e (0,tt/2), jestliže oblast SI a její absolutní stability obsahuje segment {z | tt - a < arg z < tt 4- a} komplexní roviny. BDFy metody (pro v < 6) jsou A(a)-stabilní, příslušné úhly ve stupních udává následující tabulka
v	1 2	3 4	5 6
a	90° 90°	88° 73°	52° 18°
Úhel 18° metody BDF6 je příliš malý a proto se tato metoda obvykle nepoužívá. V programu ODE15S v MATLABu jsou implementovány metody zpětného derivování řádů 1 až 5. Program ODE15S je typu VSVO, tj. volí optimální délku kroku i řád metody. Protože stabilita s růstem v klesá, uživatel může určit, že se budou používat jen metody řádů v < vmax < 5. Výběr délky kroku a řádu metody vychází z odhadu (4.65) lokální
chyby, významnou roli ale hraje také rychlost konvergence zjednodušené Newtonovy metody použité k řešení nelineární rovnice (4.64).
Ao-stabilní metoda. Metoda, jejíž interval absolutní stability je interval (—oo, 0), se nazývá Ao-stabilní. A(a)-stabilní metoda je A0-stabilní pro každé a, tj. BDFi/ metody jsou Ao-stabilní. A0-stabilní metody se používají např. pro řešení soustav ODR, které vzniknou po provedení prostorové diskretizace nestacionární úlohy vedení tepla, viz odstavec 6.2.
4.5. Tuhé problémy
Při řešení mnoha praktických úloh vznikají soustavy diferenciálních rovnic, které vykazují jev označovaný v anglicky psané literatuře jako „stiffness", samotné soustavy (problémy) se opatřují přívlastkem „stiff", tj. hovoří se o „stiff systems (problems)". Doslovný překlad termínu „stiffness" je „tuhost", budeme tedy hovořit o tuhosti a o tuhých soustavách (problémech), někdy ale použijeme „originální" označení a budeme hovořit o stiff soustavách (problémech). V česky psaných učebnicích [37], [38] se pro tuhé soustavy používá opisné označení „úlohy se silným tlumením", proto příležitostně použijme také toto označení.
Slovo „stiffness" vstoupilo do literatury pravděpodobně z teorie řízení: řídicí mechanismus, jehož chování může být popsáno stiff soustavou ODR, má tendenci přejít do ustáleného stavu a setrvat v ném jako tuhý mechanismus.
4.5.1. Co je to tuhost?
Kvalitní programy založené na explicitních Rungových-Kuttových nebo Adamsových metodách jsou obecně velmi efektivní. Jestliže však pomocí nich chceme řešit tuhý problém, jsou najednou značně neefektivní, dokonce natolik, že jejich použití je zcela nepraktické: pro pevnou délku kroku dostáváme numerické řešení, které (často) osciluje a (vždy) explozivně roste; pokud použijeme automatické řízení délky kroku, lze sice tuhý problém explicitní metodou vyřešit, ale jen za cenu toho, že budou použity extrémně malé kroky. Takové chování tuhých problémů se zdá být záhadné a matoucí. Abychom pochopili o co jde, je třeba tuhou soustavu nějak popsat. Než se k tomu dostaneme, řekněme si, co budeme chápat pod pojmem stabilní (počáteční) problém.
Stabilní počáteční problém. Zhruba řečeno, problém je stabilní vzhledem k počáteční podmínce, když malá změna v počáteční podmínce Způsobí malou změnu řešení. Podobně řekneme, že problém je stabilní vzhledem k pravé straně, když malá změna pravé strany (diferenciální rovnice) způsobí malou změnu řešení. A souhrnně řekneme, že problém je stabilní, když je stabilní jak vzhledem k počáteční podmínce tak vzhledem k pravé straně.
Abychom mohli popsat stabilitu přesněji, zaveďme si jeden nový pojem: řekneme, že funkce f(x,y) splňuje jednostrannou Lipschitzovu podmínku, jestliže platí
(f(x,u)-f(x,v),u-v)<^ř||u-v||2,
kde (•, ■) je nějaký skalární součin generující normu || ■ ||, tj. || • ||2 = (■, ■)• Je-li L Lip-schitzova konstanta, pak vždy můžeme položit iř= L > 0. Zajímavá situace však nastane, když existuje Jžř< 0, elementární příklad je funkce f(x, y) = Jŕy pro .T < 0. Diferenciální
98
99
rovnice, jejíž pravá strana f splňuje jednostranou Lipschitzovu podmínku s konstantou ,Sŕ< 0, se nazývá disipativní.
Stabilita řeší otázku, do jaké míry se řešení úlohy
u' = f (x, u),      u(a) dané,
může lišit od řešení úlohy
v' = f (x, v) + g (x), v(a)dané.
když f splňuje jednostrannou Lipschitzovu podmínku s konstantou JĚf Odpověď poskytuje následující nerovnost
||v(x)-u(x)|| < ||v(a)-u(a)||e^-">+ f ||g(í)||^(l-t) di,
J a
(4.66)
jejíž důkaz lze najít např. v [34]. Uvedený odhad je zvláště pěkný pro disipativní problémy (iř < 0), neboť v tom případě jsou exponenciály ohraničeny 1, takže jde nepochybně o stabilní problém a to na libovolně dlouhém intervalu (a, b). Speciálně pro g = O ze (4.66) plyne
||v(x2)-u(x2)|| < ||v(zi)-u(a;i)||,      kde a < xx < x2 je libovolné.
(4.67)
Stabilita numerické metody. Měli bychom si také říci, co rozumíme stabilitou numerické metody. Zhruba řečeno, numerická metoda je stabilní, jestliže, aplikována na (libovolný) stabilní počáteční problém, má tuto vlastnost: malá změna startovacích hodnot způsobí malou změnu numerického řešení (pro fc-krokovou metodu jsou startovací hodnoty yo,yi, ■ ■ • lYfc-i)- Zpřesňujících definic je celá řada, pro disipativní počáteční problém lze stabilitu numerické metody definovat (s přihlédnutím k (4.67)) např. (viz [17]) takto :
Definice. Necht {uť} a {vť} jsou dvě numerická řešení téhož (ale jinak libovolného) disi-pativního počátečního problému, spočtená zvolenou k-krokovou metodou, avšak pro různé (ale jinak libovolné) startovací hodnoty, tj. {uQ,Ui,..., xzfc—i}- / {v0lvi,..., vfc_i}. Pak řekneme, že numerická metoda je stabilní, pokud
HUí+i-Vi+ill^llUi-Vil,
k-l,k,..
■>vi+l-*J ■
kde U; = (uf, uf_1:..., uf+l_k)T, Vi = (vf, vl
V dalším se budeme podrobně zabývat stabilitou speciálního modelového problému, který popisuje následující
Věta (o stabilitě). Uvažujme diferenciální rovnici
y' = Jy + g(*), (4-68)
kde J je konstantní matice, která má úplný systém vlastních vektorů. Nechť u(x) a \(x) jsou dvě řešení soustavy (4-68) s počátečními podmínkami u(a) = r) a v (a) = r) + S. Nectí vlastní čísla matice 3 jsou {\}}. Jestliže Re(A3) > 0 pro nějaké \j, pak existuje 6 takové, že ||u(x) - v(x)|| -4 co exponenciálně pro x -4 co. Jestliže Re(Aj) < 0 pro
100
1
všechna j, pak pro každé 6 je \\u(x) - v(x)|| stejnoměrně omezená vzhledem k x. Navíc, pokud Re(Aj) < 0 pro všechna j, pak pro každé S platí \\u(x) - v(x)\\ -> 0 exponenciálně pro x —> co.
Důkaz. Jestliže matice J řádu n má úplný systém vlastních vektorů znamená to, že matice J má n lineárně nezávislých vlastních vektorů vy, j = 1,2,..., n, (příslušných k vlastním číslům Aj) a tedy matice V = (vi,v2,..., v„) je regulární. Pak V^JV = D = diag{Aj je diagonální matice. Rozdíl z(x) = v[x) - u(x) splňuje z' = Jz. V nové proměnné w = V_1z máme
w' = V"
Odtud dostaneme
V^JW^z = Dw.
10,(1) = ^"-^(0),
j = l,2,...,n.
Jestliže Re(Aj) > 0 a Wj(a) ^ 0, pak Wj(x) exponenciálně roste. Jestliže Re(Aj) < 0, pak \u>j{x)\ < \wj(a)\ pro všechna x, a když Re(AJ) < 0, pak Wj(x) -> 0 exponenciálně. Přechodem k původním proměnným z = Vw obdržíme tvrzení věty. □
Úloha (4.68) je tedy stabilní vzhledem k počáteční podmínce, právě když všechna vlastní čísla matice J mají nekladnou reálnou složku. Dá se ukázat, že když J je symetrická matice s nekladnými vlastními čísly, pak úloha (4.68) je disipativní (pro skalární součin (u,v) = uTv dostaneme Jžf= maxAj < 0).
Nyní se již můžeme pokusit o charakteristiku tuhého problému ODR. Přesná definice tuhého systému k dispozici není, proto postupně uvedeme několik tvrzení, která tento pojem spíše objasňují než definují. Dosti výstižné a přitom poměrně obecné je
Tvrzení 1. Tuhý problém je stabilní problém s Lipschitzovou konstantou L takovou, že L(6-o)»l.
Pro lineární problém (4.68) je l > g(3), kde q{3) = maxj |A,| je spektrální poloměr matice J. Skutečně, je-li w = u - v vlastní vektor příslušný vlastnímu číslu A, pak
||(Ju+g)-(Jv+g)|| = ||Jw|| = ||Aw|| = |A|||w|| < I||n-v||.
Jestliže má matice J navzájem různá vlastní čísla, můžeme řešení úlohy (4.68) zapsat ve tvaru
n
y(x) = Yicje^xvj+p{x),
3=1
kde \j jsou vlastní čísla matice J, v j jsou odpovídající vlastní vektory, p(x) je partikulární řešení nehomogenního systému a c,- jsou konstanty, které určíme z počátečních podmínek. Předpokládejme, že Re(Aj) < 0, j = 1,2,..., n. Pak pro obecné řešení w(x) homogenního systému -w7 = Jw platí
n
w(x) =      CjeXjXVj -4 0      pro x -4 oo.
3=1
101
Z toho důvodu mluvíme o funkci w(x) jako o přechodovém řešení a. o funkci p(x) jako o ustáleném řešení. Označme dále Xmax a \min taková vlastní čísla, že
|Re(Amal)j > |Re(Aj)| > |Re(Amin)|,   j = 1,2,...,n.
Je-li naším úkolem spočítat ustálené řešení, musíme řešit soustavu alespoň do toho bodu, v němž je nejpomaleji klesající exponenciála v přechodovém řešení už zanedbatelná. Cím menší je tedy |Re(Amiíl)|, tím na delším intervalu je třeba soustavu řešit, takže délka b - a je přímo úměrná číslu |l/Re(Amj„)|. Jestliže je číslo |Re(Ama3:)| velké, je velká také Lipschitzova konstanta L. Proto je (v souladu s trvzením 1) přirozené považovat podíl
S =
|Re(AmM)!
|Re(Amin)| '
nazývaný poloměr tuhosti, za měřítko tuhosti soustavy (4.68). Můžeme proto vyslovit další Tvrzení 2. Problém (4-68) je tuhý, jestliže všechna vlastní čísla matice 3 mají zápornou reálnou část a poloměr tuhosti S je velký .
Tvrzení 2 poněkud zužuje obecnější tvrzení 1. Podle tvrzení 1 lze totiž problém (4.68) považovat za tuhý také v případě, že vlastní čísla matice J jsou velká ryze imaginární čísla. Tak přicházíme k definici nového pojmu: řekneme, že problém (4.68) je oscilatoricky tuhý, jestliže vlastní čísla {Aj} matice J splňují podmínky
Re(Aj) < 0,      |Re(Aj)| < 1,      (b-a) max |A,| » 1.
Věnujme se nyní vztahu mezi tuhým problémem a numerickou metodou jeho řešení. Uvažujme homogenní problém (4.68), tj. problém
ý = Jy >    y(a) = v ■
Pak w = V_1y je řešením ekvivalentního problému w' = Dw,      w(a) = V-177,
(4.69)
(4.70)
kde (stejně jako v důkazu výše uvedené věty) V = (vi, v2,..., v„) je matice vlastních vektorů matice J a D = diag{Aj} je diagonální matice odpovídajících vlastních čísel. Nechť yť, i = 0,1,..., je numerické řešení problému (4.69) a w;, i = 0,1,..., je numerické řešení problému (4.70), přičemž obě řešení jsou získána stejným numerickým výpočtem, tj. stejnou metodou pro stejné délky kroků. Pak není těžké prověřit, že Wj = V_1yí Jak pro Rungovy-Kuttovy metody tak pro LMM.
Je-li Re(Aj) < 0, pak y(x) -> 0 pro x -y co, tj. úloha (4.69) je stabilní (vzhledem k počáteční podmínce). Stabilitu numerické metody stačí studovat na formálně jednodušším problému (4.70). Ten se rozpadá na samostatné úlohy
: XjWj ,
přičemž Wj(x) = Wj(a)ex'<-X °' 4 0 pro x přibližné řešení uij^ w Wj(xí) platilo uijj -
j = l,2,...,n,
-+ co. Je proto přirozené požadovat, aby pro -> 0 pro i —> co, a to pro každé j = 1,2,..., n.
102
1
Když si vzpomeneme na definici absolutní stability zjistíme, že náš požadavek je splněn, právě když (obecně komplexní) čísla h = hXj leží v oblasti e S&a absolutní stability užité numerické metody. Jestliže je oblast 3£A omezená, je omezení, které zmíněný požadavek klade na velikost integračního kroku tím přísnější, čím delší je interval integrace (a, b). Je-li |Re(AmM)| podstatně větší než |Re(Ami„)|, můžeme se dostat do velmi nepříjemné situace vyžadující řešit na dlouhém intervalu diferenciální rovnici a užít při tom integrační krok malé délky splňující hXj e MA Vj. Potíže, které vznikají při numerickém řešení tuhých problémů ilustruje následující
Příklad 4.2. Uvažujme počáteční úlohu
0
-1000
1
-1001
x e (0, i),
l/i(0) 2/2(0)
1
-1
(4.71)
Řešení je yi = e x, y2 = —e x. Matice soustavy
J =
0
-1000
1
-1001
-1,
má vlastní čísla Ai = -1, A2 = -1000, spektrální poloměr g(J) = 1000, v || • 1^ normě platí ||J(u — v)||i < ||J||i||u —v||i, přičemž ||J||i = 1002, takže Lipschitzovu konstantu L lze vybrat tak, že 1 000 < L < 1 002. Připomeňme si, že podle tvrzení 1 tuhost problému reprezentuje číslo LL
Úlohu (4.71) jsme řešili explicitní BS32 metodou (program ODE23 v MATLABu) a implicitní TR metodou (program ODE23T v MATLABu). Oba programy jsme použili se stejným požadavkem na přesnost (konkrétně e r = 10~3 a e a = 10"6). BS32 je (explicitní) Rungova-Kuttova metoda řádu 3, která má interval absolutní stability (—2,51; 0). Proto délka h kroku musí splňovat podmínku hX2 < 2,51, tj. h < 2,51 ■ 10~3. TR metoda je implicitní Rungova-Kuttova metoda řádu 2 a je A-stabilní, takže délka kroku z důvodu stability omezena není.
Efektivnost obou metod lze přibližně porovnat podle počtu úspěšně provedených kroků pk nebo lépe, podle počtu pf vyhodnocení pravé strany. V následující tabulce jsou uvedeny hodnoty pk/pf pro několik délek i intervalu integrace.
1	io-2	ío-1	10°	101	102
BS32	10/32	40/128	399/1211	3 982/11960	39 799/119 411
TR	10/15	U/21	16/24	67/79	86/108
Pro větší £ je problém (4.71) tuhý a to je důvod, proč k jeho efektivnímu řešení nelze použít explicitní metodu BS32 s omezeným intervalem absolutní stability.
Délka kroku BS32 metody je určena „přísnou" podmínkou stability. Proto při zvětšení intervalu integrace z i = 10 na l = 100 počet vyhodnocení pravé strany vzroste přibližně 10-krát. Jinak je tomu u A-stabilní TR metody. Tam je délka kroku řízena velikostí lokální chyby a ta se s rostoucím x zmenšuje. O tom svědčí to, že pro i = 10 se provedlo 67 kroků a pro £ = 100 celkem 86 kroků, tj. na interval (10,100) připadlo jen 19 dalších kroků. Na intervalu (10,100) lze tedy pozorovat paradoxní situaci: přesnější metoda BS32 vyžaduje
103
výrazně menší časový krok než méně přesná TR metoda. Délku kroku zde totiž omezuje stabilita metody a ne její přesnost.
Nenechte se zmást tím, že nestabilita se projevuje prostřednictvím lokální chyby, tj. pomocí nástroje pro měření přesnosti metody: jestliže délka kroku nesplňuje podmínku stability, dojde v několika málo krocích k takovému nárůstu lokální chyby, že na to zareaguje mechanismus automatického řízení délky kroku a krok patřičně zkrátí. □
Právě uvedený příklad vystihuje charakteristikou vlastnost tuhých soustav ODR. Tuto vlastnost sformulujeme jako
Tvrzení 3. Tuhost soustavy ODR se projevuje tím, že při jejím numerickém řešení délku kroku omezuje spíše stabilita metody než její přesnost.
Až dosud jsme naše úvahy o stabilitě numerické metody zakládali na analýze jejího chování při řešení lineární soustavy (4.68) s konstantní maticí J. Je jistě rozumné požadovat, aby metoda fungovala dobře na speciální třídě rovnic (4.68), je však třeba mít na paměti, že tak dostáváme jen hrubou informaci o možném chování metody při řešení problémů složitějších. Ukažme si nyní, proč se můžeme domnívat, že taková informace má vůbec nějakou vypovídací hodnotu.
Pro obecný problém y' = i(x,y) aproximujeme funkci f v blízkosti bodu (xi,y(xi)) lineární funkcí v proměnné y,
f(x,y) » f(xi,y{xi))+ô^(xi,y(xi))(x-
xi)+j-(xúy(x'))(y-y(xi)) ■
a doufáme, že chování numerického řešení původního problému lze v blízkosti bodu [xí, y(xi)) objasnit z chování numerického řešení aproximujícího problému
u' = Jiu + gi(x)
s konstantní Jacobiovou maticí     J1 = fy{xi,y(xi))
a s funkcí      g(x) - f{xuy(zť)) + íx(xhy{xi)){x - xt) - J'y (a*).
Aproximující problém je téhož typu jako problém (4.68), takže stabilitu numerické metody můžeme posoudit známým způsobem pomocí délky kroku h a vlastních čísel {Xj} konstantní Jacobiovy matice J'. Jestliže to uděláme, tak vlastně tiše předpokládáme, že stabilita metody aplikované na aproximující problém je stejná jako stabilita metody aplikované na problém obecný. To není vždy pravda, přesto však tento přístup může být pro praxi přínosem, nesmíme ale zapomenout, že případná podmínka stability omezující délku kroku má lokální charakter, tj. platí jen v okolí zvoleného bodu (xi,y(xi)).
Právě uvedený postup odpovídá klasickému postupu, který se používá v teorii diferenciálních rovnic, kdy se stabilita řešení nelineárního problému zkoumá prostřednictvím stability aproximujícího lineárního problému. V našem případě je třeba analyzovat navíc ještě stabilitu numerické metody aplikované na aproximující problém. Z teorie obyčejných diferenciálních rovnic je známo, že analýza lineární stability je užitečná, je však třeba mít na zřeteli její omezenou působnost. Totéž platí tím spíše i pro analýzu lineární stability numerických metod. Dokonalejší nástroje pro posuzování stability numerických metod při řešení nelineárních problémů vycházejí z teorie nelineární stability ODR, tím se však my zde zabývat nebudeme, základní informace k tomuto tématu lze najít např. v [17].
1
Příklad 4.3. Uvažujme počáteční problém
xe(0,2/í),      2/(0) = 5.
y' = y2-y\
(4.72)
Diferenciální rovnice má dvě konstantní řešení y = 0 a y = 1. Zvolíme-li počáteční podmínku y(Q) < 0, pak y(x) -t 0 pro x -> oo, zatímco pro 2/(0) > 0 dostaneme y(x) -> 1 pro x —> oo. Zvolíme-li ô > 0 velmi malé, pak pravá strana y2 — ys diferenciální rovnice nabývá malých kladných hodnot, tj. funkce y(x) 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(x) < 10~2 ještě pro x = 9 900, pak y(x) začíná prudce růst a pro x > 10 020 je už y(x) prakticky rovno 1. Roli Lipschitzovy konstanty L v nelineární úloze může do jisté míry reprezentovat velikost derivace pravé strany diferenciální rovnice, tj. výraz \J\ := \2y - 3y2\. Pro malé y « 0 je |J| « 0, pro y blízké 1 je však \J\ « 1. Na intervalu (0,9900) je výraz 9 900| J| (charakterizující podle tvrzení 1 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í (což se někdy nepřesně označuje jako snížená hladkost) než narůstající tuhost. Poté, co řešení nabude hodnotu rovnou přibližně 1, je délka kroku metody řízena stabilitou.
1.0001
0.9999
Obr. 4.2. Příklad 4.3 řešený DP54 metodou
Úlohu (4.72) jsme řešili explicitní metodou DP54 (program ODE45 v MATLABu) a implicitní TR metodou (program ODE23T v MATLABu). Oba programy jsme použili se stejnou přesností (er = 10-4, sA = 10~7). DP54 je explicitní Rungova-Kuttova metoda řádu 5 s intervalem absolutní stability (—3,30; 0). Délka kroku proto musí splňovat podmínku stability h\J\ < 3,3, což pro x > 10 020 znamená volit h = 3,3. TR metoda
104
105
(nebo-li lichoběžníková 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 e (9800,11200). V následující tabulce uvádíme pro intervaly (OJ), kde £ je postupně rovno 9900, 10020 a 20000, údaj pk/pf, kde pk je počet (úspěšně) provedených kroků a pf je celkový počet vyhodnocení pravé stany:
e	9900	10020	20 000
DP54	17/151	36/331	3041/20 245
TR	85/170	184/385	192/399
Vidíme, že v intervalu (10 020,20000), kde přesné řešení je prakticky rovno 1, je TR-me-toda velmi efektivní, na zdolání intervalu (10020,20 000) potřebovala jen 8 kroků. Zato metoda DP54 na tomto intervalu provedla 3 005 kroků délky h = 3,3, takže její použití rozhodně vhodné není. □
x 10
1.0001
0.9999
0.98        1 1.02      1.04      1.06      1.08       1.1 1.12
„4
x 10
Obr. 4.2. Příklad 4.3 řešený TR metodou
Výsledky příkladu 2 umožňují vyslovit praktické kritérium, podle něhož lze vcelku spolehlivě zjistit, zda uvažovaný problém je tuhý.
Tvrzení 4. Jestliže numerická metoda s omezenou oblastí absolutní stability, aplikovaná na počáteční problém, je nucena v jistém intervalu integrace používat krok, jehož délka je nepřiměřeně malá vzhledem k hladkosti přesného řešení v tomto intervalu, pak to znamená, že problém je. v tomto intervalu tuhý.
106
Velmi zjednodušený postup řešení počáteční úlohy, o níž nevíme, zda je či není tuhá, je tedy tento: zkusíme úlohu vyřešit přesnou explicitní metodou (v MATLABu třeba programem ODE113 založeným na Adamsových formulích řádů 1 až 12) a pokud řešení bude trvat příliš dlouho v důsledku volby extrémně krátkých kroku, výpočet přerušíme a zkusíme program pro řešení tuhých úloh (v MATLABu třeba program ODE15S založený na metodách zpětného derivování řádů 1 až 5). Pokud pravá strana diferenciální rovnice není příliš hladká, je vhodnější použít metody nižších řádů, třeba explicitní metodu BS32 (v MATLABu program ODE23) a v případě neúspěchu A-stabilní implicitní metodu, třeba TR-BDF2 (v MATLABu program ODE23TB).
Pokud dopředu víme, že problém je resp. není tuhý, neexperimentujeme zbytečně a přímo použijeme program určený pro řešení příslušného typu úloh.
Chceme-li efektivně využít všechny možnosti, které nám dostupné programy nabízejí, je třeba se seznámit s jejich vlastnostmi, vybrat správný program a vhodně nastavit jeho řídicí parametry (v kolekci programů v MATLABu jsou to například parametry ovlivňující řízení délky kroku, výstup výsledků, u implicitních metod je třeba určit způsob výpočtu Jacobiovy matice pravé strany, a jsou ještě další parametry, které zde neuvádíme a o nichž se lze poučit z programové dokumentace).
4.5.2. Řešení soustav nelineárních rovnic
Numerické metody pro řešení tuhých soustav jsou vždy implicitní. Jako příklad si uveďme TR metodu
yť+i = y; + ^h[f(xi+1,yi+l) +{(xí,yr)].
Rovnice pro yj+i má tvar y = hyí(y) +tp .
(4.73)
Pro větší přehlednost jsme v tomto zápisu potlačili závislost f (x, y) na x, neboť x = xi+i se při výpočtu y nemění. 7 je charakteristická konstanta metody a ip je konstantní vektor reprezentující dříve spočtená data. Konkrétně pro TR metodu
7 :
</> = yi + ™/if(xi,yi) ■
Rovnice (4.73) je vhodným modelem pro každou implicitní LMM, zejména pro BDF metody, ale také např. pro semi-implicitní Rungovy-Kuttovy metody TRX2 a TR-BDF2.
Nejdříve si objasníme, proč rovnici (4.73) nemůžeme řešit metodou prosté iterace v případě, když problém je tuhý. Předpokládejme, že y* je řešení a yM je už spočtená aproximace v s—té iteraci. Pak další aproximace je
yť*1) = hjf (y(s)) + i/>. Odečteme-li rovnici, kterou splňuje y*, dostaneme rovnici pro chybu:
y<s+1>-y' = /i7[f(y(s))-f(y*)].
107
Speciálně pro diferenciální rovnici (4.68) proto platí
y(s+1)-y* = ft7J(y(s)-y*)-
Jestliže jako počáteční aproximaci y'°) zvolíme vlastní vektor v matice J příslušný vlastnímu číslu A, pak chyba yw -y* = (hy\)'v. Konvergence tedy nastane, jen když \hyX\ < 1 pro všechna vlastní čísla matice J. Takové omezení délky kroku je však pro tuhý systém nepřijatelné.
Rovnici (4.73) budeme řešit iteračně, jinak to obecně nejde, víme však již, že nepůjde o iteraci prostou. Počáteční aproximace y(°> se obvykle určí dostatečně přesnou interpolací z hodnot yj,yi-i,____Funkci f linearizujeme okolo stávající aproximace y'5',
f(y)«f(yW) + J(y-yW), a další aproximaci y(s+1' dostaneme řešením rovnice
y(»+i) = ^ + A7f(yW) + Ä7j(y(h-i)_yM). (4.74)
Jestliže zvolíme J = fy(y(s)), dostáváme známou Newtonovu metodu. Ta vzhledem k dobré počáteční aproximaci konverguje velmi rychle. Problém je ale v tom, že takový postup je výpočetně příliš nákladný. V každé iteraci se totiž musí sestavit Jacobiova matice a řešit nová soustava lineárních rovnic. Proto všechny moderní programy pracují tak, že matici J drží co nejdéle konstantní. Aby iterace y(s> konvergovaly, měla by matice J být dobrou aproximací Jacobiovy matice v řešení y*, tj. J fa fy(y*)- Prakticky se postupuje tak, že na začátku řešení spoustavy ODR, tj. pro x = b, se za J vezme Jacobiova matice funkce f v bodě [a,rj\, a pokud se při řešení rovnice (4.73) zjistí, že iterační metoda (4.74) konverguje příliš pomalu, Jacobiova matice se přepočítá, tj. za J se vezme Jacobiova matice funkce f v bodě [:ri+i,yw]. Takto upravená Newtonova metoda se nazývá zjednodušená Newtonova metoda nebo také metoda tětiv.
Výpočet je dobré organizovat tak, že počítáme A = y(s+1> - yW řešením soustavy rovnic
a pak položíme y(s+1) = yW+A. Pro malé h je matice soustavy dobře podmíněná. Navíc, protože pravá strana „je malá" (jde to reziduum rovnice (4.73) dělené číslem hy), jsou malé také zaokrouhlovací chyby vznikající při řešení soustavy rovnic. Soustavu rovnic (4.75) můžeme řešit buďto pomocí vhodné iterační metody nebo metodou přímou. Všimněme si blíže druhého případu (použitého např. v MATLABu). Matice
soustavy (4.75) obsahuje tři členy, které se mohou měnit: h při změně délky kroku, y 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 rovnic s trojúhelníkovou maticí soustavy, tj. L<5 = b a pak UA = <5, kde b je pravá strana rovnice (4.75). Program, který má pracovat efektivně, musí mít promyšlenou strategii, podle níž rozhodne, kdy změní J, což je výpočetně nejnáročnější, a kdy jen h nebo 7. Do detailů zde zacházet nebudeme, zájemce odkazujeme na skvělou monografii [34], Neocenitelným zdrojem poučení je také studium kódů programů pro řešení ODR v MATLABu, viz [19].
Výpočet Jacobiovy matice. Časově velmi náročnou součástí algoritmu řešení rovnice (4.73) je výpočet Jacobiovy matice J. Je žádoucí, aby uživatel poskytl maximum informací pro její výpočet. Rada programů umožňuje zadávat přímo předpis pro výpočet prvků Jacobiovy matice. To ale většinou uživatel neumí a proto Jacobiovu matici program počítá numericky. Postupuje se obvykle tak, že prvky dí(x,y)/dye jejího ^-tého sloupce se spočtou přibližně užitím nejjednoduššího vzorce numerického derivování
t{x,y + Set) - t(x,y) 5
kde S > 0 je malé kladné číslo a e^ = (0,..., 0,1,0,..., 0)r je sloupcový vektor, který obsahuje jedničku v t-té pozici a ostatní složky má nulové. Numerický výpočet lze podstatně urychlit, pokud uživatel popíše strukturu nenulových prvků Jacobiovy matice (tj. pokud programu dodá indexy (j, k) všech pozic, v nichž prvky Jacobiovy matice mohou nabývat nenulových hodnot).
Poměrně častý je případ, kdy Jacobiova matice je pásová. Jestliže j-tá pravá strana fj(x, y) neobsahuje proměnné yk pro \k —j\ > s, pak Jacobiova matice může mít nenulové prvky jedině v pásu kolem hlavní diagonály. Číslo s + 1 se nazývá šíře polopásu. V tomto případě lze numerický výpočet Jacobiovy matice zorganizovat tak, aby se funkce f počítala jen 2(s + l)-krát, viz [34].
Přesnější výpočet derivace. Rada programů pro řešení tuhých problémů považuje za neznámou zj+1 = /if(a:i+1,yí+1) místo yj+1. Jestliže v rovnici (4.73) položíme z = ftf(y), dostaneme y = 7Z + i/> a odtud
z = Af(y) = M(7z + V')-Linearizace okolo yz + ip vede na iterační předpis
z(s+1> = /if(7ZM + j/>) + /i7J(z(s+1)-z(s)), a odtud dostaneme pro ô = z's+1> — z's) rovnici
(éI-,)'=^+*-^"t,)'
jejíž matice soustavy je stejná jako matice soustavy (4.75). Zřejmě z's+1) = z^+ó. Pokud zvolíme počáteční aproximaci z'°) = (y'0' — ip)jy, pak y<s+1) = 7z's+1' +t/> lze vyjádřit ve
108
109
tvaru y(s+1> = y's' + 7<5, jak lze snadno ověřit. Výpočet tedy můžeme organizovat takto: určíme y^, položíme z(°> = (yW —
z rovnice    {~l -3)6 = -f(yW) - -^-z^   vypočteme ô , (4.76) \"7       /       7 hj
určíme  z<i+l> = z« + £,   y(s+1) = y(s) + jô ■
Objasněme si přínos výpočtu podle algoritmu (4.76). Nechť z* = hf (y*). Jestliže spočteme yW podle (4.75) a položíme z'*' = /if(yM), pak pro chybu z* - z's' dostaneme
Äf(z,)-Äf(zW) = AJ(*)(y*-y
kde J(+) je Jacobiova matice funkce f, jejíž prvky jsou vyhodnoceny v (různých) bodech úsečky spojující body y* a y(*\ Pokud však z'3' počítáme podle (4.76), pak
Z--ZW = I(y._^)_I(yť.)_^) = Í(y*-yW).
7 7 7
Pro tuhý systém lze očekávat, že ||hJ(*)|| » 1, a v tom případě je přibližný výpočet derivace podle (4.76) výrazně přesnější.
Statistika o činnosti 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 poblému počínaly. Tak třeba 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í pravá strany f
pj   počet výpočtů Jacobiovy matice J
pr   počet rozkladů G = LU
ps   počet řešených soustav LU<5 = b Pro ilustraci uvádíme
Příklad 4.4. Robertsonův problém
j/Í = -0,04í/1 + 104j/22/3,
y2 = 0,04y1-104!/2y3-3-1072/22,
y3 = 3-107?y22,
ž/i(0) = l,
y2(0) = o, Ifa(O) = o
popisuje koncentrace tří příměsí v chemické reakci, tj. 0 < y\.y2,yz < 1, nezávisle proměnná í je čas, blíže viz [34], Jacobiova matice této soustavy je
/-0,04 104y3 104y2\
J =    0,04   -104 y3 - 6 ■ 107y2  -10*y2 ■ \   0 6-107y2 0 /
V čase t = 0, tj. pro yx = 1, y2 = y3 = 0, má Jacobiova matice vlastní čísla {—0,04; 0; 0}. Z fyzikálních úvah plyne, že yu y2 —► 0 a y3 -> 1 pro t -+ co. Vlastní čísla Jacobiovy
matice pro yi = y2 = 0, y3 = 1, jsou {-10000,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 í > 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 v MATLABu třemi programy určenými pro tuhé problémy: progra-mem ODE23T (metoda TR), programem ODE23TB (metoda TR-BDF2) a programem ODE15S (metody BDF1-5). Jacobiova matice se počítala přesně, délku kroku jsme řídili pomocí tolerancí eR = 10~3, sA = 10~6, činnost programu ODE15S 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
ODE23TB	140	13	630	10	93	728
ODE15S	245	15	504	11	67	458
Pro úplnost uvádíme, že MATLAB nabízí pro řešení tuhých problémů ještě program ODE23S (založený na modifikované Rosenbrockově metodě). Kromě toho, program ODE15S lze použít ve variantě založené na jisté optimalizaci BDF metod (jde o tzv. metody numerického derivování), podrobnosti viz [19], [36]. □
110
111