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. Metoda konečných prvků

Metoda konečných prvků je jedna z nejběžnějších numerických metod pro řešení diferenciálních rovnic v inženýrských aplikacích. Ukážeme si základní myšlenky této metody na jednoduchém příkladu rovnice

\[-\frac{\mathrm d^2 u}{\mathrm dx^2}=f(x)\tag{E}\]

s okrajovými podmínkami \(\displaystyle u(0)=u(1)=0.\) Tato rovnice může vystupovat například při modelování teplotního pole v tyči, ve které jsou zdroje tepla popsané funkcí \(\displaystyle f(x)\). Z praktického hlediska je snadné rovnici vyřešit přímo integrací, ale pro ilustraci metody konečných prvků je tento příklad vhodný. Metoda je obecná a použitelná i pro složitější rovnice ve vícedimenzionálním prostoru, kde analytické metody selhávají: buď jsou nepoužitelné, nebo vedou na řešení zapsané pomocí komplikovaných speciálních funkcí a nekonečných řad, s nimiž se v praxi špatně pracuje.

12.3.1. Slabá formulace rovnice

Nejprve převedeme rovnici na tzv. slabou formulaci. To je formulace, kde nevystupuje druhá derivace funkce \(\displaystyle u\), ale jenom derivace první. To umožní pracovat například s úlohami vedení tepla, kde se materiálové vlastnosti mění skokem. Kromě toho slabá formulace převádí úlohu na úlohu integrální, což je výhodné pro numerické metody. Na tuto výhodu se zaměříme a proto nebudeme rozebírat aspekty spojené s hladkostí funkcí. Budeme předpokládat, že funkce, se kterými pracujeme, jsou dostatečně hladké.

