Autonomní rovnice a systémy¶

Tento zápisník obsahuje výpočty k cvičení z Aplikované matematiky věnované autonomním systémům.

Příklady s velikostmi populací a se syntézou kolagenu jsou převzaty nebo silně inspirovány knihou Alan Garfinkel, Jane Shevtsov, Yina Guo: Modeling Life, The Mathematics of Biological Systems.
Příklad se skladováním stavebního materiálu je inspirovaný okolím vlakového nádraží, ze kterého jezdím na fakultu.

import matplotlib.pyplot as plt   # knihovna pro statické grafy, grafy budeme kreslit sem do zápisníku
%matplotlib inline               

import numpy                      # knihovna pro numerické výpočty
from scipy.integrate import solve_ivp  # řešení diferenciálních rovnic

import bokeh.io                   # knihovna pro kreslení interaktivních grafů
import bokeh.plotting             # knihovna pro kreslení interaktivních grafů
bokeh.io.output_notebook()        # grafy kreslit sem do zápisníku

colors = 10*bokeh.palettes.Greens5  # paleta barev

BokehJS 2.3.0 successfully loaded.

Skladování stavebního recyklátu¶

Hromada sypkého materiálu má tvar kužele. Úhel u vrcholu je konstantní, daný mechanickými vlastnostmi materiálu a je nezávislý na objemu.

Předpokládejme, že personál stavebnin přisypává na hromadu materiál konstantní rychlostí (v jednotkách objemu za jednotku času).
Tato hromada je však v poměrně otevřené krajině a vítr rozfoukává materiál po okolí. Je rozumné předpokládat, že rozfoukávání (opět v jednotkách objemu za jednotku času) se děje rychlostí úměrnou povrchu návětrné strany pláště.

Diferenciální rovnice popisující růst hromady je $\frac{\mathrm dV}{\mathrm dt} = R - k_0S$ a po vyjádření pomocí objemu $\frac{\mathrm dV}{\mathrm dt} = R - k V^{\frac 23},$ kde $R$ a $k$ jsou konstanty. Vhodnou volbou jednotky objemu můžeme docílit toho, že jsou tyto konstanty numericky stejné a vhodnou volbou jednotky času poté toto, že tato společná hodnota je rovna jedné. Budeme tedy modelovat rovnici pro $R=k=1$ . Na obrázku vidíme jediné stabilní řešení, které je globálně atraktivní.

tspan = numpy.linspace(0, 10, 1000)            # Definice časových značek
p = bokeh.plotting.figure(title="Skladování recyklátu", x_axis_label='čas', y_axis_label='výška hromady',plot_width=1000)

for i in range(0,15):                          # Cyklus přes počátečnmí podmínky
    sol = solve_ivp(lambda t, y: 1-y**(2/3), [tspan[0], tspan[-1]], [i/10], t_eval=tspan)  # Vyřešení rovnice pro danou počáteční podmínku
    p.line(sol.t, sol.y[0], color=colors[i], legend_label="V(0)=%s"%str(i/10), width=3)    # Vykreslení řešení

p.legend.click_policy="hide"                   # Skrývat křivky po kliknutí na legendu
bokeh.plotting.show(p)                         # Vykreslení řešení

Propeptid kolagenu¶

Kolagen je klíčový protein pojivových tkání. Jeden z kroků při syntéze kolagenu spočívá v reakci tří molekul propeptidu kolagenu, zkráceně propeptidu. Tento propeptid se formuje konstantní rychlostí a kromě toho, že je surovinou pro produkci kolagenu, se ještě rozpadá rychlostí úměrnou koncentraci. Modelem je rovnice $\frac{\mathrm dP}{\mathrm dt}=-k_1 P^3 +k_2-k_3 P,$ kde $k_i$ jsou konstanty. Rovnici je možno přepsat do tvaru $\frac{\mathrm dP}{\mathrm dt}=\Bigl(k_2-k_3 P\Bigr)-\Bigl(k_1 P^3\Bigr)$ s rozdílem klesající a rostoucí funkce na pravé straně. Tyto funkce mají jediný průsečík a proto má rovnice jediný stacionární bod. Tento bod je stabilní, protože pro vysoké hodnoty $P$ je pravá strana rovnice $-k_1P^3$ záporná (dominantní člen je ) a pro malé hodnoty je pravá strana $k_2-k_3P$ kladná (v nule je rovna konstantě $k_2$ ).

tspan = numpy.linspace(0, 10, 1000)           # Časy, ve kterých se bude počítat řešení
plt.figure(figsize=(12,6))                    # Inicializace obrázku, rozměry
for podminka in numpy.linspace(0,2,20):       # Cyklus pro počáteční podmínky
    sol = solve_ivp(lambda t, P: -0.5*P**3+1-P, [tspan[0], tspan[-1]], [podminka], t_eval=tspan)
    plt.plot(sol.t, sol.y[0])
