Numerické řešení soustav obyčejných diferenciálních rovnic v modelování a simulaci

Tomáš Vogeltanz, Roman Jašek

Univerzita Tomáše Bati ve Zlíně
Fakulta aplikované informatiky
Nad Stráněmi 4511, 760 05 Zlín
Česká republika

vogeltanz@fai.utb.cz, jasek@fai.utb.cz

Abstrakt: Tato práce pojednává o nejčastěji využívaných numerických metodách při řešení soustav obyčejných diferenciálních rovnic vmodelování a simulaci systémů. Celá řada úloh, které se vyskytují vmodelování asimulací systémů, není analyticky řešitelná nebo je nalezení přesného řešení příliš obtížné. Numerická matematika se vtěchto případech snaží nalézt alespoň řešení přibližné. Na začátku této práce je uveden popis soustav obyčejných diferenciálních rovnic a jejich numerického řešení. Dále jsou popsány Eulerovy a Runge-Kuttovy numerické metody, které se nejčastěji využívají pro řešení obyčejných diferenciálních rovnic, a to včetně uvedení odhadu chyb a stability těchto metod.

Klíčová slova: Numerické metody, Obyčejné diferenciální rovnice, Eulerova metoda, Runge-Kuttovy metody, Modelování, Simulace

Abstract: This paper discusses the most common numerical methods used to solve systems of ordinary differential equations in modeling and simulation. A number of tasks which occur in modeling and simulation cannot be solved analytically or it is too difficult to find their exact solution. In these cases, numerical mathematics tries to find at least approximate solution. At the beginning of this paper are mentioned the description of systems of ordinary differential equations and the description of their numerical solution. After that, the Euler and Runge-Kutta methods, which are the most commonly used numerical methods for finding the solution of ordinary differential equations, are described, including the error estimates and stability.

Keywords: Numerical Methods, Ordinary Differential Equations, Euler Method, Runge-Kutta Methods, Modeling, Simulation

Úvod

Většina fyzikálních zákonů je formulována ve formě diferenciálních rovnic. Jsou to bud’ parciální nebo obyčejné diferenciální rovnice. V této práci jsou popsány obyčejné diferenciální rovnice. Tyto rovnice nejčastěji popisují závislost fyzikálních veličin na čase. Hlavním představitelem těchto rovnic ve fyzice jsou rovnice pohybové, které popisují pohyb těles pod vlivem vzájemných a vnějších sil. Tyto rovnice není často možné analyticky řešit a je tedy třeba použít metod numerické matematiky. [1]

Omezíme se na diferenciální rovnice prvního řádu, ve kterých je derivace přímo vyjádřena. Hledáme tedy reálnou funkci y(x), která splňuje rovnici (1) – tato počáteční (Cauchyho) úloha bývá zapisována také jako (2). Funkcí, které splňují tuto rovnici, existuje většinou nekonečně mnoho. Jedno konkrétní řešení určíme tzv. počáteční podmínkou, které je nezbytnou součástí zadání úlohy. Počáteční podmínka má tvar (3), kde x0 a y0 jsou zadané hodnoty. Dále předpokládejme, že jsou splněny podmínky existence a jednoznačnosti řešení voblasti Ω × I, kde Ω je stavový prostor uvažovaného systému a I je nějaký interval reálných čísel. [1] [2] [3]

	(1)
	(2)
	(3)

Soustavy nelineárních rovnic ve tvaru (4) můžeme zapisovat také ve vektorovém tvaru (5). Zde se však jedná o vyjádření stacionárního (ustáleného) stavu systému na rozdíl od (1), která vyjadřuje dynamiku systému. [4]

	(4)
	(5)

Pokud Fx vyjádříme jako derivaci n vektorových funkcí podle vektoru x, tak platí (6). [5]

(6)

Pro jednoduchost uvažujme soustavu n obyčejných diferenciálních rovnic 1. řádu (7) spočátečními podmínkami . [2] [3]

(7)

Jako příklad vyjádření soustavy rovnic (7) můžeme uvést soustavu dvou rovnic pro dvě neznámé funkce y1(x), y2(x) (8). [1]

(8)

Pokud je třeba řešit diferenciální rovnici, ve které vystupuje vyšší derivace, lze zavedením další proměnné převést tuto rovnici na soustavu rovnic prvního řádu. [1]

