9. Lineární diferenciální rovnice#
Anotace.
Ve všech úvodních příkladech se budeme věnovat asi nejjednodušší diferenciální rovnici: lineární s konstantními koeficienty. Tato tovnice má jediné konstantní řešení, kdy se nic neděje. Nekonstantná řešení jsou popsány exponenciálními funkcemi.
Ve dvou úvodních příkladech se podíváme na řešení rovnice popisující radioaktivní rozpad a exponenciální pokles aktivity. V jednom případě je pokles k nule (radioaktivní prvek se časem rozpadne) a ve druhém případě k nenulové hodnotě (množství radioaktivního prvku je doplňováno difuzí a nastolí se rovnováha).
V dalších dvou příkladech si ukážeme modely vedoucí k rovnici popisující exponenciální růst, resp. exponenciální vzdalování se od konstantního řešení. Ve druhém příkladě si zkusíme numericky demonstrovat intenzitu exponenciálního růstu.
V posledním příkladě se zaměříme na vektorový případ a ukážeme, že zpětná vazba reagující na hodnoty sledované veličiny je vhodnou taktikou pro udržení konstantního stavu této veličina a podle povahy může a nemusí vést k oscilacím. Ukážeme si toto na příkladě regulace teploty ale v podstatě stejně je udržována rovnováha enzymů a proteinů v tělech rostlin a živočichů.
9.1. Radioaktivní rozpad#
Množství radioaktivního materiálu se snižuje rychlostí úměrnou tomuto množství.
Napište vhodný matematický model a vyřešte jej.
Určete dobu, za kterou se množství radioaktivního materiálu sníží na polovinu.
Řešení
Modelem je rovnice
\[\frac{\mathrm dy}{\mathrm dt}=-ky\]s jedním řešením\[y=e^{-kt}\]a obecným řešením\[y=Ce^{-kt}.\]Doba \(\displaystyle T\) za kterou klesne hodnota na polovinu se určí ze vztahu
\[\frac 12=e^{-kT},\]tj.\[kT=\ln 2\]a\[T=\frac 1k \ln 2.\]Tato hodnota se nazývá poločas rozpadu a je \(\displaystyle \ln2\) násobkem převrácené hodnoty rozpadové konstanty \(\displaystyle k\).
9.2. Radioaktivní rozpad a radon v budovách#
Radon je po kouření druhou nejčastější příčinou rakoviny plic.
Obr. 9.7 V případě malého rizika stačí izolace proti radonu jako na obrázku. V oblastech s vyšším rizikem je nutná sofistikovanější ochrana. Zdroj: https://www.estav.cz/cz/5599.ochrana-proti-pronikani-radonu-do-objektu-protiradonova-izolace#
V zemích s přirozenou radioaktivitou je nebezpečím radon. Na rozdíl od ostatních radioaktivních prvků je plynný a může prostupovat do budov, kde se hromadí zejména v nižších patrech, protože má poměrně velkou hustotu.
Předpokládejme, že radon se do budovy dostává konstantní rychlostí.
V budově se rozpadá v souladu s obecným zákonem radioaktivního rozpadu, tj. rychlostí úměrnou množství.
Kromě toho může unikat z budovy difuzí do vnějšího prostředí. Je realistické předpokládat, že rychlost snižování množství radonu tímto efektem je úměrná množství radonu.
Napište matematický model pro množství radonu v budově.
Najděte stacionární (konstantní) řešení modelu.
Najděte obecné řešení modelu. Ukažte, že rozdíl aktuálního stavu od stacionárního řešení exponenciálně klesá.
Obranou proti radonu je izolace. Obranou proti radonu v nezaizolovaných objektech je časté větrání. Které parametry modelu jsou těmito systémy obrany proti radonu ovlivněny a jak? Ukažte, že opravdu vedou ke snížení množství radonu ve stacionárním stavu.
Více viz Radon in the environment.
Řešení
Je-li \(\displaystyle q\) rychlost s jakou radon proniká do budovy, \(\displaystyle k\) rozpadová konstanta a \(\displaystyle d\) konstanta úměrnosti charakterizující intenzitu s jakou plyn uniká difuzí z budovy, je modelem rovnice
\[\frac{\mathrm dy}{\mathrm dt}=q-ky-dy,\]nebo po úpravě\[\frac{\mathrm dy}{\mathrm dt}=q-(k+d)y.\]Stacionární řešení je
\[y=\frac{q}{k+d},\]což zjistíme, pokud do modelu za derivaci dosadíme nulu a vyřešíme vzhledem k \(\displaystyle y\).Rovnice je lineární, což vidíme z tvaru
\[\frac{\mathrm dy}{\mathrm dt}+(k+d)y=q.\]Asociovaná homogenní rovnice\[\frac{\mathrm dy}{\mathrm dt}+(k+d)y=0\]má obecné řešení\[y=Ce^{-(k+d)t}\]a řešením modelu je tedy\[ y = \frac {q}{k+d} + Ce^{-(k+d)t}.\]Rozdíl aktuálního stavu od stacionárního řešení je dán vztahem\[ \frac {q}{k+d} - y = -Ce^{-(k+d)t}\]a exponenciálně klesá. V prakticky smysluplných aplikacích je koncentrace pod stacionárním řešením a roste k němu a pravá strana je tedy kladná. Aby to bylo patrnější, označíme \(\displaystyle K=-C\) a řešení je možno psát jako\[ \frac {q}{k+d} - y = Ke^{-(k+d)t}\]a\[y=\frac {q}{k+d} - Ke^{-(k+d)t}.\]Izolace sníží rychlost s jakou radon prostupuje do budovy a sníží tedy konstantu \(\displaystyle q\). Větrání zvýší rychlost difuze ven a zvýší koeficient \(\displaystyle d\). Každý z těchto efektů nezávisle sníží čitatel resp. zvýší jmenovatel podílu
\[\frac {q}{k+d},\]což vede k menší hodnotě. Hodnota ve stacionárním stavu klesne.
9.3. Časový rozestup mezi trolejbusy#
Obr. 9.8 Zdroj: vlastní#
Uvažujme dva trolejbusy jedoucí za sebou po stejné trati. Označme \(\displaystyle x(t)\) jejich časový odstup. Pokud první trolejbus zastaví na určité zastávce v čase \(\displaystyle t\), druhý trolejbus na tuto zastávku dorazí v čase \(\displaystyle x(t)\). Naším úkolem je zjistit, jak se \(\displaystyle x(t)\) mění s rostoucím \(\displaystyle t\).
Předpokládejme, že
pokud žádní pasažéři nečekají na druhý vůz, druhý vůz se pohybuje rychleji než první vůz a oba vozy se „sjedou“, tj. \(\displaystyle x(t)\) klesá konstantní rychlostí, pokud na druhý vůz nečekají žádní pasažéři
rychlost druhého vozu klesá (a rozestup roste) s rostoucím počtem pasažérů čekajících na zastávce
a počet pasažérů kteří čekají na zastávce roste s rostoucím intervalem mezi oběma vozy.
Navrhněte model pro rozestup trolejbusů, najděte stacionární řešení a posuďte kvalitativní chování nekonstantních řešení.
Podle knihy P. Blanchard, R. L. Devaney, G. R. Hall: Differential equations, Cengage Learning (2006), 828 pp.
Řešení
Situaci je možno modelovat diferenciální rovnicí
Není těžké uhodnout, že v úloze je konstantní řešení. Protože derivace konstanty je nulová, musí toto řešení splňovat vztah
Ve tvaru s lineárním operátorem na levé straně má rovnice tvar
Protože výraz \(\displaystyle e^{\beta t}\) je rostoucí funkcí \(\displaystyle t\), znamená to, že vzdálenost řešení od stacionárního řešení s časem roste. V našem modelu to znamená, že trolejbus, který nabral zpoždění, bude nabírat zpoždění stále větší. Naopak trolejbus, který nedodržel odstup a zkrátil interval od předchozího bude tento předchozí stále více dohánět. To souvisí s tím, že se změnou intervalu se změní počet pasažérů čekajících na zastávce a toto má vliv na rychlost trolejbusu.
Vzhledem k nestabilitě stacionárního řešení nemůžeme nechat řidiče veřejné dopravy jezdit „jak jim to vyjde“. Situace by směřovala k tomu, že cestující budou nejprve dlouho čekat na trolejbus a nakonec přijede několik trolejbusů těsně za sebou.
9.4. Pásová brzda#
Obr. 9.9 Pásová brzda v provedení pro majitele jachet. Zdroj: https://www.sailspar.co.uk/yacht-products/scott-boomlock/#
Obr. 9.10 Nechtěný princip pásové brzdy v tomto případě znesnadňuje ovládání rolety. Zdroj: wikipedia, autor: Laelele#
Lano je částečně obtočeno okolo dřevěného kůlu. Fyzikální rozbor ukazuje, že třecí sílou mezi lanem a kůlem se mění tah \(\displaystyle T\) v laně a tím je způsobeno to, že tah na jednom konci lana se přes kůl sníží a na druhém konci lana pociťujeme tah menší. Je tak možno omotáním okolo kůlu o vhodný úhel vyrovnat s malým úsilím i působení velkých sil. Tah ve vlákně je funkcí úhlu \(\displaystyle \theta\) udávajícího míru obtočení a splňuje diferenciální rovnici
Vyřešte diferenciální rovnici. Ukažte, že brzdné účinky rostou exponenciálně s úhlem. (Brzdné účinky definujeme jako podíl tahů na obou koncích.)
Demonstrujte brzdné účinky na konkrétních číslech. Předpokládejme otočení lana o čtvrtinu otáčky. Tedy táhneme za lano, které je omotáno okolo dřevěného kůlu a zatáčí do pravého úhlu. Předpokládejme, že v tomto uspořádání se tahová síla sníží o 25 procent. Jak se sníží tahová síla, pokud je lano obtočeno o půl otáčky, o celou otáčku, o dvě otáčky?
Předpokládejme, že lano je obtočeno okolo dvou kůlů. Okolo prvního o úhel \(\displaystyle \theta_1\), okolo druhého o úhel \(\displaystyle \theta_2\). Jak se sníží tahová síla v tomto uspořádání? Ukažte, že úhly se sčítají. Toto může být modelová situace nevhodného vedení ovládací šňůry k roletě na obrázku. Optimální a v tomto případě i logické uspořádání je znázorněno modrou čarou.
Poznámky.
Zařízení zkoumané v tomto příkladě se nazývá pásová brzda. Uplatnění najde ve strojírenství, v praktických úlohách každodenního života (zejména života lezců či námořníků), ale například i při posouzení soudržnosti splétaných lan.
Exponenciální růst využívají arboristé u spouštěcího bubnu nebo námořníci u ráhnové brzdy. Na principu více brzd spojených za sebou pracuje slaňovací osma a její zdokonalené varianty. V tomto případě je materiálem kov s nízkým koeficientem tření, které se vyrovná několika otáčkami. Počtem otoček je možné regulovat účinnost.
Řešení
Jedná se o klasickou homogenní lineární diferenciální rovnici s konstantním koeficientem. Jedno řešení je
\[T=e^{\mu\theta}\]a obecné řešení\[T=Ce^{\mu\theta}.\]Pro \(\displaystyle \theta = 0\) dostáváme \(\displaystyle T=C\) a je-li \(\displaystyle T_0\) tah na jenom konci (méně namáhaném), je tah na druhém konci\[T=T_0 e^{\mu\theta}.\]Tah se tedy zvýší multiplikativním faktorem \(\displaystyle e^{\mu\theta}.\) Pokud brzdné účinky definujeme jako podíl tahů na obou koncích, dostáváme\[\frac{T(\theta)}{T_0}=e^{\mu\theta}\]a tento podíl roste exponenciálně s úhlem \(\displaystyle \theta\).Podle zadání platí
\[\frac {1}{0.75}=e^{\mu \frac \pi 2}\]a chceme najít \(\displaystyle e^{\mu \pi }\), \(\displaystyle e^{\mu 2\pi }\) a \(\displaystyle e^{\mu 4\pi }.\) Umocněním získáme\[ e^{\mu \pi } =\left(e^{\mu \frac \pi 2}\right)^2=\frac {1}{0.75^2}=\frac{1}{0.5625}\]a podobně\[ e^{\mu 2\pi } =\left(e^{\mu \frac \pi 2}\right)^4=\frac {1}{0.75^4}=\frac{1}{0.3164}\]a\[ e^{\mu 4\pi } =\left(e^{\mu \frac \pi 2}\right)^8=\frac {1}{0.75^8}=\frac{1}{0.1001}.\]Postupně tedy poklesne tahová síla o 44 procent (půlotáčka), o 70 procent (celá otáčka) a o 90 procent tj. na desetinu (dvě otáčky).Pokud se první brzdy se tah mění násobkem \(\displaystyle e^{\mu\theta_1}\) a u druhé násobkem \(\displaystyle e^{\mu\theta_2}\), vynásobením těchto faktorů s pomocí vzorce pro součin exponenciálních výrazů se stejným základem (mocniny o stejném základu se násobí tak, že se sčítají exponenty) dostáváme
\[e^{\mu\theta_1}e^{\mu\theta_2}=e^{\mu(\theta_1+\theta_2)}.\]Tah se tedy mění s faktorem odpovídajícím součtu úhlů, což bylo dokazováno.
9.5. Regulace stacionárního stavu a jeho robustnost#
Následující příklad je z knihy Uri Alon, An introduction to systems biology, Design principles of biological circuits. Kniha popisuje svět syntéz bílkovin a dalších sloučenin nutných pro život i tak jednoduchých organismů, jako je buňka. Aby byly buňky a živé organismy schopny reagovat na vnější podněty, musí syntéza potřebných bílkovin probíhat včas, rychle a stacionární stavy musí být robustní (nezávislé na změnách parametrů).
Obr. 9.11 V chemických procesech probíhajících v živých organizmech je jistá setrvačnost a pro životaschopnost organismů je nutné udržovat rovnováhu důležitých enzymů a dalších látek. Trošku by to mohlo připomínat vytápění krbem, kde není zapínání a vypínání tak pružné, jako u plynového kotle a proto je potřeba mít pro řízení teploty vhodný plán. V případě buněk takový plán vyvinula evoluce. Zdroj: pixabay.com, autor joseclaudioguima#
Pro jednoduchost si naformulujeme problematiku pro jednoduše představitelný případ regulace teploty v místnosti. Podobné mechanismy se však dlouhodobou evolucí vyvinuly a optimalizovaly i v buňkách a fungují tedy v tělech rostlin a živočichů.
Uvažujme místnost, ve které je umístěno topné zařízení, další zdroje tepla a teplo může unikat z místnosti stěnami.
Příkon zařízení je \(\displaystyle P\) a toto zařízení zvyšuje teplotu místnosti rychlostí úměrnou příkonu. Konstanta úměrnosti je \(\displaystyle a\).
Další zdroje tepla zvyšují teplotu konstantní rychlostí \(\displaystyle S\). V tomto jsou započítáni například lidé.
Místnost se ochlazuje rychlostí úměrnou teplotě \(\displaystyle T\). Konstanta úměrnosti je \(\displaystyle b\). (Popsaného stavu můžeme dosáhnout tak, že teplotní stupnici volíme tak, aby teplota okolí byla nulová.)
Napište model pro teplotu v místnosti a za předpokladu konstantního \(\displaystyle P\) najděte stacionární stav. Všimněte si, že tento stacionární stav závisí na parametrech systému.
Chceme udržovat systém na teplotě \(\displaystyle T_0\). Přidejte k modelu předpoklad, že výkon \(\displaystyle P\) je regulován. Tato regulace je realizována vztahem
\[\frac{\mathrm dP}{\mathrm dt}=-k(T-T_0).\]Okomentujte slovně, v čem regulace spočívá. Poté přidejte tento vztah v rovnici pro teplotu, čímž vznikne lineární autonomní systém. Ukažte, že stacionární teplotou je teplota \(\displaystyle T_0\). Všimněte si, že tato hodnota je robustní: nezávisí na konstantách \(\displaystyle a\), \(\displaystyle b\), \(\displaystyle k\) nebo \(\displaystyle S\).Předpokládejme nejobecnější lineární systém regulace výkonu
\[\frac{\mathrm dP}{\mathrm dt}=\alpha P-k T + \gamma.\]Ukažte, že stacionární stav je robustní (roven \(\displaystyle T_0=\frac\gamma k\) a nezávislý na ostatních parametrech) tehdy a jen tehdy, když platí \(\displaystyle \alpha =0\).Najděte vlastní čísla autonomního systému vytvořeného v tomto modelu a okomentujte, jak se budou chovat řešení. Ukažte, že parametry systému je možné shluknout do jednoho bezrozměrného parametru, jehož hodnota dokáže odlišit kvalitativně různá chování systému. Všimněte si, že situace je skoro stejná jako u soustavy popisující kmity tlumeného oscilátoru v přednášce.
Řešení
Modelem je rovnice
\[\frac{\mathrm dT}{\mathrm dt}=aP+S-bT\]a stacionární řešení pro konstantní výkon \(\displaystyle P\) získáme řešením rovnice\[0=aP+S-bT,\]tj.\[T=\frac{aP+S}{b}.\]Evidentně toto řešení závisí na parametrech systému.Vztah
\[\frac{\mathrm dP}{\mathrm dt}=-k(T-T_0)\]modeluje regulaci výkonu topného zařízení tak, aby se výkon snižoval při teplotě větší než \(\displaystyle T_0\) a zvyšoval při teplotě menší. Rychlost změny výkonu je úměrná odchylce aktuální teploty od teploty \(\displaystyle T_0\). Spojením rovnic získáme systém\[\begin{split}\begin{aligned} \frac{\mathrm dT}{\mathrm dt}&=aP+S-bT\\ \frac{\mathrm dP}{\mathrm dt}&=-k(T-T_0) \end{aligned} \end{split}\]se stacionárním řešením daným soustavou\[\begin{split}\begin{aligned} 0&=aP+S-bT\\ 0&=-k(T-T_0). \end{aligned} \end{split}\]Druhá rovnice okamžitě dává \(\displaystyle T=T_0\) a toto řešení evidentně nezávisí na parametrech jako \(\displaystyle a\), \(\displaystyle b\), \(\displaystyle k\) nebo \(\displaystyle S\).Systém
\[\begin{split}\begin{aligned} \frac{\mathrm dT}{\mathrm dt}&=aP+S-bT\\ \frac{\mathrm dP}{\mathrm dt}&=\alpha P-k T+\gamma \end{aligned} \end{split}\]má stacionární řešení dané soustavou\[\begin{split}\begin{aligned} 0&=aP+S-bT\\ 0&=\alpha P-k T+\gamma . \end{aligned} \end{split}\]Pokud první rovnici vynásobíme faktorem \(\displaystyle -\alpha\), druhou rovnici faktorem \(\displaystyle a\) a sečteme, dostáváme\[0=-\alpha S +\alpha b T - a k T+\gamma a\]a odsud\[T=\frac{\alpha S-a\gamma}{\alpha b-ak}.\]Aby řešení nezáviselo na \(\displaystyle S\), musí být \(\displaystyle \alpha =0.\) Potom dostáváme\[T=\frac{\gamma}{k}.\]V označení podle předchozího bodu máme \(\displaystyle \gamma=kT_0\), tj.\[T=T_0\]a řešení nezávisí na parametrech \(\displaystyle a\), \(\displaystyle b\), \(\displaystyle k\) nebo \(\displaystyle S\). Tomu říkáme, že stacionární stav je robustní.Model
\[\begin{split}\begin{aligned} \frac{\mathrm dT}{\mathrm dt}&=aP+S-bT\\ \frac{\mathrm dP}{\mathrm dt}&=-k(T-T_0) \end{aligned}\end{split}\]má maticový zápis\[\begin{split}\frac{\mathrm d}{\mathrm dt}\begin{pmatrix}T\\P\end{pmatrix}=\begin{pmatrix}-b&a\\-k&0\end{pmatrix}\begin{pmatrix}T\\P\end{pmatrix}+\begin{pmatrix}S\\kT_0\end{pmatrix}.\end{split}\]Charakteristický polynom v proměnné \(\displaystyle \lambda\) je\[\begin{split}0=\begin{vmatrix}-b-\lambda&a\\-k&-\lambda\end{vmatrix}=\lambda^2+b\lambda+ak.\end{split}\]a kořeny jsou\[\lambda_{1,2}=\frac{-b\pm\sqrt{b^2-4ak}}{2}.\]Pokud je \(\displaystyle b^2-4ak>0\), jsou oba kořeny reálné záporné a ve stacionárním bodě je stabilní uzel. Systém přejde do stacionárního bodu bez oscilací. Pokud je \(\displaystyle b^2-4ak<0\), jsou oba kořeny komplexně sdružené se zápornou reálnou částí a systém přejde do stabilního stavu oscilacemi se snižující se amplitudou. Z podmínek vidíme, že o chování rozhoduje to, zda je výraz \(\displaystyle \frac{b^2}{ak}\) větší či menší v porovnání s hodnotou \(\displaystyle 4\). Pokud větší, oscilace se nevyskytnou. Pokud menší, systém osciluje k rovnováze. Tyto oscilace se tedy objeví, pokud je \(\displaystyle a\) nebo \(\displaystyle k\) velké. Pro systém to znamená, že oscilace se objeví, pokud příkon \(\displaystyle P\) výrazně reaguje na odchylku od rovnovážné teploty (velké \(\displaystyle k\)), nebo teplota výrazně reaguje na změny v příkonu (velké \(\displaystyle a\)). Že se v takových situacích může rovnovážná teplota „přestřelit“ a budeme se ní blížit zase z opačného směru není proto nijak překvapující. Přirozenou analýzou vyplynulo, že kvalitativní chování systému (oscilace versus neoscilace) nezávisí na jednotlivých parametrech samostatně, ale na čísle \(\displaystyle \frac{b^2}{ak}\). Pokud prohodíme pořadí proměnných, má soustava tvar\[\begin{split}\frac{\mathrm d}{\mathrm dt}\begin{pmatrix}P\\T\end{pmatrix}=\begin{pmatrix}0&-k\\a&-b\end{pmatrix}\begin{pmatrix}P\\T\end{pmatrix}+\begin{pmatrix}kT_0\\S\end{pmatrix}.\end{split}\]Až na konstantní člen na pravé straně je systém stejný jako systém pro mechanický oscilátor. Přítomnost konstantního člene způsobí, že stacionární bod není v počátku, ale je posunutý.