7. Parciální derivace#

Co se dozvíte v tomto textu

V této části se naučíme sledovat rychlosti růstu funkce dvou a více proměnných.

Mezi aplikace patří ochrana ohrožených druhů, kdy při snaze podpořit růst populace musíme vybrat ze škály parametrů ovlivňujících vývoj ten, který má na růst populace nejpodstatnější vliv. Postup je možné použít i naopak a při snaze o kontrolu škůdců vytipovat parametr, který vývoj populace škůdce nejvíce postihne.

Další aplikací je model šíření invazních druhů, kdy se zastoupení živočišného druhu mění v čase i v prostoru tím, jak se druh množí a rozšiřuje území, na kterém žije.

Foto: Ondatra se jako invazní druh velmi rychle rozšířila z Dobříšského panství do celé Evropy. Tato invaze byla popsána v jednom z prvním článků věnovaných invazím živočišných druhů a díky tomu se naše země dostala do většiny ekologických učebnic. D. Gordon E. Robertson, https://commons.wikimedia.org/

Pomocí derivace jsme se naučili sledovat rychlost růstu funkce jedné proměnné. Tato veličina je užitečná protože rychlost růstu je často determinována měřitelnými veličinami a je k dispozici zákon vyjadřující souvislost rychlosti růstu sledované veličiny s ostatními parametry systému. Často je mezi těmito parametry i ona sledovaná veličina, což při sestavení modelu vede na diferenciální rovnice.

Někdy není možné si vystačit pouze s funkcemi jedné proměnné. Proto si nyní koncept derivace rozšíříme i na funkce více proměnných.

Navážeme na příklad s růstem populace tuleňů. Je-li \(\lambda\) největší reálné vlastní číslo matice projekce a \(v\) příslušný vlastní vektor, potom velikost populace roste ve všech věkových kategoriích geoemtrickou řadou s kvocientem \(\lambda\) a poměr velikostí populace v jednotlivých věkových kategoriích je stejný, jako poměr velikostí komponenty vlastního vektoru \(v\). Změna hodnot v matici, která by vyjařovala změnu v pravděpodobnosti přežívání nebo změnu v počtu potomků, se zcela jistě projeví i na dynamice růstu. Zde však je nutné si uvědomit, že parametrů v matici je několik a situace je tedy odlišná od úloh, jaké jsme řešili doposud. Dosud pomocí derivací umíme popsat rekci závislé proměnné, závisející na jedné nezávislé proměnné. Touto nezávislou proměnnou byl zpravidla čas. Nyní však je nezávislých proměnných hned několik. Rychlost růstu se v tkaových případech modeluje pomocí parciální derivace.

7.1. Zavedení parciální derivace#

Změna funkce více proměnných může být způsobena změnou libovolné nezávislé proměnné. Pokud sledujeme například rychlost růstu populace popsané Leslieho maticí, může nás zajímat, jak tato rychlost růstu souvisí s jednotlivými prvky projekční matice, tj. s plodností jednotlivých věkových tříd a s jejich pravděpodobností přežití. Zdá se býti rozumné oddělit jednotlivé změny a studovat každou samostatně. Buď v daném bodě fixovat všechny parametry až na jeden a sledovat například vliv pravděpodobnosti přežití u nejmladší věkové kategorie na celkový růst populace.

Následující definice výše uvedenou myšlenku odděleného sledování změny funkce (závislé veličiny) jako reakce na změnu jedné jediné vstupní informace (jedné nezávislé veličiny) uvádí v život. Definice je stejná jako u derivace funkce jedné proměnné, změna je pouze v tom, že je přítomna i další proměnná.

Definice (Parciální derivace funkce dvou proměnných)

Buď \(f\colon \mathbb R^2\to\mathbb R\) funkce dvou proměnných, \(x\) a \(y\), tj. \(f(x,y)\). Výraz