plt.title("Časový vývoj množství propeptidu kolagenu")    # 
plt.xlabel('čas')
plt.ylabel('množství propeptidu')
None

      
    

_images/Autonomni_rovnice_a_systemy_5_0.png

Jelen a los¶

Uvažujme populaci jelenů a losů. Tyto populace spolu soupeří o potravu. Bez konkurence by populace jelena rostla rychlostí $3$ a populace losa rychlostí $2$ na jeden kus. Vnitrodruhová konkurence se projevuje v obou populacích stejně a je rovna druhé mocnině příslušné velikosti populace. Mezidruhová konkurence je vyjádřena členem rovným součinu velikosti populací a tato konkurence se projeví s koeficientem $0.5$ v populaci losa a s koeficientem $1$ v populaci jelena.

Sestavte matematický model a otestujte jej numerickým experimentem na stabilitu stacionárních bodů. Poté zdvojnásobte parametry mezidruhové konkurence a sledujte změnu odezvy.

Modelem je následující autonomní systém. $\begin{aligned} \frac{\mathrm dx}{\mathrm dt}&=3x-xy-x^2 \cr \frac{\mathrm dy}{\mathrm dt}&=2y-0.5xy-y^2 \end{aligned}$

def jelen_a_los(t, z, parametry):   # Definice interakce mezi jeleny a losy
    x, y = z
    a, b, c, d, e, f = parametry
    return [a*x - b*x*y - c*x**2, d*y - e*x*y - f*y**2]

sol=[None,None]                                        # pole pro reseni
par = [ [3, 1, 1, 2, 0.5, 1] , [3, 2, 1, 2, 1  , 1] ]  # pole pro nastavení parametrů
# vyřešení rovnice pro jednotlivé hodnoty parametrů, u každého nastavení řešíme dvě počáteční úlohy a budeme kreslit dvě trajektorie
sol[0] = [ solve_ivp(jelen_a_los, [0, 20], pocatek, args=([par[0]]), dense_output=True) for pocatek in [[2,2],[1,0.5]] ]
sol[1] = [ solve_ivp(jelen_a_los, [0, 20], pocatek, args=([par[1]]), dense_output=True) for pocatek in [[2,2],[1,0.5]] ]

      
    

import warnings
warnings.filterwarnings('ignore')              # ignorovat chyby při dělení nulou při dělení délkou vektoru

fig, ax = plt.subplots(4, 2, figsize=(15,15))  # nastavenní velikosti obrázku a počtu podobrázků
t = numpy.linspace(0, 20, 1000)                # nastavení časových značek

colors=['red','orange','green']                # pole pro střídání barev

for j in [0,1]:                                # cyklus pro nastavení parametrů, vždy kreslí do dvou řad
    for i in [0,1]:                            # cyklus přes počáteční podmínky s daným nastavením parametrů
        reseni=sol[j][i].sol(t)                # uložení odpovídajícího řešení
        ax[2*j,0].plot(t,reseni[0], color=colors[i])  # vykreslení časovového průběhu první komponenty řešení, jeleni, obrázek vlevo nahoře
        ax[2*j,1].plot(t,reseni[1], color=colors[i])  # vykreslení časovového průběhu druhé komponenty řešení, losi, obrázek vpravo nahoře
        ax[2*j+1,1].plot(reseni[0],reseni[1], color=colors[i], linewidth=4) # vykreslení trajektorie do obou obrázků v dolní řadě čtveřice
        ax[2*j+1,0].plot(reseni[0],reseni[1], color=colors[i], linewidth=4) 

    y, x = numpy.mgrid[0:3:20j, 0:3:20j]                       # mřížka pro vykreslení směrového pole a trajektorií
    U,V = jelen_a_los(t,[x,y],par[j])                          # směrové pole
    ax[2*j+1,0].streamplot(x, y, U, V, color='blue')           # trajektorie do obrázku vlevo dole
    U, V =U/numpy.sqrt(U**2+V**2), V/numpy.sqrt(U**2+V**2)     # směrové pole s normalizovanou délkou vektoru
    ax[2*j+1,1].quiver(x, y, U, V, color='blue', lw=.01)       # vykreslení směrového pole

    ax[2*j,0].set_title("Model %s. Populace jelenů"%str(j+1))  # nadpis v grafu nahoře vlevo
    ax[2*j,1].set_title("Model %s. Populace losů"%str(j+1))    # nadpis v grafu nahoře vpravo
    ax[2*j,0].set_ylim([0,None])                               # svislou osu zacinat od nuly v grafu nahoře vlevo
    ax[2*j,1].set_ylim([0,None])                               # svislou osu zacinat od nuly v grafu nahoře vpravo
    ax[2*j+1,0].set_title("Model %s. Trajektorie v prostoru jelenů a losů"%str(j+1))
    ax[2*j+1,1].set_title("Model %s. Fázový portrét v prostoru jelenů a losů"%str(j+1))
