Obyčejné diferenciální rovnice (ODE)

Obyčejné diferenciální rovnice (ODE) Obyčejné diferenciální rovnice  N ­tého řádu převádíme na soustavy  N diferenciálních rovnic prvního řádu. V rovnici  f  x , y , y ', y '',, y N  =g x  se substituují  y '=z1 , y ''=z 2 , ... ,  y N −1=z N −1 a máme tedy soustavu rovnic  y' z1 ' z2 '  f  x , y , z1 , z 2 , , z N −1 , z 'N −1 

= = = = =

z1  z 2  z3  ,  g x 

poslední rovnici obvykle lze psát ve tvaru  z 'N −1= g  x , y , z1 , z 2 ,  , z N −1  a soustavu lze celkově zapsat jako y z0 ' z1 ' z2 '  zN−1 '

= = = = = =

z0 z1 z2 . z3  gx , z0 , z1 , z2 , ,zN−1 

K jednoznačnému řešení musí být známo  N podmínek. Podle toho, ve kterých bodech jsou zadány se řešení ODE dělí na: ●

●

 Počáteční problém – všech  N podmínek je zadáno v jednom bodu, sleduje se přímo řešení vycházející z tohoto bodu, příkladem může být rovnice harmonického oscilátoru a k ní zadané počáteční podmínky v čase  t=0 ,  d 2 y t  =y t  , y 0=0 , y '0=1 d t2  Okrajový problém – ne všechny podmínky jsou zadány v jednom bodě, často jsou podmínky zadány na dvou koncích intervalu, na kterém řešíme, například rovnice popisující pohyb hmotného bodu v gravitačním poli a podmínky  y  x 0  =0 , y  x n =0 a  v 0=v 0 popisující úlohu, kterou řeší střelec z děla, když hledá elevační úhel   takový, aby se z bodu  x 0 trefil střelou s počáteční rychlostí  v 0 do bodu  x n

v0 x0



xN

Runge-Kuttovy medoty pro řešení počátečního problému Budeme se zabývat odvozeními pouze pro 1 diferenciální rovnici 1. řádu (vzorce zobecněné na systém diferenciálních rovnic 1. řádu budou uvedeny na konci kapitol).

Eulerova metoda Řešíme tedy rovnici 

dy = f  x , y s počáteční podmínkou  y  x 0  =y 0 . Přibližné dx nalezneme pomocí rozvoje funkce  y  x do Taylorovy řady v

řešení v bodě  x n1 bodě  x n   y  x n1 =y  xn h≈y  x n h y ' x n =y  x n h f  x n , y n 

Nejnižší zanedbaný člen určuje odhad chyby metody, chyba v každém kroku je tedy úměrná  h2 , a protože počet kroků na daném intervalu je nepřímo úměrný  h , je celková chyba metody úměrná  h . Eulerova metoda je tedy 1. řádu přesnosti, tedy poměrně málo přesná, a proto se v praxi nepoužívá (k dobré přesnosti je třeba velkého množství kroků). Postup při použití Eulerovy metody je znázorněn na následujících obrázcích. 

lim y  x n h−y  x n  K Eulerově metodě lze dospět i ze vztahu  d y = h  0 = f  x , y . dx h Informace o konvergenci metody ve slidech, postačující podmínka konvergence je, aby funkce  f měla na oblasti, kde hledáme řešení omezenou parciální derivaci podle  y . Metody vyššího řádu, které by využívaly Taylorova rozvoje se v praxi nepoužívají, protože by  byly potřeba vyšší derivace  y a tedy parciální derivace  f podle  x a y . Taylorův rozvoj je ale důležitý i pro odvození dalších metod.

Princip Runge-Kuttových metod Jsou v praxi často používané, a protože pro výpočet  y n1 využívají pouze předchozí bod   x n , y n  a nevyužívají body   x k , y k  , kde  k n , říká se jim metody jednokrokové. Runge­Kuttovy metody jsou založené na postupném zpřesňování hodnot ležících mezi  x n a  x n1 . Obecný výpočetní vzorec metody má tvar  y n1 =y nhRK  xn , y n , h ,  kde   x n , y n , h= p1 k 1 h p 2 k 2 h p r k r h . Vzorec je v podstatě hodně podobný Eulerově metodě, pouze místo  f  x n , y n  používáme funkci   RK  x n , y n , h , která je defakto lineární kombinací funkčních hodnot funkce  f . Pro  k 1 h , , k r h totiž platí vztahy k 1 h k 2 h  k r h

