TIP: Přetáhni ikonu na hlavní panel pro připnutí webu
Článek volně navazuje na předchozí díly a rozebírá nejdůležitější vlastnosti Eulerovy metody řešení Cauchyho počáteční úlohy, chyby a nestability výpočtu.
Od píky
Vezměme si tedy nejjednodušší model, který popisuje například základní vlastnosti (odborně řečeno lineární model prvního řádu). A protože lenost zasáhla i mou osobu, odvolám se na své předcházející články o simulaci. První oblast zájmů – abychom měli nějakou hodnotnou vztažnou informaci – bude analytické řešení dané úlohy. Jen připomenu, že jde o diferenciální rovnici se speciální (konstantní) pravou stranou, jejíž řešení jsem popsal v článku simulace3. Výsledkem je exponenciální funkce proudu, která je závislá na proměnné času (proud se mění s časem). Tento výsledek považujeme za přesný, protože jsme jej dostali analytickým postupem, který vychází z přesně definovaných operací nad množinou komplexních čísel – žádné zanedbávání drobných čísel, žádné zaokrouhlování, jen exaktní matematika (pozn.: Těm z vás, kteří neovládají vyšší matematiku, musí stačit výsledek – poslední řádek odvození).
Tuto časovou funkci (průběh proudu cívkou v závislosti na čase) si vykreslíme do grafu pro trochu vyšší počet hodnot. Je to taková příprava na porovnání s numerickým řešením, kde budeme potřebovat větší množství kroků.
Nyní přistoupíme k výpočtu numerickému pomocí základní Eulerovy metody. Je to metoda tzv. jednokroková, což znamená, že k výpočtu následující hodnoty nám postačí jeden vzorec. Opět jsme ho tu už měli, takže jen pro zopakování:
Problém, jak zapsat tento vzorec do buňky v Excelu, jsem vyřešil za vás, ale pro vysvětlení: odpor R, indukčnost L a napětí U jsou z pohledu jednoho kroku řešení konstantní. Proto jsem uvedl absolutní odkazy (se znaky $). Jedinou proměnnou je stavová veličina – proud, který získáme v každém kroku jako výsledek kroku předcházejícího. Ještě nesmíme zapomenout celý výsledek vynásobit krokem a přičíst jej k předcházející hodnotě výsledku. Potom stačí už jen vzorec překopírovat na celý sloupec a výsledné hodnoty vykreslit.
Pod lupou
Abychom mohli nějak danou metodu ohodnotit, zavedeme pojem odchylky numerického řešení (chyby výpočtu). Bude to rozdíl mezi hodnotou přesnou a hodnotou získanou numerickým počtem. Odchylku jsem převedl na procenta správné hodnoty.
Pokud se podíváme na oba průběhy proudu, zjistíme, že proud získaný numerickým řešením je větší (pro názornější ukázku jsem zvolil krok dt = 0,01).
S rostoucím krokem se zvětšuje i chyba výpočtu, přičemž platí předchozí vztah mezi chybou a velikostí kroku – pohybují se ve stejném řádu. Pokud si s velikostí kroku zahrajeme až moc (v našem případě dt = 0,2), výsledek se nám odporoučí do neskutečných hodnot.
Proč? Odpověď musíme najít v několika prvních krocích výpočtu. První krok je totiž počítán pouze jako součin budící veličiny (napětí) a kroku simulace. Pokud tedy bude krok velký, bude relativně velký i výsledek po prvním kroku, resp. chyba hned v prvním kroku bude velká. Ve druhém kroku se celá záležitost může a nemusí vyřešit – vše záleží na velikosti kroku. V případě, že rovnice modelu ve druhém kroku nabude hodnoty kladné, je tu naděje, že se budoucí výsledek ustálí na správné hodnotě (s tím vědomím, že přechodný děj neodpovídá skutečnosti). Pokud ale hodnota modelu dosáhne záporných hodnot, výsledek začne kmitat z kladných hodnot do záporných a jeho absolutní hodnota poroste. Můžeme to nazývat rozkmitáním systému nebo divergencí výsledku. Při běhu simulace ale často nemůžeme kontrolovat chybu (neznáme analytické řešení) a kmitání v systémech vyššího řádu je zcela normální. Obecně platí, že velikost kroku by měla být o mnoho menší než nejmenší časová konstanta systému. V případě systému prvního řádu je časová konstanta jen jedna a je rovna kladně vzaté převrácené hodnotě řešení charakteristické rovnice. Následující zápis vychází z výše uvedeného analytického řešení.
V našem případě je tedy volba dt = 0,01 postačující.
Eulerova metoda tedy není za všech okolností stabilní. V případě, že vhodně zvolíme krok, se nemusíme tomuto stavu ani přiblížit. Některé složitější systémy vyšších řádů ale mohou dosáhnout těchto stavů v nečekané okamžiky a v tom případě se můžeme jen divit, jak nám to ta simulace vychází. Nemusí to ale vždy být problém systému, naopak, dost často to bývá problém metody výpočtu. Proto jsou tu metody další, složitější, náročnější, které tento problém řeší. O jiných metodách si ale povíme příště.
Krok-sun-krok
Tento problém kroku lze řešit celkem elegantně, ale už musíme programovat. Můžeme totiž stanovit maximální odchylku mezi dvěma po sobě jdoucími kroky (rozdíl mezi i(n) a i(n+1)). Pokud bude odchylka větší než maximální, program zopakuje výpočet, ale zmenší podle nějakého pravidla krok výpočtu (zpravidla vydělí dvěma). Podobně je zkontrolován i výsledek po změně kroku. Rozhodování ilustrují následující dva obrázky:
Hodnota kroku výpočtu se zvětší v případě, že odchylka mezi dvěma kroky je malá (obvykle o řád nebo dva řády menší než maximální odchylka). Tímto postupem sice vnášíme do algoritmu další výpočty, ale v případě, že se stavová veličina ustálila, je přesnější výpočet (s malým krokem) zbytečný. Celková doba výpočtu může být (a většinou bývá) potom dokonce i kratší. Tento způsob používá i známý a často používaný program MATLAB a jeho nadstavba Simulink, kde je označován jako variable-step (proměnný krok), ale jeho implementace je daleko složitější a obsahuje i prvky proti rozkmitání kroku. Uvedený základní princip lze bez problémů naprogramovat pomocí podmínek.
Když je problémů málo, určitě se najdou další
V některých případech na vás může vybafnout ještě jeden záludný problém. Vyskytuje se ale jen u systémů, jejichž parametry jsou řádově značně odlišné. Je to problém interpretace čísel v počítači, jinak nazývaný zaokrouhlovací chybou. Číslo je totiž interpretováno dvěma způsoby – pevnou a plovoucí řádovou čárkou. V pevné řádové čárce osobní počítače (a tedy i simulace) obvykle nepracují, neboť mají implementovánu variantu druhou, lepší. Plovoucí řádová čárka ale skrývá jeden neduh. Číslo je uloženo jen na omezený počet platných číslic a doplněno kódem exponentu (řádu). Takto je zvětšen rozsah hodnot na neskutečných 10^(+/−300) (přibližně). Přesto ale dokáže uchovat jen omezený počet platných číslic. Upozorňuji na to proto, že ve výpočtech se může objevit kombinace postupných násobení dvou čísel řádově značně odlišných (např. 11e−6 × 1.12), přičemž např. jeden z operandů je stavovou (a tedy měnící se) veličinou. Zaokrouhlovací chyba se neprojeví jen jednou, ale postupně se zvětšuje a může dost ovlivnit výsledek, odborně se to nazývá kumulativní chyba. Nejnázornější by to mohlo být v některém z programovacích jazyků při výpočtech ve float. Proto je lepší (ale časově náročnější) používat typy objemnější, např. double. V Excelu se zaokrouhlovací chyba ukázat nedá, protože čísla v něm jsou interpretována ve velkém rozsahu (tj. přesnosti), větším než double formát.
Pro vás, kteří byste si chtěli se simulací v Excelu vyhrát, tu mám mou skromnou verzi ke stažení: num_derivace.xls (380 kB)
O čem to vlastně bylo?
Chtěl jsem vám v tomto dílu ukázat nejen možnost simulovat v tabulkovém editoru, ale i střízlivý přístup k výsledku simulace a jedno z možných řešení problému s krokem. V žádném případě se nejedná o vysoce odborné záležitosti. Moderní matematik by se nad tímto článkem pousmál. Přesto je ale spousta nematematických oborů, kde se využívají simulace i v této jednoduché podobě. V některém z příštích dílů se podíváme na zpětnou Eulerovu metodu, ve které otázka stability odpadá, ale vyvstávají další problémy (jak to tak chodí).