\[\frac{\partial f}{\partial x}:=\lim_{h\to 0}\frac{f(x+h,y)-f(x,y)}h\]
se nazývá parciální derivace funkce \(f\) podle \(x\). Podobně,
\[\frac{\partial f}{\partial y}:=\lim_{h\to 0}\frac{f(x,y+h)-f(x,y)}h\]
je parciální derivace funkce \(f\) podle \(y\).

Podobně můžeme definovat parciální derivaci pro funkce libovolného konečného počtu proměnných. V těchto parciálních derivacích vlastně sledujeme, jak reaguje veličina \(f\) na změny jenom v jedné proměnné. Proměnná, přes kterou se nederivuje, má vlastně roli parametru, nijak se nemění.

Kromě \(\frac{\partial f}{\partial x}\) se běžně parciální derivace označuje symboly \(f'_x\), \(f_x\), \(\partial_x f\) nebo \(D_x f\).

Parciální derivaci \(\frac{\partial f}{\partial x}\) je tedy možné chápat jako rychlost s jakou se mění veličina \(f\) jako funkce proměnné \(x\). Ekvivalentně jako změnu veličiny \(f\) při změně veličiny \(x\) o jednotku. Tato hodnota bude malá, pokud se funkční hodnoty při změnách \(x\) moc nemění. To by znamenalo, že rozkolísání hodnot veličiny \(x\) má malý vliv na veličinu \(f\). V tomto smyslu se jedná o jakousi míru citlivosti (sensitivity) a parciální derivace se používají pro citlivostní analýzu.

Vektor

\[\nabla f=\left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial x}\right)\]
se nazývá gradient funkce \(f\). Gradient míří směrem, ve kterém roste funkce \(f\) nejrychleji a jeho délka je rovna nárůstu funkční hodnoty na jednotce délky. V přírodních vědách často pracujeme se záporně vzatým gradientem, který udává směr a intenzitu maximálního poklesu funkčních hodnot funkce \(f\).

7.2. Sensitivita a elasticita projekční matice#

Zkoumejme model

\[x(n+1)=Ax(n)\]
věkově nebo jinak strukturované populace, kde \(A=(a_{ij})\) je projekční matice zobrazující současný stav na stav v následujícím časovém okamžiku. Bude nás zajímat, jak se bude měnit dynamika modelu při změnách parametrů. To je důležité pro nastavení správných postupů při ochraně živočišných druhů nebo při nastavení dlouhodobě udržitelné intenzity lovu či sklizně.

Rychlost růstu populace je dána dominantní vlastní hodnotou matice \(A\). Derivace \(\frac{\partial \lambda}{\partial a_{ij}}\) dominantní vlastní hodnoty \(\lambda\) podle prvku \(a_{ij}\) udává, jak rychle se mění tato vlastní hodnota při změnách prvku \(a_{ij}\). Jedná se o citlivost (sensitivitu) dominatní vlastní hodnoty na komponentu \(a_{ij}\) projekční matice \(A\). Podle [5] platí

\[ \frac{\partial \lambda}{\partial a_{ij}}=v_iw_j, \]
kde \(w\) a \(v\) jsou pravý a levý vlastní vektory příslušné vlastní hodnotě \(\lambda\) a jejichž skalární součin je roven jedné.

Změna \(\Delta \lambda\) vlastní hodnoty způsobená změnou parametru \(a_{ij}\) o \(\Delta a_{ij}\) je

\[\Delta \lambda = \frac{\partial \lambda}{\partial a_{ij}} \Delta a_{ij}.\]
Matice \(\left(\frac{\partial \lambda}{\partial a_{ij}}\right)\) se nazývá matice sensitivity nebo matice citlivosti dominantní vlastní hodnoty na složky projekční matice. Matice ukazuje, jaké mají jednotlivé parametry populace vliv na celkový růst. Zpravidla ji počítáme pomocí sloupcového pravého vlastního vektoru \(\vec w\) a řádkového levého vlastního vektoru \(\vec v\) pomocí vztahu
\[\left(\frac{\partial \lambda}{\partial a_{ij}}\right)=\vec v^T \vec w^T,\]
přičemž tyto vektory splňují \(\vec v \vec w = 1\) (skalární součin je roven jedné). Podmínka na skalární součin bývá splněna automaticky, pokud pokud levý vlastní vektor najdeme pomocí inverzní matice složené z pravých vlastních vektorů.