= = = =

f  x n , yn f  x n 2 h , y n21 h k 1  .  f [ x n r h , y nh r1 k 1  r2 k 2 rr−1 k r−1 ]

Volbou  r =1 dostaneme přímo Eulerovu metodu. Pro  r≤4 se řád metody může rovnat  r . Pro konstrukci medoty 5. řádu je už třeba  r=6 . Konstrukce RK metody 2. řád u Vyjdeme opět jako u Eulerovy metody z Taylorova rozvoje funkce  y  x h2 h2 d f y  x = y  x h= y  x h y '  x  y ' '  x = y h f  x , y   x , y = n1 n n n n n n n   2 2 dx n n h2 ∂ f ∂ f ∂y = y nh f  x n , y n   x n , y n  x , y  = 2 ∂x ∂y ∂x n n h2 =y nh f  xn , y n   f xf f y   xn , yn  2





Tentokrát jsme použili o jeden člen rozvoje více, abychom byli schopni dopočítat konstanty  p1, p 2,  2, 21 používané v RK metodě . h2 V Taylorově rozvoji je tedy  y n1=y nh f x n , y n   f xf f y  xn , y n  2 y n1 =y nhRK xn , y n , h= v RK metodě  . =yn h  p1 f xn , yn p2 f xn 2 h , y n21 h f xn , y n   Abychom dostali v RK metodě rovněž členy s  h 2 , provedeme také Taylorův rozvoj funkce   RK  x n , y n , h= RK  x n , y n ,0h  RK ' x n , y n ,0 . Vyjádřením tohoto rozvoje dostaneme  p1  f  x n , y n  p2  f  x n , y n h p 2    2  f x 21 f f y   x n , y n  . Porovnáním s Taylorovým rozvojem funkce  y pro jednotlivé mocniny  h pak dostaneme soustavu rovnic pro neznámé koeficienty RK metody: h 2  y n h f   f x  f f y  2  y n h  p1   p 2  f h 2  p2   2  f x 21 f f y  h

:

h 2  f x : h 2  f y :

p1  p2 =1 1  p 2   2 = 2  1  p2  21= 2

Tato soustava má nekonečně mnoho řešení, v RK metodách se však používají jen dvě z nich: 1  ●    2  =1,  21 =1,  p1  = p 2  = , je analogická lichoběžníkové metodě integrace 2 h k 1  = f  x n , y n  , k 2 = f  x n h , y n h f  x n , y n  , y n1=y n   k 1  k 2   2  ●

1  p =0,  p 2 =1, je analogická obdélníkové metodě integrace 2 ,  1  h h k 2 = f x n  y n  f  x n , y n  , y n1=y n h k 2 2 ,  2     2 =21=





Runge­Kuttova metoda 4. ř ádu Nejčastěji se při výpočtech používá RK metoda 4. řádu, v jednom kroku se používá h h k 1  = f  xn , yn k 2  = f x n  y n  k 1  2 ,  2  k 3 

=



f x n

h h y n  k 2  2 ,  2 





k 4 

=

f  x n h , y n h k 3 



y n1=y n 

h k 2  k 2  2  k 3  k 4. 6  1 

Chyba v jednom kroku je řádu  h 5 , chyba metody v zadaném intervalu se zvětšuje lineárně s počtem kroků, který je nepřímo úměrný  h , a tedy celková chyba na intervalu je řádu  h 4 . Odhad chyby a automatická volba kroku Existují 2 metody pro odhad chyby, který se používá při automatické volbě kroku. ● ●

 Srovnání výsledků 2 kroků o délce  h s výsledkem jednoho kroku o délce  2 h  Srovnání výsledků 2 RK metod různého řádu

Různá  délka kroku Odhad chyby je  =y h −y 2 h . Postup při automatické volbě kroku je následující:   požadovaná maximální lokální chyba je  0   je­li  ∣∣≤∣0 ∣ , krok přijmeme a v následujícím kroku zvětšíme  h , h '=S h  

∣ 5

∣

0  , kde  S≃0.9 ,  h ' obvykle omezíme např.  h '≤4  h   

∣

4  je­li  ∣∣∣0 ∣ , krok nepřijmeme, ale opakujeme ho s  h '=S h

0  

