7    Numerické metody

počáteční úloha
diskretizace proměnných
uzly sítě
diferenční metody
krok metody
jednokrokové metody
Eulerova metoda
řád metody
metody Runge-Kutta
vícekrokové metody
implicitní a explicitní metody
samostartovací metody
Adams-Bashfothova a Adams-Moultonova metoda
prediktor, korektor
krokové a globální chyby výpočtu

7.1    Numerická metoda Runge-Kutta
7.2    Stoermerova formule a metody GBS
7.3    Vícekrokové Adamsovy metody
7.4    Algoritmus prediktor-korektor
7.5    Odhad chyby výpočtu

Při simulaci dynamických systémů reprezentovaných soustavou diferenciálních rovnic (2.2) se používá řada numerických metod, vyvíjených od počátku tohoto století. Metody lze rozčlenit podle mnoha kriterií do různých tříd a studovat přitom jejich efektivitu z hlediska přesnosti a numerické náročnosti. Šíře teorie v této oblasti často odrazuje od použití numerických metod v běžné praxi, proto se zde seznámíme pouze s jednoduššími postupy, které postačí na simulaci demonstrovaných soustav (je třeba poznamenat, že pro chaotické systémy stejně každá metoda dříve nebo později selže v tom smyslu, že pro každou přesnost existuje konečný čas, po kterém je požadovaná tolerance překročena).

Uvažujme nejprve nejjednodušší systém s jedinou stavovou proměnnou (fázovým prostorem je přímka R):

(7.1)

s počáteční podmínkou x(t₀)=x₀, jehož řešením je trajektorie x=x(t). Protože se jedná o nalezení trajektorie splňující počáteční podmínku zadávanou obvykle definicí stavu systému v čase t₀, nazývá se uvedený problém počáteční úloha.

U většiny praktických úloh nelze nalézt analytické řešení, proto je třeba trajektorii hledat numerickou metodou. Základním principem numerického řešení je obvykle diskretizace proměnných, kdy řešení (spojitou trajektorii) nahradíme posloupností diskrétních bodů

x₀(t₀), x₁(t₁), x₂(t₂), ... , kde t₀ < t₁ < t₂ < ...

Hodnoty nezávisle proměnné t tvoří tzv. uzly sítě. Pokud budou body posloupnosti dostatečně husté, lze je považovat za přibližnou reprezentaci hledané trajektorie. Uvedené metody se nazývají diferenční. Vzdálenost mezi dvěma sousedními body sítě

(7.2)

se nazývá krok metody a nemusí být obecně konstantní (síť není ekvidistantní), avšak my budeme používat metody výhradně s konstantním krokem (h_i= h). Nezávisle proměnnou v našich systémech je totiž vždy čas a při simulaci (tj. nalezení a zobrazení nového stavu systému) potřebujeme, aby byl čas homogenní, neboť pak zobrazíme správně i rychlost vývoje systému.

Nový stav systému v čase t_i vždy počítáme na základě předcházejích stavů systému v čase t_i-1, t_i-2 , ... , t_i-k. Pokud je k =1 hovoříme o jednokrokové metodě, tj. nový stav systému odvodíme z předcházejícího stavu.

Nejjednodušším představitelem jednokrokových metod je Eulerova metoda, kde nový stav lze stanovit z předcházejícího na základě vztahu

(7.3)

Uvedený vztah lze jednoduše odvodit dvěma způsoby:

1. Nahrazením vztahu (7.1) diferenčním vztahem

(7.4)

odkud ihned plyne Eulerův vztah

2. Taylorův rozvoj v okolí bodu x_n= x(t_n) lze vyjádřit řadou

(7.5)

což pro ekvidistantní síť dává po dosazení z (7.1)

(7.6)

Vidíme, že Eulerův vztah je lineární aproximací Taylorova rozvoje.

Protože je funkce f(x,t) ohraničená včetně všech svých derivací (systém by v žádném stavu neměl mít nekonečnou rychlost, zrychlení apod.), je Eulerův vztah přibližný s chybou úměrnou druhé mocnině kroku h (tzv. diskretizační chyba metody je O(h²), což znamená, že pro dostatečně malý krok je chyba odhadu nového stavu |e_n|< Ch², kde C je konstanta). Metoda s chybou O(h²) je tzv. metoda prvního řádu (obecně pro chybu kroku O(h^p) je metoda řádu p-1).

Eulerova metoda je v praxi obvykle nepoužitelná. Vypočtený nový stav odpovídá skutečnosti jenom při konstantní rychlosti systému (viz obr.7.1), při každém zakřivení trajektorie (zrychlení), přejdeme při Eulerově metodě na jinou blízkou trajektorii.

Obr.7.1 Geometrická interpretace Eulerovy metody. Nový bod x_i+1 leží ve směru tečny v bodě x_i ke skutečné trajektorii x(t). Po dvou krocích způsobí chyba metody přechod z trajektorie x_a(t)na x_b(t).