Často je vhodné pracovat s relativními změnami \(\frac{\Delta \lambda}{\lambda}\) a \(\frac{\Delta a_{ij}}{a_{ij}}\), pro které plyne vztah

\[\frac{\Delta \lambda}{\lambda} = \frac{\partial \lambda}{\partial a_{ij}} \frac{a_{ij}}{\lambda}\frac{\Delta a_{ij}}{a_{ij}},\]
kde výraz
\[\frac{\partial \lambda}{\partial a_{ij}} \frac{a_{ij}}{\lambda}\]
se nazývá elasticita dominantní vlastní hodnoty vzhledem ke složce \(a_{ij}\). Udává, o kolik procent se změní maximální vlastní číslo, pokud se koeficient \(a_{ij}\) změní o jedno procento. Součet všech elasticit je roven jedné a proto je možno elasticitu chápat jako míru, s jakou se jednotlivé komponenty matice projekce podílejí na růstu populace. Bývá výhodné elasticity uspořádat do matice
\[\left(\frac{\partial \lambda}{\partial a_{ij}} \frac{a_{ij}}{\lambda}\right),\]
která se nazývá matice elasticity.

Prozkoumejme matici elasticity na příkladu s populací tuleňů z minulé přednášky.

Hide code cell source
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
np.set_printoptions(suppress=True)

L = np.matrix([
    0, 1.26, 2.0,
    0.614, 0, 0,
    0, 0.808, 0.808
]).reshape(3,3)
v,P = np.linalg.eig(L)
l = max(v).real  # maximální vlastní hodnota
l_index = np.argmax(v) # pozice maximální vlastní hodnoty

W = P[:,l_index].real # pravý sloupcový vlastni vektor příslušný maximální vlastní hodnotě
V = (P**(-1))[l_index:].real # levý řádkový vlastní vektor příslušný maximální vlastní hodnotě
matice_sensitivity = V.T*W.T
matice_elasticity = np.matrix(np.array(matice_sensitivity)*np.array(L)/l)
matice_elasticity
matrix([[0.        , 0.10157076, 0.19054952],
        [0.29212028, 0.        , 0.        ],
        [0.        , 0.19054952, 0.22520993]])

Vidíme, že největší vliv na růst populace mají parametry spojené s přežíváním nejstarší věkové kategorie a s dospíváním mláďat, toto jsou dvě jediné hodnoty nad 20 procent. Pokud populace vymírá, snažíme se při ochraně tohoto druhu navýšit parametry s nejvyšší elasticitou. Změna těchto parametrů má totiž největší efekt. Fertilita jednotlivých věkových tříd (čísla v prvním řádku) nemá tak výrazný vliv a navíc ji zpravidla nemůžeme tak efektivně ovlivnit, jako koeficienty přežití.

Poznámka (Matice elasticity v ochraně přírody)

Podle [16] (str. 127) by se v ochraně přírody při ochraně ohrožených druhů měly navyšovat ty koeficienty, které jsou spojeny s největšími elasticitami. Při redukci populace škůdců se také věnujeme koeficientů s největšími elasticitami, ale tyto koeficienty se snažíme snižovat. Pokud populaci hospodářsky využíváme, poté lov či sklizeň zaměříme naopak na kategorie, kde jsou koeficienty s nejnižšími elasticitami. V těchto případech totiž vnější zásahy ovlivňují populaci relativně nejméně. Vždy je však nutné přihlédnout k tomu, že některé koeficienty ovlivníme relativně snadno, některé jenom s velkými ekonomickými náklady a některé prakticky ovlivnit moc nemůžeme.