∣

 při přijetí výsledku je možné ho zpřesnit použitím vztahu 16  1  y  x2 h= y h  x2 h− y  x2 h , tím získáme metodu 5. řádu 15  15  2 h přesnosti, tedy o řád vyšší

Srovnání výsledků dvou RK metod (tzv. Embedded RK metody) Někdy se jim také říká Runge­Kutta­Fehlbergovy. Pro metodu 5. řádu přesnosti např. máme vztah  y n1=y n c1  k 1c 2  k 2 c 3  k 3 c 4  k 4 c 5  k 5 c 6  k. Jednotlivá  c lze volit 6 i tak, že dostaneme metodu 4. řádu y 'n1 =y n a 1  k 1 a 2  k 2 a 3  k 3  a 4  k 4  a 5  k 5  a 6  k. 6 6

Pak je odhad chyby dán jako  =y n1 −y 'n1=∑ c i −a i  k i . i=1

Spojité RK metody se používají pokud potřebujeme poměrně přesně znát hodnoty  y i v bodech mezi  y n a  y n1 a nepostačuje nám obyčejná interpolace (více informací ve slidech a v Numerical recipies). Věty o lokálních a globálních chybách RK metod viz. slidy.

Vlastnosti RK metod ●  samostartující ●  robustní a odolné ●  dostupné v numerických knihovnách ●  mnoho vyčíslení funkce  f (někdy pomalé) ●  nehodí se pro stiff rovnice Pokud je funkce  f nespojitá, můžeme nespojitost ignorovat a spolehnout se na automatickou volbu kroku, nebo nespojitost nalézt, výpočet v tomto bodě ukončit a začít opět za ním. O konzistentnosti, regulárnosti a konvergenci obecně viz. slidy. Zaokrouhlovací chyby jsou větší při volbě menšího  h , a proto existuje pro volbu  h optimum. Pokud použijeme příliš malé  h může se i stát, že  y n  y n1 −y n =y n , ačkoliv to ve skutečnosti neplatí. Příklady v PASCALU ODE1.PAS , ODE2.PAS

Bulirsch-Stoerova metoda Je obdobná Rombergově integraci a hodí se pokud je funkce  f dostatečně hladká a zadaná rovnice nemá singulární bod. Postup je následující:   provede se výpočet pro několik  h , používá se sudá metoda (chyba je úměrná sudé mocnině  h ) 2   výsledek se interpoluje racionální lomenou funkcí proměnné  h a provede se extrapolace na  h=0   výpočet se provádí postupně s počtem kroků n=2,4 ,6 ,8 ,12 ,16 ,24 ,32 ,48,64 ,96 , více kroků než  96 se nepoužívá   pro interpolaci se používá maximálně 7 z vypočtených hodnot pro různá  n Pro jednotlivé výpočty se zadaným  h se používá zpravidla nějaká jednoduchá metoda, např. metoda středního bodu uvedená ve slidech. Opět je snaha o co nejméně nutných vyčíslení funkce  f .  S Bulirsch­Stoerovou metodou lze dosáhnout přesnosti 14. řádu, metoda se opět nehodí pro řešení stiff rovnic. Příklady v PASCALU BS41.PAS.

Vícekrokové metody Používají se v jednom kroku hodnoty funkce  y a  f ve více bodech (ne pouze

k

k

j=1

j=0

x n , y n ). Obecně pro ně platí vzorec  y n1 =∑ a j y n1− j h ∑ b j f i1− j . Pokud platí  b 0  =0 nazýváme metody obecně explicitní, v opačném případě implicitní. Vícekrokové metody nejsou samostartující. Adamsovy metody explicitní metoda 4. řádu Adams­Bashforthova y n1=y n 

h 55  f  x n , y n −59  f  x n−1 , y n−1 37  f  x n−2 , y n−2 −9  f  x n−3 , y n−3  24  implicitní metoda 4. řádu Adams­Moultonova

h 9  f  x n1 , y n1 19  f  x n , y n −5  f  x n−1 , y n−1  f  x n−2 , y n−2  24  pokud je  f lineární funkce  y dají se všechny  y n1 převést na jednu stranu rovnice a rovnice se dá jednoduše vyřešit. V opačném případě by bylo třeba řešit nelineární rovnici h g y=y n  9  f  x n1 , y19  f  x n , y n −5  f  x n−1 , y n−1  f  x n−2 , y n−2 −y=0 24  y n1=y n 

