12. Metody řešení diferenciálních rovnic 1

Anotace

Ukážeme si některé metody pro numerické řešení obyčejných i parciálních diferenciálních rovnic. Probereme metody jako Eulerova metoda, Runge-Kuttovy metody, metody konečných diferencí, metodu konečných prvků a další. Jako doplněk si ukážeme řešení parciálních diferenciálních rovnic Fourierovou metodou separace proměnných, která umožňuje odseparovat časovou a prostorovou složku řešení a tím lépe pochopit některé vlastnosti řešení, zejména u rovnic popisujících vlnění.

12.1. Metody pro obyčejné diferenciální rovnice

12.1.1. Numerické řešení IVP

Numerické řešení diferenciální rovnice je zpravidla založeno na aproximaci derivace konečnou diferencí a postupným prodlužováním řešení od počáteční podmínky směrem dopředu nebo dozadu v čase. Otevřít prezentaci Spustit video

Obr. 12.1 Eulerova metoda s velmi dlouhým krokem (modrou barvou) zaostává za přesným řešením (šedou barvou). Pro lepší výsledek můžeme zmenšit krok nebo vylepšit metodu.

Obr. 12.2 Metoda Runge Kutta s velmi dlouhým krokem (modrou barvou, jde jasně vidět aproximace lomenou čarou). Přesné řešení je nakresleno šedou barvou.

Numerické řešení diferenciálních rovnic je základním nástrojem pro ukázku průběhu simulací pro dané hodnoty parametrů a počátečních podmínek. Jedná se o velice užitečnou a široce používanou činnost při inženýrských simulacích. Neprofesionálům často musí stačit použít hotové postupy, procedury a nástroje. Například Python je jednou z nejvhodnějších voleb.

Řešení počáteční úlohy lze numericky aproximovat poměrně snadno: začneme v bodě zadaném počáteční podmínkou a v okolí tohoto bodu nahradíme integrální křivku její tečnou. Tím se dostaneme do dalšího bodu, odkud opět integrální křivku aproximujeme tečnou. Směrnici tečny zjistíme z diferenciální rovnice, buď přímo z derivace (Eulerova metoda).

Vyjdeme-li z počáteční úlohy