Použijeme základní poznatky z integrálního počtu a derivaci součinu funkcí \(\displaystyle u'\) a \(\displaystyle v\)

\[(u'v)'=u''v+u'v',\]

která má po integraci na intervalu \(\displaystyle [a,b]\) podobu

\[u'(b)v(b)-u'(a)v(a)=\int_a^b u''v\,\mathrm dx +\int_a^b u'v'\,\mathrm dx.\]

Pokud funkce \(\displaystyle v\) splňuje okrajové podmínky \(\displaystyle v(a)=v(b)=0\), dostaneme

\[-\int_a^b u''v\,\mathrm dx = \int_a^b u'v'\,\mathrm dx.\tag{*}\]

Nyní přistoupíme k převodu rovnice (E) na obecnější formulaci. Rovnici napíšeme pro jednoduchost s derivacemi vyjádřenými čárkami a s vynechanou závislostí na \(\displaystyle x\)

\[-u''=f\]

a tuto rovnici vynásobíme funkcí \(\displaystyle v\), která splňuje okrajové podmínky \(\displaystyle v(0)=v(1)=0\). Výslednou rovnici

\[-u''v=fv\]

integrujeme přes interval \(\displaystyle [0,1]\). Tím dostaneme

\[-\int_0^1 u''v\,\mathrm dx = \int_0^1 fv\,\mathrm dx.\]

Ze vztahu (*) poté plyne, že získanou rovnost můžeme přeformulovat na tvar

\[\int_0^1 u'v'\,\mathrm dx = \int_0^1 fv\,\mathrm dx. \tag{W}\]

Pokud nějaká funkce \(\displaystyle u\) splňuje tuto rovnost pro každou hladkou funkci \(\displaystyle v\) s okrajovými podmínkami \(\displaystyle v(0)=v(1)=0\), říkáme, že \(\displaystyle u\) je slabým řešením rovnice (E) (s uvažovanými okrajovými podmínkami). Rovnice (W) se nazývá slabá formulace (angl. weak form) rovnice (E). Slabá formulace rovnice úzce souvisí s hledáním minima energie nebo obecněji nějaké vhodné veličiny související s úlohou a proto ji řadíme mezi variační metody. Proto se slabá formulace někdy nazývá variační formulace.

12.3.2. Galerkinova metoda

Vyjdeme ze slabé variační formulace (W) a budeme hledat přibližné řešení ve tvaru

\[u(x)=\sum_{i=1}^n u_i \varphi_i(x),\]

kde funkce \(\displaystyle \varphi_i(x)\) jsou zvolené hladké funkce splňující okrajové podmínky \(\displaystyle \varphi_i(0)=\varphi_i(1)=0\) a \(\displaystyle u_i\) jsou neznámé koeficienty, které budeme určovat. To vlastně znamená, že hledáme řešení v podprostoru generovaném funkcemi \(\displaystyle \varphi_i(x)\). Funkce \(\displaystyle \varphi_i(x)\) se proto nazývají bázové funkce. Pro snadné odlišení funkcí a koeficientů už nebudeme vynechávat závislost na \(\displaystyle x\).

Derivace funkce \(\displaystyle u(x)\) je

\[u'(x)=\sum_{j=1}^n u_j \varphi_j'(x).\]

Dosadíme-li toto vyjádření do slabé formulace (W), dostaneme

\[\int_0^1 \Bigl(\sum_{j=1}^n u_j \varphi_j'(x)\Bigr) v'(x)\,\mathrm dx = \int_0^1 f(x)v(x)\,\mathrm dx.\]

Využitím linearity integrálu je možné rovnici přepsat do tvaru

\[\sum_{j=1}^n u_j \int_0^1 \varphi_j'(x) v'(x)\,\mathrm dx = \int_0^1 f(x)v(x)\,\mathrm dx.\]

Galerkinova metoda spočívá v tom, že za funkci \(\displaystyle v(x)\) volíme postupně jednotlivé bázové funkce, tedy \(\displaystyle v(x)=\varphi_i(x)\) pro \(\displaystyle i=1,2,\ldots,n\). Tím dostaneme soustavu rovnic

\[\sum_{j=1}^n u_j \int_0^1 \varphi_j'(x) \varphi_i'(x)\,\mathrm dx = \int_0^1 f(x)\varphi_i(x)\,\mathrm dx, \quad i=1,2,\ldots,n\]

nebo po přeznačení

\[a_{ij}=\int_0^1 \varphi_i'(x) \varphi_j'(x)\,\mathrm dx\]

a

\[b_j=\int_0^1 f(x)\varphi_j(x)\,\mathrm dx\]

ve tvaru

\[\sum_{j=1}^n a_{ij} u_j = b_i, \quad i=1,2,\ldots,n.\]

Nyní už je patrné, že se jedná o soustavu lineárních rovnic a po zavedení matice \(\displaystyle A=(a_{ij})\) a vektoru \(\displaystyle \vec b=(b_i)\) můžeme tuto soustavu psát v maticovém tvaru

\[A\vec u=\vec b.\]

Matice \(\displaystyle A\) se nazývá (z historických důvodů) matice tuhosti a vektor \(\displaystyle \vec b\) se nazývá vektor zatížení.

Obr. 12.5 Bázové funkce pro metodu konečných prvků na intervalu [0,1] s dělením na deset dílčích intervalů. V obrázku jsou vybrány tři bázové funkce a lineární kombinace bázových funkcí aproximující parabolu \(\displaystyle y=2x(1-x)\).

Je účelné volit bázové funkce tak, aby na jednu stranu generovaly dostatečně širokou škálu funkcí, ale také aby byla matice \(\displaystyle A\) vhodná pro numerické řešení — například, aby byla řídká. Vhodnou volbou jsou trojúhelníkové funkce dle připojeného obrázku. Pomocí lineární kombinace těchto funkcí je možné vyjádřit libovolnou po částech lomenou funkci splňující nulové okrajové podmínky. Při takto zvolených funkcích je \(\displaystyle a_{ij}\) nulové, pokud se \(\displaystyle i\) a \(\displaystyle j\) liší více než o jedničku a matice \(\displaystyle A\) má nenulové prvky jenom na hlavní diagonále a na dvou přilehlých diagonálách (je tridiagonální). Výpočet integrálů pro \(\displaystyle a_{ij}\) vede k následující matici.

\[\begin{split} A = \frac 1h\begin{pmatrix} 2 & -1 & 0 & 0 & \cdots & 0 \\ -1 & 2 & -1 & 0 & \cdots & 0 \\ 0 & -1 & 2 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & -1 & 2 \end{pmatrix} \end{split}\]

12.3.3. Metoda konečných prvků (FEM)

Obr. 12.6 Ukázka programu pro modelování pomocí metody konečných prvků (FEM). V tomto případě je modelován tříbodově podepřený nosník zatížený silou uprostřed. Program umožňuje měnit parametry nosníku a zatížení a okamžitě vidět výsledky a identifikovat kritická místa konstrukce.

Výše uvedený postup tvoří jádro metody konečných prvků (angl. finite element method, FEM). Problém je nejprve naformulován obecněji než v původní diferenciální rovnici (slabá variační formulace) a poté je hledáno přibližné řešení v podprostoru generovaném zvolenými bázovými funkcemi. Dosazením tohoto tvaru řešení do slabé formulace a volbou testovacích funkcí je získána soustava lineárních rovnic pro neznámé koeficienty v rozvoji řešení podle bázových funkcí. Tuto soustavu je poté možné vyřešit běžnými numerickými metodami pro řešení soustav lineárních rovnic.

Uvedený postup je možné zobecnit na složitější rovnice a úlohy. Na rozdíl od analytických metod si metoda konečných poradí i s komplikovanějšími geometrickými tvary a nespojitými materiálovými vlastnostmi, kdy na sebe navazují dva odlišné materiály. Je možné metodu použít i pro nelineární rovnice, i když v tomto případě je výpočet numericky náročnější.

Výhodou metody konečných prvků oproti metodě konečných diferencí je větší volnost diskretizaci. Metoda konečných diferencí vyžaduje pravidelnou síť bodů, kdežto metoda konečných prvků umožňuje využít nepravidelnou síť a přizpůsobit hustotu bodů místním potřebám. To je výhodné například při zjemnění diskretizace v místech, kde se očekávají velké změny řešení.

Poznámka.

Metodu podrobněji rozebereme a zařadíme do širšího kontextu numerických metod při studiu numerických metod pro parciální diferenciální rovnice v poslední přednášce semestru.