Metody prediktor­korektor Postup je následující:  n1 explicitní metodou (např. Adams­Bashforthovou)   prediktor – odhad  y y ' = f  x n1 , y n1    výpočet   n1  'n1 určíme  y n1   korektor – implicitní metodou s použitím  y   výpočet  y ' n1 = f  x n1 , y n1  Je složitější automatická volba kroku, ale v knihovnách implementována je. Výhodou metody je rychlost, nevýhodou je, že není samostartující, je citlivá na vlastnosti funkce f a nehodí se pro stiff rovnice. Stabilita, konvergence a konzistence viz. slidy.

Špatně podmíněné úlohy Přík lad: Řešení diferenciální rovnice  y ' ' = y S počátečními podmínkami  y 0=1 a  y ' 0=−1 je analytické řešení  y=e−x . S počátečními podmínkami  y 0=1  a  y ' 0=−1 je analytické řešení y=e−x  e x . Pro libovolně malé   existuje  x takové, že složka řešení   e x bude převažovat nad složkou  e−x . A protože se chybě při numerickém řešení prostě nevyhneme, objeví se v řešení parazitní složka, která poroste a po určité době převáží řešení. Takovým úlohám se říká špatně podmíněné. Příklad v PASCALU ODE3.PAS .

Stiff rovnice (rovnice se silným tlumením) Obsahují v sobě útlum s charakteristickým časem   mnohem menším, než jsou zbylé charakteristické časy úlohy. V čase větším než   je rychle se tlumící složka zanedbatelně malá, proto je u dosud popsaných metod nutné volit krok  h . (Nemusí se samozřejmě jednat jen o čas, ale ve fyzice se takové rovnice často vyskytují a jsou většinou v nich jako proměnná vystupuje právě čas.) Přík lad 1. stiff rovnice  y ' ' 101 y ' 100 y=0 má řešení  y=c 1 exp−100 x c 2 exp−x  Na následujícím obrázku je řešení rovnice s konstantami  c 1=c 2=1 a je vidět, že první složka řešení pro větší  x úplně vymizí, naopak pro malá  x je mnohem důležitější, než složka druhá.

Ještě názornější je příklad řešení rovnice  y ' =−100 y100 s počáteční podmínkou y 0= y 0 . Analytické řešení je  y= y 0 −1exp−100 x1 tedy pro veklá  x v podstatě  1 . Řeší­li se rovnice Eulerovou metodou y n1= y nh−100 y n100 , z čehož lze vyjádřit  y n= y 0 −11−100 hn 1 . Je vidět, že pokud se použije krok  h0.02 , první člen v absolutní hodnotě roste a řešení je zcela chybné. Pro stiff problémy jsou tedy mnohem vhodnější implicitní metody, protože dovolují podstatně delší krok  h . Příklad v PASCALU ODE41.PAS . Implicitní Eulerova metoda má tvar  y n1 = y nh f  x n1 , y n1  , a pokud ji aplikujeme na předchozí rovnici y 0 −1 dostaneme  y n =1  n , tedy řešení konvergující k  1 pro libovolné  h . 1 100 h Metody, které dovolují u stiff rovnic použít delší krok se nazývají stiff stabilní a zabýval se jimi C.W.Gear. Stabilita s konečným krokem viz. slidy. Příklad v PASCALU IME41.PAS . Semiimplicitní Eulerova metoda Implicitní metody jsou vhodné pro lineární diferenciální rovnice, kde se při výpočtu

y n1 řeší lineární rovnice. Řešení nelineárních rovnic je obtížnější a u systémů rovnic se to ještě víc komplikuje. Proto se implicitní rovnice předělá tak, aby byla lineární v y n1 a takové metodě se říká semiimplicitní. Linearizace se provádí rozvojem funkce ∂f x , y  . f do Taylorovy řady  f  x n1 , y n1 = f  x n1 , y n  y n1− y n  ∂ y n1 n Dosazením do implicitní Eulerovy metody dostaneme  −1 ∂f y n1= y n h 1 −h  x , y  f  x n1 , y n  . ∂ y n1 n

[

]

Pro řešení stiff problémů se používají: ●  Rosenbrockovy metody (semiimplicitní zobecnění Runge­Kuttových metod) ●  Semiimplicitní zobecnění Bulirsch­Stoerovy metody ●  Vícekrokové Gearovy metody (semiimplicitní zobecnění metod prediktor­ korektor)

Okrajová úloha