Pro kontrolu můžeme vypočítat matici elasticity po složkách pomocí numerické aproximace derivace centrální diferencí. S každou komponentou matice posuneme o zvolený krok nahoru a dolů, v každém z těchto případů vypočteme dominantní vlastní hodnotu a určíme derivaci (rychlost změny) pomocí poměru změn a poté vypočteme elasticitu jako relativní rychlost změny. Přitom není nutné počítat derivace pro nulové komponenty matice, protože tyto derivace budou při přepočtu na relativní rychlost změny násobeny nulou a neuplatní se.

Hide code cell source
h = 0.01
matice = np.zeros([3,3])
matice
for f1 in range(3):
    for f2 in range(3):
        if L[f1,f2] == 0:
            matice[f1,f2] = 0
        else:
            L2 = L.copy()
            L3 = L.copy()
            L2[f1,f2] = L2[f1,f2] + h
            L3[f1,f2] = L3[f1,f2] - h
            v2 = np.linalg.eigvals(L2)
            v3 = np.linalg.eigvals(L3)
            matice[f1,f2] = ((max(v2)-max(v3))/(2*h)).real*L[f1,f2]/l
matice = np.matrix(matice)
matice
matrix([[0.        , 0.10157068, 0.19055025],
        [0.29212876, 0.        , 0.        ],
        [0.        , 0.19055404, 0.22521132]])

Součet komponent matice elasticity by měl být v rámci zaokrouhlovací chyby roven jedné.

Hide code cell source
np.sum(matice_elasticity)
1.0

Pokud potřebujeme regulovat rychlost růstu populace, má smysl soustředit se na ty komponenty matice elasticity, které jsou numericky největší. Jejich změna se na růstu populace projeví nejvíce. (Některé koeficienty projekční matice je ovšem z principu těžké nebo nemožné regulovat, potom se soustředíme na ty, které regulovat můžeme.)

Koeficienty matice elasticity můžeme navzájem sčítat a součty interpretovat jako agregovaný podíl na růstu populace. Například součty po sloupcích odhalí, jak se jednotlivé věkové kategorie podílejí na celkovém růstu.

Hide code cell source
np.sum(matice, axis=0)
matrix([[0.29212876, 0.29212473, 0.41576157]])

7.3. Lineární aproximace#

Jedna z úspěšných inženýrských metod je nahrazení nelineární funkce pomocí funkce lineární. Linearizace nám objasní, proč je v naprosté většině fyzikálních zákonů přímá úměrnost resp. přímá úměrnost nějaké mocnině. Není v tom nic složitějšího než prostá selská logika představená v následující poznámce.

Poznámka (Lineární aproximace selskou logikou)

Uvažujme modelový případ populace jelenů, která má v roce 2020 celkem 1023 členů a roste rychlostí 12 jedinců za rok. (Čistý růst, tedy přírůstek snížený o úmrtí a lov.) Odhadneme velikost populace v budoucích letech. V roce \(x\) je populace celkem \(\Delta x = x-2020\) let od roku 2020. V tomto roce bude velikost \(N(x)\) populace dána součtem velikosti v roce 2020 a přírůstku za daný počet let. Tento přírůstek je roven přírůstku za rok vynásobenému počtem let, tedy

\[\Delta N = 12 \Delta x.\]
Velikost populace v roce \(x\) je
\[N(x)=N(2020)+\Delta N = 1023+12\Delta x = 1023 + 12(x-2020).\]

V předchozím příkladě je rychlost růstu vlastně derivace velikosti populace podle času. Stejným postupem můžeme určit časový vývoj libovolné veličiny, u níž je znám výchozí stav a rychlost růstu. Má to pouze jednu vadu na kráse. Málokdy se totiž setkáme s tím, že rychlost růstu je konstantní. Proto uvedená myšlenka bude platit jenom přibližně a jenom lokálně v okolí bodu, z nějž vycházíme. To proto, že na krátkém intervalu se rychlost růstu funkce zpravidla nestihne podstatně změnit. Dostáváme následující přibližné vzorce pro změnu funkce a pro funkční předpis v okolí bodu aproximace.

\[\Delta f\approx \frac{\mathrm {d}f}{\mathrm {d}x}(x_0)\Delta x \]
\[f(x)\approx f(x_0) + \frac{\mathrm {d}f}{\mathrm {d}x}(x_0)(x-x_0)\]