None

      
    

_images/Autonomni_rovnice_a_systemy_8_0.png

Puštík obecný¶

Puštík obecný se téměř výhradně živí malými hlodavci. Předpokládejme následující vztahy.

Populace hlodavců má porodnost $0.1$ na jedince a úmrtnost $0.025$ na jedince za jednotku času.
Rychlost s jakou jeden puštík konzumuje hlodavce je úměrná počtu hlodavců s konstantou úměrnosti $0.01$ . Porodnost v populaci puštíka je úměrná množství zkonzumované potravy s konstantou úměrnosti $0.05$ .
Úmrtnost v populaci puštíka je $0.1$ na jedince za jednotku času.

Uvedné vztahy je možné modelovat následující soustavou diferenciálních rovnic. $\begin{aligned} \frac{\mathrm dx}{\mathrm dt}&=0.1 x-0.025x-0.01xy\cr \frac{\mathrm dy}{\mathrm dt}&=-0.1y+0.05xy \end{aligned}$

def pustik_a_hlodavec(t, z):
    x, y = z
    a, b, c, d, e = 0.1, 0.025, 0.01, 0.1, 0.05
    return [a*x - b*x - c*x*y, -d*y + e*x*y]

xmax, ymax=5, 20

y, x = numpy.mgrid[0:ymax:20j, 0:xmax:20j]
U,V = pustik_a_hlodavec(t,[x,y])

pocatek = [2,10]
sol = solve_ivp(pustik_a_hlodavec, [0, 200], pocatek, dense_output=True)  
t = numpy.linspace(0, 200, 1000)
reseni=sol.sol(t)

      
    

fig = plt.figure(figsize=(10,10))
plt.streamplot(x, y, U, V, color='blue')
plt.plot(reseni[0],reseni[1], color='red', linewidth=2)
plt.title("Trajektorie v prostoru hrabošů a puštíků")
plt.xlabel('Populace hraboše')
plt.ylabel('Populace puštíka')
None

      
    

_images/Autonomni_rovnice_a_systemy_11_0.png

Vidíme oscilatorický průběh. Ten je dán tím, že puštík sám bez sebe bez potravy nepřežije. V období hojnosti puštík má dostatek potravy, jeho populace roste. Tím však spotřebovává stále více hlodavců a od určitého limitu začne hlodavců ubývat. Tento vývoj trvá, dokud není populace puštíka zredukována natolik, že již nebrání hlodavcům v růstu jejich populace. Tím se populace hlodavců namnoží a začne období hojnosti pro puštíka, čímž se uzavře cyklus.

fig = plt.figure(figsize=(10,10))                 # definice obrazku
ax = fig.add_subplot(1, 1, 1)                     # zpristupneni manipulace s osami
plt.plot(t,reseni[1], color='red', linewidth=2)   # nakresleni druhe komponenty reseni (Python indexuje od nuly)
plt.title("Časový vývoj populace puštíka")        # nadpis
ax.set_ylim([0,None])                             # osa y začne na nule
ax.set_xlabel("čas")                              # popisek na ose x
ax.set_ylabel("populace puštíka")                 # popisek na ose y
None  

      
    

_images/Autonomni_rovnice_a_systemy_13_0.png

Kůň Převalského¶

Kůň Převalského je divoký kůň ze střední Asie, jediný druh koně, který nebyl domestikován. V divočině jsou tyto koně loveni vlky. Napište matematický model založený na následujících předpokladech. Porodnost v populaci koní je $0.15$ na jedince. Úmrtnost v populaci koní je $0.01$ na jedince. Vlci se živí i jinou potravou, mají tedy kladnou porodnost. Ta je $0.1$ na jedince. Vlci mají konstantní úmrtnost $0.05$ na jedince. Pravděpodobnost s jakou je kůň uloven vlkem je úměrná počtu vlků s konstantou úměrnosti $0.02$ .

Matematickým modelem je následující dvojice rovnic. $\begin{aligned} \frac{\mathrm dx}{\mathrm dt}&=0.15 x-0.01x-0.05xy\cr \frac{\mathrm dy}{\mathrm dt}&=0.1y-0.05y \end{aligned}$