Numerickým řešením počáteční úlohy se rozumí posloupnost hodnot {yi}, tj. (9), které odpovídají hodnotám xi, i = 1,2,...,k, nezávisle proměnné x (zpravidla času) v intervalu <x0, xk>. [2] [3]

(9)

Hodnota výrazu přitom udává velikost integračního kroku. Exaktní řešení uvažované počáteční úlohy značíme jako . [3]

1 Numerické metody

Cílem numerických metod je nalezení numerického řešení. Přitom se požaduje, aby posloupnost {yi} konvergovala pro h → 0 k exaktnímu řešení Y(xi). [2] [3]

Při numerickém řešení diferenciální rovnice vyjdeme z bodu x0 ve kterém máme zadánu počáteční podmínku. Poté se budeme snažit najít řešení v dalších bodech xi. V těchto bodech získáme odhady přesného řešení yi. [1]

Pokud k výpočtu nové hodnoty použijeme pouze jednu předchozí hodnotu, tj. (10), nazývá se tato metoda jednokroková. Pokud použijeme n předchozích bodů (11) nazývá se tato metoda mnohokroková, přesněji n-kroková. Vzdálenost dvou sousedních bodů (neboli krok) se dá vyjádřit jako (12). [1]

Krok může být i proměnlivý, ovšem podle praktických zkušeností je lepší využít konstantní krok, protože algoritmy adaptace velikosti kroku mohou svou nepřesností „přehlédnout“ kmity (a tudíž je neuvažovat) nebo je proložit přímkou (aproximovat), což ovlivní průběh i výsledek řešení.

	(10)
	(11)
	(12)

Kvalita použité numerické metody se hodnotí podle těchto základních kritérií: [2] [3]

• Přesnost metody (velikost lokální chyby a prostředky pro její odhad) - Metody vyšších řádů jsou obecně přesnější (mají menší lokální chybu) než metody nižších řádů téhož typu.
• Stabilita metody
• Časová náročnost výpočtu - Rychlost řešení je pro daný integrační krok h nepřímo úměrná řádu metody. Vícekrokové metody při daném integračním kroku jsou rychlejší než metody jednokrokové.
• Nároky na operační paměť počítače - Nároky, které klade metoda na paměť, jsou přímo úměrné řádu metody. Pokud uvažujeme stejný řád metod a stejný počet iterací v jednom kroku, tak lze říci, že metody vícekrokové mají větší nároky na paměť, než metody jednokrokové. U vícekrokových metod mají implicitní metody větší nárok na paměť, než metody explicitní. Dále je zřejmé, že metody s proměnným krokem integrace mají větší nároky na paměť, než metody s pevným krokem integrace.
• Využitelnosti

1.1 Chyby

Vypočítané hodnoty řešení yi nejsou zcela přesné, ale jsou zatíženy chybou. Přesnost výsledku vyjadřuje tzv. akumulovaná diskretizační chyba. Tato chyba udává konečný rozdíl numerického řešení yi a přesného řešení Y(xi) po n krocích (13) – vněkterých literaturách je také uveden výpočet (14), který je svou podstatou totožný, jen bude ve výsledku opačné znaménko. [1] [3]

	(13)
	(14)

Akumulovaná chyba vzniká jako součet dílčích chyb v jednotlivých krocích řešení. Chyba, která vznikne v jednom kroku, se nazývá lokální diskretizační chyba. Tuto chybu definujeme jako chybu našeho přibližného vzorce po dosazení přesného řešení. Tato chyba je součtem dvou chyb:

· chyby metody (truncation error - eT), která je způsobena respektováním jen konečného počtu členů Taylorova rozvoje funkce f

· zaokrouhlovací chyby (rounding error - eR), jež souvisí somezenou délkou slova v paměti počítače (dnes se již ovšem nejedná o až tak velký problém)

Chyba jednoho kroku pak přirozeně ovlivňuje výsledky kroků následujících. Lokální diskretizační chybu většinou upravujeme rozvojem do Taylorovy řady. Pokud rozvoj chyby začíná mocninou hp+1, tak říkáme, že metoda má řád p. [1] [3]

1.2 Eulerova metoda

Nejjednodušší metodou na řešení diferenciálních rovnic je Eulerova metoda. Pokud vnaší diferenciální rovnici (15) nahradíme derivaci přibližným vzorcem, dostaneme rovnici (16). [1] [6]

	(15)
	(16)

Když vyjádříme yi+1 dostaneme Eulerovu metodu:

.

(17)

Tento vzorec nám umožňuje při znalosti řešení v xi vypočítat řešení v xi+1. Jedná se tedy o jednokrokovou metodu. [1]

Taylorův rozvoj odhadu lokální diskretizační chyby pro Eulerovu metodu začíná (18). V tomto výrazu je druhá mocnina h a Eulerova metoda je tedy metodou prvního řádu. [1]

(18)

Eulerova metoda je velice jednoduchá, ale bohužel nepřesná [1]. Výrok o nepřesnosti ovšem platí jen tehdy, pokud uvažujeme použití stejného kroku pro jinou metodu. Když použijeme u Eulerovy metody krok 0.001 a u Runge-Kuttovy 0.1, tak můžeme dostat lepší výsledky u Eulerovy metody.

Hlavní nepřesnost metody je způsobena tím, že při kroku z xi do xi+1 používáme na celém intervalu konstantní hodnotu derivace . Tato hodnota by se však správně měla vprůběhu kroku měnit. Tuto chybu se snaží odstranit modifikace Eulerovy metody - jedná se opět o jednokrokové metody. [1]

První modifikace se snaží vystihnout změnu derivace tím, že místo derivace na začátku intervalu použije derivaci uprostřed intervalu. Doprostřed intervalu se dostaneme obyčejnou Eulerovou metodou s poloviční délkou kroku. Výsledný vzorec je:

.

(19)

Druhá modifikace používá průměr z derivace na začátku a na konci intervalu:

.

(20)

Lze dokázat, že obě tyto modifikace jsou metody druhého řádu. [1]

1.3 Runge-Kuttovy metody

Dalšími modifikacemi Eulerovy metody získáme tzv. Runge-Kuttovy metody. Při těchto metodách získáme několik odhadů derivace ki vs různých bodech. Pro celý krok potom použijeme vážený průměr těchto derivací s váhami wi tj. (21). Odhady derivací se vypočítají pomocí vztahů (22). [3] [6]

	(21)
	(22)

Pokud použijeme jen jeden bod (tj. s=1) dostaneme metodu prvního řádu – Eulerovu metodu. Pokud použijeme dva body, můžeme dostat modifikované Eulerovy metody. První modifikaci odpovídá volba , druhé modifikaci odpovídá , . S rostoucím počtem bodů dostáváme metody vyšších řádů. [1] [3]

Při dvou odhadech derivace dostáváme metodu druhého řádu (23). Je-li f(x,y) pouze funkcí x, tak jde o aplikaci lichoběžníkového pravidla. [3]

(23)

Při třech odhadech derivace můžeme dostat metodu třetího řádu (24). Je-li f(x,y) pouze funkcí x, tak jde v podstatě o použití Simpsonova pravidla. [3]

(24)

Při čtyřech odhadech derivace můžeme dostat následující metodu čtvrtého řádu (25) [1] [6]. Samozřejmě, že lze sestavit více takových rovnic, které budou různé, ovšem principiálně stejné. Např. pro k3 je možné počítat ve vzorci i sk1, stejně tak i pro k4. Pro daný případ níže jsou u předchozích odhadů derivací kuvedeny nulové parametry, takže se neuvažují.

(25)

U dosud uvedených Runge-Kuttových metod se řád metody vždy rovnal použitému počtu odhadů derivace. Tak to ale vždy není, např. při pěti odhadech derivace se dá sestrojit pouze metoda čtvrtého řádu. I z toho důvodu patří uvedená metoda čtvrtého řádu k velice oblíbeným. [1]

K odhadování přesnosti Runge-Kuttových metod se používá tzv. metoda polovičního kroku (nebo také Richardsonův princip), tj. dvojí výpočet: jednak s krokem h, jednak s krokem 2h. Pro odhad chyby při kroku 2h dostaneme (26), kde (2)y resp. (1)y značí aproximaci řešení s krokem 2h resp. h a p je řád metody. Z hlediska rychlosti výpočtu je použití tohoto vztahu značnou ztrátou - vyžaduje provedení alespoň o třetinu více strojových operací, než je třeba k vlastnímu výpočtu. Tento odhad lze použít i pro vícekrokové metody. [2] [3]

(26)

