Miroslav Bartušek, Univerzita J. E. Purkynš - 1 - 1. CAUCI1YHO ÚLOHA PRO OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE. 1.1. Formulace úlohy. Budeme se zabývat řešením počáteční úlohy pro diferenciální rovnici U> y' = f(x.y). y(a) • y0 . Budeme předpokládat, že 1° funkce f(x,y) je definována a spojitá v páse G : aáxíb - es 0 taková, ie pro každé x« [a,b] a libovolnou dvojici y,z platí lf(x,y) - f(x,z)| á L ly - z I . Za těchto podmínek existuje pravé jedna funkce, která je spojitá na intervalu [a,b], má tam spojitou derivaci a vyhovuje diferenciální rovnici (1). Podmínka 2* se nazývá Lipschitzova podmínka. Uá-li funkce f spojitou a ohraničenou parciální derivaci fy v G : lfy| í L •< oe , pak podmínka 2 je splněna, nebol z věty o střední hodnotě plyne If (x,y) - f (x,z)l = lfy(x,ij) I . ly - z I é L ly - z I . Obecněji můžeme rovnici (1) pokládat za vektorový zápis systému N-rovnic prvního řádu, kde y, yo, f jsou vektory s N složkami. Mnoho z toho, co odvodíme, bude platit nejen pro jednu rovnici, nýbrž i pro systém. l'ři odvozeni metod pro numerické řešení diferenciálních rovnic bude důležité zabývat se těmito otázkami: 1) Jak velké chyby se dopouštíme v každém kroku výpočtu (a£ jde o chybu metody nebo o chybu zaokrouhlovací), jak tato chyba ovlivni výsledky v následujících krocích. Tento jev se studuje obvykle pod hlavičkou stability, přičemž stabilní metodou se rozumí taková metoda, ve které chyby, kterých se dopustíme v jednom kroku, nemají tendenci se zvětšovat v dalších krocích. 2) S problémem chyb a jejich šířením souvisí problém, jak stanovit chybu v dané fázi výpočtu jako funkci vypočítaných výsledků. 3) Otázka začátku řešení. Rovnice (1) obsahuje počáteční podmínku v bodě x = %o. Ale mnoho metod potřebuje pro výpočet hodnoty funkce y v daném bodě znalost hodnot funkce y v několika jiných bodech. Často je tedy potřeba pomocných prostředků pro začátek řešení. 4) Rychlost metody. Doba pro výpočet velkých soustav (N 3 10) může být značná i na velice rychlých počítačích. Proto při volbě konkrétní - 2 - metody musíme přihlížet i k rychlosti výpočtu. Jev stability se nesmí zaměňovat s vnitřní nestabilitou, která je vlastnosti diferenciálních rovnic samotných. Např. systém který má obecné řešení yl - + a2e" , y,, « a1ex - a2e~* je vnitřně nestabilní vzhledem k počátečním podmínkám y^O) = -J^'0' " v tomto případě i O, «a * 1. Ale libovolné malá porucha počátečních podmínek se projeví tím, že se stane nenulovým a to se projeví i v řešeni, kde člen e* převládne. ITotože chyba způsobená libovolnou numerickou metodou má za následek, že numerické řešení v každém bodě x± splňuje diferenciální rovnici s porušenými počátečními podmínkami, není možné získat přesné numerické řešení této soustavy s uvedenými počátečními podmínkami . 1.2. Metody numerické integrace. V následujících odstavcích si všimneme numerických metod řešení diferenciální rovnice y' = f(x,y). pomocí kterých hledáme řešení jen na zvolené množine hodnot nezávisle proměnné x z intervalu [a.bj" : a = převedeme na integrální integrací v mezích od x n x_ do (2) y(xn+1) = y(xn) ♦ n =0,1.....N - 1. Wl ) f(t,y(t))dt, y<*J Většina uvedených numerických metod se liší různým přístupem k výpočtu Integrálu v (2 ). Nejdříve zavedeme některé označení. Nechí Y(x) je přesné řešení a y(x) přibližné řešení naší úlohy. i>ak položme: i íl- dx ^i y,r fix. Interval h mezi uzly zůstává obvykle konstantní (v tom případě říkáme, ze uzly jsou ekvidístantní) a nazývá se krokem numerické integrace. Lze jej však změnit, jestliže se z rozboru chyb ukáže, že je to žádoucí. Je nutno si také dobře ujasnit otázku chyb. U interpolace se z přesných funkčních hodnot počítá přibližná hodnota v jiném bodé. U diferenciálních rovnic se z přibližných funkčních hodnot počítá další přibližná hodnota. Dochází zde teny ke kumulaci chyb. Nejjednodušší a Klasickou metodou integrace rovnice Y' = f(x,Y), Y(x ) = Y je Eulerova metoda: - 3 - <3> yn+i - yn + h ■ f(JW- *0 ■ Yo • Tuto metodu, na které si osvětlíme některé základní pojmy, obdržíme z integrální rovnice (2) přibližným nahrazením integrandu f(t,Y) polynomem nultého stupně na intervalu [x^ x„+1] . tj- hodnotou f(xn,yn). Definice. Funkce f(x) se nazývá řádu g(x)>0, jestliže existuje konstanta M>0 taková, že na průniku definičních oborů funkci g,f platí lf (x)I í H • g(x) . Označení f(x) ■ 0(g(x)). Při numerické realizaci kulérový metody vzhledem k zaokrouhlováni nesplníme rovnici (3), nýbrž rovnici W " *n + hf<*n-?n> + rn • kde rn je zaokrouhlovací chyba. Definice. Veličině Yn - yn se říká chyba metody, veličině Yn - yn akumulovaná (celková) chyba. Jestliže platí pro libovolné x £ [a,b] n.h ■ x - x. , ;im ♦ Yn - yn = 0 h-» o pak se příslušná metoda nazývá konvergentní. Vztah Yn dává rychlost konvergence. yn=0(hP) nám Věta 1. Je-li y řešení rekurentního vztahu (3) a Yn přesné řešeni rovnice (1), pak platí odhad iy LU -a) n - y„l * w(h) w(d) = sup Iy'(x.) - y'(xa)l Ix.-xJád 1 2 "1_Ä2 je modul spojitosti první derivace, L>0 Je Lipschitzova konstanta. Poznámka. Ze spojitosti funkce y'(x) plyne, že lim w(h) = 0 h-»o a tedy Eulerova metoda je konvergentní. Důkaz. Podle věty o střední hodnotě platí (n*l) Yn+1 ' Yn + h Y'(x„ e. h>, Oznaííme-11 en = Yn - y„, pak odečteme-li tuto rovnici od (3) obdržíme: en+l " en + h (Y'Un+ Ô'h) " '^n" •n+l = en + V " f ««n'7n> > ' h"1 Yn " Y' (*n + 0h »" Odtud {za využiti Lipschitzovy podmínky) - 4 - - 5 - Ip , ! á )• ! ♦ h.l. !e ! + h.w(h) l*n + 1 I ŕ (1 + hl.) Ien I ♦ h w (h ) . Tvrzení vrty dokážeme indukci, ľro n - (I věta piati. neboC z pnčáteč-nich podmínek pjyne \eQ \ = lYo - yol = 0. Obecný krok: t poslední nerovnosti a indukčního předpokladu plyne L<*n-> , lentllř(l + hl.).w(h).--—------ ♦ h«(h) * ,. , f hl. 1 (xn~a> 1 1 h . ] s w(h) .je e .j- - p - h + hj = -- w(h) .*- -1 čímž Je věta dokázána. Kulérová metoda je ne jjednodušši metodou jedné třídy numerických metod, které mi nazývají Adamsovy metody. Jejich zakladeni je následující úvaha. Předpokládejme, že jsme nějakým způsobem vypočítali přibližné hodnoty yo.....y() řešeni y rovnice (1), které odpovídají ekvidistant- n í ni uzlům k ...... xn. Přibližnou hodnotu ľI)tj iíe vypočítat ze vzorce U) ii.l kde funkci f nahradíme interpolačním polynomem, který v bodech x = . j - n -p. n - p * I.....n nabývá hodnot f. - f(Xj,y^). ('oznámeno jme. že tento lnterpolar.nl polynom se jen přibližně shoduje s interpolačním po-lyiionieiu funkce l"(x,V(x)). protože obecně y^ č \y Použijme pro aproximaci funkce f NewtonAv interpolační polynom vzad ru.y,x,) . rn ♦ . ... . H-1)^ r„.p . x = xn . t 1. . uosadlwe toto vyjádřeni do (4) a provedeme v integrálu substituci x = x„ . th: y . = y •■ » í U < t ac (t+,i-1)Ap rn „ ] ot ■ ' 11«1 Jn j() Ln n-1 P n-p J »nil - V' h[fn< ''lArn-3 " — * bp*P'n-pJ' 1 J h " í ( ' '"I"1) 1,1 ~ i'' í Mt ' 1)... (t «■ j -1 >dt . Te.ly kile rento vzorec se nazývá Adamsův extrapolační vzorec. Koeficienty bj nezávisí n>! na funkci I ani ru. kroku h a lze je jednou provždy spočítat ;í j sou talii- lovány. Hodnoty prvních přjti jsou následující: k 1 2 | 3 4 5 \ i/a 5/12 1 0/8 251/720 95/288 Vzorec se nazývá extrapolační, protože jej používáme pro výpočet hodnoty ve vnějším bodě interpolace. Aproxitnujerae-li funkci f Lagrangeovým interpolačním polynomem, lze stejným postupem ukázat, že AdamsAv extrapolační vzorec má tvar í5) *«1 " ín * h . £ bpj Vj PÍ j!(p i_ f t(t+n...(t+P) dt -j>! Jo t + j Nahradíme-li funkci f ve vzorci (4) interpolačním polynomem s uzly xn-p' xn~p*l'''",xn*xn+l ' iZř odvodit stejným způsobem Adamsův interpolační vzorec: *n+l ' ^ + h jinak bychom mohli zmenšit číslo p. Jestliže = 0, pak se metoda (7) nazývá explicitní a hod- notu yn+1 lze přímo vypočítat. Jestliže b_t ŕ 0, pak metodu (7) nazýváme implicitní a je to implicitní rovnice pro neznámou Yn+j> kterou je nutno řešit iterací. Koeficienty a^. b^ určíme tak, aby formule (7) byla přesná pro všechny.polynomy až do určitého stupně. Nebudeme však často určovat tyto koeficienty tak, abychom dosáhli přesnosti formule pro polynomy maximálního stupně. Necháme některé volné parametry a ty určíme tak, aby např, chyba byla co nejmenší nebo aby některé koeficienty byly nulové nebo aby vlastnosti šířeni chyb byly co najvýhodnejší. Pozadujeme-li, aby metoda (7) byla přesná pro polynomy do stupně r vřetni, pak koeficienty jsou vzájemně vázány r+i podmínkami, které obdržíme tak, ít- do vzorce (7) dosadíme funkce y(x) * x-1, j = 0,1,...,r, t j. podmínkami (8) / i *i i . - 7- iai i»o £ (-i)Ja j =o i * i=-l (-1) 1, 2,3, Jsou-li splněny první dvě rovnice (přesnost pro lineární polynomy) formule se nazývá konzistentní. Jsou-li splněny všechny rovnice, formule se nazývá řádu r. V tomto případě je vázáno 2p +3 koeficientu r+í podmiň kami . 1.3. Lokální chyba. Přesné řešení Y splňuje rovnici " / "i Ví + 11 í , bi Vi + Tn • (9) n + 1 i=-l kde Tn značí lokální chybu v kroku od xn do xn+j ■ 7 Taylorova rozvoje (v bodě xn) plyne platnost vztahu i2h2 i => Y - ihY' + Y" + n-i n n jji n (-l)rirhr (r) r! n 7Í í 'xn-i-8' r v(r+l) (s)ds Vi = yn -ihy;+ •• í n-i <7=*> J U xn n-i + (-Dr-1ir-1h1-1 „(r) (r - 1)! n - s)1"-1 Y(r+1,(s)ds . kde r je řád formule (7). Dosadíme-li tyto vztahy do (9), pak za použití podmínek (8) pro koeficienty a±, b± obdržíme vyjádření pro lokálni chybu. L5 (xn+1-8)ry(s)d8-í a. j (x^-L xn i'O *n .y(s)ds - rh| ft JD~\*n í-D'-'t"*11!.)*!] i=-l xn (Začátek Taylorovýcb. rozvoju pro dosazení vymizí, protože formule je přesná pro polynomy stupně r ). Tento vzorec lze přepsat na tvar: T = -n r Tn = FT J G<8) Ylr+1) , pak podle věty o střední hodnotě n-p' n+3 T_ = ,(r+l), ^ĽD 71G(s)ds. x„.p<^< sn+l * n-p Za tohoto předpokladu je chyba tvaru (10) t = C.hr+1y(r+1)(T> . což je vidět z následujícího vyjádření integrálu (substituce s * xnp + "roWŮB = hr+1 ["j {(p+l-t)ľ-rb í{p*í-t)r*i] át * xn-p p p + £ } /a.(p-i-t)r + rbi(p-i-t)r"1} dt] . i=o p-i Koeficient C v (10) můžeme určit také tak, že do formule (7) dosadíme funkci y(x) = (x-Xjj)1-*1, neboi pro tuto funkci známe hodnotu (r+1) derivace y(r+1)(|) = (r+1)! Pokud ovšem účinková funkce mění znaménko v intervalu [Xn.p'^+iJ • pak větu o střední hodnotě užít nelze a chybu nelze vyjádřit ve tvaru (10). Platí ovšem odhad |Tn,šJi!!{mi xf,G(8)i de kde xn_p^ é xn+i Je 1,0,1 ■ ve kterém funkce lY*r'| nabývá maximální hodnoty na intervalu [xn , xn+j ] • 1.4. Lineární diferenční rovnice s konstantními koeficienty. V tomto odstavci se stručně ztulnime o řešení lineárních diferenčních rovnic, tj; rovnic tvaru (11) y„ * Mn) ■ Zde Bj jsou známé konstanty, ak / 0, k se nazývá řád diferenční rov-nice. Uvažujme k rovnici (11) příslušnou homogenní rovnici (12) a, . y_„ +. Jestliže známe počáteční hodnoty yo,yd' spočítat rekurentně y pro libovolné *'yk-l pak mužeme podle (11) Lehce se ověří, že platí: l) lineární Kombinace řešení rovnice (12) je opět řešení této rov- 2) je-li dáno k řešení rovnice (12) y^ j>7^ g' 1,. . .) takových, že ,yi ,k (i = 0, (13) yi,i 'i,k yk,i 'k,k * 0 (tzv. fundamentální systém řešeni), pak libovolné řešení zn rovnice (12) lze vyjádřit ve tvaru zn " ¥n,d + °2yn,2 +"-+ Vn,k ' kde Cj jsou vhodné konstanty 3) obecné řešení rovnice (11) obdržíme jako součet partikulárního řešení «n rovnice (12) a obecného řešení rovnice (11). * + c.y .+...+ c, y Ľ n l/n,i kJn,k Cj libovolné konstanty, fundamentálni systém řešení rovnice (12) lze získat i v analytickém tvaru. Hledejme řešení ve tvaru yn = zn, kde z je neznámé číslo. Po dosazení do rovnice (12) obdržíme, že z musí být kořenem tzv, charakteristické rovnice (14) k-1 "k-1* Jsou-li kořeny této rovnice různé, pak fundamentální systém je tvaru: z" , .... z" . Je-li z = f>(cos ý + i sin V*' komplexní číslo, pak rovnice (14) má i kořen komplexně sdružený a existují dvě reálná řešení tvaru pncos n , pn sin n y> . Je-li z. kořen násobnosti r, pak vezmeme funkce z" , nzn. n2zn i ' i i r-1 n pro z^ reálné a p"cos n ip . p" sin n \p , npncosny>, .... n1^1^" cos n y>, nr"^>nstany> pro z^ = p (cos (f + i sin tyj). Ve všech případech obdržíme fundamentální systém řešení rovnice (12). 1.5. Stabilita, konvergence^ Nyní si všimneme, jak dobře řešení diferenční rovnice (7) aproximuje přesné řešení rovnice (1). Budeme se zabývat pouze konzistentními vlce-krokovými metodami. Prozkoumáme situaci nejprve pro konkrétní diferenciální rovnici (15) ~Ky . y(x0) = y0 řešením Y(x) = y . exp {-K(x-xo)J- . Pro tuto rovnici přejde rovnice (7) do (16) yn+1(l-ehKb^) = J_ (a. - hKbt )yn_1 s obecným řešením i=o (17) kde Cj jsou konstanty a r^ řešeni charakteristické rovnice (1+hKb_1)rp'fl = £ (at -hKbj^)^"1 . í=o (18) Přitom předpokládáme, že její kořeny jsou reálné různé. Jsou-li násobné, ovlivní to sice detaily, ale ne podstatu. Nejprve ukážeme, že jeden kořen, označíme jej r (jinak přehodíme pořadí členů v (17)), mé asymptotické vyjádření (16) rQ = d - Kh + 0(h2) (všimněte si, že tato vlastnost je vnitřní vlastnosti, diferenční rovnice 116) a vůbec nezávisí na diferenciální rovnici (15)! ). Z první rovnice konzistence (8) obdržíme, že rovnice (18) má pro h = 0 kořen r0 = 1. 9/9 1 ledy hledejme pro h = 0 odpovídající kořen ve tvaru r = 1+ Z. PHli . 1=1 Dosadíme toto vyjádření do (18): (1+ hKb .)(!+£ (3. h1)"*1 = f (a. - hK^ )(! + I /3.hl)P-' "a i = l 1 i=o 1 1 i=l 1 Provedením naznačených operací a porovnáním koeficientů u stejných mocnin h obdržíme podmínky pro koeficienty ■ - 10 - p i=o p p Or p K I/t + X i-8! *(\ [p + i-p.£ M i=o i=o J (je to první rovnice konzistence) p 0r p 1 = 0 Odtud, za využití druhé rovnice konzistence obdržíme (3^ = -K , což dokazuje vztah (19). -Kh P Protože 1 - Kh = e Nn + 0(h ) , platí (20) (1- Kh+0(h2))k = e"khK + 0(h2) = e~K + 0(h2). A z vyjádření (17) můžeme určit pomocí (p+1) počátečních c je dáno pomocí Konstanty c podmínek, dosazených do (17) (n = 0,1.....p). Napr. Cramerova pravidla 1 ... 1 / 1 ... 1 r» ... r_ / r r_ Počáteční podmínky jsou zatíženy jistou chybou, ale je rozumné předpoklá- dat, že se blíží k přesným hodnotám, tj. k y pro h ■* 0 , i -■ 0, ,p. Odtud za využití (19) vidíme, že lim c h*o+ což spolu 8 (20) značí, že výraz c0-r0 aproximuje přesné řešení diferenciální rovnice (15) pro malé hodnoty h. Ostatní členy v (17) se nazývají parazitní ře- äení a lze stejným způsobem ukázat, že lim c- «0, i ŕ 0. Tato řešení íHo+ 1 vznikají tím, že řád (p+1) diferenční rovnice je větší než řád i diferenciální rovnice. Úkolem je zvolit metodu (7) tak, aby parazitní řešení měla co možná nejmenší vliv na řešení diferenční rovnice (17). Definice. Necht počáteční podmínky yk = 3^'"'« k * 0,1,, při řešení rovnice (7) jsou takové, že (ai) ,p užité lim+yk(h) = , h-»o k = 0. Pak se metoda (7) nazývá konvergentní, jestliže pro libovolnou počáteční úlohu (1) má (7) vlastnost (2:j) lira y„ = Y(x), -n h = x - a h+o + " (n-»eo) pro všechna x € [a.t>J • Zde Y je přesné řešení rovnice (1). De f inice ■ linkncw, že nmtoda (7) je stabilní podle Dahlquista, jestliže všechny kořeny r charakteristického polynomu 183) J r'"1 - .£ a^P"1 11 jsou takové, že lij | é 1 a ty, pro které I rj I »1 jsou jednoduché. Veta 2. Metoda (7) je konvergentní pravé tehdy, když je konzistentní a stabilní podle Dahlquista. Důkaz. Dokážeme vetu pouze z jedné strany - konzistence a stabilita je nutnou podmínkou pro konvergenci. Necht tedy je metoda (7) konvergentní. Uvažujme úlohu y' ■ 0, y(0) = 0 s přesným řešením V s 0. Předpokládejme, že kořeny rg,...,r^ charakteristické rovnice (23) jsou reálné a jednoduché (pro komplexní a násobné kořeny je důkaz obdobný, ale složitější). Ukážeme, že pro tuto diferenciální rovnici jsou funkce i h.r? i=o k x p+l.p+2,, -P+1 Odtud (násobením výrazem rH"P) řešením rovnice (7). Ze vztahu (23) platí p i=o J i=o 1 J Sečtěme tyto výrazy přes j = 0,...,p a vynásobme h: P P P p . k+1 J=0 1=0 j = 0 •r Tedy y^ je skutečně řešení rovnice (7) pro y's 0. Navíc z vyjádřeni yfe, k = 0,...,p vidíme, že je splněna podmínka (21), nebof Koeficienty a. a tedy ani r. nezávisí na h. Aby platila konvergence lim y » 1 1 h*o* n = Y(x) 0, h.n * x, tj. P lim £ h r" = lim £ h r.5 «0 ll->0 + 1=0 tl-*0 + iso 1 musí být absolutní hodnota každého r, menší nebo rovna Jedné, tj. metoda musí být stabilní podle Dahlquista. Uvažujme nynJ rovnici y'= 0, y(0) = 1 s přesným řešením Y(x) £ 1. f'ro tuto rovnici má vicekroková metoda tvar y , = \ ajy„ \ • Necht po- i=o čáteční hodnoty jsou přesné, tj. jsou rovny jedné. Z definice konvergence plyne, že lim y » 1, což dosazeno do předchozí rovnice dává přímo h+o+ 1 orvní podmínku konzistence. Uvažujme nakonec rovnici y = 1, y(0) « 0 s přesným řešením Y(x) = x. Pak rovnice (7) přejde do tvaru (24) 'n+t t Vn-i * n í 1 = 0 Uvažujme posloupnost (25) yn = n.h.A. A V P , JI b..(l+ Za,.i)" , n = 0.1,2,... - 12 - lato posloupnost splňuje požadavky na počáteční podmínky (21) a vyhovuje rovnici (24): (n *l)h A = ^ a. (n-i)h A + h £ b. l=o 1 i = -l 1 P P P (n + 1) = J a. (n-i) + i + J i«. = n2a.+l = n + l. i=o 1 i=o 1 i=o 1 Protože však z předpokladu konvergence máme lim y = Y(x) * jc = n.h. h-»o+ n plyne odtud porovnáním s (25), že A = 1, což značí platnost druhé rovnice konzistence a důkaz je ukončen. Konvergenci jsme definovali pomoci limitního přechodu pro h -* O. V praxi se ale zajímáme o to, co se děje pro nenulové hodnoty h," které musíme používat při konkrétních výpočtech. Budeme se nyní zabývat akumulovanou chybou. Definujme £n = ?n - 7n- Platí P P (26) V* = IoVn-i + h ii-l^ň-1 + Tn P P (27) yn+i ■ X Vn-i 1 o + h J.jVn-i + «» 13 - (29) eB+1(i- "^VW = .|0 * f. i = o hbiK) £n_1 Je ji řešení je, jak lze snadno ověřit E (31) P Z i=o hK Y b. i=-l 1 kde Jsou kořeny charakteristické rovnice (32) (1 thKb^írP*1 = ^ (at- hKb. )rp_i . kůe Tn Je ohyba "etody a 1^ je lokální zaokrouhlovacl chyba, neboi yn+l nePOíítame přesně. Odhad pro 1^, jsou-li y„_i. y^_j zaokrouhleny na d-desetinných míst je 1«! «5 . 10-'d+1)(Z LJ ♦ h £ IbJ ). i=o i=-l Ve většině případů je chyba metody větší než zaokrouhlovacl chyba. Avšak v některých aplikacích (např. let řízených střel apod.) tomu může být naopak. Odečtením rovnic (26).(27) obdržíme: (28) n+1 ." .ž »iÉn-l + b í . bi*n-i + En • i=o i=-i kde En je lokální chyba v kroku od xn do xn+1 , která zahrnuje chybu metody i zaokrouhlení. Podle věty o střední hodnotě £n-i " Yn-i " K-i - f<*n-i'Vi> " "Vi^n-i' = ■ lYn-i - W "Vn-i•?»-!> - £n-ify(xn-i"?n-l) ■ kde ^n_1 leží mezi hodnotami y„_j. V(J_i . Tedy po dosazení do (28) Zde pro jednoduchost předpokládáme, že kořeny jsou jednoduché. Definice. Nechf metoda (7) je konzistentní 0. r^, i ■ 0,1.....p jsou kořeny charakteristické rovnice (32)j kde ro je kořen, pře který platí lim r * 1, Řekneme, že metoda je stabilní na intervalu h-K>+ ° který musí obsahovat nulu, když pro všechna hK na tomto intervalu je r. I Ä1, a když pro je r. i - 1.....P jednoduchý kořen. Poznámka. Existence kořene ro plyne z dokázaného vztahu (19) pro charakteristickou rovnici (18), resp. (32). Budo-li metoda (7) stabilní na intervalu fx./?]. P*ť parazitní řešení v (31) budou potlačena řešením d0r" a celkovou chybu lze přibližně odhadnout vztahem ř «d f" < -1- . <-n o o p hK V b. ( i = -l i Koeficient dn lze opět určit z počátečních hodnot. Označme pro jedno- ,y jako e. Pak z Crame- duchost chybu v počátečních podmínkách yg rova pravidla, aplikovaného na (31) pro n ,p plyne: - 14 d = o .p o - • -1 Foložíme-li přibližné r *I (viděli jsme v (18), že r je blízko jed- né), pak vztahem. (33) d0«* e + f a celková akumulovaná chyba je dána přibližným E A h E En*(« hK L=-l 1 hKÍ bi i«-l Všimněme si, že stabilita metody závisí na dané diferenciální rovnici, nebo£ za -K musíme vzit hodnotu charakterizující fy. Tím. že požadujeme, aby interval [P-.fil obsahoval nulu, máme zaručeno, že pro libevol-né K můžeme učinit řešení stabilní tím, že zvolíme dostatečně malé h. Navíc porovnáním takte definované stability na intervalu a stability podle Dahlquista zjistíme, že stabilita podle Dahlqulsta, a tedy i konvergence, je nutnou podmínkou pre platnost stability na intervalu. Poznamenejme ještě, že jestliže Il> 1 pro malé h, což odpovídá případu -K > 0, je řešení rostoucí exponenciálni funkce. Pak nelze očekávat, že udržíme chybu omezenou, a proto nepřekvapuje, že člen dor" ve vyjádření chyby (33) je neohraničený (vzhledem k n). To ale nevadí, neboř nám stačí udržet chybu malou vzhledem k přesnému řešeni, a te pravé vyjadřuje pojem stability na intervalu. Je však jasné, že nám chyba nesmí přerůst řešeni, tj. číslo d0 musí být podstatně menší než co. 1.6. Metody prediktor - korektor Uvažujme dvě numerické metody (R značí příslušné chyby metody). 'n+1 £ h(y" . + y') 2 Jn+1 3n' 1/12 hJY" <7i> yn+i ■ yn-2 + 3/211'yň+ yn-i> H = 3/4hJV- ly2) Obě jsou druhého řádu. První rovnice je obecně asi 9x přesnější. Ukazuje to na obecné pravidlo, že je vhodné užívat implicitních metod i přes obtíže, které při tom vznikájí. Obecně implicitní rovnice musíme ovšem řešit iteraci. Nějakým způsobem určíme první aproximaci 7^*1.' vyP0CÍtaBle f'xn+1. yn+í' a dosadíme do pravé strany (7) a provádíme iteraci, lze přepsat na tvar: yU+1) = r (a v i =o Přesné řešení splňuje Rovnici (7) + hbiyn- .) * hb.i(ynJ») - 15 yn+i ■ £0 + hb-iyň+i Odečtením těchto rovnic obdržíme yn+i-yn£4) ■ bb-i [yn+i-'3 !ki.i-Vvi' ^^t-útb- leží mezi body y„+1. y„+j • Jestliže platí |fy I ^ K v určitém okolí bodu \*n+i> yn+ll • ktere obsahujfi body Dw y„ií], J - O'1-2.....pflk Iv v'J+1)l á hb K Iv - v(j) I |yn+l yn+l I * ao-iK lyn+i yn+l1 . a odtud indukcí i 34) v v'J+1)l < íhb KlJ1"1 Iv - vvo' yn+l yn+l 1 - *nD-4*'' • ''n+1 yn+i ,<«>l Tedy iterace konverguje, je-li hb_4K-^l a rychlost iterace je dána vztahem (34). Důležitá věc je určení počáteční aproximace. Ze vztahu (34) vidíme, že člm přesnější bude počáteční aproximace, tím rychleji budou iterace konvergovat. Nejlepší pro určení poč. aproximace je použit některou explicitní formuli. Takto použitý systém, formuli se nazývá metoda prediktor-korektor (předpověá- oprava ). Příkladem soustavy prediktor - korektor 4.řádu je Milnova metoda. Irediktor: y<°> = yn_3 + \ h <2yn - y^ + 2yn_2) (v(0)) f(x, n + 1 T(o), yn+l' Korektor: y<> = ♦ £h [(yU>)' + 4y„ ♦y^j , j - 0,1.2. Chyby jsou j| h5 Y(v>(7|1) resp. Lh5Y(v) «?2>' (Korektor je Simpsonovo pravidlo). Je důležité, jak uvidíme z rozboru chyby, aby prediktor a korektor byl stejného řádu. Jako prediktory mohou být použity např. Adamsovy extrapolačnl vzorce aj. V praxi se nejčastěji používají metody prediktor-korektor 4.řádu. Prediktory nejsou tak důležité pro výpočet a při Jejich výběru se řídíme malým koeficientem u ohyby. Na Milnové metodě ukážeme, jak lze rozboru chyb užit ke zlepšeni metody. U jiných mptod to lze udělat úplně stejně. Prediktor lze zapsat takto: y^í - V3 + Th(2*n -*„-l+2Yn-2> - 6n-3 " - 16 - - 4/3h(2Én - + 2f n-2> * Yn+1 " £n-3 " - 4/3h(2€; - e;.! ♦ 2€ ;.2) -£f h5Y<*><7l) Zde £„ je chyba. f„ = Y„ - y„. Podobne korektor (vynecháme indexy) kde y je hodnota, kterou bychom získali, kdybychom rozřešili korek-n+1 tor přesné. Odečtením těchto rovnic získáme: yn+1-Vnľ^M^(5) ♦ £„.3 - Cn.4 . Předpokládárae-li nyni, že se chyba £n od jednoho kroku ke druhému mění jen pomalu (pak £^»0) a že Yť5' se příliš nemění mezi hodnotami Tl' ^2 ' pa't mu*eBie P8a* tedy , , h5Y(5),_, fyj 90 , _ (o). Odtud lze určit přibližně chyby metody. Korektor: T„ = h5Y(5) ty « -|g (y„+1 - y„°i> Prediktor: = |J h5Y(5>ty«* §§ • Jelikož podle našeho předpokladu se rozdíl mezi předpověděnou a opravenou chybou mění pomalu, mažeme položit: T(0)s» rnc (V y<0)). *n " 25 vyn yn ' To nám dovoluje zlepšit Milného metodu: Prediktor: y„^ = yn.3 + 4/3h (2y^ - y^ ♦ 2yn.2) Modifikátor: y<;> « yn?> ♦ §| (y„ - y' ■ f<"n+l- yníi> • Korektor: y^*1' - y„.1 ♦ */3h [' + * yn.4 J • Pro nejlepši určení korektoru je nejdůležitějěí hledisko stabilita. Vzhledem k tomu, že metody 4.řádu jsou najpoužívanejší, nastíníme způsob, jak odvodit korektor co nejlepších vlastností. Zvolíme tvar korektoru: 17 - (35) yn+l " aoyn +Vn-l + a2yn-2 + h K dosažení 4. řádu je třeba 5 koeficientů. My máme 6,.proto zbývající koeficient použijeme k tomu, aby korektor byl stabilní. Požadavek na přesnost formule vede na: (36) o a. s a 1/8 (9 -9aj) , 1 ~ "i a2 « -1/8 (1-Bj) b-l " 2T <9-«i> bo 'h <9 + 7-i) bj «= 1/24 (-8+ 17a4) Zde a4 Je volný parametr. Účinková funkce ve výraze pro lokální ohybu Je G(s) = (xn+1 - s)4 ♦ a^x^ - s)* ♦ a2(«n.2 - ■>* + + 4h [-b.i^+i - »>3 ♦ ^(vT^)3] • Zde nám jde o to. pro které a^ tato funkce nemění znaménko. Lze zjistit, že to je pro a. £ <-0,6; 1,0> . Pak chyba Je tvaru 5 (5) R ■ C.h Y' 'ty. Charakteristická rovnice (32) z definice stability na intervalu je tvaru ■'• (l + hKb.jtr3 = (a0- hKb0)r2 + (aj-hKbjJr + aa Hledáme, pro které hodnoty hK platí Pro h ■ 0 (tj. sta- bilita podle Dahlquista) lze zjistit, 0 že kořeny charakteristické rovnice mají průběh znázorněný na obr. 1. r. -1 -0,13 o m 1. Vidíme tedy, že stabilita může nastat Jen pro -0,6 á a4 ŕ 1,0. Nejlepší poměr |p^| je okolo hodnoty a^ « 0. Volme tedy dále .'■ 0; Na obr. 2 je nakreslen průběh kořenů charakteristické rovnice pro různá v 1 L8 -i 1 t 1 1 1--r* -1 ■-- -i. -1- * i i Je vidět, že pro hK S 0.69 je podmínka stability splněna. Koeficient b_4 = 3/8 a tedy pro konvergenci iteraci korektoru je třeba lhKl<£8/3 (viz (34)). Tedy hodnota |hKláO,69 není příliš omezující pro rychlost konvergence iteraci. Závěr: Ponecháme-li Milneův prediktor a použijeme-li tento korektor, obdržíme Hammingovu metodu: IhK I *0,69 Prediktor: Modifikátor: Korektor: y(0) /n+i yn+l „(o) yn+l 112 rlJll) - k <9yn- wf* HtW*^-yU]- Vícekrokové metody (7) maji některé nevýhody. Pro své použití vyžadují znalost (p+i) počátečních podmínek (diferenciální rovnice (1) nám poskytuje pouze jednu) a nelze u nich vzhledem k ekvidistantním uzlům měnit krok. Jak překonat tyto obtíže, uvidíme později. 1.7. Metody Runge-Kutta. Tyto metody jsou jednou z nejrozéířenějších metod v praxi. Nevyžadují dodatečných počátečních hodnot a lze Je snadno naprogramovat. Nevýhodou je obtížné určení chyby a malá rychlost výpočtu, neboí se musí v jednom kroku počítat mnoho funkčních hodnot f. Základem všech R - K metod je vyjádření rozdílů mezi hodnotami v bodech x n+11 ve tvaru - 19 (37) W " 'n " £,»1*1 kde W± jsou konstanty a i-i (38) h " *n+l - V Äl * °-Naším cílem bude určit konstanty (U i( 0t. tak, aby (37) dostatečně aproximovala řešení diferenciální rovnice (1). Konkrétně Je určíme tak, aby koeficienty * lir t Taylorově rozvoji pravé a levé strany (37) v bodě [xn.yn] souhlasily pro r « 1,2,... N. V ton případě výsledný vzorec nazýváme metoda Runge-Kutta řádu N. Nelze dosáhnout lepšího výsledku než N = m. Odvozeni koeficientů naznačíme pro Jednoduchost pro m * 3. ľro ostatní m je to obdobné, ale výpočty Jsou příliš složité. Pro jednoduchost nebudeme rozlisovat označení mezi přesným a přibližným řešením. Budeme tedy uvažovat metodu (39) 'n+i Wik! + W2k2 + W3k3 + R K * h flxn'y„» k2 = h flkl> kde R je chyba metody. Rozvoj levé strany (39) je yn+l t=l '(xn,yn) kde = f = f. h 2f .f + f .ľ xy yy V» * fy Zavedeme-li oznaiení D a ^ + f.j^ , pak (40) w - yn - [»■' + ŕDf + í «2h. yn * /321h.f(xn,yn)) - = h.f * h2D2f + j£ ľ|f * 3T D2f In * 0 • Dosadlme-li tyto vztahy a (40) do rovnic© (39) dostaneme na obou stranách rozvoje podle mocnin h. Porovnejme koeficienty u h,h2,h3: Wl +W2 + W3 " 1 (41) (42) W2D2f *W3D3f - §f 02f °2f D,f „ ■ , «2 ~T .* W3 "T" + W3/VyD2f 1 3T } (D3f>/332ry . • V + 6^32VD3Vln V . ; Tím máme ve skutečnosti ne 3, ale 6 rovnic pro určeni koeficientů, nebol Wj, cti( musí být nezávislé na funkci t, má-lií být metoda užitečná. RozepiSme poslední dvě rovnice (41) do parciálních derivaci a porovnejme koeficienty u stejných derivaci: *2 < Ajl + r332> W3 = £ W2Ä2 + UgOtg * 1/3 "zPti * '.'" 1/3 ' ; •4 «3 ' Až " 1/6 ^ ' : i ^Ami • ^zi = 1/6 Z posledních dvou rovnic plyne oíg = /321 a z prvních dvou pak otj = /^31 * fls2' Celkov6 tedy máme pro určeni koeficientů vztahy:. yy" • 'Vy' Vy = = 1 . funkci y(x) = (x-xn)* 6. Dokažte, že jestliže účinková funkce G(s) mění znaménko v integračním intervalu [a,b], existuje spojitá funkce y(s) tak, že G(s)yU ) ds * y(q) 1 G(s)ds pro každé <0 [a.b] 24 - 7. Uvažujte numerickou metodu yn+iB <1-»>Jn + ayn-i*Tí [(5-a)yn+i + <8+8a>y„+<5»-i>yn.i]. Ukažte, že je tato metoda 3.řádu a zjistěte, pro která a účinková funkce G(s) nesáni znaménko. Určete chybu. 8. Pro a = \ určete hodnoty ;hk, pro něž je metoda z předcházejícího přikladu stabilní. 9. Vysetřte stabilitu nejjednodušších Adamsových metod. la Určete, jaké hodnoty h je možno volit pro konkrétní výpočet, rešíme-li rovnici y' * sin (x,y), x e [o.looj , y(o) =1 , " metodou w ?n * h (5yňn + 8yn-yn-i'- 11. Určete, jak velký může být krok h, aby iterační prooes pro Adamsovu interpolační formuli 4.řádu konvergoval, aplikujeme-li Ji na rovnici: y' = -10y + 15, y(0) = 1. 12. Nalezněte přibližné řešeni rovnice y' + 2xy = 2x3, y(0) = 0 v bodech x » 0,1; 0,2; 0,3 pomocí nSkteré Runge-Kuttovy metody 4.řádu. - 25 - OKRAJOVÉ PROBLÉMY DIFERENCIÁLNÍCH ROVNIC 2. VARIAČNÍ METODY. 2.1. Okrajové problémy. Okrajové problémy obyčejných a parciálních diferenciálních rovnic, popř. systémů rovnic se řeši různými metodami. Nejpoužívanější jsou tyto* a) Variační metody - do této kategorie patří i metoda konečných prvků, která je v poslední době velice moderní. b) Metoda sítí - je to nejrozéířenějíl metoda, v poslední době však ustupuje modernějším metodám. c) Metody, které úlohu ve více dimenzích převádějí na posloupnost úloh jednodimenzionálních. 2.2, Hilbertův prostor. V této kapitole shrneme některé základní pojmy a věty z funkcionální analýzy, které budeme dále potřebovat. Důkazy a podrobnou teorii lze nalézt např. v knize (jY] • Říkáme, že M je lineál (lineární prostor), jestliže a) pro libovolné prvky u e M, v e M a libovolné reálné číslo a je definován soufet u + v € M a součin a.u € M s vlastnostmi u + v = v + u a (u + v) = au + av a(bu) = (ab)u u + (v + z) = (u + v) ♦ z (a + b)u = au + bv 1 . u = u b) existuje jediný prvek 0, který má vlastnost u + 0 = u, u€M c) ke každému prvku u e M existuje jediný prvek v « M s vlastností Prvky u,,... vislé, když platí *2»2 + u + v ■ 0, lineárního prostoru M se nazývají lineárně nezá- Vi A„u„ n n ... ■ A ■ o. V opačném případě se nazývají lineárně závislé. Prvky množiny {%}l M nazývají lineárně nezávislé, jsou-li lineárně nezávislé prvky libovolné její konečné podmnožiny. Nechf P = {u } jo nejvýše spočetná množina, P C M. Množina všech k . prvků x tvaru x = "£ A.u. , A. libovolná reálná čísla, k libovol- 1=1 1 1 1 né přirozené číslo, se nazývá lineární obal množiny P. ťííkárae, že na lincálu M je definován skalární součin, je-li každé