Geometricky uvedená aproximace znamená, že křivku definovanou grafem funkce nahrazujeme přímkou. Tímto jsou dány i meze použitelnosti: metoda výborně funguje pro funkce, jejichž derivace se mění pomalu a graf připomíná lineární závislost. Kvantitativně to je možno vyjádřit tak, že derivace má malou derivaci, tj. že druhá derivace je malá (ve smyslu numericky blízko k nule, tj. že absolutní hodnota druhé derivace je malá).

Pokud funkce závisí na dvou nebo více proměnných, jednotlivé příspěvky ke změně sečteme.

\[\Delta f\approx \frac{\partial f}{\partial x}(x_0,y_0)\Delta x + \frac{\partial f}{\partial x}(x_0,y_0)\Delta y \]

\[f(x,y)\approx f(x_0,y_0) + \frac{\partial f}{\partial x}(x_0,y_0)(x-x_0) + \frac{\partial f}{\partial x}(x_0,y_0)(y-y_0) \]

Toto je možno zapsat vektorově pomocí gradientu a součinu řádkové a sloupcové matice ve tvaru

\[f(x,y)\approx f(x_0,y_0) + \nabla f(x_0,y_0)\begin{pmatrix}x-x_0\\y-y_0\end{pmatrix}.\]

Je-li \(\vec F(x,y)=(f_1(x,y),f_2(x,y))\) vektorová funkce dvou proměnných, je možné napsat linearizaci podle výše uvedeného vzorce pro každou komponentu samostatně a v maticovém zápisu máme

\[\vec F(x,y)=\vec F(x_0,y_0) + J(x_0,y_0)\begin{pmatrix}x-x_0\\y-y_0\end{pmatrix} ,\]
kde
\[ J(x,y)=\begin{pmatrix}\frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y}\\\frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y}\end{pmatrix}\]
je Jacobiho matice zobrazení \(\vec F\)

Jacobiho matice má tedy v řádcích gradienty jednotlivých komponent.

V praxi je často vstupem funkce nějaký podnět a výstupem odezva na tento podnět. Funkční závislosti tohoto typu se nazývají konstituční vztahy. Protože v těchto typech závislostí zpravidla platí, že bez podnětu není odezva, studované funkční závislosti procházejí bodem \([0,0]\), tj. počátkem. Lineární aproximace v tomto bodě je poté lineární funkce procházející počátkem a tedy přímá úměrnost. To je i důvod, proč se v konstitučních vztazích tak často přímá úměrnost vyskytuje.

7.4. Rovnice vedení tepla#

Rovnice vedení tepla popisuje vedení tepla v jednorozměrném objektu, například v tyči. Teplo se předává z teplejšího místa do místa o nižší teplotě a tím se v čase teplejší místo ochlazuje a chladnější otepluje.

Vedení tepla je proces, který si většina z nás umí snadno představit a má s ním praktické zkušenosti. Model je však plně univerzální a hodí se pro libovolný transportní děj v přírodě. To zahrnuje nejenom transport energie ve formě tepla, ale i transport látek, například pomocí difuze nebo prouděním v pórovitém prostředí. Tím je podchyceno například šíření nečistot ze zdroje do okolí, proudění podzemní vody, sušení dřeva ale i rozšiřování chorob či invazních druhů v ekosystému. To si ukážeme později (Fisherova–KPP rovnice).

Shrňme si nejprve základní fyzikální představy o vedení tepla.

  • Rychlost růstu teploty v čase je úměrná rychlosti, s jakou je do daného místa dodáváno teplo. Pokud do daného místa dodáváme teplo, teplota roste tím rychleji, čím více tepla dodáváme za jednotku času.

  • Rychlost, s jakou je do daného místa dodáváno teplo, se určí jako pokles toku tepla podél tyče v daném místě.

  • Tok tepla je úměrný rychlosti, s jakou klesá teplo podél tyče.