Vychází opět z Taylorova rozvoje a bere do úvahy i členy vyšších řádů. Potřebné derivace funkce f(t,x) počítá složitější diferenční metodou pomocí dalších pomocných bodů mezi sousedními uzly v síti. Metod Runge-Kutta je více, nejoblíbenější je tzv. klasická metoda 4. řádu, která pro naše potřeby postačí. Vyjádření nového stavu systému pomocí předcházejícího stavu je dáno vztahy:

(7.7)

Pomocné hodnoty k_i představují derivace stavu (rychlosti) systému ve speciálních bodech na začátku, konci a uprostřed intervalu <t_n , t_n+1>.

Poznámka:
Diskretizační chyba metody je O(h⁵), tedy pro krok h=0.1 je chyba řádu 10^-5. Pro dosažení stejné přesnosti bychom museli v rozmezí kroku h vykonat 1000 mezikroků Eulerovou metodou. Takové úvahy je ovšem třeba brát s rezervou, neboť skutečná chyba závisí na složitosti reálné trajektorie a odhad chyby nic neříká o ohraničující konstantě C (viz výše).

V případě fázového prostoru vyšší dimenze m (m>1) je stav systému popsán vektorem x a dynamický systém vektorovou funkcí f(t,x). Algoritmus Runge-Kutta je pak třeba psát ve vektorové formě, tj. pomocné hodnoty k_i budou nyní m-složkové vektory, takže např. vztah v (7.7) k₁= f (t_n , x_n) tvoří soustavu rovnic:

Klasická metoda Runge-Kutta nepatří k nejmodernějším, ale dosud se často používá, i když se pro přesnější výpočet příliš nehodí. Populární je také adaptivní metoda Runge-Kutta, kdy se vypočítá nový stav pomocí jednoho a následně podruhé pomocí více kroků, a pokud diference výsledků překračuje zadanou toleranci, dělení sítě se zjemní. Při animaci však tuto metodu použít nemůžeme, neboť potřebujeme pro výpočet nového stavu vždy stejný čas výpočtu pro znázornění odpovídající rychlosti.

Pro srovnání zde připomeneme velmi starou numerickou metodu (1907), kterou lze použít pro řešení počátečních úloh pro konzervativní dynamické systémy definované diferenciální rovnicí druhého řádu:

(7.8)

s počátečními podmínkami

Vztah (7.8) je obvykle přímo pohybová rovnice simulovaného systému. Místo obvykle doporučovaného přeformulování rovnice (7.8) do soustavy rovnic prvního řádu (př.2.1), lze použít Stoermerovu formuli, založenou na rozdělení diskretizačního intervalu H na m ekvidistantních podintervalů h. Zde použijeme Henricim modifikovaný tvar formule pro zmenšení zaokrouhlovací chyby:

Pro i = 1,2,..,m-1

Výsledný nový stav po kroku H je pak dán dvojicí x_m,v_m, kde

(7.9)

Formule tvoří metodu 2. řádu a lze ji použít přímo na pohybové rovnice systému. Metoda není přesná, ale je velmi stabilní, a někdy se proto používá pro simulaci vícečásticových soustav. Pro fázový prostor vyšší dimenze n (n>1) je opět nutné přeformulovat algoritmus do vektorové podoby (srov. metodu Runge-Kutta), kde pomocné hodnoty D_i jsou vektory a funkce f(t,x) je rovněž vektorová.

Obdobné metody dělení základního kroku se nazývají souhrnně modifikované obdélníkové pravidlo a používají se v souvislosti s velmi efektivní metodou GBS (Gragg, Bulirsch, Stoer). Podstatou této metody je výpočet nového stavu v několika verzích pro různé dělení základního kroku H. Výsledná posloupnost stavů pro různá m se pak extrapoluje (např. Richardsonovou extrapolační metodou) pro , což odpovídá nekonečně jemnému dělení základního intervalu.

Vícekrokové metody používají k vyjádření nového stavu systému znalost k předcházejících stavů. Obecně nelze říci, že vícekroková metoda je nutně efektivnější než jednokroková. Avšak za předpokladu spojitého systému, který se nevyvíjí příliš dramaticky, je pravděpodobné, že odhad nového stavu bude přesnější. Pokud odhlédneme od paměťových nároků metody (nutnost uložení více stavů), obecně neroste výpočetní složitost s počtem kroků k. Nový stav systému vyjádříme obecně pomocí vztahu (tzv. lineární k-kroková metoda):

(7.10)

