8. Diferenciální rovnice#
Anotace.
V přednášce se seznámíme s rovnicemi obsahujícími derivace neznámé funkce. Jejich využití je všude tam, kde rychlosti změn veličin jsou dány hodnotami těchto veličin.
Typickým příkladem je radioaktivita, protože množství rozpadlých atomů je dáno množstvím nestabilních atomů. Aplikace potom můžeme najít například při ochraně budov před radioaktivním radonem.
Jiným typickým příkladem jsou populační modely, kdy přírůstek populace je dán počtem jedinců schopných reprodukce a ten zpětně souvisí s velikostí populace. Využití je při návrhu trvale udržitelného hospodaření s přírodními zdroji při lovu.
Technicky významným příkladem je i model tepelná výměny, kdy se rychlost změny teploty při tepelné výměně mění podle intenzity toku a ta se mění s teplotním rozdílem.
Řada diferenciálních rovnic má speciální vlastnosti, které můžeme využít při prozkoumávání řešení. Dokonce můžeme například popsat, jak vypadají všechna řešení, aniž bychom je museli počítat. Některé z těchto taktik se naučíme v přednáškách v dalších týdnech věnovaných lineárním rovnicím (následující přednáška) a autonomním systémům (přednáška následující po přednášce o lineárních rovnicích). Teď to zmiňujeme proto, aby šlo vidět, že v případě diferenciálních rovnic nejsou dovednosti spojené s výpočtem jejich řešení tak důležité, jak jsme zvyklí u jiných druhů rovnic. Proto jsou v následujícím seznamu dovedností až na konci.
Důležité dovednosti, které se naučíme v souvislosti s diferenciálními rovnicemi, jsou zejména * schopnost naformulovat diferenciální rovnici podle slovního popisu mechanismu modelovaného děje, * dovednost posoudit existenci a jednoznačnost řešení, * dovednost snížit transformací počet parametrů rovnice * a až v poslední řadě najít řešení numericky nebo analytickou cestou.
Prerekvizity.
Diferenciální rovnice souvisí s derivacemi. Pro úspěšné rozhodnutí, zda se úloha dá modelovat pomocí diferenciální rovnice nutně potřebujeme spolehlivě znát využití derivace. V podstatě s jistotou všude tam, kde se mluví o rychlostech, ale aplikace jsou i jinde.
Pro nalezení analytického řešení diferenciální rovnice je třeba ovládat integrál funkce jedné proměnné.
8.1. Modely založené na rychlostech (derivacích)#
8.1.1. Tepelná výměna, káva v hrnku#
Obr. 8.1 Tepelná výměna probíhá intenzivněji při velkém rozdílu teplot, https://pixabay.com#
Z fyziky víme, že rychlost tepelné výměny mezi dvěma tělesy je úměrná rozdílu jejich teplot (Newtonův zákon tepelné výměny). Rychlostí tepelné výměny můžeme rozumět například rychlost růstu teploty studeného tělesa v teplém prostředí. (Nebo rychlost poklesu teploty horkého tělesa umístěného v chladnějším prostředí.)
Rychlost růstu teploty \(\displaystyle T\) tělesa s časem je derivace teploty podle času. Pokud potřebujeme pracovat s poklesem, uvažujeme záporně vzatou derivaci. Úměrnost matematicky vyjádříme násobením konstantou a teplotní rozdíl může být například při umístění horkého tělesa o teplotě \(\displaystyle T\) v chladné místnosti o teplotě \(\displaystyle T_0\) vyjádřen rozdílem \(\displaystyle T-T_0\).
Proces tepelné výměny probíhající podle Newtonova zákona je tedy možno modelovat vztahem
\[ \frac{\mathrm dT}{\mathrm dt}=-k(T-T_0). \]K rovnici v ideálním případě dodáváme materiálovou charakteristiku (konstantu úměrnosti \(\displaystyle k\)) a počáteční teplotu. Řešením je funkce udávající závislost teploty na čase. Chceme-li znát teplotu za určitý čas, není nutné provádět pokus a čekat na uplynutí požadované doby. Můžeme teplotu přímo vypočítat.
Někdy může být vhodné nesledovat teplotu \(\displaystyle T\), ale rozdíl oproti okolní teplotě, \(\displaystyle \tau=T-T_0\). Model se potom zjednoduší na
\[ \frac{\mathrm d\tau}{\mathrm dt}=-k\tau,\]tedy na model, kdy rychlost změny je úměrná funkční hodnotě.
8.1.2. Radioaktivní rozpad, radon ve sklepě#
Obr. 8.2 Model úbytku rychlostí úměrnou množství modeluje radioaktivní rozpad. V České republice nás zajímá radioaktivní rozpad například z hlediska nežádoucího hromadění radonu v obytných budovách. Zdroj: pixabay.com, rabedirkwennigsen#
Radioaktivní prvky se rozpadají rychlostí, která je úměrná množství dosud nerozpadnutého materiálu. Rychlost, s jakou se mění množství (a tedy i koncentrace \(\displaystyle y\) v daném vzorku) nerozpadnutého radioaktivního materiálu, je tedy popsána matematickým modelem
\[\frac{\mathrm dy}{\mathrm dt}=-\lambda y,\]kde \(\displaystyle \lambda\) je konstanta úměrnosti. Tato rovnice je přirozeným důsledkem toho, že pro daný nestabilní izotop mají všechny atomy stejnou pravděpodobnost, že u nich dojde k rozpadu a tato pravděpodobnost se s časem nemění.Nejznámější aplikací této rovnice je datování archeologických vzorků pomocí radioaktivního uhlíku \(\displaystyle ^{14}C\). V tomto případě se sleduje vzájemná relace mezi množstvím tohoto nestabilního uhlíku a množstvím stabilního \(\displaystyle ^{12}C\). Počáteční podmínka je známa (předpokládáme stejný poměr zastoupení jako relativně nedávno, před průmyslovou revolucí) a díky tomu můžeme najít funkci udávající, jak s časem klesá zastoupení radioaktivního uhlíku. Obsah radioaktivního i stabilního uhlíku je možné změřit a tím získáme odhad, kolik procent radioaktivního uhlíku se rozpadlo. Řešení počáteční úlohy poté použijeme pro odhad doby, kdy organismus přestal spotřebovávat uhlík z atmosféry, tj. odhad stáří vzorku.
Při pokusu o datování kostí dinosaurů klesne množství radioaktivního uhlíku pod měřitelnou úroveň. Proto se v tomto případě používají látky s delším poločasem rozpadu.
Optikou běžného života je nejzajímavější aplikací této rovnice model rozpadu v radioaktivní řadě uranu. V tomto případě vzniká plynný radioaktivní radon, který se může hromadit ve stavbách a je po kouření druhou nejčastější příčinou vzniku rakoviny plic. Proto je v rizikových lokalitách nutné hromadění radonu eliminovat. Toho dosahujeme buď vhodnými konstrukčními přístupy nebo aktivními zařízeními na lapání a odvětrávání radonu.
8.1.3. Samočištění jezer, kontaminace v jezeře#
Obr. 8.3 Jezero, ve kterém se přirozeně obměňuje znečištěná voda za čistou, se dokáže samo zotavit ze znečištění. Rychlost vyplavování nečistot je úměrná míře znečištění. https://pixabay.com#
Nechť veličina \(\displaystyle y\) udává množství látky, která znečišťuje vodu v jezeře o objemu \(\displaystyle V\).
Předpokládejme, že do jezera přitéká čistá voda a stejnou rychlostí odtéká voda s nečistotami (hladina se nemění, je v ustáleném stavu). Nechť veličina \(\displaystyle r\) udává, jaký objem vody se v jezeře takto vymění za jeden den. Předpokládejme dále (poněkud nerealisticky), že rozdělení znečišťujících částic v jezeře je rovnoměrné.
Úbytek hmotnosti nečistot za časovou jednotku je dán derivací \(\displaystyle \frac{\mathrm dy}{\mathrm dt}\).
Protože koncentrace nečistot v jezeře a v odtékající vodě je \(\displaystyle \frac yV\), je úbytek znečištění možno vyjádřit též ve tvaru \(\displaystyle \frac rVy\). Podíl \(\displaystyle \frac rV\) je pro dané jezero kladná konstanta udávající, jak velká část z celkového množství vody se v jezeře vymění za časovou jednotku. Označíme-li tuto konstantu symbolem \(\displaystyle k\), je proces úbytku nečistot v jezeře popsán vztahem
\[ \frac{\mathrm dy}{\mathrm dt} =-ky. \]Výše uvedený model se nazývá rovnice samočištění jezer, ale tento název je čistě formální. Jedná se vlastně o stejnou rovnici, která popisuje radioaktivní rozpad nebo změnu rozdílu mezi teplotou horkého nápoje a místnosti při chladnutí nápoje.
Stejnou rovnicí je možné popsat nejenom odbourávání nečistot z životního prostředí, ale i odbourávání léků nebo drog z těla. Považujme krevní oběh za jezero a lék nebo drogu za znečišťující látku. V případě, že rychlost odbourávání je úměrná koncentraci (platí pro farmakokinetiku prvního řádu, toto splňuje většina léčiv za běžných koncentrací), řídí se proces odbourávání stejnou diferenciální rovnicí.
8.1.4. Akutní normovolemická hemodiluce aneb o krvinky neradi přicházíme#
Obr. 8.4 Při operaci ztrácí pacient krvinky rychlostí úměrnou koncentraci krvinek. Pokud je tato koncentrace malá, pacient ztratí krvinek málo. Zdroj: https://pixabay.com#
Při chirurgické operaci dochází ke krvácení. Pacient ztrácí krev a s ní i krvinky. Při konstantní intenzitě krvácení to znamená, že pacient ztrácí krvinky rychlostí úměrnou počtu krvinek. Formálně na krvinky v krvi můžeme pohlížet stejně jako na znečištění jezera. Jedná se o stejný proces vyplavování látek obsažených v tekutině, jenom měníme interpretaci veličin.
Pokud očekáváme takový průběh operace, že i po uvedeném poklesu bude pořád množství krvinek nad minimální přípustnou hodnotou, je možné před operací toto množství snížit tím, že se část krve odebere a krev se poté doplní vhodnými roztoky.
Protože pacient bude mít po výše uvedeném zákroku už od začátku operace menší počet krvinek, ztrácí tyto krvinky pomaleji a celkový úbytek během operace je menší. Na konci operace se pacientovi vrátí dříve odebraná krev. Výsledkem je, že po operaci v jeho těle koluje více krvinek, než pokud by byl operován s „původní krví“.
Aby metoda fungovala, je nutné odhadnout ztrátu krve během operace. Modelování pomocí matematických metod dokáže předpovědět, kolik krve odebrat na začátku tak, aby i po plánované době operace zůstaly krevní hodnoty pacienta v bezpečných mezích. Pokud na začátku operace část krve dáme bokem a poté tekutiny doplňujeme fyziologickým roztokem (s tím, že vlastní krev vrátíme po skončení operace), jedná se o stejný proces a stejnou rovnici jako samočištění jezer. Pokud krev doplňujeme během operace z krve dopředu odebrané, dokážeme model samočištění jezer modifikovat pro daný proces.
Metoda akutní normovolemické hemodiluce nachází v současné praxi široké využití v řadě operačních oborů. Poskytuje totiž možnost vyhnout se podání alogenní krevní transfuze a tím eliminovat rizika z ní vyplývající. Současně je tato metoda výrazně finančně levnější a její přínos je tak i ekonomický. (Podle https://zdravi.euro.cz/)
8.1.5. RC obvod a chytré stěny ve dřevostavbách.#
Obr. 8.5 Senzor pro sledování vlhkosti dřeva vyvinutý na ÚNOD LDF MENDELU a zabudovaný do dřevostavby. Zdroj: R. Slávik et. al., A Nondestructive Indirect Approach to Long-Term Wood Moisture Monitoring Based on Electrical Methods (2019)#
Při nabíjení kondenzátoru o kapacitě \(\displaystyle C\) přes odpor o velikosti \(\displaystyle R\) roste napětí na kondenzátoru, tím se mění nabíjecí proud. Tím se mění i rychlost nabíjení. Pomocí zákonů elektrotechniky je možno ukázat, že nabíjecí proud \(\displaystyle i\) kondenzátoru se řídí vztahem
Rovnice je tedy stejná jako rovnice radioaktivního rozpadu a rovnice samočištění jezer. Vhodnou manipulací s parametry součástek je možno měnit koeficient u této rovnice a vhodným spojováním těchto obvodů dokážeme podobně simulovat i složitější rovnice. To bylo základem analogových počítačů, které nepracovaly s čísly, ale s napětími. Tyto počítače sehrály svou roli v době, kdy číslicové počítače byly nedostupné, pomalé a nespolehlivé. Tím byla historická úloha analogových počítačů splněna a již se nepoužívají.
RC obvod jako takový má však důležité místo i dnes. Dokáže například filtrovat signály podle frekvence. Výpočet jeho charakteristiky (tj. vyřešení rovnice) a sledování napětí na kondenzátoru umožní měření elektrického odporu tam, kde není vhodné odpor určovat z proudu a napětí pomocí Ohmova zákona. Typickým příkladem je odpor dřeva a jeho vodivost, tj. převrácená hodnota odporu. Tato veličina se používá k rychlému stanovení vlhkosti dřeva, nebo je možno ji dlouhodobě sledovat pomocí senzorů zabudovaných do dřevostavby.
Ve skutečnosti žádná elektronická součástka nemá ideální vlastnosti, a proto se v obvodu projevují i nežádoucí parazitní charakteristiky. Pokud by toto bylo limitující, je možné obvod nahradit podobně se chovajícím zapojením s operačním zesilovačem (odkazovaná stránka pracuje s rovnicí v integrálním tvaru).
8.1.6. Vývoj populace a její ekologický lov#
Obr. 8.6 Při intenzivním lovu může dojít ke zničení populace https://pixabay.com#
Zkoumejme velikost \(\displaystyle y\) určité populace v prostředí s nosnou kapacitou \(\displaystyle K\).
Realistickým předpokladem dodaným biologickými vědami je, že v prostředí s omezenými úživnými vlastnostmi specifická míra růstu populace (rychlost, s jakou se velikost populace zvětšuje vztažená na jednotkové množství populace) klesá s tím, jak se velikost populace přibližuje k nosné kapacitě, a rychlost růstu populace je modelována funkcí \(\displaystyle ry\left(1-\frac yK\right)\). Podle velkosti koeficientů v této funkci dělíme živočichy na r-stratégy a K-stratégy a toto dělení odráží, jak se snaží druh vyrovnávat se změnami prostředí.
Za uvedených předpokladů je možno vývoj populace popsat modelem
\[\frac{\mathrm dy}{\mathrm dt}=ry\left(1-\frac yK\right),\]který se nazývá logistická rovnice.Pokud lovem snížíme přírůstky populace, můžeme tento proces popsat modelem
\[\frac{\mathrm dy}{\mathrm dt}=ry\left(1-\frac yK\right)-h(y),\]kde \(\displaystyle h(y)\) je intenzita lovu populace o velikosti \(\displaystyle y\). Modelování tohoto procesu umožní nalezení trvale udržitelné strategie lovu.
8.1.7. Lovci meteoritů z ČSSR a ČR#
Obr. 8.7 Tři dosud nalezené meteority Benešov. foto: Pavel Spurný, převzato z https://dvojka.rozhlas.cz/#
Česká republika je na světové špičce ve oblasti propočítávání dráhy meteoritů ze světelné stopy zachycené sítí bolidových kamer. Vědcům z Astronomického ústavu se podařilo
jako prvním na světě najít pozůstatky meteoritu propočítáním jeho dráhy ze snímků zachycených speciálními kamerami a zpětně propočítat, odkud meteorit přiletěl (meteorit Příbram, 1959, první „meteorit s rodokmenem“, tj. s doloženým původem),
jako prvním na světě najít pozůstatky meteoritu 20 let po dopadu použitím analýz, které v době dopadu meteoritu nebyly k dispozici (meteorit Benešov, dopad 1991, nalezen 2011),
propočítat a najít (mimo jiné i na dně jezera!) zbytky meteoritu Čeljabinsk z roku 2013.
Meteority s vystopovaným původem jsou extrémně vzácné (do roku 2000 jenom 5 meteoritů, do roku 2016 pouze 31 meteoritů) a tým založený Zdeňkem Ceplechou a nyní vedený Pavlem Spurným se podílel na výpočtu drah většiny z nich. Použité metody jsou popsány například v článku Ceplecha, Revelle: Fragmentation model of meteoroid motion, mass loss, and radiation in the atmosphere, Meteoritics & Planetary Science 40, Nr 1, 35–54 (2005). Například ztráta rychlosti třením v atmosféře je modelována rovnicí
8.2. Obyčejná diferenciální rovnice prvního řádu#
Obyčejná diferenciální rovnice je rovnice, kde vystupuje neznámá funkce a její derivace. Setkáváme se s ní například všude tam, kde rychlost růstu nebo poklesu veličiny souvisí s její velikostí. Například rychlost změny teploty horkého tělesa je funkcí teploty samotné. Rychlost tepelné výměny mezi dvěma tělesy je totiž úměrná rozdílu jejich teplot (Newtonův zákon). Takto se přirozeně diferenciální rovnice objevují v modelech nejrůznějších dějů a jevů. Podstatu děje, který modelujeme, musí dodat fyzika, biologie nebo jiná aplikovaná věda. To v matematice obsaženo není. Matematika poté poslouží k analýze, jaké jsou pozorovatelné důsledky a tím se ověří, jestli příslušná aplikovaná věda správně vystihuje podstatu modelovaného děje.
Definice (diferenciální rovnice)
Obyčejnou diferenciální rovnicí prvního řádu rozřešenou vzhledem k derivaci (stručněji též diferenciální rovnicí, DR) s neznámou \(\displaystyle y\) rozumíme rovnici tvaru
kde \(\displaystyle \varphi\) je funkce dvou proměnných.
(anglicky ordinary differential equation, ODE)
Další formy zápisu rovnice (8.1) jsou
Příklad. Najděte všechny funkce splňující \(\displaystyle y'=2xy\). (Naučíme se řešit později.)
Diferenciální rovnice udává scénář vývoje systému. K jednoznačnému předpovězení budoucího stavu je ovšem nutno znát nejenom, jaký mechanismus ovlivňuje vývoj systému, ale také stav současný.
Definice (počáteční podmínka, Cauchyova úloha)
Nechť \(\displaystyle x_0\), \(\displaystyle y_0\) jsou reálná čísla. Úloha najít
řešení rovnice
se nazývá počáteční (též Cauchyova) úloha.
Řešení Cauchyovy úlohy nazýváme též partikulárním řešením rovnice. Graf libovolného partikulárního řešení se nazývá integrální křivka.
(anglicky initial condition, IC, initial value problem, IVP)
Příklad. Najděte všechny funkce splňující \(\displaystyle y'=2xy\) a \(\displaystyle y(0)=3\). (Naučíme se řešit později.)
Věta (existence a jednoznačnost řešení Cauchyovy úlohy)
Má-li funkce \(\displaystyle \varphi (x,y)\) ohraničenou parciální derivaci \(\displaystyle \frac{\partial \varphi}{\partial y}\) v okolí počáteční podmínky, potom má počáteční úloha (8.1)-(8.2) právě jedno řešení definované v nějakém okolí počáteční podmínky.
Příklad. Rovnice
má řešení \(\displaystyle y=e^x\), což nahlédneme snadno, protože exponenciální funkce se nemění derivováním. Dosazením je možné ukázat, že má dokonce řešení
kde \(\displaystyle C\) je libovolné číslo.
Příklad. Řešení počáteční úlohy
8.2.1. Obecné a partikulární řešení#
Řešení diferenciální rovnice je nekonečně mnoho. Zpravidla je dokážeme zapsat pomocí jediného vzorce, který obsahuje nějakou (alespoň do jisté míry libovolnou) konstantu \(\displaystyle C\). Takový vzorec se nazývá obecné řešení rovnice. Pokud není zadána počáteční podmínka a mluvíme o partikulárním řešení, máme tím na mysli jednu libovolnou funkci splňující diferenciální rovnici.
Příklad: Obecným řešením diferenciální rovnice
Online řešiče ODE (symbolicky):
8.3. Úvod do problematiky numerického řešení diferenciálních rovnic#
Nejprve si naznačíme možnosti numerického řešení. To vychází z grafické interpretace diferenciální rovnice a odpovídá v podstatě modelování, kdy postupně prodlužujeme řešení od zadané počáteční podmínky dopředu či dozadu v čase. Přitom musíme řešit situaci vždy pro konkrétní numerické hodnoty počáteční podmínky a všech parametrů. Naštěstí se dá vhodnou transformací (resp. vhodnou volbou jednotek) počet parametrů zredukovat a tím se zvýší obecná použitelnost numerického výpočtu.
8.3.1. Geometrická interpretace ODE#
Obr. 8.8 Směrové pole diferenciální rovnice, integrální křivky, isokliny#
Protože derivace funkce v bodě udává směrnici tečny ke grafu funkce v tomto bodě, lze rovnici
Počáteční podmínka \(\displaystyle y(x_0)=y_0\) geometricky vyjadřuje skutečnost, že graf příslušného řešení prochází v rovině bodem \(\displaystyle [x_0,y_0]\). Má-li tato počáteční úloha jediné řešení, neprochází bodem \(\displaystyle [x_0,y_0]\) žádná další křivka. Má-li každá počáteční úloha jediné řešení (což bude pro nás velice častý případ), znamená to, že integrální křivky se nikde neprotínají.
Křivky s konstantní hodnotou \(\displaystyle \varphi(x,y)\) mají tu vlastnost, že je všechna řešení protínají pod stejným úhlem, měřeným od kladné části osy \(\displaystyle x\). Například v bodech kde platí \(\displaystyle \varphi(x,y)=0\) míří všechny integrální křivky vodorovně. Proto se křivky, kde je \(\displaystyle \varphi(x,y)\) konstantní, nazývají izokliny.
8.3.2. Numerické řešení IVP#
Obr. 8.9 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. 8.10 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.
Vyjdeme-li z počáteční úlohy
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 následujícím textu.
Online řešiče ODE (numericky):
8.3.3. Transformace diferenciální rovnice#
Obr. 8.11 Letecký snímek údolí Vajont krátce po katastrofě. Video ukazuje, že při modelování procesu ve zmenšeném měřítku je nutné transformovat ostatní veličiny, například čas. Pro nás klíčová slova v čase 3706 dokumentu jsou „tým techniků odhaduje nejvyšší možnou reálnou rychlost sesuvu půdy na jednu minutu, kterou pro simulaci přepočítají na čtyři sekundy“. Čas ve zmenšeném modelu ubíhá jinou rychlostí než čas v reálném ději. Foto: Wikipedia.#
Naučíme se vyjadřovat diferenciální rovnici v jiných proměnných tak, aby bylo možné snížit počet parametrů v této rovnici. Pro jednoduchost budeme uvažovat jenom případ, kdy nová proměnná je lineární funkcí původní proměnné.
Uvažujme funkci \(\displaystyle y\) proměnné \(\displaystyle x\). Připomeneme si vzorce pro derivaci součtu, derivaci konstantního násobku a derivaci složené funkce, ale uvedeme si je v kontextu vhodném pro studium diferenciálních rovnic.
Z derivace součtu a z derivace konstanty plyne pro funkci \(\displaystyle y\) a konstantu \(\displaystyle y_0\) vztah
\[ \frac{\mathrm d (y\pm y_0)}{\mathrm dx} = \frac{\mathrm d y}{\mathrm dx} \pm \frac{\mathrm d y_0}{\mathrm dx} = \frac{\mathrm d y}{\mathrm dx} \pm 0= \frac{\mathrm d y}{\mathrm dx}.\]Z derivace konstantního násobku funkce plyne pro funkci \(\displaystyle y\) a konstantu \(\displaystyle k\) vztah
\[ \frac{\mathrm d (ky)}{\mathrm dx} = k\frac{\mathrm d y}{\mathrm dx}.\]Z derivace složené funkce plyne pro konstantu \(\displaystyle k\) a veličinu \(\displaystyle X = kx\) vztah
\[ \frac{\mathrm d y}{\mathrm d x} = \frac{\mathrm d y}{\mathrm dX} \frac{\mathrm d X}{\mathrm dx} = \frac{\mathrm d y}{\mathrm d X} k \]tj.\[ \frac{\mathrm d y}{\mathrm d (kx)} = \frac{\mathrm d y}{\mathrm d X} = \frac 1k \frac{\mathrm d y}{\mathrm d x}.\]
Výše uvedené výpočty je možno shrnout do pravidla v následující poznámce.
Poznámka (transformace diferenciální rovnice do jiných jednotek)
Pro \(\displaystyle Y=k_1(y-y_0)\) a \(\displaystyle X=k_2 x\) platí
Navíc vzorec z poznámky silně připomíná klasické počítání se zlomky. Proto máme Leibnizův tvar zápisu derivací \(\displaystyle \frac{\mathrm dy}{\mathrm dx}\) při studiu diferenciálních rovnic více v oblibě, než zápis Lagrangeův, \(\displaystyle y'\).
Příklad. Diferenciální rovnice tepelné výměny
obsahuje tři parametry: teplotu okolního prostředí \(\displaystyle T_\infty\), počáteční teplotu \(\displaystyle T_0\) a konstantu \(\displaystyle k\) související s fyzikálními vlastnostmi prostředí. Postupně můžeme posunout teplotní stupnici tak, aby teplota okolí byla nula a počáteční teplota jedna, tj. hodnotu \(\displaystyle T\) snížíme o \(\displaystyle T_\infty\) a upravíme dílek stupnice \(\displaystyle (T_0-T_\infty)\)-krát
Nová rovnice (**) neobsahuje žádné parametry a proto je pro studium jednodušší. Přesto je v ní obsažena veškerá informace obsažená v rovnici (8.5). Tuto informaci je však nutno interpretovat v kontextu definice nových proměnných. Například to, že všechna řešení rovnice (8.5) konvergují k nule znamená, že všechna řešení rovnice (8.5) konvergují k \(\displaystyle T_0\). To, že řešení rovnice (8.6) klesne na poloviční hodnotu za čas \(\displaystyle \ln 2\) znamená, že vzdálenost řešení rovnice (8.5) od rovnovážného stavu se na polovinu zmenší za čas \(\displaystyle \frac 1k \ln 2\).
Poznámka (nondimenzinalizace, rozměrová analýza)
Proces eliminace parametrů z modelu popsaného diferenciální rovnicí se nazývá nondimenzionalizace nebo rozměrová analýza modelu, protože eliminaci parametrů je vhodné provádět tak, aby výsledné nové veličiny vycházely bez fyzikálních jednotek. K tomu se provádí rozbor jednotek jednotlivých veličin. V jednoduchých případech však stačí primitivní postup popsaný v odstavcích výše a ukázaný na příkladu. V tomto příkladě veličina \(\displaystyle x\) nemá fyzikální jednotku, protože je součinem konstanty \(\displaystyle k\) (s jednotkou \(\displaystyle \mathrm s^{-1}\)) a času \(\displaystyle t\) (s jednotkou \(\displaystyle \mathrm s\)). Je možné ji považovat za bezrozměrný čas. Veličina \(\displaystyle y\) také nemá fyzikální jednotku, protože je podílem dvou teplot a je možné ji považovat za bezrozměrnou teplotu.
V této úloze bylo zavedení nových veličin přirozené. I u méně zřejmých úloh zkušenosti ukazují, že je vhodné volit transformaci tak, aby vznikly veličiny bezrozměrné, které nemají fyzikální jednotku. Například v Horáček, Fyzikální a mechanické vlastnosti dřeva I je zavedena bezrozměrná vlhkost, bezrozměrný čas a bezrozměrná vzdálenost na straně 61 pro rovnici popisující difuzi a charakteristická délka, Biotovo číslo (bezrozměrná tepelná vodivost) a bezrozměrná teplota, bezrozměrný čas a bezrozměrná vzdálenost pro rovnici popisující vedení tepla na stranách 88 a 89.
8.4. ODE tvaru \(\displaystyle \frac{\mathrm dy}{\mathrm dx}=f(x)g(y)\) (rovnice se separovanými proměnnými)#
Najít řešení obecné diferenciální rovnice je nemožné, ani však takové ambice mít nemusíme. V praxi se setkáváme s poměrně speciálními druhy diferenciálních rovnic a pro ně jsou metody řešení k dispozici. Jeden takový jednoduše řešitelný druh diferenciální rovnice je představen v následujícím textu.
Definice (ODE se separovanými proměnnými)
Diferenciální rovnice tvaru
Příklad: Rovnice
8.5. Řešení ODE se separovanými proměnnými#
Má-li algebraická rovnice \(\displaystyle g(y)=0\) řešení \(\displaystyle k_1\), \(\displaystyle k_2\), …, \(\displaystyle k_n\), jsou konstantní funkce \(\displaystyle y\equiv k_1\), \(\displaystyle y\equiv k_2\), …, \(\displaystyle y\equiv k_n\) řešeními rovnice.
Pracujme na intervalech, kde \(\displaystyle g(y)\neq 0\) a odseparujeme proměnné.
\[ \frac{\mathrm{d}y}{g(y)}=f(x)\mathrm{d}x\]Získanou rovnost integrujeme. Tím získáme obecné řešení v implicitním tvaru.
\[ \int \frac{\mathrm{d}y}{g(y)}=\int f(x)\mathrm{d}x+C\]Pokud je zadána počáteční podmínka, je možné ji na tomto místě dosadit do obecného řešení a určit hodnotu konstanty \(\displaystyle C\). Tuto hodnotu poté dosadíme zpět do obecného řešení a obdržíme řešení partikulární.
Pokud je to možné, převedeme řešení (obecné nebo partikulární) do explicitního tvaru (vyjádříme odsud \(\displaystyle y\)).
Poslední krok (převod do explicitního tvaru) je volitelný, zpravidla záleží na tom, co dalšího hodláme s řešením dělat. Pro většinu výpočtů je však explicitní tvar vhodnější než tvar implicitní, a proto se o něj vždy snažíme.
Poznámka (zápis partikulárního řešení pomocí určitého integrálu)
V případě počáteční podmínky \(\displaystyle y(x_0) = y_0\) je možné spojit třetí a čtvrtý krok a použít určitý integrál
Počáteční úloha má jediné řešení, pokud má pravá strana ohraničenou parciální derivace podle \(\displaystyle y\), jak je zmíněno v úvodu přednášky. Nicméně pro diferenciální rovnici se separovanými proměnnými je možné vyslovit následující mnohem jednodušší postačující podmínku pro jednoznačnost řešení.
Věta (existence a jednoznačnost řešení Cauchyovy úlohy pro rovnici se separovanými proměnnými)
Je-li \(\displaystyle g(y_0)\neq 0\), má počáteční úloha
Poznámka (existence a jednoznačnost konstantního řešení)
Je-li \(\displaystyle g(y_0)=0\), potom má počáteční úloha
8.6. Redukce parciální diferenciální rovnice na obyčejnou#
V předchozích týdnech jsme se seznámili s modely založenými na parciálních derivacích, zejména s difuzní rovnicí. V případě, kdy hledaná stavová veličina je funkcí jenom jedné proměnné, se parciální derivace redukují na obyčejné derivace a můžeme takové modely řešit v rámci obyčejných diferenciálních rovnic.
8.6.1. Jednorozměrný případ#
Ukážeme si , že parciální diferenciální rovnice popisující tok tepla nebo tok podzemní vody se ve speciálních případech redukují na diferenciální rovnice, jaké jsme se právě naučili řešit.
Uvažujme tok tepla stěnou o tloušťce \(\displaystyle d\), která odděluje dvě prostředí o teplotách \(\displaystyle T_1\) a \(\displaystyle T_2\).
Stacionární tok tepla v jedné dimenzi je dán rovnicí
8.6.1.1. Lineární materiálové vztahy, tj. konstantní materiálová charakteristika#
Je-li \(\displaystyle k\) konstantní, dostáváme
\[\frac{\mathrm dT}{\mathrm dx}=\frac {C_1}k\]a integrací dostáváme\[T=\frac {C_1}k x + C_2.\]Konstanty \(\displaystyle C_1\) a \(\displaystyle C_2\) určíme z podmínek na teplotu na jednotlivých stranách stěny. Vidíme, že teplota ve stěně klesá lineárně.Stejná rovnice a stejné řešení vychází i pro piezometrickou hladinu při rovinném ustáleném proudění podzemní vody v případě, že materiálová charakteristika je konstantní, tj. při proudění s napjatou hladinou (podzemní kolektor s nepropustným stropem a pod tlakem).
8.6.1.2. Nelineární materiálové vztahy, tj. nekonstantní materiálová charakteristika#
Zopakujme předchozí výpočet pro materiál s nelineární materiálovou odezvou, kdy Fourierův (Darcyho v případě podzemní vody) zákon není lineární, tj. \(\displaystyle k\) závisí na teplotě. Nejjednodušší zobecnění je případ, kdy \(\displaystyle k(T)\) je lineární, tj. platí
\[k(T)=aT+b.\]Poté má rovnice tvar\[({aT+b})\frac{\mathrm dT}{\mathrm dx}={C_1}\]a po separaci proměnných dostáváme\[({aT+b}){\mathrm dT}={C_1}{\mathrm dx}\]a\[\frac 12 aT^2+bT=C_1x+C_2.\]Teplotní profil není lineární, ale parabolický s parabolou otočenou naležato. Kterou polovinu paraboly vybrat poznáme z toho, že teplota uvnitř stěny je mezi teplotami na okrajích.Stejný výpočet pro \(\displaystyle b=0\) odpovídá proudění podzemní vody s volnou hladinou. Toto je jiným způsobem (přímé odvození rovnice z Darcyho zákona) odvozeno v textu Dana Říhová a Jana Marková, Poznámky k přednáškám z Hydrauliky, přednáška č. 9. Hladina podzemní vody tedy klesá jako ležatá parabola.
8.6.2. Dvourozměrný radiálně symetrický případ#
Jiný případ, kdy je možno redukovat složitost problému na jednu dimenzi, je stacionární děj v rovině, kdy je situace radiálně symetrická. K tomu je nutno transformovat divergenci a gradient do polárních souřadnic. Příslušné vzorce nebudeme odvozovat, dodá je Wikipedie.
Obr. 8.12 Radiální proudění směrem k čerpanému vrtu. Zdroj: http://ecoursesonline.iasri.res.in.#
Uvažujme například horkou trubku ochlazovanou zvenčí a proudění tepla radiálně směrem od středu. Teplota \(\displaystyle T\) je funkcí vzdálenosti \(\displaystyle r\) od středu a po transformaci gradientu a divergence do polárních souřadnic se stacionární bezzdrojová rovnice vedení tepla
\[0=\nabla\cdot(k\nabla T)\]redukuje na\[\frac 1r \frac{\partial}{\partial r}\left(kr\frac{\partial T}{\partial r}\right)=0.\]Parciální derivace se opět redukují pro funkci jedné proměnné na obyčejné derivace a stejně jako v předchozím případě můžeme integrovat na\[kr\frac{\mathrm d T}{\mathrm d r}=C_1.\]Odsud\[k\mathrm d T=\frac{C_1}r {\mathrm d r}\]a\[kT={C_1}\ln(r)+C_2.\]Konstanty \(\displaystyle C_1\) a \(\displaystyle C_2\) se určí z teplot na vnitřním a vnějším povrchu trubky.Stejný vzorec platí pro analogické radiální proudění podzemní vody při proudění s napjatou hladinou. Toho se využívá při čerpacích zkouškách nebo při umělém snižování hladiny spodní vody. Po dosazení relevantních veličin a výpočtu konstant se odvozený vzorec uvádí ve tvaru
\[h-h_0=\frac{Q}{2\pi T}\ln \frac r{r_0}\]a nazývá Thiemova rovnice.Předchozí postup můžeme modifikovat i pro radiální proudění s volnou hladinou, tj. proudění modelované rovnicí
\[\nabla \cdot(K\nabla h)=0,\]kde \(\displaystyle K=kh\) je materiálová konstanta pro proudění s volnou hladinou. Jako v předchozím případě přejdeme do proměnné \(\displaystyle r\) a dostáváme\[khr\frac{\mathrm d h}{\mathrm d r}=C_1.\]Odsud\[h{\mathrm d h}=\frac{C_1}{kr}\mathrm dr\]a\[\frac 12 h^2=\frac{C_1}{k}\ln(r)+C_2.\]Zpravidla se tato rovnice používá pro stanovení vydatnosti čerpané studny a konstanty \(\displaystyle C_1\) a \(\displaystyle C_2\) určíme z výšky hladiny ve studni a z výšky hladiny v kontrolním vrtu nedaleko studny. Tento vztah umožňuje například navrhnout průměr studny, odhadnout vydatnost studny, nebo pomocí odčerpávaného vrtu a menších pomocných vrtů sledujících pokles hladiny v okolí odčerpávaného vrtu stanovit filtrační součinitel \(\displaystyle k\). Využití vzorce je však mnohem rozmanitější, umožňuje vypočítat poměry ve stavebních jámách a v jejich okolí. To je užitečné například při odhadu, kolik vody se hromadí ve výkopu. Další využití je, že dokážeme odhadnout vliv stavební jámy na hydrologické poměry v okolí a tyto poměry dokážeme měnit a přizpůsobovat našim potřebám. Častou aplikací je také hydraulická clona (soustava prvků rozmístěných a provozovaných tak, aby nedocházelo k šíření kontaminace z chemické výroby do vodárensky využívaných vod). V tomto případě je však situace komplikovanější, protože je nutné zkombinovat dostředivé proudění k čerpanému vrtu s rovinným prouděním podzemní vody. Toto se naučíme v příští přednášce a využijeme linearitu.
8.7. Diferenciální rovnice růstu vodní kapky#
Obr. 8.13 Londýnská mlha. Dnes už to není jako za časů Sherloka Holmese. Poslední velká mlha (Pea soup fog) byla v roce 1952. Zdroj: Wikipedia.#
Modelujme růst kulové kapky. Ta roste tak, že na povrchu kondenzují vodní páry. Kapka proto roste tak, že její objem se zvětšuje rychlostí úměrnou povrchu. Povrch je zase úměrný druhé mocnině poloměru a poloměr je úměrný třetí odmocnině objemu. Platí tedy (po sloučení všech konstant úměrnosti do jedné)
Tato rovnice má konstantní řešení \(\displaystyle V=0\). Nekonstantní řešení dostaneme po úpravě
Všimněte si, že počáteční úloha s počáteční podmínkou \(\displaystyle V(0)=0\) má konstantní nulové řešení