Rychlost růstu je vždy vyjádřena derivací, rychlost poklesu záporně vzatou derivací. Pokud sledujeme, jak se veličina mění v čase, jedná se o derivaci podle času, pokud sledujeme změnu podél tyče, jedná se o derivaci podle polohy. Tímto jsme si připravili vše potřebné k formulaci matematického modelu.

Teplota \(T\) souvisí s polohou a s časem, jedná se tedy o funkci dvou proměnných a pro vyjádření rychlosti změny musíme použít parciální derivace. Parciální derivace teploty podle času \(\frac{\partial T}{\partial t}\) udává, jak rychle s daném místě teplota roste s časem a souvisí s rychlostí dodávání tepla. Rychlost dodávání tepla je možné měřit tím, jak v daném místě slábne tok tepla \(q\), tj. jaká je záporně vzatá derivace toku tepla podle polohy. Konstanty úměrnosti dodá fyzika, jedná se o hustotu \(\varrho\) a měrnou tepelnou kapacitu \(c\), tj.

\[ \varrho c \frac{\partial T}{\partial t}=-\frac{\partial q}{\partial x}.\]

Teplo teče z místa s větší teplotou do místa s menší teplotou a teče intenzivněji, pokud je mezi těmito místy velký teplotní rozdíl. Pokles tepla podél tyče je dán záporně vzatou derivací teploty podle polohy a tok je tomuto poklesu úměrný (s konstantou úměrnosti, která se nazývá specifická vodivost a značí \(\lambda\)). Tedy

\[ q = -\lambda \frac{\partial T}{\partial x}.\]

Spojením obou vztahů dostáváme rovnici vedení tepla ve tvaru

\[ \varrho c \frac{\partial T}{\partial t}=\frac{\partial }{\partial x} \left(\lambda \frac{\partial T}{\partial x}\right). \]
Je-li \(\lambda\) konstantní (nezávislá na teplotě a poloze), potom je možné rovnici upravit na
\[ \varrho c \frac{\partial T}{\partial t}=\lambda\frac{\partial^2 T}{\partial x^2}, \]
kde \(\frac{\partial^2 T}{\partial x^2} = \frac{\partial }{\partial x} \left(\frac{\partial T}{\partial x}\right)\) je druhá derivace teploty podle polohy.

7.5. Fisherova–KPP rovnice#

Analogický postup jako u vedení tepla funguje pro libovolný transportní děj. Používá se například pro model šíření živočišného druhu v životním prostředí nebo šíření genu v populaci. Ukážeme si zde lineární verzi pro jednorozměrné životní prostředí, například podél řeky. Vícerozměrné zobecnění si uvedeme na konci semestru.

Tento model se nazývá Fisherova–KPP rovnice a jedná se vlastně o rovnici vedení tepla, do které jsou doplněny zdroje.

Fisherova–KPP rovnice popisuje populaci, která se může šířit v prostoru. Kromě časové závislosti tedy musíme uvažovat i závislost na poloze. Tato rovnice je vyjádřena pro funkci vyjadřující hustotu populace. Hustotou populace může být například množství jedinců na jednotku délky, v případě zobecnění na více dimenzí množství jedinců na jednotku plochy. Procesy vedoucí k migraci a tedy i změně hustoty populace jsou stejné jako u vedení tepla. Populace má tendenci migrovat z míst, kde je vysoká hustota do míst, kde je hustota menší a v místě, kde převažuje imigrace nad emigrací velikost populace roste. Protože se populace může rozmnožovat, je do rovnice navíc započítán člen modelující vlastní reprodukci. Zpravidla je uvažován logistický růst a rovnice má poté tvar

\[\frac{\partial u}{\partial t}=D\frac{\partial^2 u}{\partial x^2}+ru\left(1-\frac uK\right).\]
Tato rovnice byla původně navržena jako model šíření výhodného genu populací. To potvrzuje, že v rovnici máme opravdu nástroj pro šíření nebo transport různých objektů, od energie, přes molekuly nebo částice v látce či geny v populaci až po jedince invazního druhu v ekosystému.