Eulerova metoda je explicitní jednokrokový řešič prvního řádu (s chybou metody O(h^1)). Použitá metoda Runge-Kutta je potom explicitní jednokrokový řešič řádu čtvrtého (chyba metody O(h^4)).
Explicitní metody řešení Cauchyho úlohy pro SODE jsou vhodné pro řešení takzvaných non-stiff systémů. Pro pochopení jejich nedostatků při řešení stiff systémů (systémy, kde mají jednotlivé členy řádově rozdílné hodnoty) je potřeba pochopit koncept A-stability numerických metod.
Jednokrokové metody jsou vhodné pro řešení malých soustav nebo soustav pro které hrozí oscilace řešení.
Námi programované řešiče nedisponují funkcí adaptivní délky kroku. Délka kroku, H, musí být zadána jako jeden ze vstupních parametrů. Pokročilejší řešiče jsou schopny na základě a posteriori odhadu chyby řešení upravovat délku kroku v průběhu výpočtu tak, aby bylo dosaženo požadované přesnosti řešení. Příklad metody s adaptivní délkou kroku je možné najít například na wikipedii (poté, co si přečtete popis příkladu, zkuste se zamyslet, proč se algoritmy v matlabu jmenují ode45, ode23 atd., atp.).
K souboru se kterým jsme pracovali na cvičení jsem přidal problem3. Jedná se o soustavu dravec-kořist. Na cvičení jsme vždy zobrazovali takzvané integrální křivky jednotlivých funkcí řešení, x_i = x_i(t), i=1,2,...,n.
V případě problem3 jsem nezobrazoval integrální křivky, ale trajektorie řešení ve stavovém prostoru pro různé počáteční podmínky. Tomuto zobrazení se říká také fázový portrét, což je výraz, který doporučuji si zapamatovat.
V případě problem3 je stavový prostor R2, takže fázový portrét vykresluji do roviny.
Upravený kód naleznete zde. Obdobný kód, ale implementovaný v MAPLE je ke stáhnutí zde.
Na cvičení SODR III (zaměřeném na okrajové úlohy) bude nutné umět řešit soustavy nelineárních algebraických rovnic. Základním pro tento účel používaným algoritmem je Newtonova metoda.
Zkuste si projít přiložený kód. Proč Newtonova metoda selhává v případě problem2?
Případ problem3 je zase převzat z modelu dravec kořist (viz. výše). Tentokrát soustavu diferenciálních rovnic neřešíme, ale hledáme takový bod, kde se všechny derivace (rychlosti změny jednotlivých souřadnic) rovnají nule. Jinak řečeno takový bod, který když zvolíme jako počáteční podmínku, tak jej za libovolný časový interval (kladný i záporný) neopustíme. Jedná se tedy o jednobodovou trajektorii ve stavovém prostoru.
Takovému bodu potom říkáme stacionární/rovnovážný stav (steady state) soustavy. Která z trajektorií vykreslená souborem cv4_1.m odpovídá tomuto bodu?
Kód Newtonovy metody včetně příkladů naleznete zde.
V souboru výše je možné nalézt kód pro řešení SODE pomocí vícekrokového explicitního řešiče. Konkrétně pomocí Adams-Bashforthovy metody 4. řádu a tedy s chybou O(h^4). Součástí tohoto řešiče je i vestavěná metoda Runge-Kutta 4. řádu použitá k nastartování výpočtu.
Schéma řešení počáteční úlohy SODE Adams-Bashforthovou metodou k-tého řádu je přibližně:
Budete-li studovat kód, zkuste si zodpovědět následující otázky.
Jak se od sebe liší řešiče ab4 a ab4O? Kterou z těchto implementací byste použili? V čem jsou výhodné vícekrokové metody oproti jednokrokovým metodám vyšších řádů? V čem jsou naopak lepší metody jednokrokové?
Vysvětlení kódu metody prediktor-korektorV případě, že vyžadujeme vyšší přesnost řešení a nejsme ochotni/nemůžeme integrovat s kratším krokem (například případy, kdy integrace s krátkým krokem už interferuje s přesností počítače v některé z prováděných operací), je vhodné použít k řešení SODE některou z metod prediktor-korektor.
Tyto metody jsou obvykle (rozumnějme prakticky vždy) kombinací explicitní a implicitní metody.
Jenom pro pořádek, explicitní metoda určuje hodnotu řešení y_N na základě hodnoty předchozího kroku, y_N-1, a vyčíslení pravých stran (vektorového pole) v obecně k předchozích krocích. Implicitní metoda potřebuje vyčíslení vektorového pole i v kroku aktuálním (N-tém).
Prediktor je v této kombinaci tedy explicitní metoda, ze které získáme odhad řešení v N-tém kroku.
Korektor je metoda implicitní. Tedy abychom nalezli y_N, musíme řešit soustavu nelineárních algebraických rovnic
(1) 0 = y_N - y_N-1 - MAP(y_N,knowns),
kde MAP je nějaká lineární kombinace hodnot vektorového pole v k předchozích bodech (knowns) a neznámém y_N.
Jako počáteční nástřel pro numerický řešič této soustavy použijeme hodnotu dopočtenou v prediktoru. Samotný řešič můžeme použít libovolný. V souboru, který máte k dispozici je použita metoda prostých iterací, která konverguje pro dostatečně malý krok metody a dostatečně dobrý odhad získaný z prediktora. Další možností je použít například Newtonovu metodu, která konverguje výrazně rychleji, ale je citlivější na počáteční nástřel. Robustní volbou by bylo startovat řešič pomocí metody půlení intervalu, takto zpřesnit počáteční odhad prediktoru a až tento zpřesněný odhad použít jako počáteční nástřel Newtonovy metody.
V souboru, který máte k dispozici je jako prediktor použita Adams-Bashforthsova metoda 4. řádu (ta, kterou jsme procházeli posledně na cvičení). Jako korektor potom metoda Adams-Moultonova, také 4. řádu. Jako řešič soustavy (1) je (jak již bylo uvedeno výše) použita metoda prostých iterací
Jelikož je hlavní integrátor vícekrokový, je třeba jej startovat. Zde je k nastartování metody použita metoda Runge-Kutta 4. řádu, bez korektoru.