Implementace Runge-Kuttových metod na počítači je velmi jednoduchá, krok lze libovolně měnit. Mají přibližně stejnou, často i vyšší přesnost než metody prediktor-korektor stejného řádu. Pro řešení na celém intervalu jsou někdy považovány za méně vhodné než metody vícekrokové - potřebují totiž v každém kroku tolik výpočtů funkční hodnoty f(x,y), jaký je řád metody. [2] [3]

1.4 Stabilita numerického řešení

Metoda je absolutně stabilní pro daný krok h a danou soustavu diferenciálních rovnic, jestliže chyba vzniklá při výpočtu hodnoty yn se nezvětšuje při výpočtu následujících hodnot yk, k > n. K vyšetřování absolutní stability se používá testovací rovnice (27), kde λ je konstanta (obecně komplexní číslo). Množina hodnot (28), pro které je metoda absolutně stabilní, se nazývá obor absolutní stability příslušné soustavy. [3]

	(27)
	(28)

Všechny Runge-Kuttovy metody mají ohraničený obor absolutní stability, definovaný nerovností (29), kde p je řád metody a je definováno v (28). Oblast absolutní stability se zvětšuje s rostoucím řádem metody. [2] [3]

(29)

Na Obr. 1 jsou zobrazeny oblasti absolutní stability nejpoužívanějších metod. [2]

Obr. 1 Oblasti absolutní stability

Závěr

Vtéto práci byl uveden popis soustav obyčejných diferenciálních rovnic včetně různých forem zápisu a jejich numerického řešení. Dále byly popsány Eulerovy a Runge-Kuttovy metody pro výpočet numerického řešení obyčejných diferenciálních rovnic. Byl zde uveden nejen postup těchto metod, ale i odhad chyby (resp. přesnosti) řešení a stability. Byl uveden graf se zobrazením oblasti absolutní stability nejpoužívanějších metod (včetně metody prediktor-korektor a Adamsovy explicitní metody).

Zgrafu jasně vyplývá, že největší oblast absolutní stability má vdaném případě Runge-Kuttova metoda čtvrtého řádu. Dále bylo konstatováno, že se, také i vzhledem kvysoké přesnosti, řadí tato metoda kvelice oblíbeným. Při zvyšování řádu Runge-Kuttovy metody však vzniká nutnost počítat stále více a více odhadů derivací, což může vést kvýraznému zpomalení výpočtu.

Eulerova metoda je na rozdíl od Runge-Kuttovy metody jednodušší, rychlejší a méně přesná, ovšem jen pokud uvažujeme stejný krok. Při výrazně menším kroku může být Eulerova metoda přesnější. Stejně tak je to i srychlostí výpočtu těchto metod – obecně je rychlejší metoda Eulerova (resp. v podstatě Runge-Kuttova metoda prvního řádu). Ovšem při použití menšího kroku může být Eulerova metoda pomalejší a zároveň přesnější než Runge-Kuttova metoda čtvrtého řádu.

Poděkování

Tento článek vznikl za podpory Interní Grantové Agentury IGA/FAI/2014/006.

Seznam použité literatury

[1] VICHER, Miroslav. Numerická matematika. UJEP. Ústí n/L. 7. 1. 2003. ISBN 80-7044-516-5.
[2] RÁBOVÁ, Zdeňka, Milan ČEŠKA, Jaroslav ZENDULKA, Petr PERINGER a Vladimír JANOUŠEK. Modelování a simulace. Brno: VUT, 2005, 227 s. ISBN 80-214-0480-9.
[3] KŘIVÝ, Ivan a Evžen KINDLER. Simulace a modelování. Vyd. 1. Ostrava: Ostravská univerzita Ostrava, 2001, 146 s. ISBN 80-704-2809-0.
[4] JEŘÁBEK, Tomáš. Metody pro řešení soustav nelineárních rovnic [online]. 2009 [cit. 2013-11-27]. Diplomová práce. Masarykova univerzita, Přírodovědecká fakulta. Vedoucí práce: Ivanka Horová. Dostupné z: <http://is.muni.cz/th/77840/prif_m/>.
[5] ŠTECHA, Jan a Vladimír HAVLENA. Teorie dynamických systémů. Praha: Vydavatelství ČVUT, 2005, 248 s. ISBN 80-010-1971-3.
[6] KUBÍČEK, Milan, Miroslava DUBCOVÁ a Drahoslava JANOVSKÁ. Numerické metody a algoritmy. 2. oprav. vyd. Praha: Vydavatelství VŠCHT Praha, 2005, 188 s. ISBN 80-708-0558-7.