def kun_a_vlk(t, z):
    x, y = z
    return [0.15*x - 0.01*x - 0.05*x*y, 0.1*y -0.05*y]

xmax, ymax=20, 20
y, x = numpy.mgrid[0:ymax:20j, 0:xmax:20j]
U, V = kun_a_vlk(t,[x,y])

fig = plt.figure(figsize=(10,10))
plt.streamplot(x, y, U, V, color='blue')
plt.title("Trajektorie v prostoru koní a vlků")
plt.xlabel('Populace koní')
plt.ylabel('Populace vlků')
None

      
    

_images/Autonomni_rovnice_a_systemy_15_0.png

Populace vlků stále roste. Populace koní má při malých hodnotách vlků prostor pro růst (v dolní části obrázku trajektorie směřují doprava), ale tento růst je časově omezený a populace koní tedy koexistenci s populací vlků nepřežije, pokud se vzájemná interakce bude řídit uvedenými předpoklady. Podle Wikipedie kůň Převalského přežil jenom díky péči zoologických zahrad a z rodokmenu je zřejmé, že 70 procent jedinců tohoto druhu má původní předky ze zoologické zahrady v Praze.

Systém, vlastní čísla a vlastní směry¶

Úloha ze cvičení. Nakreslíme si do obrázku vlastní vektory a trajektorie v okolí stacionárního bodu. Stacionární bod je nestabilním uzlem. $\begin{aligned} \frac{\mathrm dx}{\mathrm dt}&= 4x^2y+y^3-5\cr \frac{\mathrm dy}{\mathrm dt}&= 3xy^2-3y \end{aligned}$

def system(t, z):
    x, y = z
    return [4*x**2*y+y**3-5, 3*x*y**2-3*y]

def Jacobi(z):
    x, y = z
    return [[8*x*y, 4*x**2 + 3*y**2],[3*y**2, 6*x*y-3]]    

      
    

from numpy import linalg as LA
x0, y0 = 1, 1
J=Jacobi([x0,y0])
vlastni_cisla, vlastni_smery = LA.eig(J)
e1=vlastni_smery[:,0]
e2=vlastni_smery[:,1]

      
    

Do jednoho obrázku nakreslíme trajektorie a vlastní vektory Jacobiho matice. Je vidět, že v okolí stacionárního bodu se trajektorie řídí vypočtenými vlastními směry. Vlastně jenom jedním, protože vlastní hodnoty jsou $10$ a $0.3$ a v exponenciální závislosti jasně dominuje vyšší hodnota. Proto trajektorie “drží” směr vlastního vektoru příslušného vyšší vlastní hodnotě, tj. červený směr.

fig = plt.figure(figsize=(10,10))

xmin, xmax, ymin, ymax = 0.9, 1.1, 0.9, 1.1
y, x = numpy.mgrid[ymin:ymax:20j, xmin:xmax:20j]
U, V = system(t,[x,y])
plt.streamplot(x, y, U, V, color='blue')

krok = 0.04
plt.title("Trajektorie a vlastní směry")
plt.plot([x0-krok*e1[0],x0+krok*e1[0]],[y0-krok*e1[1],y0+krok*e1[1]], color='red', lw=3, label=str(vlastni_cisla[0]))
plt.plot([x0-krok*e2[0],x0+krok*e2[0]],[y0-krok*e2[1],y0+krok*e2[1]], color='green', lw=3, label=str(vlastni_cisla[1]))
plt.legend()

# Ukazka recyklace kodu. Az na nastaveni xmin, xmax, ... je vse stejne jako vyse
# Neni to uplne optimalni postup, pro kod u ktereho se predpoklada kratka zivotnost 
# (kod vznikly behem uceni nebo pro jednorazove vyuziti) to muze byt akceptovatelne

fig = plt.figure(figsize=(10,10))

xmin, xmax, ymin, ymax = 0, 5, 0, 5
y, x = numpy.mgrid[ymin:ymax:20j, xmin:xmax:20j]
U, V = system(t,[x,y])
plt.streamplot(x, y, U, V, color='blue')

krok = 0.04
plt.title("Trajektorie a vlastní směry")
plt.plot([x0-krok*e1[0],x0+krok*e1[0]],[y0-krok*e1[1],y0+krok*e1[1]], color='red', lw=3, label=str(vlastni_cisla[0]))
plt.plot([x0-krok*e2[0],x0+krok*e2[0]],[y0-krok*e2[1],y0+krok*e2[1]], color='green', lw=3, label=str(vlastni_cisla[1]))
plt.legend()


None

      
    

_images/Autonomni_rovnice_a_systemy_21_0.png

_images/Autonomni_rovnice_a_systemy_21_1.png

A-Ro-Ma-Math