\[y'=\varphi(x,y), \quad y(x_0)=y_0,\]

má lineární aproximace řešení v bodě \(\displaystyle [x_0,y_0]\) tvar

\[y=y_0+\varphi(x_0,y_0)(x-x_0).\]

Funkční hodnotu v bodě \(\displaystyle x=x_1\) označíme \(\displaystyle y_1\) a tento bod bude dalším body lomené čáry, tj.

\[y_1=y_0+\varphi(x_0,y_0)(x_1-x_0).\]

Hodnota \(\displaystyle x_1-x_0\) je krok Eulerovy metody označovaný \(\displaystyle h\). Tento postup opakujeme s počáteční podmínkou \(\displaystyle y(x_1)=y_1\). Iterační formule Eulerovy metody má potom následující tvar.

\[\begin{split}\begin{aligned}x_{n+1}&=x_n+h, \\ y_{n+1}&=y_n+\varphi(x_n,y_n)h.\end{aligned}\end{split}\]

Stačí tedy mít zvolen krok numerické metody (délku intervalu, na kterém aproximaci tečnou použijeme) a výstupem metody bude aproximace integrální křivky pomocí lomené čáry.

Vylepšení

• Pro přesnější aproximaci je možné zjemnit krok \(\displaystyle h\) (buď všude, nebo jenom tam, kde „je to potřeba“).

• Pro přesnější aproximaci je možné použít místo \(\displaystyle \varphi(x_n,y_n)\) lepší směrnici, která dokáže zohlednit, jestli se růst zrychluje nebo zpomaluje (metoda Runge Kutta druhého nebo čtvrtého řádu, …).

• Modely obsahující diferenciální rovnice obsahují zpravidla sadu parametrů charakterizujících fyzikální vlastnosti studovaných objektů. Pro numerické řešení musíme těmto parametrům dát konkrétní hodnoty a přicházíme tak o cennou informaci, jak řešení závisí na těchto parametrech. Vhodnou úpravou rovnice dokážeme počet parametrů eliminovat. Jednoduchým a často dostatečným způsobem je volba jednotek, obecnější metodou je transformace diferenciální rovnice uvedená v úvodní kapitole věnované diferenciálním rovnicím.

Online řešiče ODE (numericky):

• Sage

• Python

12.1.2. Malá odbočka - zaokrouhlovací chyby v numerických výpočtech

Obr. 12.3 Součást protiraketového systému Patriot. Raketu Scud vystřelenou 25.2.1991 systém nesestřelil vinou zaokrouhlovací chyby. Zdroj: U.S. Army.

Uvedli jsme, že počáteční úlohu umíme vyřešit numericky. Ukázali jsme si základní algoritmus (Eulerův) a řekli, že existují algoritmy pokročilejší. Na tomto místě upozorníme na záludnosti skryté v numerických výpočtech. Je iluzorní se domnívat, že zjemněním kroku při numerickém řešení diferenciální rovnice vždy dostaneme přesnější řešení. Toto platí jenom dokud se nedostaneme ke kritické hodnotě kroku, kdy další snižování vede k tomu, že zpřesnění díky kratšímu kroku se přebije akumulovanou chybou z velkého množství výpočtů nutně zatížených zaokrouhlováním a dalším zjemňováním přesnost ztrácíme.

Zajímavá léčka je v samé podstatě výpočtů na počítači a to v reprezentaci desetinných čísel ve dvojkové soustavě. Například číslo 0.1 je ve dvojkové soustavě periodické! Proto desetinásobným sečtením tohoto čísla ve dvojkové soustavě nedostaneme (překvapivě) jedničku! Je to podobné, jako bychom v námi běžně používané desítkové soustavě třikrát sečetli jednu třetinu v desetinném tvaru reprezentovaném konečným počtem desetinných míst, tj. například třikrát sečetli číslo \(\displaystyle 0.33333333\). Nedostaneme přesně jedničku.

Tento efekt měl i tragický důsledek. Software protiraketového systému Patriot počítal čas postupným přičítáním desetiny sekundy. Protože systém byl vytvořen a testován na mobilním zařízení, které se často restartovalo a běželo krátkou dobu, ničemu to nevadilo. Dlouhodobé nasazení systému Patriot ve válečné operaci na osvobození Kuvajtu však odhalilo závažný nedostatek. Při ostrém nasazení systém běžel dlouho a zaokrouhlovací chyba se kumulovala například 100 hodin. I když za tu dobu chyba dosáhla pouze zlomku sekundy, raketa letící vysokou rychlostí již byla jinde, než systém Patriot propočítal. Dne 25.2.1991 systém Patriot během operace Pouštní bouře na osvobození Kuvajtu od irácké okupace nesestřelil útočící raketu Scud a ta zabila 28 vojáků osvobozující armády a okolo 100 osob zranila.

S chybami plynoucími ze zaokrouhlování se setkáme i při výpočtech mimo modelování diferenciálních rovnic. Viz například Floating-point arithmetic may give inaccurate results in Excel.

12.2. Metoda konečných diferencí

12.2.1. Diskretizace rovnice vedení tepla pomocí konečných diferencí

Rovnice obsahující parciální derivace jsou přirozeným jazykem, kterým modelujeme fyzikální děje. To jsme viděli na rovnici vedení tepla výše a setkáme se s tím i dále. Bohužel tyto rovnice umíme ručně vyřešit jenom v poměrně speciálních případech a ani v těchto případech to není snadná práce. Proto v inženýrské praxi dáváme přednost numerickému řešení rovnice. To je založeno na numerické aproximaci derivací a převádí řešení rovnic s parciálními derivacemi na řešení lineárních rovnic. Možnosti si naznačíme na příkladu rovnice vedení tepla. Tato ukázka je důležitá pro pochopení, co nám z rovnic vlastně může vyplývat a jaké jsou zhruba požadavky na výpočetní prostředky.

Obr. 12.4 Konečné diference umožňují převést parciální diferenciální rovnici na soustavu lineárních rovnic. Červený rámeček označuje neznámé hodnoty v dalším časovém kroku. U explicitní metody je tato hodnota jediná a je snadné ji vypočítat. U implicitní metody jsou neznámé hodnoty tři, každá z nich figuruje ve třech rovnicích a je nutné řešit soustavu rovnic s tridiagonální maticí.

Po převedení derivací z rovnice vedení tepla

\[\rho c\frac{\partial T}{\partial t}=k \frac{\partial ^2 T}{\partial x^2}\]

na konečné diference bychom dostali

\[\rho c\frac{T(x,t+\Delta t)-T(x,t)}{\Delta t}= k\frac{T(x-\Delta x,t)-2T(x,t)+T(x+\Delta x,t)}{\Delta x^2},\]

kde \(\displaystyle \Delta x\) a \(\displaystyle \Delta t\) jsou intervaly oddělující body a časy, ve kterých aproximujeme teplotu. Odsud

\[T(x,t+\Delta t)=T(x,t)+\frac{k\Delta t}{\rho c (\Delta x)^2}\Bigl[T(x-\Delta x,t)-2T(x,t)+T(x+\Delta x,t)\Bigr]\]

a teplotu \(\displaystyle T(x,t+\Delta t)\) v následujícím časovém okamžiku v libovolném bodě \(\displaystyle x\) dokážeme vypočítat ze současné teploty v tomto bodě a z teploty v sousedních bodech \(\displaystyle x+\Delta x\) a \(\displaystyle x-\Delta x\). Konkrétní tvar vzorce není v této chvíli až tak důležitý, podstatné je, že teplotu v dalším časovém okamžiku dokážeme vypočítat z teplot v současném čase. Proto se tato metoda nazývá explicitní metoda.

Explicitní metodu je snadné implementovat programovým kódem. Pokud teploty v čase \(\displaystyle t\) uspořádáme do sloupcového vektoru \(\displaystyle \vec T(t)\), je dokonce možno předchozí vztah zapsat pro všechny body současně jedinou maticovou rovnicí

\[\vec T(t+\Delta t)=\vec T(t)+\frac{k \Delta t}{\rho c (\Delta x)^2} A \vec T(t),\]

kde \(\displaystyle A\) je matice, která má v hlavní diagonále čísla \(\displaystyle -2\), podél diagonály má čísla \(\displaystyle 1\) a jinak nuly s výjimkou prvního a posledního řádku, které jsou nulové. Viz výsledný kód, kde je jenom jeden cyklus pro posun v čase a namísto cyklu přes všechny body v tyči je zde maticové násobení.

Ještě existuje metoda implicitní metoda řešení založená na zpětné diferenci v čase namísto dopředné, tj.

\[\frac{\partial T(x,t)}{\partial t}=\frac{T(x,t)-T(x,t-\Delta t)}{\Delta t}\]

a odsud

\[T(x,t) = T(x,t-\Delta t) +\frac{k\Delta t}{\rho c (\Delta x)^2}\Bigl[T(x-\Delta x,t)-2T(x,t)+T(x+\Delta x,t)\Bigr].\]

Toto vztah umožňující výpočet teplot v čase \(\displaystyle t\) z teplot v čase \(\displaystyle t-\Delta t\). Programová realizace je založena na řešení rovnice a může vypadat následovně. Obsahuje pro každý časový krok řešení soustavy lineárních rovnic s tridiagonální maticí (nenulová čísla jsou na hlavní diagonále a vedle ní).

Rozdíl mezi implicitní a explicitní metodou je v tom, že u explicitní metody máme v každé rovnici jednu neznámou a tuto neznámou je snadné určit. Formálně metoda vede na soustavu rovnic, ale řešení této soustavy je triviální. Oproti tomu u implicitní metody máme v každém vztahu tři neznámé a řešení soustavy rovnic je komplikovanější. Zdá se tedy, že explicitní metoda je výhodnější. Bohužel v praxi explicitní metoda vyžaduje dostatečně jemný časový krok, což může vést k nutnosti použít velmi jemnou časovou diskretizaci a tato skutečnost navyšuje výpočetní náročnost jak z hlediska času, tak i z hlediska paměťových nároků. Implicitní metoda je sice komplikovanější na výpočet, ale dovoluje použít větší časové kroky a v praxi se ukazuje jako výhodnější. Často se pro zvýšení přesnosti používá i kombinace obou metod, kdy je v diferenčním schematu použita dopředná i zpětná diference pro derivaci podle času a obě jsou vhodně zkombinovány (Crank-Nicolsonova metoda).

12.3. Variační metody

12.3.1. Slabá formulace parciálních diferenciálních rovnic

12.3.2. Souvislost řešení rovnice a minima funkcionálu

12.3.3. Metoda konečných prvků