kde a_i a b_i jsou vhodné konstanty. Je-li konstanta b₀ą0, je nový stav systému závislý na neznámé hodnotě funkce f_n+1= f ( t_n+1, x_n+1) a vztah (7.10) je třeba řešit iteračně. Metoda se v takovém případě nazývá implicitní a její výpočetní složitost závisí na rychlosti konvergence iteračního procesu. Pro b₀=0 se metoda nazývá explicitní, její výpočet je jednodušší, avšak oblast stability metody je vždy menší než v případě implicitních metod. Řešení explicitních metod se provádí numerickou integrací soustavy (7.1), přičemž nahradíme předpokládanou trajektorii polynomem stupně k a využijeme znalost k předcházejících stavů. Uvedená metoda se nazývá Adams-Bashforthova a nový stav systému vyjadřuje vztahem

(7.11)

kde konstanty b_i závisí na počtu kroků k a pro k<6 jsou uvedeny v tabulce tab.7.1. Řád metody odpovídá počtu kroků, k výpočtu nového stavu je třeba mít k dispozici k předcházejících stavů. Při startu metody se tyto počáteční stavy zjišťují např. metodou Runge-Kutta (z tohoto důvodu se nazývají jednokrokové metody též samostartovací).

Tab. 7.1    Koeficienty Adams-Bashforthovy metody

Stejným postupem nahrazení trajektorie polynomem, avšak s připuštěním neznámé hodnoty f_n+1, dostáváme po numerické integraci (7.1) tzv. Adams-Moultonův vztah

(7.12)

kde konstanty rovněž závisí na počtu kroků k a pro k<6 jsou uvedeny v tabulce tab.7.2. Vidíme, že pro k kroků je metoda řádu k+1, její důležitou předností je kromě toho vyšší stabilita, avšak problémem je nalezení neznámé hodnoty f_n+1.

Tab. 7.2   Koeficienty Adams-Moultonovy metody

Problém výpočtu neznámé hodnoty f_n+1 ve vztahu (7.12) iteračním postupem lze vyjádřit vztahem

(7.13)

kde hodnota ^s+1x_n+1 vyjadřuje s+1 iteraci a člen g je v iteračním vztahu konstantní (reprezentuje explicitní část metody).

Lze ukázat, že pro spojitou a omezenou funkci f, uvedená metoda postupných aproximací vždy konverguje. Účinnost metody závisí mimo jiné na výběru prvního odhadu stavu ⁰x_n+1. Algoritmus prediktor-korektor používá pro první odhad Adams-Bashforthovu metodu pokud možno stejného řádu, jako je následně použitá implicitní metoda. První odhad explicitní metodou se nazývá prediktor (P), jeden krok iterace se skládá z výpočtu funkce f  (evaluace - E) a následné opravy nového stavu podle vztahu (7.13), tzv. korektor (C). Pro n iterací se obecně metoda značí P(EC)ⁿ. Teorie ukazuje, že z hlediska efektivity algoritmu obvykle postačuje n<3. Stabilitu metody lze zvýšit dalším vyčíslením funkce f při konečně stanovené hodnotě x_n+1. Tento nejčastěji použitý tvar algoritmu se označuje P(EC)ⁿE. Pro simulaci systémů, kde jsou jednokrokové metody Runge-Kutta již zjevně neefektivní, budeme používat algoritmus typu PECE, který má největší oblast stability.

Celkový algoritmus má tedy následující kroky:

1. pro k krokovou metodu vypočítáme např. klasickou metodou Runge-Kutta k hodnot stavů x_i a odpovídajících hodnot funkcí f_i (nastartování metody)
2. podle vztahu (7.11) vyjádříme odhad nového stavu x_n+1 (prediktor)
3. vypočítáme hodnotu funkce f_n+1 (evaluace)
4. podle vztahu (7.12) opravíme hodnotu x_n+1 a definujeme ji jako nový stav systému (korektor).
5. vyjádříme novou hodnotu funkce f_n+1 (další evaluace)
6. přejdeme na krok 2 pro počátek výpočtu stavu x_n+2.

Chyby numerických metod se zásadně dělí na krokové a globální (chyby vznikající kumulací chyb v jednotlivých krocích). Obecně platí, že je-li diskretizační chyba řádu O(h^p), pak globální chyba je O(h^p-1), detailní teorie je však značně složitá. Chybu lze v době výpočtu sledovat např. srovnáváním výsledků při zvýšení hustoty diskretizačních bodů, musíme však mít jistotu, že použitá metoda konverguje. Jinou možností je srovnávat výsledek s očekávaným chováním vyplývajícím ze známých zákonitostí simulovaného systému (např. symetrie). Při simulaci konzervativních dynamických systémů máme jednoduchou možnost přímého sledování chyby výpočtu. Chyba, které se dopustíme při výpočtu nového stavu totiž způsobí, že se v novém stavu přesuneme na jinou blízkou trajektorii, odpovídající jiným počátečním podmínkám (obr.7.1). Protože však v systému existují integrály pohybu (např. celková energie, případně další integrály - hybnost, různé momenty apod.), lze během výpočtu neustále ověřovat jejich zachování, neboť na jiných trajektoriích budou mít integrály pohybu patrně jinou hodnotu.